SphinxBase 0.6
|
00001 /* 00002 NOTE: This is generated code. Look in README.python for information on 00003 remaking this file. 00004 */ 00005 #include "sphinxbase/f2c.h" 00006 00007 #ifdef HAVE_CONFIG 00008 #include "config.h" 00009 #else 00010 extern doublereal slamch_(char *); 00011 #define EPSILON slamch_("Epsilon") 00012 #define SAFEMINIMUM slamch_("Safe minimum") 00013 #define PRECISION slamch_("Precision") 00014 #define BASE slamch_("Base") 00015 #endif 00016 00017 00018 extern doublereal slapy2_(real *, real *); 00019 00020 00021 00022 /* Table of constant values */ 00023 00024 static integer c__0 = 0; 00025 static real c_b163 = 0.f; 00026 static real c_b164 = 1.f; 00027 static integer c__1 = 1; 00028 static real c_b181 = -1.f; 00029 static integer c_n1 = -1; 00030 00031 integer ieeeck_(integer *ispec, real *zero, real *one) 00032 { 00033 /* System generated locals */ 00034 integer ret_val; 00035 00036 /* Local variables */ 00037 static real nan1, nan2, nan3, nan4, nan5, nan6, neginf, posinf, negzro, 00038 newzro; 00039 00040 00041 /* 00042 -- LAPACK auxiliary routine (version 3.0) -- 00043 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00044 Courant Institute, Argonne National Lab, and Rice University 00045 June 30, 1998 00046 00047 00048 Purpose 00049 ======= 00050 00051 IEEECK is called from the ILAENV to verify that Infinity and 00052 possibly NaN arithmetic is safe (i.e. will not trap). 00053 00054 Arguments 00055 ========= 00056 00057 ISPEC (input) INTEGER 00058 Specifies whether to test just for inifinity arithmetic 00059 or whether to test for infinity and NaN arithmetic. 00060 = 0: Verify infinity arithmetic only. 00061 = 1: Verify infinity and NaN arithmetic. 00062 00063 ZERO (input) REAL 00064 Must contain the value 0.0 00065 This is passed to prevent the compiler from optimizing 00066 away this code. 00067 00068 ONE (input) REAL 00069 Must contain the value 1.0 00070 This is passed to prevent the compiler from optimizing 00071 away this code. 00072 00073 RETURN VALUE: INTEGER 00074 = 0: Arithmetic failed to produce the correct answers 00075 = 1: Arithmetic produced the correct answers 00076 */ 00077 00078 ret_val = 1; 00079 00080 posinf = *one / *zero; 00081 if (posinf <= *one) { 00082 ret_val = 0; 00083 return ret_val; 00084 } 00085 00086 neginf = -(*one) / *zero; 00087 if (neginf >= *zero) { 00088 ret_val = 0; 00089 return ret_val; 00090 } 00091 00092 negzro = *one / (neginf + *one); 00093 if (negzro != *zero) { 00094 ret_val = 0; 00095 return ret_val; 00096 } 00097 00098 neginf = *one / negzro; 00099 if (neginf >= *zero) { 00100 ret_val = 0; 00101 return ret_val; 00102 } 00103 00104 newzro = negzro + *zero; 00105 if (newzro != *zero) { 00106 ret_val = 0; 00107 return ret_val; 00108 } 00109 00110 posinf = *one / newzro; 00111 if (posinf <= *one) { 00112 ret_val = 0; 00113 return ret_val; 00114 } 00115 00116 neginf *= posinf; 00117 if (neginf >= *zero) { 00118 ret_val = 0; 00119 return ret_val; 00120 } 00121 00122 posinf *= posinf; 00123 if (posinf <= *one) { 00124 ret_val = 0; 00125 return ret_val; 00126 } 00127 00128 00129 /* Return if we were only asked to check infinity arithmetic */ 00130 00131 if (*ispec == 0) { 00132 return ret_val; 00133 } 00134 00135 nan1 = posinf + neginf; 00136 00137 nan2 = posinf / neginf; 00138 00139 nan3 = posinf / posinf; 00140 00141 nan4 = posinf * *zero; 00142 00143 nan5 = neginf * negzro; 00144 00145 nan6 = nan5 * 0.f; 00146 00147 if (nan1 == nan1) { 00148 ret_val = 0; 00149 return ret_val; 00150 } 00151 00152 if (nan2 == nan2) { 00153 ret_val = 0; 00154 return ret_val; 00155 } 00156 00157 if (nan3 == nan3) { 00158 ret_val = 0; 00159 return ret_val; 00160 } 00161 00162 if (nan4 == nan4) { 00163 ret_val = 0; 00164 return ret_val; 00165 } 00166 00167 if (nan5 == nan5) { 00168 ret_val = 0; 00169 return ret_val; 00170 } 00171 00172 if (nan6 == nan6) { 00173 ret_val = 0; 00174 return ret_val; 00175 } 00176 00177 return ret_val; 00178 } /* ieeeck_ */ 00179 00180 integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1, 00181 integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen 00182 opts_len) 00183 { 00184 /* System generated locals */ 00185 integer ret_val; 00186 00187 /* Builtin functions */ 00188 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen); 00189 integer s_cmp(char *, char *, ftnlen, ftnlen); 00190 00191 /* Local variables */ 00192 static integer i__; 00193 static char c1[1], c2[2], c3[3], c4[2]; 00194 static integer ic, nb, iz, nx; 00195 static logical cname, sname; 00196 static integer nbmin; 00197 extern integer ieeeck_(integer *, real *, real *); 00198 static char subnam[6]; 00199 00200 00201 /* 00202 -- LAPACK auxiliary routine (version 3.0) -- 00203 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00204 Courant Institute, Argonne National Lab, and Rice University 00205 June 30, 1999 00206 00207 00208 Purpose 00209 ======= 00210 00211 ILAENV is called from the LAPACK routines to choose problem-dependent 00212 parameters for the local environment. See ISPEC for a description of 00213 the parameters. 00214 00215 This version provides a set of parameters which should give good, 00216 but not optimal, performance on many of the currently available 00217 computers. Users are encouraged to modify this subroutine to set 00218 the tuning parameters for their particular machine using the option 00219 and problem size information in the arguments. 00220 00221 This routine will not function correctly if it is converted to all 00222 lower case. Converting it to all upper case is allowed. 00223 00224 Arguments 00225 ========= 00226 00227 ISPEC (input) INTEGER 00228 Specifies the parameter to be returned as the value of 00229 ILAENV. 00230 = 1: the optimal blocksize; if this value is 1, an unblocked 00231 algorithm will give the best performance. 00232 = 2: the minimum block size for which the block routine 00233 should be used; if the usable block size is less than 00234 this value, an unblocked routine should be used. 00235 = 3: the crossover point (in a block routine, for N less 00236 than this value, an unblocked routine should be used) 00237 = 4: the number of shifts, used in the nonsymmetric 00238 eigenvalue routines 00239 = 5: the minimum column dimension for blocking to be used; 00240 rectangular blocks must have dimension at least k by m, 00241 where k is given by ILAENV(2,...) and m by ILAENV(5,...) 00242 = 6: the crossover point for the SVD (when reducing an m by n 00243 matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds 00244 this value, a QR factorization is used first to reduce 00245 the matrix to a triangular form.) 00246 = 7: the number of processors 00247 = 8: the crossover point for the multishift QR and QZ methods 00248 for nonsymmetric eigenvalue problems. 00249 = 9: maximum size of the subproblems at the bottom of the 00250 computation tree in the divide-and-conquer algorithm 00251 (used by xGELSD and xGESDD) 00252 =10: ieee NaN arithmetic can be trusted not to trap 00253 =11: infinity arithmetic can be trusted not to trap 00254 00255 NAME (input) CHARACTER*(*) 00256 The name of the calling subroutine, in either upper case or 00257 lower case. 00258 00259 OPTS (input) CHARACTER*(*) 00260 The character options to the subroutine NAME, concatenated 00261 into a single character string. For example, UPLO = 'U', 00262 TRANS = 'T', and DIAG = 'N' for a triangular routine would 00263 be specified as OPTS = 'UTN'. 00264 00265 N1 (input) INTEGER 00266 N2 (input) INTEGER 00267 N3 (input) INTEGER 00268 N4 (input) INTEGER 00269 Problem dimensions for the subroutine NAME; these may not all 00270 be required. 00271 00272 (ILAENV) (output) INTEGER 00273 >= 0: the value of the parameter specified by ISPEC 00274 < 0: if ILAENV = -k, the k-th argument had an illegal value. 00275 00276 Further Details 00277 =============== 00278 00279 The following conventions have been used when calling ILAENV from the 00280 LAPACK routines: 00281 1) OPTS is a concatenation of all of the character options to 00282 subroutine NAME, in the same order that they appear in the 00283 argument list for NAME, even if they are not used in determining 00284 the value of the parameter specified by ISPEC. 00285 2) The problem dimensions N1, N2, N3, N4 are specified in the order 00286 that they appear in the argument list for NAME. N1 is used 00287 first, N2 second, and so on, and unused problem dimensions are 00288 passed a value of -1. 00289 3) The parameter value returned by ILAENV is checked for validity in 00290 the calling subroutine. For example, ILAENV is used to retrieve 00291 the optimal blocksize for STRTRI as follows: 00292 00293 NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) 00294 IF( NB.LE.1 ) NB = MAX( 1, N ) 00295 00296 ===================================================================== 00297 */ 00298 00299 00300 switch (*ispec) { 00301 case 1: goto L100; 00302 case 2: goto L100; 00303 case 3: goto L100; 00304 case 4: goto L400; 00305 case 5: goto L500; 00306 case 6: goto L600; 00307 case 7: goto L700; 00308 case 8: goto L800; 00309 case 9: goto L900; 00310 case 10: goto L1000; 00311 case 11: goto L1100; 00312 } 00313 00314 /* Invalid value for ISPEC */ 00315 00316 ret_val = -1; 00317 return ret_val; 00318 00319 L100: 00320 00321 /* Convert NAME to upper case if the first character is lower case. */ 00322 00323 ret_val = 1; 00324 s_copy(subnam, name__, (ftnlen)6, name_len); 00325 ic = *(unsigned char *)subnam; 00326 iz = 'Z'; 00327 if (iz == 90 || iz == 122) { 00328 00329 /* ASCII character set */ 00330 00331 if (ic >= 97 && ic <= 122) { 00332 *(unsigned char *)subnam = (char) (ic - 32); 00333 for (i__ = 2; i__ <= 6; ++i__) { 00334 ic = *(unsigned char *)&subnam[i__ - 1]; 00335 if (ic >= 97 && ic <= 122) { 00336 *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32); 00337 } 00338 /* L10: */ 00339 } 00340 } 00341 00342 } else if (iz == 233 || iz == 169) { 00343 00344 /* EBCDIC character set */ 00345 00346 if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 162 && 00347 ic <= 169) { 00348 *(unsigned char *)subnam = (char) (ic + 64); 00349 for (i__ = 2; i__ <= 6; ++i__) { 00350 ic = *(unsigned char *)&subnam[i__ - 1]; 00351 if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 00352 162 && ic <= 169) { 00353 *(unsigned char *)&subnam[i__ - 1] = (char) (ic + 64); 00354 } 00355 /* L20: */ 00356 } 00357 } 00358 00359 } else if (iz == 218 || iz == 250) { 00360 00361 /* Prime machines: ASCII+128 */ 00362 00363 if (ic >= 225 && ic <= 250) { 00364 *(unsigned char *)subnam = (char) (ic - 32); 00365 for (i__ = 2; i__ <= 6; ++i__) { 00366 ic = *(unsigned char *)&subnam[i__ - 1]; 00367 if (ic >= 225 && ic <= 250) { 00368 *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32); 00369 } 00370 /* L30: */ 00371 } 00372 } 00373 } 00374 00375 *(unsigned char *)c1 = *(unsigned char *)subnam; 00376 sname = *(unsigned char *)c1 == 'S' || *(unsigned char *)c1 == 'D'; 00377 cname = *(unsigned char *)c1 == 'C' || *(unsigned char *)c1 == 'Z'; 00378 if (! (cname || sname)) { 00379 return ret_val; 00380 } 00381 s_copy(c2, subnam + 1, (ftnlen)2, (ftnlen)2); 00382 s_copy(c3, subnam + 3, (ftnlen)3, (ftnlen)3); 00383 s_copy(c4, c3 + 1, (ftnlen)2, (ftnlen)2); 00384 00385 switch (*ispec) { 00386 case 1: goto L110; 00387 case 2: goto L200; 00388 case 3: goto L300; 00389 } 00390 00391 L110: 00392 00393 /* 00394 ISPEC = 1: block size 00395 00396 In these examples, separate code is provided for setting NB for 00397 real and complex. We assume that NB will take the same value in 00398 single or double precision. 00399 */ 00400 00401 nb = 1; 00402 00403 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) { 00404 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00405 if (sname) { 00406 nb = 64; 00407 } else { 00408 nb = 64; 00409 } 00410 } else if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, 00411 "RQF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen) 00412 3, (ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) 00413 == 0) { 00414 if (sname) { 00415 nb = 32; 00416 } else { 00417 nb = 32; 00418 } 00419 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) { 00420 if (sname) { 00421 nb = 32; 00422 } else { 00423 nb = 32; 00424 } 00425 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) { 00426 if (sname) { 00427 nb = 32; 00428 } else { 00429 nb = 32; 00430 } 00431 } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) { 00432 if (sname) { 00433 nb = 64; 00434 } else { 00435 nb = 64; 00436 } 00437 } 00438 } else if (s_cmp(c2, "PO", (ftnlen)2, (ftnlen)2) == 0) { 00439 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00440 if (sname) { 00441 nb = 64; 00442 } else { 00443 nb = 64; 00444 } 00445 } 00446 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) { 00447 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00448 if (sname) { 00449 nb = 64; 00450 } else { 00451 nb = 64; 00452 } 00453 } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00454 nb = 32; 00455 } else if (sname && s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) { 00456 nb = 64; 00457 } 00458 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) { 00459 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00460 nb = 64; 00461 } else if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00462 nb = 32; 00463 } else if (s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) { 00464 nb = 64; 00465 } 00466 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) { 00467 if (*(unsigned char *)c3 == 'G') { 00468 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00469 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00470 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00471 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00472 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00473 ftnlen)2, (ftnlen)2) == 0) { 00474 nb = 32; 00475 } 00476 } else if (*(unsigned char *)c3 == 'M') { 00477 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00478 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00479 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00480 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00481 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00482 ftnlen)2, (ftnlen)2) == 0) { 00483 nb = 32; 00484 } 00485 } 00486 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) { 00487 if (*(unsigned char *)c3 == 'G') { 00488 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00489 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00490 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00491 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00492 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00493 ftnlen)2, (ftnlen)2) == 0) { 00494 nb = 32; 00495 } 00496 } else if (*(unsigned char *)c3 == 'M') { 00497 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00498 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00499 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00500 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00501 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00502 ftnlen)2, (ftnlen)2) == 0) { 00503 nb = 32; 00504 } 00505 } 00506 } else if (s_cmp(c2, "GB", (ftnlen)2, (ftnlen)2) == 0) { 00507 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00508 if (sname) { 00509 if (*n4 <= 64) { 00510 nb = 1; 00511 } else { 00512 nb = 32; 00513 } 00514 } else { 00515 if (*n4 <= 64) { 00516 nb = 1; 00517 } else { 00518 nb = 32; 00519 } 00520 } 00521 } 00522 } else if (s_cmp(c2, "PB", (ftnlen)2, (ftnlen)2) == 0) { 00523 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00524 if (sname) { 00525 if (*n2 <= 64) { 00526 nb = 1; 00527 } else { 00528 nb = 32; 00529 } 00530 } else { 00531 if (*n2 <= 64) { 00532 nb = 1; 00533 } else { 00534 nb = 32; 00535 } 00536 } 00537 } 00538 } else if (s_cmp(c2, "TR", (ftnlen)2, (ftnlen)2) == 0) { 00539 if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) { 00540 if (sname) { 00541 nb = 64; 00542 } else { 00543 nb = 64; 00544 } 00545 } 00546 } else if (s_cmp(c2, "LA", (ftnlen)2, (ftnlen)2) == 0) { 00547 if (s_cmp(c3, "UUM", (ftnlen)3, (ftnlen)3) == 0) { 00548 if (sname) { 00549 nb = 64; 00550 } else { 00551 nb = 64; 00552 } 00553 } 00554 } else if (sname && s_cmp(c2, "ST", (ftnlen)2, (ftnlen)2) == 0) { 00555 if (s_cmp(c3, "EBZ", (ftnlen)3, (ftnlen)3) == 0) { 00556 nb = 1; 00557 } 00558 } 00559 ret_val = nb; 00560 return ret_val; 00561 00562 L200: 00563 00564 /* ISPEC = 2: minimum block size */ 00565 00566 nbmin = 2; 00567 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) { 00568 if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", ( 00569 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, ( 00570 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0) 00571 { 00572 if (sname) { 00573 nbmin = 2; 00574 } else { 00575 nbmin = 2; 00576 } 00577 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) { 00578 if (sname) { 00579 nbmin = 2; 00580 } else { 00581 nbmin = 2; 00582 } 00583 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) { 00584 if (sname) { 00585 nbmin = 2; 00586 } else { 00587 nbmin = 2; 00588 } 00589 } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) { 00590 if (sname) { 00591 nbmin = 2; 00592 } else { 00593 nbmin = 2; 00594 } 00595 } 00596 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) { 00597 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00598 if (sname) { 00599 nbmin = 8; 00600 } else { 00601 nbmin = 8; 00602 } 00603 } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00604 nbmin = 2; 00605 } 00606 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) { 00607 if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00608 nbmin = 2; 00609 } 00610 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) { 00611 if (*(unsigned char *)c3 == 'G') { 00612 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00613 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00614 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00615 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00616 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00617 ftnlen)2, (ftnlen)2) == 0) { 00618 nbmin = 2; 00619 } 00620 } else if (*(unsigned char *)c3 == 'M') { 00621 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00622 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00623 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00624 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00625 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00626 ftnlen)2, (ftnlen)2) == 0) { 00627 nbmin = 2; 00628 } 00629 } 00630 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) { 00631 if (*(unsigned char *)c3 == 'G') { 00632 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00633 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00634 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00635 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00636 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00637 ftnlen)2, (ftnlen)2) == 0) { 00638 nbmin = 2; 00639 } 00640 } else if (*(unsigned char *)c3 == 'M') { 00641 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00642 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00643 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00644 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00645 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00646 ftnlen)2, (ftnlen)2) == 0) { 00647 nbmin = 2; 00648 } 00649 } 00650 } 00651 ret_val = nbmin; 00652 return ret_val; 00653 00654 L300: 00655 00656 /* ISPEC = 3: crossover point */ 00657 00658 nx = 0; 00659 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) { 00660 if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", ( 00661 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, ( 00662 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0) 00663 { 00664 if (sname) { 00665 nx = 128; 00666 } else { 00667 nx = 128; 00668 } 00669 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) { 00670 if (sname) { 00671 nx = 128; 00672 } else { 00673 nx = 128; 00674 } 00675 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) { 00676 if (sname) { 00677 nx = 128; 00678 } else { 00679 nx = 128; 00680 } 00681 } 00682 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) { 00683 if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00684 nx = 32; 00685 } 00686 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) { 00687 if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00688 nx = 32; 00689 } 00690 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) { 00691 if (*(unsigned char *)c3 == 'G') { 00692 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00693 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00694 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00695 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00696 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00697 ftnlen)2, (ftnlen)2) == 0) { 00698 nx = 128; 00699 } 00700 } 00701 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) { 00702 if (*(unsigned char *)c3 == 'G') { 00703 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00704 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00705 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00706 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00707 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00708 ftnlen)2, (ftnlen)2) == 0) { 00709 nx = 128; 00710 } 00711 } 00712 } 00713 ret_val = nx; 00714 return ret_val; 00715 00716 L400: 00717 00718 /* ISPEC = 4: number of shifts (used by xHSEQR) */ 00719 00720 ret_val = 6; 00721 return ret_val; 00722 00723 L500: 00724 00725 /* ISPEC = 5: minimum column dimension (not used) */ 00726 00727 ret_val = 2; 00728 return ret_val; 00729 00730 L600: 00731 00732 /* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) */ 00733 00734 ret_val = (integer) ((real) min(*n1,*n2) * 1.6f); 00735 return ret_val; 00736 00737 L700: 00738 00739 /* ISPEC = 7: number of processors (not used) */ 00740 00741 ret_val = 1; 00742 return ret_val; 00743 00744 L800: 00745 00746 /* ISPEC = 8: crossover point for multishift (used by xHSEQR) */ 00747 00748 ret_val = 50; 00749 return ret_val; 00750 00751 L900: 00752 00753 /* 00754 ISPEC = 9: maximum size of the subproblems at the bottom of the 00755 computation tree in the divide-and-conquer algorithm 00756 (used by xGELSD and xGESDD) 00757 */ 00758 00759 ret_val = 25; 00760 return ret_val; 00761 00762 L1000: 00763 00764 /* 00765 ISPEC = 10: ieee NaN arithmetic can be trusted not to trap 00766 00767 ILAENV = 0 00768 */ 00769 ret_val = 1; 00770 if (ret_val == 1) { 00771 ret_val = ieeeck_(&c__0, &c_b163, &c_b164); 00772 } 00773 return ret_val; 00774 00775 L1100: 00776 00777 /* 00778 ISPEC = 11: infinity arithmetic can be trusted not to trap 00779 00780 ILAENV = 0 00781 */ 00782 ret_val = 1; 00783 if (ret_val == 1) { 00784 ret_val = ieeeck_(&c__1, &c_b163, &c_b164); 00785 } 00786 return ret_val; 00787 00788 /* End of ILAENV */ 00789 00790 } /* ilaenv_ */ 00791 00792 /* Subroutine */ int sposv_(char *uplo, integer *n, integer *nrhs, real *a, 00793 integer *lda, real *b, integer *ldb, integer *info) 00794 { 00795 /* System generated locals */ 00796 integer a_dim1, a_offset, b_dim1, b_offset, i__1; 00797 00798 /* Local variables */ 00799 extern logical lsame_(char *, char *); 00800 extern /* Subroutine */ int xerbla_(char *, integer *), spotrf_( 00801 char *, integer *, real *, integer *, integer *), spotrs_( 00802 char *, integer *, integer *, real *, integer *, real *, integer * 00803 , integer *); 00804 00805 00806 /* 00807 -- LAPACK driver routine (version 3.0) -- 00808 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00809 Courant Institute, Argonne National Lab, and Rice University 00810 March 31, 1993 00811 00812 00813 Purpose 00814 ======= 00815 00816 SPOSV computes the solution to a real system of linear equations 00817 A * X = B, 00818 where A is an N-by-N symmetric positive definite matrix and X and B 00819 are N-by-NRHS matrices. 00820 00821 The Cholesky decomposition is used to factor A as 00822 A = U**T* U, if UPLO = 'U', or 00823 A = L * L**T, if UPLO = 'L', 00824 where U is an upper triangular matrix and L is a lower triangular 00825 matrix. The factored form of A is then used to solve the system of 00826 equations A * X = B. 00827 00828 Arguments 00829 ========= 00830 00831 UPLO (input) CHARACTER*1 00832 = 'U': Upper triangle of A is stored; 00833 = 'L': Lower triangle of A is stored. 00834 00835 N (input) INTEGER 00836 The number of linear equations, i.e., the order of the 00837 matrix A. N >= 0. 00838 00839 NRHS (input) INTEGER 00840 The number of right hand sides, i.e., the number of columns 00841 of the matrix B. NRHS >= 0. 00842 00843 A (input/output) REAL array, dimension (LDA,N) 00844 On entry, the symmetric matrix A. If UPLO = 'U', the leading 00845 N-by-N upper triangular part of A contains the upper 00846 triangular part of the matrix A, and the strictly lower 00847 triangular part of A is not referenced. If UPLO = 'L', the 00848 leading N-by-N lower triangular part of A contains the lower 00849 triangular part of the matrix A, and the strictly upper 00850 triangular part of A is not referenced. 00851 00852 On exit, if INFO = 0, the factor U or L from the Cholesky 00853 factorization A = U**T*U or A = L*L**T. 00854 00855 LDA (input) INTEGER 00856 The leading dimension of the array A. LDA >= max(1,N). 00857 00858 B (input/output) REAL array, dimension (LDB,NRHS) 00859 On entry, the N-by-NRHS right hand side matrix B. 00860 On exit, if INFO = 0, the N-by-NRHS solution matrix X. 00861 00862 LDB (input) INTEGER 00863 The leading dimension of the array B. LDB >= max(1,N). 00864 00865 INFO (output) INTEGER 00866 = 0: successful exit 00867 < 0: if INFO = -i, the i-th argument had an illegal value 00868 > 0: if INFO = i, the leading minor of order i of A is not 00869 positive definite, so the factorization could not be 00870 completed, and the solution has not been computed. 00871 00872 ===================================================================== 00873 00874 00875 Test the input parameters. 00876 */ 00877 00878 /* Parameter adjustments */ 00879 a_dim1 = *lda; 00880 a_offset = 1 + a_dim1; 00881 a -= a_offset; 00882 b_dim1 = *ldb; 00883 b_offset = 1 + b_dim1; 00884 b -= b_offset; 00885 00886 /* Function Body */ 00887 *info = 0; 00888 if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) { 00889 *info = -1; 00890 } else if (*n < 0) { 00891 *info = -2; 00892 } else if (*nrhs < 0) { 00893 *info = -3; 00894 } else if (*lda < max(1,*n)) { 00895 *info = -5; 00896 } else if (*ldb < max(1,*n)) { 00897 *info = -7; 00898 } 00899 if (*info != 0) { 00900 i__1 = -(*info); 00901 xerbla_("SPOSV ", &i__1); 00902 return 0; 00903 } 00904 00905 /* Compute the Cholesky factorization A = U'*U or A = L*L'. */ 00906 00907 spotrf_(uplo, n, &a[a_offset], lda, info); 00908 if (*info == 0) { 00909 00910 /* Solve the system A*X = B, overwriting B with X. */ 00911 00912 spotrs_(uplo, n, nrhs, &a[a_offset], lda, &b[b_offset], ldb, info); 00913 00914 } 00915 return 0; 00916 00917 /* End of SPOSV */ 00918 00919 } /* sposv_ */ 00920 00921 /* Subroutine */ int spotf2_(char *uplo, integer *n, real *a, integer *lda, 00922 integer *info) 00923 { 00924 /* System generated locals */ 00925 integer a_dim1, a_offset, i__1, i__2, i__3; 00926 real r__1; 00927 00928 /* Builtin functions */ 00929 double sqrt(doublereal); 00930 00931 /* Local variables */ 00932 static integer j; 00933 static real ajj; 00934 extern doublereal sdot_(integer *, real *, integer *, real *, integer *); 00935 extern logical lsame_(char *, char *); 00936 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 00937 sgemv_(char *, integer *, integer *, real *, real *, integer *, 00938 real *, integer *, real *, real *, integer *); 00939 static logical upper; 00940 extern /* Subroutine */ int xerbla_(char *, integer *); 00941 00942 00943 /* 00944 -- LAPACK routine (version 3.0) -- 00945 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00946 Courant Institute, Argonne National Lab, and Rice University 00947 February 29, 1992 00948 00949 00950 Purpose 00951 ======= 00952 00953 SPOTF2 computes the Cholesky factorization of a real symmetric 00954 positive definite matrix A. 00955 00956 The factorization has the form 00957 A = U' * U , if UPLO = 'U', or 00958 A = L * L', if UPLO = 'L', 00959 where U is an upper triangular matrix and L is lower triangular. 00960 00961 This is the unblocked version of the algorithm, calling Level 2 BLAS. 00962 00963 Arguments 00964 ========= 00965 00966 UPLO (input) CHARACTER*1 00967 Specifies whether the upper or lower triangular part of the 00968 symmetric matrix A is stored. 00969 = 'U': Upper triangular 00970 = 'L': Lower triangular 00971 00972 N (input) INTEGER 00973 The order of the matrix A. N >= 0. 00974 00975 A (input/output) REAL array, dimension (LDA,N) 00976 On entry, the symmetric matrix A. If UPLO = 'U', the leading 00977 n by n upper triangular part of A contains the upper 00978 triangular part of the matrix A, and the strictly lower 00979 triangular part of A is not referenced. If UPLO = 'L', the 00980 leading n by n lower triangular part of A contains the lower 00981 triangular part of the matrix A, and the strictly upper 00982 triangular part of A is not referenced. 00983 00984 On exit, if INFO = 0, the factor U or L from the Cholesky 00985 factorization A = U'*U or A = L*L'. 00986 00987 LDA (input) INTEGER 00988 The leading dimension of the array A. LDA >= max(1,N). 00989 00990 INFO (output) INTEGER 00991 = 0: successful exit 00992 < 0: if INFO = -k, the k-th argument had an illegal value 00993 > 0: if INFO = k, the leading minor of order k is not 00994 positive definite, and the factorization could not be 00995 completed. 00996 00997 ===================================================================== 00998 00999 01000 Test the input parameters. 01001 */ 01002 01003 /* Parameter adjustments */ 01004 a_dim1 = *lda; 01005 a_offset = 1 + a_dim1; 01006 a -= a_offset; 01007 01008 /* Function Body */ 01009 *info = 0; 01010 upper = lsame_(uplo, "U"); 01011 if (! upper && ! lsame_(uplo, "L")) { 01012 *info = -1; 01013 } else if (*n < 0) { 01014 *info = -2; 01015 } else if (*lda < max(1,*n)) { 01016 *info = -4; 01017 } 01018 if (*info != 0) { 01019 i__1 = -(*info); 01020 xerbla_("SPOTF2", &i__1); 01021 return 0; 01022 } 01023 01024 /* Quick return if possible */ 01025 01026 if (*n == 0) { 01027 return 0; 01028 } 01029 01030 if (upper) { 01031 01032 /* Compute the Cholesky factorization A = U'*U. */ 01033 01034 i__1 = *n; 01035 for (j = 1; j <= i__1; ++j) { 01036 01037 /* Compute U(J,J) and test for non-positive-definiteness. */ 01038 01039 i__2 = j - 1; 01040 ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j * a_dim1 + 1], &c__1, 01041 &a[j * a_dim1 + 1], &c__1); 01042 if (ajj <= 0.f) { 01043 a[j + j * a_dim1] = ajj; 01044 goto L30; 01045 } 01046 ajj = sqrt(ajj); 01047 a[j + j * a_dim1] = ajj; 01048 01049 /* Compute elements J+1:N of row J. */ 01050 01051 if (j < *n) { 01052 i__2 = j - 1; 01053 i__3 = *n - j; 01054 sgemv_("Transpose", &i__2, &i__3, &c_b181, &a[(j + 1) * 01055 a_dim1 + 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b164, 01056 &a[j + (j + 1) * a_dim1], lda); 01057 i__2 = *n - j; 01058 r__1 = 1.f / ajj; 01059 sscal_(&i__2, &r__1, &a[j + (j + 1) * a_dim1], lda); 01060 } 01061 /* L10: */ 01062 } 01063 } else { 01064 01065 /* Compute the Cholesky factorization A = L*L'. */ 01066 01067 i__1 = *n; 01068 for (j = 1; j <= i__1; ++j) { 01069 01070 /* Compute L(J,J) and test for non-positive-definiteness. */ 01071 01072 i__2 = j - 1; 01073 ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j + a_dim1], lda, &a[j 01074 + a_dim1], lda); 01075 if (ajj <= 0.f) { 01076 a[j + j * a_dim1] = ajj; 01077 goto L30; 01078 } 01079 ajj = sqrt(ajj); 01080 a[j + j * a_dim1] = ajj; 01081 01082 /* Compute elements J+1:N of column J. */ 01083 01084 if (j < *n) { 01085 i__2 = *n - j; 01086 i__3 = j - 1; 01087 sgemv_("No transpose", &i__2, &i__3, &c_b181, &a[j + 1 + 01088 a_dim1], lda, &a[j + a_dim1], lda, &c_b164, &a[j + 1 01089 + j * a_dim1], &c__1); 01090 i__2 = *n - j; 01091 r__1 = 1.f / ajj; 01092 sscal_(&i__2, &r__1, &a[j + 1 + j * a_dim1], &c__1); 01093 } 01094 /* L20: */ 01095 } 01096 } 01097 goto L40; 01098 01099 L30: 01100 *info = j; 01101 01102 L40: 01103 return 0; 01104 01105 /* End of SPOTF2 */ 01106 01107 } /* spotf2_ */ 01108 01109 /* Subroutine */ int spotrf_(char *uplo, integer *n, real *a, integer *lda, 01110 integer *info) 01111 { 01112 /* System generated locals */ 01113 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 01114 01115 /* Local variables */ 01116 static integer j, jb, nb; 01117 extern logical lsame_(char *, char *); 01118 extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 01119 integer *, real *, real *, integer *, real *, integer *, real *, 01120 real *, integer *); 01121 static logical upper; 01122 extern /* Subroutine */ int strsm_(char *, char *, char *, char *, 01123 integer *, integer *, real *, real *, integer *, real *, integer * 01124 ), ssyrk_(char *, char *, integer 01125 *, integer *, real *, real *, integer *, real *, real *, integer * 01126 ), spotf2_(char *, integer *, real *, integer *, 01127 integer *), xerbla_(char *, integer *); 01128 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 01129 integer *, integer *, ftnlen, ftnlen); 01130 01131 01132 /* 01133 -- LAPACK routine (version 3.0) -- 01134 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 01135 Courant Institute, Argonne National Lab, and Rice University 01136 March 31, 1993 01137 01138 01139 Purpose 01140 ======= 01141 01142 SPOTRF computes the Cholesky factorization of a real symmetric 01143 positive definite matrix A. 01144 01145 The factorization has the form 01146 A = U**T * U, if UPLO = 'U', or 01147 A = L * L**T, if UPLO = 'L', 01148 where U is an upper triangular matrix and L is lower triangular. 01149 01150 This is the block version of the algorithm, calling Level 3 BLAS. 01151 01152 Arguments 01153 ========= 01154 01155 UPLO (input) CHARACTER*1 01156 = 'U': Upper triangle of A is stored; 01157 = 'L': Lower triangle of A is stored. 01158 01159 N (input) INTEGER 01160 The order of the matrix A. N >= 0. 01161 01162 A (input/output) REAL array, dimension (LDA,N) 01163 On entry, the symmetric matrix A. If UPLO = 'U', the leading 01164 N-by-N upper triangular part of A contains the upper 01165 triangular part of the matrix A, and the strictly lower 01166 triangular part of A is not referenced. If UPLO = 'L', the 01167 leading N-by-N lower triangular part of A contains the lower 01168 triangular part of the matrix A, and the strictly upper 01169 triangular part of A is not referenced. 01170 01171 On exit, if INFO = 0, the factor U or L from the Cholesky 01172 factorization A = U**T*U or A = L*L**T. 01173 01174 LDA (input) INTEGER 01175 The leading dimension of the array A. LDA >= max(1,N). 01176 01177 INFO (output) INTEGER 01178 = 0: successful exit 01179 < 0: if INFO = -i, the i-th argument had an illegal value 01180 > 0: if INFO = i, the leading minor of order i is not 01181 positive definite, and the factorization could not be 01182 completed. 01183 01184 ===================================================================== 01185 01186 01187 Test the input parameters. 01188 */ 01189 01190 /* Parameter adjustments */ 01191 a_dim1 = *lda; 01192 a_offset = 1 + a_dim1; 01193 a -= a_offset; 01194 01195 /* Function Body */ 01196 *info = 0; 01197 upper = lsame_(uplo, "U"); 01198 if (! upper && ! lsame_(uplo, "L")) { 01199 *info = -1; 01200 } else if (*n < 0) { 01201 *info = -2; 01202 } else if (*lda < max(1,*n)) { 01203 *info = -4; 01204 } 01205 if (*info != 0) { 01206 i__1 = -(*info); 01207 xerbla_("SPOTRF", &i__1); 01208 return 0; 01209 } 01210 01211 /* Quick return if possible */ 01212 01213 if (*n == 0) { 01214 return 0; 01215 } 01216 01217 /* Determine the block size for this environment. */ 01218 01219 nb = ilaenv_(&c__1, "SPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( 01220 ftnlen)1); 01221 if (nb <= 1 || nb >= *n) { 01222 01223 /* Use unblocked code. */ 01224 01225 spotf2_(uplo, n, &a[a_offset], lda, info); 01226 } else { 01227 01228 /* Use blocked code. */ 01229 01230 if (upper) { 01231 01232 /* Compute the Cholesky factorization A = U'*U. */ 01233 01234 i__1 = *n; 01235 i__2 = nb; 01236 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { 01237 01238 /* 01239 Update and factorize the current diagonal block and test 01240 for non-positive-definiteness. 01241 01242 Computing MIN 01243 */ 01244 i__3 = nb, i__4 = *n - j + 1; 01245 jb = min(i__3,i__4); 01246 i__3 = j - 1; 01247 ssyrk_("Upper", "Transpose", &jb, &i__3, &c_b181, &a[j * 01248 a_dim1 + 1], lda, &c_b164, &a[j + j * a_dim1], lda); 01249 spotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info); 01250 if (*info != 0) { 01251 goto L30; 01252 } 01253 if (j + jb <= *n) { 01254 01255 /* Compute the current block row. */ 01256 01257 i__3 = *n - j - jb + 1; 01258 i__4 = j - 1; 01259 sgemm_("Transpose", "No transpose", &jb, &i__3, &i__4, & 01260 c_b181, &a[j * a_dim1 + 1], lda, &a[(j + jb) * 01261 a_dim1 + 1], lda, &c_b164, &a[j + (j + jb) * 01262 a_dim1], lda); 01263 i__3 = *n - j - jb + 1; 01264 strsm_("Left", "Upper", "Transpose", "Non-unit", &jb, & 01265 i__3, &c_b164, &a[j + j * a_dim1], lda, &a[j + (j 01266 + jb) * a_dim1], lda); 01267 } 01268 /* L10: */ 01269 } 01270 01271 } else { 01272 01273 /* Compute the Cholesky factorization A = L*L'. */ 01274 01275 i__2 = *n; 01276 i__1 = nb; 01277 for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { 01278 01279 /* 01280 Update and factorize the current diagonal block and test 01281 for non-positive-definiteness. 01282 01283 Computing MIN 01284 */ 01285 i__3 = nb, i__4 = *n - j + 1; 01286 jb = min(i__3,i__4); 01287 i__3 = j - 1; 01288 ssyrk_("Lower", "No transpose", &jb, &i__3, &c_b181, &a[j + 01289 a_dim1], lda, &c_b164, &a[j + j * a_dim1], lda); 01290 spotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info); 01291 if (*info != 0) { 01292 goto L30; 01293 } 01294 if (j + jb <= *n) { 01295 01296 /* Compute the current block column. */ 01297 01298 i__3 = *n - j - jb + 1; 01299 i__4 = j - 1; 01300 sgemm_("No transpose", "Transpose", &i__3, &jb, &i__4, & 01301 c_b181, &a[j + jb + a_dim1], lda, &a[j + a_dim1], 01302 lda, &c_b164, &a[j + jb + j * a_dim1], lda); 01303 i__3 = *n - j - jb + 1; 01304 strsm_("Right", "Lower", "Transpose", "Non-unit", &i__3, & 01305 jb, &c_b164, &a[j + j * a_dim1], lda, &a[j + jb + 01306 j * a_dim1], lda); 01307 } 01308 /* L20: */ 01309 } 01310 } 01311 } 01312 goto L40; 01313 01314 L30: 01315 *info = *info + j - 1; 01316 01317 L40: 01318 return 0; 01319 01320 /* End of SPOTRF */ 01321 01322 } /* spotrf_ */ 01323 01324 /* Subroutine */ int spotrs_(char *uplo, integer *n, integer *nrhs, real *a, 01325 integer *lda, real *b, integer *ldb, integer *info) 01326 { 01327 /* System generated locals */ 01328 integer a_dim1, a_offset, b_dim1, b_offset, i__1; 01329 01330 /* Local variables */ 01331 extern logical lsame_(char *, char *); 01332 static logical upper; 01333 extern /* Subroutine */ int strsm_(char *, char *, char *, char *, 01334 integer *, integer *, real *, real *, integer *, real *, integer * 01335 ), xerbla_(char *, integer *); 01336 01337 01338 /* 01339 -- LAPACK routine (version 3.0) -- 01340 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 01341 Courant Institute, Argonne National Lab, and Rice University 01342 March 31, 1993 01343 01344 01345 Purpose 01346 ======= 01347 01348 SPOTRS solves a system of linear equations A*X = B with a symmetric 01349 positive definite matrix A using the Cholesky factorization 01350 A = U**T*U or A = L*L**T computed by SPOTRF. 01351 01352 Arguments 01353 ========= 01354 01355 UPLO (input) CHARACTER*1 01356 = 'U': Upper triangle of A is stored; 01357 = 'L': Lower triangle of A is stored. 01358 01359 N (input) INTEGER 01360 The order of the matrix A. N >= 0. 01361 01362 NRHS (input) INTEGER 01363 The number of right hand sides, i.e., the number of columns 01364 of the matrix B. NRHS >= 0. 01365 01366 A (input) REAL array, dimension (LDA,N) 01367 The triangular factor U or L from the Cholesky factorization 01368 A = U**T*U or A = L*L**T, as computed by SPOTRF. 01369 01370 LDA (input) INTEGER 01371 The leading dimension of the array A. LDA >= max(1,N). 01372 01373 B (input/output) REAL array, dimension (LDB,NRHS) 01374 On entry, the right hand side matrix B. 01375 On exit, the solution matrix X. 01376 01377 LDB (input) INTEGER 01378 The leading dimension of the array B. LDB >= max(1,N). 01379 01380 INFO (output) INTEGER 01381 = 0: successful exit 01382 < 0: if INFO = -i, the i-th argument had an illegal value 01383 01384 ===================================================================== 01385 01386 01387 Test the input parameters. 01388 */ 01389 01390 /* Parameter adjustments */ 01391 a_dim1 = *lda; 01392 a_offset = 1 + a_dim1; 01393 a -= a_offset; 01394 b_dim1 = *ldb; 01395 b_offset = 1 + b_dim1; 01396 b -= b_offset; 01397 01398 /* Function Body */ 01399 *info = 0; 01400 upper = lsame_(uplo, "U"); 01401 if (! upper && ! lsame_(uplo, "L")) { 01402 *info = -1; 01403 } else if (*n < 0) { 01404 *info = -2; 01405 } else if (*nrhs < 0) { 01406 *info = -3; 01407 } else if (*lda < max(1,*n)) { 01408 *info = -5; 01409 } else if (*ldb < max(1,*n)) { 01410 *info = -7; 01411 } 01412 if (*info != 0) { 01413 i__1 = -(*info); 01414 xerbla_("SPOTRS", &i__1); 01415 return 0; 01416 } 01417 01418 /* Quick return if possible */ 01419 01420 if (*n == 0 || *nrhs == 0) { 01421 return 0; 01422 } 01423 01424 if (upper) { 01425 01426 /* 01427 Solve A*X = B where A = U'*U. 01428 01429 Solve U'*X = B, overwriting B with X. 01430 */ 01431 01432 strsm_("Left", "Upper", "Transpose", "Non-unit", n, nrhs, &c_b164, &a[ 01433 a_offset], lda, &b[b_offset], ldb); 01434 01435 /* Solve U*X = B, overwriting B with X. */ 01436 01437 strsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b164, 01438 &a[a_offset], lda, &b[b_offset], ldb); 01439 } else { 01440 01441 /* 01442 Solve A*X = B where A = L*L'. 01443 01444 Solve L*X = B, overwriting B with X. 01445 */ 01446 01447 strsm_("Left", "Lower", "No transpose", "Non-unit", n, nrhs, &c_b164, 01448 &a[a_offset], lda, &b[b_offset], ldb); 01449 01450 /* Solve L'*X = B, overwriting B with X. */ 01451 01452 strsm_("Left", "Lower", "Transpose", "Non-unit", n, nrhs, &c_b164, &a[ 01453 a_offset], lda, &b[b_offset], ldb); 01454 } 01455 01456 return 0; 01457 01458 /* End of SPOTRS */ 01459 01460 } /* spotrs_ */ 01461