SphinxBase 0.6

src/libsphinxbase/util/slapack_lite.c

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