1 /* -------------------------------------------------------------
2
3 This file is a component of SDPA
4 Copyright (C) 2004-2012 SDPA Project
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 ------------------------------------------------------------- */
21 /*-----------------------------------------------
22 rsdpa_dpotrf.cpp
23 modification of ATL_dpotrfL
24 for dealing with numerical error
25 in diagonal elements.
26
27 int rATL_dpotrfL(int N, double *A,int lda)
28
29 modified by Makoto Yamshita 2002.07.11
30 -----------------------------------------------*/
31 #define POTRF_NONZERO (1.0e-14)
32 #define POTRF_ASSIGN (1.0e+100)
33 #define POTRF_LIMIT (-1.0e-6)
34
35 /*
36 * Automatically Tuned Linear Algebra Software v3.4.0
37 * (C) Copyright 1999 R. Clint Whaley
38 *
39 * Redistribution and use in source and binary forms, with or without
40 * modification, are permitted provided that the following conditions
41 * are met:
42 * 1. Redistributions of source code must retain the above copyright
43 * notice, this list of conditions and the following disclaimer.
44 * 2. Redistributions in binary form must reproduce the above copyright
45 * notice, this list of conditions, and the following disclaimer in the
46 * documentation and/or other materials provided with the distribution.
47 * 3. The name of the ATLAS group or the names of its contributers may
48 * not be used to endorse or promote products derived from this
49 * software without specific written permission.
50 *
51 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
52 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
53 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
54 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
55 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
56 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
57 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
58 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
59 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
60 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
61 * POSSIBILITY OF SUCH DAMAGE.
62 *
63 */
64
65
66
67 #include "sdpa_include.h"
68 #include "sdpa_algebra.h"
69 #if 0
70 #define CHOLESKY_ADJUST(val) rMessage("Choleksy adjust from " << val << " to " << POTRF_NONZERO);
71 #else
72 #define CHOLESKY_ADJUST(val) ;
73 #endif
74
75 namespace sdpa {
76
77 extern "C" {
potrf4(double * A,const int n)78 static int potrf4(double* A,const int n)
79 {
80 double* A1 = A+n+1;
81 double* A2 = A1+n+1;
82 double* A3 = A2+n+1;
83 double L11 = *A;
84 double L21 = A[1], L22 = *A1;
85 double L31 = A[2], L32 = A1[1], L33 = *A2;
86 double L41 = A[3], L42 = A1[2], L43 = A2[1], L44 = *A3;
87
88 if (L11 < POTRF_LIMIT) {
89 return 1;
90 }
91 if (L11 < POTRF_NONZERO) {
92 CHOLESKY_ADJUST(L11);
93 L11 = POTRF_ASSIGN;
94 }
95
96 *A = L11 = sqrt(L11);
97 L11 = 1.0/L11;
98 L21 *= L11;
99 L31 *= L11;
100 L41 *= L11;
101
102 L22 -= L21*L21;
103 if (L22 < POTRF_LIMIT) {
104 return 2;
105 }
106 if (L22 < POTRF_NONZERO) {
107 CHOLESKY_ADJUST(L22);
108 L22 = POTRF_ASSIGN;
109 }
110 *A1 = L22 = sqrt(L22);
111 L22 = 1.0/L22;
112 L32 = (L32 - L31*L21)*L22;
113 L42 = (L42 - L41*L21)*L22;
114
115 L33 -= L31*L31 + L32*L32;
116 if (L33 < POTRF_LIMIT) {
117 return 3;
118 }
119 if (L33 < POTRF_NONZERO) {
120 CHOLESKY_ADJUST(L33);
121 L33 = POTRF_ASSIGN;
122 }
123 *A2 = L33 = sqrt(L33);
124 L43 = (L43-L41*L31-L42*L32)/L33;
125 L44 -= L41*L41 + L42*L42 + L43*L43;
126 if (L44 < POTRF_LIMIT) {
127 return 4;
128 }
129 if (L44 < POTRF_NONZERO) {
130 CHOLESKY_ADJUST(L44);
131 L44 = POTRF_ASSIGN;
132 }
133 *A3 = sqrt(L44);
134
135 A[1] = L21;
136 A[2] = L31; A1[1] = L32;
137 A[3] = L41; A1[2] = L42; A2[1] = L43;
138 return 0;
139 }
140
potrf3(double * A,const int n)141 static int potrf3(double* A,const int n)
142 {
143 double* A1 = A+n+1;
144 double* A2 = A1+n+1;
145 double L11 = *A;
146 double L21 = A[1], L22 = *A1;
147 double L31 = A[2], L32 = A1[1], L33 = *A2;
148
149 if (L11 < POTRF_LIMIT) {
150 return 1;
151 }
152 if (L11 < POTRF_NONZERO) {
153 CHOLESKY_ADJUST(L11);
154 L11 = POTRF_ASSIGN;
155 }
156
157 *A = L11 = sqrt(L11);
158 L11 = 1.0/L11;
159 L21 *= L11;
160 L31 *= L11;
161
162 L22 -= L21*L21;
163 if (L22 < POTRF_LIMIT) {
164 return 2;
165 }
166 if (L22 < POTRF_NONZERO) {
167 CHOLESKY_ADJUST(L22);
168 L22 = POTRF_ASSIGN;
169 }
170 L22 = sqrt(L22);
171 L32 = (L32 - L31*L21)/L22;
172
173 L33 -= L31*L31 + L32*L32;
174 if (L33 < POTRF_LIMIT) {
175 return 3;
176 }
177 if (L33 < POTRF_NONZERO) {
178 CHOLESKY_ADJUST(L33);
179 L33 = POTRF_ASSIGN;
180 }
181 *A2 = sqrt(L33);
182
183 A[1] = L21; *A1 = L22;
184 A[2] = L31; A1[1] = L32;
185 return 0;
186 }
187
potrf2(double * A,const int n)188 static int potrf2(double* A,const int n)
189 {
190 double* A1 = A+n+1;
191 double L11 = *A;
192 double L21 = A[1], L22 = *A1;
193
194 if (L11 < POTRF_LIMIT) {
195 return 1;
196 }
197 if (L11 < POTRF_NONZERO) {
198 CHOLESKY_ADJUST(L11);
199 L11 = POTRF_ASSIGN;
200 }
201
202 *A = L11 = sqrt(L11);
203 L21 /= L11;
204 L22 -= L21*L21;
205 if (L22 < POTRF_LIMIT) {
206 return 2;
207 }
208 if (L22 < POTRF_NONZERO) {
209 CHOLESKY_ADJUST(L22);
210 L22 = POTRF_ASSIGN;
211 }
212 *A = L11;
213 A[1] = L21; *A1 = sqrt(L22);
214
215 return 0;
216 }
217
rATL_dpotrfL(int N,double * A,int lda)218 int rATL_dpotrfL(int N, double *A,int lda)
219 {
220 double *An, *Ar;
221 int Nleft, Nright, ierr;
222
223 if (N > 4) {
224 Nleft = N >> 1;
225 #if 0
226 int nb = ilaenv_fc(&IONE, "DPOTRF", "L", &N,
227 &IMONE,&IONE, &IMONE, strlen("DPOTRF"), strlen("L"));
228 if (Nleft > nb<<1) Nleft = (Nleft/nb)*nb;
229 #endif
230 #if 0
231 if (Nleft > 64) {
232 Nleft = 64;
233 }
234 #endif
235 Nright = N - Nleft;
236 ierr = rATL_dpotrfL(Nleft, A,lda);
237 if (!ierr) {
238 Ar = A + Nleft;
239 An = Ar + lda * Nleft;
240 dtrsm_fc ((char *)"R",(char *)"L",(char *)"T",(char *)"N",
241 &Nright,&Nleft,&DONE,A,&lda,
242 Ar, &lda, strlen("R"),strlen("L"),
243 strlen("T"),strlen("N"));
244 dsyrk_fc ((char *)"L",(char *)"N",&Nright,&Nleft,&DMONE,
245 Ar, &lda, &DONE,An,&lda,strlen("L"),strlen("N"));
246 ierr = rATL_dpotrfL(Nright, An,lda);
247 if (ierr) return(ierr+Nleft);
248 }
249 else return(ierr);
250 }
251 else if (N==4) return(potrf4(A,lda));
252 else if (N==3) return(potrf3(A,lda));
253 else if (N==2) return(potrf2(A,lda));
254 else if (N==1) {
255 if (*A < POTRF_LIMIT) {
256 return 1;
257 }
258 if (*A < POTRF_NONZERO) {
259 CHOLESKY_ADJUST(*A);
260 *A = POTRF_ASSIGN;
261 }
262 *A = sqrt(*A);
263 }
264 return(0);
265 }
266
rdpotrfl_(int * N,double * A,int * lda,int * info)267 void rdpotrfl_(int* N, double *A,int* lda,int* info)
268 {
269 *info = rATL_dpotrfL(*N,A,*lda);
270 }
271
272 }; // end of extern "C"
273
274 } // end of namespace 'sdpa'
275