src/libsphinxbase/util/slapack_lite.c

   1 /*
   2 NOTE: This is generated code. Look in README.python for information on
   3       remaking this file.
   4 */
   5 #include "sphinxbase/f2c.h"
   6
   7 #ifdef HAVE_CONFIG
   8 #include "config.h"
   9 #else
  10 extern doublereal slamch_(char *);
  11 #define EPSILON slamch_("Epsilon")
  12 #define SAFEMINIMUM slamch_("Safe minimum")
  13 #define PRECISION slamch_("Precision")
  14 #define BASE slamch_("Base")
  15 #endif
  16
  17
  18 extern doublereal slapy2_(real *, real *);
  19
  20
  21
  22 /* Table of constant values */
  23
  24 static integer c__0 = 0;
  25 static real c_b163 = 0.f;
  26 static real c_b164 = 1.f;
  27 static integer c__1 = 1;
  28 static real c_b181 = -1.f;
  29 static integer c_n1 = -1;
  30
  31 integer ieeeck_(integer *ispec, real *zero, real *one)
  32 {
  33     /* System generated locals */
  34     integer ret_val;
  35
  36     /* Local variables */
  37     static real nan1, nan2, nan3, nan4, nan5, nan6, neginf, posinf, negzro,
  38             newzro;
  39
  40
  41 /*
  42     -- LAPACK auxiliary routine (version 3.0) --
  43        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  44        Courant Institute, Argonne National Lab, and Rice University
  45        June 30, 1998
  46
  47
  48     Purpose
  49     =======
  50
  51     IEEECK is called from the ILAENV to verify that Infinity and
  52     possibly NaN arithmetic is safe (i.e. will not trap).
  53
  54     Arguments
  55     =========
  56
  57     ISPEC   (input) INTEGER
  58             Specifies whether to test just for inifinity arithmetic
  59             or whether to test for infinity and NaN arithmetic.
  60             = 0: Verify infinity arithmetic only.
  61             = 1: Verify infinity and NaN arithmetic.
  62
  63     ZERO    (input) REAL
  64             Must contain the value 0.0
  65             This is passed to prevent the compiler from optimizing
  66             away this code.
  67
  68     ONE     (input) REAL
  69             Must contain the value 1.0
  70             This is passed to prevent the compiler from optimizing
  71             away this code.
  72
  73     RETURN VALUE:  INTEGER
  74             = 0:  Arithmetic failed to produce the correct answers
  75             = 1:  Arithmetic produced the correct answers
  76 */
  77
  78     ret_val = 1;
  79
  80     posinf = *one / *zero;
  81     if (posinf <= *one) {
  82         ret_val = 0;
  83         return ret_val;
  84     }
  85
  86     neginf = -(*one) / *zero;
  87     if (neginf >= *zero) {
  88         ret_val = 0;
  89         return ret_val;
  90     }
  91
  92     negzro = *one / (neginf + *one);
  93     if (negzro != *zero) {
  94         ret_val = 0;
  95         return ret_val;
  96     }
  97
  98     neginf = *one / negzro;
  99     if (neginf >= *zero) {
 100         ret_val = 0;
 101         return ret_val;
 102     }
 103
 104     newzro = negzro + *zero;
 105     if (newzro != *zero) {
 106         ret_val = 0;
 107         return ret_val;
 108     }
 109
 110     posinf = *one / newzro;
 111     if (posinf <= *one) {
 112         ret_val = 0;
 113         return ret_val;
 114     }
 115
 116     neginf *= posinf;
 117     if (neginf >= *zero) {
 118         ret_val = 0;
 119         return ret_val;
 120     }
 121
 122     posinf *= posinf;
 123     if (posinf <= *one) {
 124         ret_val = 0;
 125         return ret_val;
 126     }
 127
 128
 129 /*     Return if we were only asked to check infinity arithmetic */
 130
 131     if (*ispec == 0) {
 132         return ret_val;
 133     }
 134
 135     nan1 = posinf + neginf;
 136
 137     nan2 = posinf / neginf;
 138
 139     nan3 = posinf / posinf;
 140
 141     nan4 = posinf * *zero;
 142
 143     nan5 = neginf * negzro;
 144
 145     nan6 = nan5 * 0.f;
 146
 147     if (nan1 == nan1) {
 148         ret_val = 0;
 149         return ret_val;
 150     }
 151
 152     if (nan2 == nan2) {
 153         ret_val = 0;
 154         return ret_val;
 155     }
 156
 157     if (nan3 == nan3) {
 158         ret_val = 0;
 159         return ret_val;
 160     }
 161
 162     if (nan4 == nan4) {
 163         ret_val = 0;
 164         return ret_val;
 165     }
 166
 167     if (nan5 == nan5) {
 168         ret_val = 0;
 169         return ret_val;
 170     }
 171
 172     if (nan6 == nan6) {
 173         ret_val = 0;
 174         return ret_val;
 175     }
 176
 177     return ret_val;
 178 } /* ieeeck_ */
 179
 180 integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
 181         integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen
 182         opts_len)
 183 {
 184     /* System generated locals */
 185     integer ret_val;
 186
 187     /* Builtin functions */
 188     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
 189     integer s_cmp(char *, char *, ftnlen, ftnlen);
 190
 191     /* Local variables */
 192     static integer i__;
 193     static char c1[1], c2[2], c3[3], c4[2];
 194     static integer ic, nb, iz, nx;
 195     static logical cname, sname;
 196     static integer nbmin;
 197     extern integer ieeeck_(integer *, real *, real *);
 198     static char subnam[6];
 199
 200
 201 /*
 202     -- LAPACK auxiliary routine (version 3.0) --
 203        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 204        Courant Institute, Argonne National Lab, and Rice University
 205        June 30, 1999
 206
 207
 208     Purpose
 209     =======
 210
 211     ILAENV is called from the LAPACK routines to choose problem-dependent
 212     parameters for the local environment.  See ISPEC for a description of
 213     the parameters.
 214
 215     This version provides a set of parameters which should give good,
 216     but not optimal, performance on many of the currently available
 217     computers.  Users are encouraged to modify this subroutine to set
 218     the tuning parameters for their particular machine using the option
 219     and problem size information in the arguments.
 220
 221     This routine will not function correctly if it is converted to all
 222     lower case.  Converting it to all upper case is allowed.
 223
 224     Arguments
 225     =========
 226
 227     ISPEC   (input) INTEGER
 228             Specifies the parameter to be returned as the value of
 229             ILAENV.
 230             = 1: the optimal blocksize; if this value is 1, an unblocked
 231                  algorithm will give the best performance.
 232             = 2: the minimum block size for which the block routine
 233                  should be used; if the usable block size is less than
 234                  this value, an unblocked routine should be used.
 235             = 3: the crossover point (in a block routine, for N less
 236                  than this value, an unblocked routine should be used)
 237             = 4: the number of shifts, used in the nonsymmetric
 238                  eigenvalue routines
 239             = 5: the minimum column dimension for blocking to be used;
 240                  rectangular blocks must have dimension at least k by m,
 241                  where k is given by ILAENV(2,...) and m by ILAENV(5,...)
 242             = 6: the crossover point for the SVD (when reducing an m by n
 243                  matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
 244                  this value, a QR factorization is used first to reduce
 245                  the matrix to a triangular form.)
 246             = 7: the number of processors
 247             = 8: the crossover point for the multishift QR and QZ methods
 248                  for nonsymmetric eigenvalue problems.
 249             = 9: maximum size of the subproblems at the bottom of the
 250                  computation tree in the divide-and-conquer algorithm
 251                  (used by xGELSD and xGESDD)
 252             =10: ieee NaN arithmetic can be trusted not to trap
 253             =11: infinity arithmetic can be trusted not to trap
 254
 255     NAME    (input) CHARACTER*(*)
 256             The name of the calling subroutine, in either upper case or
 257             lower case.
 258
 259     OPTS    (input) CHARACTER*(*)
 260             The character options to the subroutine NAME, concatenated
 261             into a single character string.  For example, UPLO = 'U',
 262             TRANS = 'T', and DIAG = 'N' for a triangular routine would
 263             be specified as OPTS = 'UTN'.
 264
 265     N1      (input) INTEGER
 266     N2      (input) INTEGER
 267     N3      (input) INTEGER
 268     N4      (input) INTEGER
 269             Problem dimensions for the subroutine NAME; these may not all
 270             be required.
 271
 272    (ILAENV) (output) INTEGER
 273             >= 0: the value of the parameter specified by ISPEC
 274             < 0:  if ILAENV = -k, the k-th argument had an illegal value.
 275
 276     Further Details
 277     ===============
 278
 279     The following conventions have been used when calling ILAENV from the
 280     LAPACK routines:
 281     1)  OPTS is a concatenation of all of the character options to
 282         subroutine NAME, in the same order that they appear in the
 283         argument list for NAME, even if they are not used in determining
 284         the value of the parameter specified by ISPEC.
 285     2)  The problem dimensions N1, N2, N3, N4 are specified in the order
 286         that they appear in the argument list for NAME.  N1 is used
 287         first, N2 second, and so on, and unused problem dimensions are
 288         passed a value of -1.
 289     3)  The parameter value returned by ILAENV is checked for validity in
 290         the calling subroutine.  For example, ILAENV is used to retrieve
 291         the optimal blocksize for STRTRI as follows:
 292
 293         NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
 294         IF( NB.LE.1 ) NB = MAX( 1, N )
 295
 296     =====================================================================
 297 */
 298
 299
 300     switch (*ispec) {
 301         case 1:  goto L100;
 302         case 2:  goto L100;
 303         case 3:  goto L100;
 304         case 4:  goto L400;
 305         case 5:  goto L500;
 306         case 6:  goto L600;
 307         case 7:  goto L700;
 308         case 8:  goto L800;
 309         case 9:  goto L900;
 310         case 10:  goto L1000;
 311         case 11:  goto L1100;
 312     }
 313
 314 /*     Invalid value for ISPEC */
 315
 316     ret_val = -1;
 317     return ret_val;
 318
 319 L100:
 320
 321 /*     Convert NAME to upper case if the first character is lower case. */
 322
 323     ret_val = 1;
 324     s_copy(subnam, name__, (ftnlen)6, name_len);
 325     ic = *(unsigned char *)subnam;
 326     iz = 'Z';
 327     if (iz == 90 || iz == 122) {
 328
 329 /*        ASCII character set */
 330
 331         if (ic >= 97 && ic <= 122) {
 332             *(unsigned char *)subnam = (char) (ic - 32);
 333             for (i__ = 2; i__ <= 6; ++i__) {
 334                 ic = *(unsigned char *)&subnam[i__ - 1];
 335                 if (ic >= 97 && ic <= 122) {
 336                     *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
 337                 }
 338 /* L10: */
 339             }
 340         }
 341
 342     } else if (iz == 233 || iz == 169) {
 343
 344 /*        EBCDIC character set */
 345
 346         if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 162 &&
 347                 ic <= 169) {
 348             *(unsigned char *)subnam = (char) (ic + 64);
 349             for (i__ = 2; i__ <= 6; ++i__) {
 350                 ic = *(unsigned char *)&subnam[i__ - 1];
 351                 if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >=
 352                         162 && ic <= 169) {
 353                     *(unsigned char *)&subnam[i__ - 1] = (char) (ic + 64);
 354                 }
 355 /* L20: */
 356             }
 357         }
 358
 359     } else if (iz == 218 || iz == 250) {
 360
 361 /*        Prime machines:  ASCII+128 */
 362
 363         if (ic >= 225 && ic <= 250) {
 364             *(unsigned char *)subnam = (char) (ic - 32);
 365             for (i__ = 2; i__ <= 6; ++i__) {
 366                 ic = *(unsigned char *)&subnam[i__ - 1];
 367                 if (ic >= 225 && ic <= 250) {
 368                     *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
 369                 }
 370 /* L30: */
 371             }
 372         }
 373     }
 374
 375     *(unsigned char *)c1 = *(unsigned char *)subnam;
 376     sname = *(unsigned char *)c1 == 'S' || *(unsigned char *)c1 == 'D';
 377     cname = *(unsigned char *)c1 == 'C' || *(unsigned char *)c1 == 'Z';
 378     if (! (cname || sname)) {
 379         return ret_val;
 380     }
 381     s_copy(c2, subnam + 1, (ftnlen)2, (ftnlen)2);
 382     s_copy(c3, subnam + 3, (ftnlen)3, (ftnlen)3);
 383     s_copy(c4, c3 + 1, (ftnlen)2, (ftnlen)2);
 384
 385     switch (*ispec) {
 386         case 1:  goto L110;
 387         case 2:  goto L200;
 388         case 3:  goto L300;
 389     }
 390
 391 L110:
 392
 393 /*
 394        ISPEC = 1:  block size
 395
 396        In these examples, separate code is provided for setting NB for
 397        real and complex.  We assume that NB will take the same value in
 398        single or double precision.
 399 */
 400
 401     nb = 1;
 402
 403     if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
 404         if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
 405             if (sname) {
 406                 nb = 64;
 407             } else {
 408                 nb = 64;
 409             }
 410         } else if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3,
 411                 "RQF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)
 412                 3, (ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3)
 413                 == 0) {
 414             if (sname) {
 415                 nb = 32;
 416             } else {
 417                 nb = 32;
 418             }
 419         } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
 420             if (sname) {
 421                 nb = 32;
 422             } else {
 423                 nb = 32;
 424             }
 425         } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
 426             if (sname) {
 427                 nb = 32;
 428             } else {
 429                 nb = 32;
 430             }
 431         } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
 432             if (sname) {
 433                 nb = 64;
 434             } else {
 435                 nb = 64;
 436             }
 437         }
 438     } else if (s_cmp(c2, "PO", (ftnlen)2, (ftnlen)2) == 0) {
 439         if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
 440             if (sname) {
 441                 nb = 64;
 442             } else {
 443                 nb = 64;
 444             }
 445         }
 446     } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
 447         if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
 448             if (sname) {
 449                 nb = 64;
 450             } else {
 451                 nb = 64;
 452             }
 453         } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
 454             nb = 32;
 455         } else if (sname && s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) {
 456             nb = 64;
 457         }
 458     } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
 459         if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
 460             nb = 64;
 461         } else if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
 462             nb = 32;
 463         } else if (s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) {
 464             nb = 64;
 465         }
 466     } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
 467         if (*(unsigned char *)c3 == 'G') {
 468             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 469                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 470                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 471                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 472                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 473                     ftnlen)2, (ftnlen)2) == 0) {
 474                 nb = 32;
 475             }
 476         } else if (*(unsigned char *)c3 == 'M') {
 477             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 478                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 479                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 480                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 481                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 482                     ftnlen)2, (ftnlen)2) == 0) {
 483                 nb = 32;
 484             }
 485         }
 486     } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
 487         if (*(unsigned char *)c3 == 'G') {
 488             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 489                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 490                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 491                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 492                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 493                     ftnlen)2, (ftnlen)2) == 0) {
 494                 nb = 32;
 495             }
 496         } else if (*(unsigned char *)c3 == 'M') {
 497             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 498                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 499                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 500                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 501                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 502                     ftnlen)2, (ftnlen)2) == 0) {
 503                 nb = 32;
 504             }
 505         }
 506     } else if (s_cmp(c2, "GB", (ftnlen)2, (ftnlen)2) == 0) {
 507         if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
 508             if (sname) {
 509                 if (*n4 <= 64) {
 510                     nb = 1;
 511                 } else {
 512                     nb = 32;
 513                 }
 514             } else {
 515                 if (*n4 <= 64) {
 516                     nb = 1;
 517                 } else {
 518                     nb = 32;
 519                 }
 520             }
 521         }
 522     } else if (s_cmp(c2, "PB", (ftnlen)2, (ftnlen)2) == 0) {
 523         if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
 524             if (sname) {
 525                 if (*n2 <= 64) {
 526                     nb = 1;
 527                 } else {
 528                     nb = 32;
 529                 }
 530             } else {
 531                 if (*n2 <= 64) {
 532                     nb = 1;
 533                 } else {
 534                     nb = 32;
 535                 }
 536             }
 537         }
 538     } else if (s_cmp(c2, "TR", (ftnlen)2, (ftnlen)2) == 0) {
 539         if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
 540             if (sname) {
 541                 nb = 64;
 542             } else {
 543                 nb = 64;
 544             }
 545         }
 546     } else if (s_cmp(c2, "LA", (ftnlen)2, (ftnlen)2) == 0) {
 547         if (s_cmp(c3, "UUM", (ftnlen)3, (ftnlen)3) == 0) {
 548             if (sname) {
 549                 nb = 64;
 550             } else {
 551                 nb = 64;
 552             }
 553         }
 554     } else if (sname && s_cmp(c2, "ST", (ftnlen)2, (ftnlen)2) == 0) {
 555         if (s_cmp(c3, "EBZ", (ftnlen)3, (ftnlen)3) == 0) {
 556             nb = 1;
 557         }
 558     }
 559     ret_val = nb;
 560     return ret_val;
 561
 562 L200:
 563
 564 /*     ISPEC = 2:  minimum block size */
 565
 566     nbmin = 2;
 567     if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
 568         if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
 569                 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, (
 570                 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0)
 571                  {
 572             if (sname) {
 573                 nbmin = 2;
 574             } else {
 575                 nbmin = 2;
 576             }
 577         } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
 578             if (sname) {
 579                 nbmin = 2;
 580             } else {
 581                 nbmin = 2;
 582             }
 583         } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
 584             if (sname) {
 585                 nbmin = 2;
 586             } else {
 587                 nbmin = 2;
 588             }
 589         } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
 590             if (sname) {
 591                 nbmin = 2;
 592             } else {
 593                 nbmin = 2;
 594             }
 595         }
 596     } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
 597         if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
 598             if (sname) {
 599                 nbmin = 8;
 600             } else {
 601                 nbmin = 8;
 602             }
 603         } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
 604             nbmin = 2;
 605         }
 606     } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
 607         if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
 608             nbmin = 2;
 609         }
 610     } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
 611         if (*(unsigned char *)c3 == 'G') {
 612             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 613                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 614                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 615                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 616                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 617                     ftnlen)2, (ftnlen)2) == 0) {
 618                 nbmin = 2;
 619             }
 620         } else if (*(unsigned char *)c3 == 'M') {
 621             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 622                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 623                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 624                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 625                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 626                     ftnlen)2, (ftnlen)2) == 0) {
 627                 nbmin = 2;
 628             }
 629         }
 630     } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
 631         if (*(unsigned char *)c3 == 'G') {
 632             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 633                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 634                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 635                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 636                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 637                     ftnlen)2, (ftnlen)2) == 0) {
 638                 nbmin = 2;
 639             }
 640         } else if (*(unsigned char *)c3 == 'M') {
 641             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 642                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 643                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 644                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 645                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 646                     ftnlen)2, (ftnlen)2) == 0) {
 647                 nbmin = 2;
 648             }
 649         }
 650     }
 651     ret_val = nbmin;
 652     return ret_val;
 653
 654 L300:
 655
 656 /*     ISPEC = 3:  crossover point */
 657
 658     nx = 0;
 659     if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
 660         if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
 661                 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, (
 662                 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0)
 663                  {
 664             if (sname) {
 665                 nx = 128;
 666             } else {
 667                 nx = 128;
 668             }
 669         } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
 670             if (sname) {
 671                 nx = 128;
 672             } else {
 673                 nx = 128;
 674             }
 675         } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
 676             if (sname) {
 677                 nx = 128;
 678             } else {
 679                 nx = 128;
 680             }
 681         }
 682     } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
 683         if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
 684             nx = 32;
 685         }
 686     } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
 687         if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
 688             nx = 32;
 689         }
 690     } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
 691         if (*(unsigned char *)c3 == 'G') {
 692             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 693                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 694                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 695                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 696                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 697                     ftnlen)2, (ftnlen)2) == 0) {
 698                 nx = 128;
 699             }
 700         }
 701     } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
 702         if (*(unsigned char *)c3 == 'G') {
 703             if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
 704                     (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
 705                     ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
 706                      0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
 707                     c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
 708                     ftnlen)2, (ftnlen)2) == 0) {
 709                 nx = 128;
 710             }
 711         }
 712     }
 713     ret_val = nx;
 714     return ret_val;
 715
 716 L400:
 717
 718 /*     ISPEC = 4:  number of shifts (used by xHSEQR) */
 719
 720     ret_val = 6;
 721     return ret_val;
 722
 723 L500:
 724
 725 /*     ISPEC = 5:  minimum column dimension (not used) */
 726
 727     ret_val = 2;
 728     return ret_val;
 729
 730 L600:
 731
 732 /*     ISPEC = 6:  crossover point for SVD (used by xGELSS and xGESVD) */
 733
 734     ret_val = (integer) ((real) min(*n1,*n2) * 1.6f);
 735     return ret_val;
 736
 737 L700:
 738
 739 /*     ISPEC = 7:  number of processors (not used) */
 740
 741     ret_val = 1;
 742     return ret_val;
 743
 744 L800:
 745
 746 /*     ISPEC = 8:  crossover point for multishift (used by xHSEQR) */
 747
 748     ret_val = 50;
 749     return ret_val;
 750
 751 L900:
 752
 753 /*
 754        ISPEC = 9:  maximum size of the subproblems at the bottom of the
 755                    computation tree in the divide-and-conquer algorithm
 756                    (used by xGELSD and xGESDD)
 757 */
 758
 759     ret_val = 25;
 760     return ret_val;
 761
 762 L1000:
 763
 764 /*
 765        ISPEC = 10: ieee NaN arithmetic can be trusted not to trap
 766
 767        ILAENV = 0
 768 */
 769     ret_val = 1;
 770     if (ret_val == 1) {
 771         ret_val = ieeeck_(&c__0, &c_b163, &c_b164);
 772     }
 773     return ret_val;
 774
 775 L1100:
 776
 777 /*
 778        ISPEC = 11: infinity arithmetic can be trusted not to trap
 779
 780        ILAENV = 0
 781 */
 782     ret_val = 1;
 783     if (ret_val == 1) {
 784         ret_val = ieeeck_(&c__1, &c_b163, &c_b164);
 785     }
 786     return ret_val;
 787
 788 /*     End of ILAENV */
 789
 790 } /* ilaenv_ */
 791
 792 /* Subroutine */ int sposv_(char *uplo, integer *n, integer *nrhs, real *a,
 793         integer *lda, real *b, integer *ldb, integer *info)
 794 {
 795     /* System generated locals */
 796     integer a_dim1, a_offset, b_dim1, b_offset, i__1;
 797
 798     /* Local variables */
 799     extern logical lsame_(char *, char *);
 800     extern /* Subroutine */ int xerbla_(char *, integer *), spotrf_(
 801             char *, integer *, real *, integer *, integer *), spotrs_(
 802             char *, integer *, integer *, real *, integer *, real *, integer *
 803             , integer *);
 804
 805
 806 /*
 807     -- LAPACK driver routine (version 3.0) --
 808        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 809        Courant Institute, Argonne National Lab, and Rice University
 810        March 31, 1993
 811
 812
 813     Purpose
 814     =======
 815
 816     SPOSV computes the solution to a real system of linear equations
 817        A * X = B,
 818     where A is an N-by-N symmetric positive definite matrix and X and B
 819     are N-by-NRHS matrices.
 820
 821     The Cholesky decomposition is used to factor A as
 822        A = U**T* U,  if UPLO = 'U', or
 823        A = L * L**T,  if UPLO = 'L',
 824     where U is an upper triangular matrix and L is a lower triangular
 825     matrix.  The factored form of A is then used to solve the system of
 826     equations A * X = B.
 827
 828     Arguments
 829     =========
 830
 831     UPLO    (input) CHARACTER*1
 832             = 'U':  Upper triangle of A is stored;
 833             = 'L':  Lower triangle of A is stored.
 834
 835     N       (input) INTEGER
 836             The number of linear equations, i.e., the order of the
 837             matrix A.  N >= 0.
 838
 839     NRHS    (input) INTEGER
 840             The number of right hand sides, i.e., the number of columns
 841             of the matrix B.  NRHS >= 0.
 842
 843     A       (input/output) REAL array, dimension (LDA,N)
 844             On entry, the symmetric matrix A.  If UPLO = 'U', the leading
 845             N-by-N upper triangular part of A contains the upper
 846             triangular part of the matrix A, and the strictly lower
 847             triangular part of A is not referenced.  If UPLO = 'L', the
 848             leading N-by-N lower triangular part of A contains the lower
 849             triangular part of the matrix A, and the strictly upper
 850             triangular part of A is not referenced.
 851
 852             On exit, if INFO = 0, the factor U or L from the Cholesky
 853             factorization A = U**T*U or A = L*L**T.
 854
 855     LDA     (input) INTEGER
 856             The leading dimension of the array A.  LDA >= max(1,N).
 857
 858     B       (input/output) REAL array, dimension (LDB,NRHS)
 859             On entry, the N-by-NRHS right hand side matrix B.
 860             On exit, if INFO = 0, the N-by-NRHS solution matrix X.
 861
 862     LDB     (input) INTEGER
 863             The leading dimension of the array B.  LDB >= max(1,N).
 864
 865     INFO    (output) INTEGER
 866             = 0:  successful exit
 867             < 0:  if INFO = -i, the i-th argument had an illegal value
 868             > 0:  if INFO = i, the leading minor of order i of A is not
 869                   positive definite, so the factorization could not be
 870                   completed, and the solution has not been computed.
 871
 872     =====================================================================
 873
 874
 875        Test the input parameters.
 876 */
 877
 878     /* Parameter adjustments */
 879     a_dim1 = *lda;
 880     a_offset = 1 + a_dim1;
 881     a -= a_offset;
 882     b_dim1 = *ldb;
 883     b_offset = 1 + b_dim1;
 884     b -= b_offset;
 885
 886     /* Function Body */
 887     *info = 0;
 888     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
 889         *info = -1;
 890     } else if (*n < 0) {
 891         *info = -2;
 892     } else if (*nrhs < 0) {
 893         *info = -3;
 894     } else if (*lda < max(1,*n)) {
 895         *info = -5;
 896     } else if (*ldb < max(1,*n)) {
 897         *info = -7;
 898     }
 899     if (*info != 0) {
 900         i__1 = -(*info);
 901         xerbla_("SPOSV ", &i__1);
 902         return 0;
 903     }
 904
 905 /*     Compute the Cholesky factorization A = U'*U or A = L*L'. */
 906
 907     spotrf_(uplo, n, &a[a_offset], lda, info);
 908     if (*info == 0) {
 909
 910 /*        Solve the system A*X = B, overwriting B with X. */
 911
 912         spotrs_(uplo, n, nrhs, &a[a_offset], lda, &b[b_offset], ldb, info);
 913
 914     }
 915     return 0;
 916
 917 /*     End of SPOSV */
 918
 919 } /* sposv_ */
 920
 921 /* Subroutine */ int spotf2_(char *uplo, integer *n, real *a, integer *lda,
 922         integer *info)
 923 {
 924     /* System generated locals */
 925     integer a_dim1, a_offset, i__1, i__2, i__3;
 926     real r__1;
 927
 928     /* Builtin functions */
 929     double sqrt(doublereal);
 930
 931     /* Local variables */
 932     static integer j;
 933     static real ajj;
 934     extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
 935     extern logical lsame_(char *, char *);
 936     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),
 937             sgemv_(char *, integer *, integer *, real *, real *, integer *,
 938             real *, integer *, real *, real *, integer *);
 939     static logical upper;
 940     extern /* Subroutine */ int xerbla_(char *, integer *);
 941
 942
 943 /*
 944     -- LAPACK routine (version 3.0) --
 945        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 946        Courant Institute, Argonne National Lab, and Rice University
 947        February 29, 1992
 948
 949
 950     Purpose
 951     =======
 952
 953     SPOTF2 computes the Cholesky factorization of a real symmetric
 954     positive definite matrix A.
 955
 956     The factorization has the form
 957        A = U' * U ,  if UPLO = 'U', or
 958        A = L  * L',  if UPLO = 'L',
 959     where U is an upper triangular matrix and L is lower triangular.
 960
 961     This is the unblocked version of the algorithm, calling Level 2 BLAS.
 962
 963     Arguments
 964     =========
 965
 966     UPLO    (input) CHARACTER*1
 967             Specifies whether the upper or lower triangular part of the
 968             symmetric matrix A is stored.
 969             = 'U':  Upper triangular
 970             = 'L':  Lower triangular
 971
 972     N       (input) INTEGER
 973             The order of the matrix A.  N >= 0.
 974
 975     A       (input/output) REAL array, dimension (LDA,N)
 976             On entry, the symmetric matrix A.  If UPLO = 'U', the leading
 977             n by n upper triangular part of A contains the upper
 978             triangular part of the matrix A, and the strictly lower
 979             triangular part of A is not referenced.  If UPLO = 'L', the
 980             leading n by n lower triangular part of A contains the lower
 981             triangular part of the matrix A, and the strictly upper
 982             triangular part of A is not referenced.
 983
 984             On exit, if INFO = 0, the factor U or L from the Cholesky
 985             factorization A = U'*U  or A = L*L'.
 986
 987     LDA     (input) INTEGER
 988             The leading dimension of the array A.  LDA >= max(1,N).
 989
 990     INFO    (output) INTEGER
 991             = 0: successful exit
 992             < 0: if INFO = -k, the k-th argument had an illegal value
 993             > 0: if INFO = k, the leading minor of order k is not
 994                  positive definite, and the factorization could not be
 995                  completed.
 996
 997     =====================================================================
 998
 999
1000        Test the input parameters.
1001 */
1002
1003     /* Parameter adjustments */
1004     a_dim1 = *lda;
1005     a_offset = 1 + a_dim1;
1006     a -= a_offset;
1007
1008     /* Function Body */
1009     *info = 0;
1010     upper = lsame_(uplo, "U");
1011     if (! upper && ! lsame_(uplo, "L")) {
1012         *info = -1;
1013     } else if (*n < 0) {
1014         *info = -2;
1015     } else if (*lda < max(1,*n)) {
1016         *info = -4;
1017     }
1018     if (*info != 0) {
1019         i__1 = -(*info);
1020         xerbla_("SPOTF2", &i__1);
1021         return 0;
1022     }
1023
1024 /*     Quick return if possible */
1025
1026     if (*n == 0) {
1027         return 0;
1028     }
1029
1030     if (upper) {
1031
1032 /*        Compute the Cholesky factorization A = U'*U. */
1033
1034         i__1 = *n;
1035         for (j = 1; j <= i__1; ++j) {
1036
1037 /*           Compute U(J,J) and test for non-positive-definiteness. */
1038
1039             i__2 = j - 1;
1040             ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j * a_dim1 + 1], &c__1,
1041                     &a[j * a_dim1 + 1], &c__1);
1042             if (ajj <= 0.f) {
1043                 a[j + j * a_dim1] = ajj;
1044                 goto L30;
1045             }
1046             ajj = sqrt(ajj);
1047             a[j + j * a_dim1] = ajj;
1048
1049 /*           Compute elements J+1:N of row J. */
1050
1051             if (j < *n) {
1052                 i__2 = j - 1;
1053                 i__3 = *n - j;
1054                 sgemv_("Transpose", &i__2, &i__3, &c_b181, &a[(j + 1) *
1055                         a_dim1 + 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b164,
1056                         &a[j + (j + 1) * a_dim1], lda);
1057                 i__2 = *n - j;
1058                 r__1 = 1.f / ajj;
1059                 sscal_(&i__2, &r__1, &a[j + (j + 1) * a_dim1], lda);
1060             }
1061 /* L10: */
1062         }
1063     } else {
1064
1065 /*        Compute the Cholesky factorization A = L*L'. */
1066
1067         i__1 = *n;
1068         for (j = 1; j <= i__1; ++j) {
1069
1070 /*           Compute L(J,J) and test for non-positive-definiteness. */
1071
1072             i__2 = j - 1;
1073             ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j + a_dim1], lda, &a[j
1074                     + a_dim1], lda);
1075             if (ajj <= 0.f) {
1076                 a[j + j * a_dim1] = ajj;
1077                 goto L30;
1078             }
1079             ajj = sqrt(ajj);
1080             a[j + j * a_dim1] = ajj;
1081
1082 /*           Compute elements J+1:N of column J. */
1083
1084             if (j < *n) {
1085                 i__2 = *n - j;
1086                 i__3 = j - 1;
1087                 sgemv_("No transpose", &i__2, &i__3, &c_b181, &a[j + 1 +
1088                         a_dim1], lda, &a[j + a_dim1], lda, &c_b164, &a[j + 1
1089                         + j * a_dim1], &c__1);
1090                 i__2 = *n - j;
1091                 r__1 = 1.f / ajj;
1092                 sscal_(&i__2, &r__1, &a[j + 1 + j * a_dim1], &c__1);
1093             }
1094 /* L20: */
1095         }
1096     }
1097     goto L40;
1098
1099 L30:
1100     *info = j;
1101
1102 L40:
1103     return 0;
1104
1105 /*     End of SPOTF2 */
1106
1107 } /* spotf2_ */
1108
1109 /* Subroutine */ int spotrf_(char *uplo, integer *n, real *a, integer *lda,
1110         integer *info)
1111 {
1112     /* System generated locals */
1113     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
1114
1115     /* Local variables */
1116     static integer j, jb, nb;
1117     extern logical lsame_(char *, char *);
1118     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
1119             integer *, real *, real *, integer *, real *, integer *, real *,
1120             real *, integer *);
1121     static logical upper;
1122     extern /* Subroutine */ int strsm_(char *, char *, char *, char *,
1123             integer *, integer *, real *, real *, integer *, real *, integer *
1124             ), ssyrk_(char *, char *, integer
1125             *, integer *, real *, real *, integer *, real *, real *, integer *
1126             ), spotf2_(char *, integer *, real *, integer *,
1127             integer *), xerbla_(char *, integer *);
1128     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
1129             integer *, integer *, ftnlen, ftnlen);
1130
1131
1132 /*
1133     -- LAPACK routine (version 3.0) --
1134        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1135        Courant Institute, Argonne National Lab, and Rice University
1136        March 31, 1993
1137
1138
1139     Purpose
1140     =======
1141
1142     SPOTRF computes the Cholesky factorization of a real symmetric
1143     positive definite matrix A.
1144
1145     The factorization has the form
1146        A = U**T * U,  if UPLO = 'U', or
1147        A = L  * L**T,  if UPLO = 'L',
1148     where U is an upper triangular matrix and L is lower triangular.
1149
1150     This is the block version of the algorithm, calling Level 3 BLAS.
1151
1152     Arguments
1153     =========
1154
1155     UPLO    (input) CHARACTER*1
1156             = 'U':  Upper triangle of A is stored;
1157             = 'L':  Lower triangle of A is stored.
1158
1159     N       (input) INTEGER
1160             The order of the matrix A.  N >= 0.
1161
1162     A       (input/output) REAL array, dimension (LDA,N)
1163             On entry, the symmetric matrix A.  If UPLO = 'U', the leading
1164             N-by-N upper triangular part of A contains the upper
1165             triangular part of the matrix A, and the strictly lower
1166             triangular part of A is not referenced.  If UPLO = 'L', the
1167             leading N-by-N lower triangular part of A contains the lower
1168             triangular part of the matrix A, and the strictly upper
1169             triangular part of A is not referenced.
1170
1171             On exit, if INFO = 0, the factor U or L from the Cholesky
1172             factorization A = U**T*U or A = L*L**T.
1173
1174     LDA     (input) INTEGER
1175             The leading dimension of the array A.  LDA >= max(1,N).
1176
1177     INFO    (output) INTEGER
1178             = 0:  successful exit
1179             < 0:  if INFO = -i, the i-th argument had an illegal value
1180             > 0:  if INFO = i, the leading minor of order i is not
1181                   positive definite, and the factorization could not be
1182                   completed.
1183
1184     =====================================================================
1185
1186
1187        Test the input parameters.
1188 */
1189
1190     /* Parameter adjustments */
1191     a_dim1 = *lda;
1192     a_offset = 1 + a_dim1;
1193     a -= a_offset;
1194
1195     /* Function Body */
1196     *info = 0;
1197     upper = lsame_(uplo, "U");
1198     if (! upper && ! lsame_(uplo, "L")) {
1199         *info = -1;
1200     } else if (*n < 0) {
1201         *info = -2;
1202     } else if (*lda < max(1,*n)) {
1203         *info = -4;
1204     }
1205     if (*info != 0) {
1206         i__1 = -(*info);
1207         xerbla_("SPOTRF", &i__1);
1208         return 0;
1209     }
1210
1211 /*     Quick return if possible */
1212
1213     if (*n == 0) {
1214         return 0;
1215     }
1216
1217 /*     Determine the block size for this environment. */
1218
1219     nb = ilaenv_(&c__1, "SPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
1220             ftnlen)1);
1221     if (nb <= 1 || nb >= *n) {
1222
1223 /*        Use unblocked code. */
1224
1225         spotf2_(uplo, n, &a[a_offset], lda, info);
1226     } else {
1227
1228 /*        Use blocked code. */
1229
1230         if (upper) {
1231
1232 /*           Compute the Cholesky factorization A = U'*U. */
1233
1234             i__1 = *n;
1235             i__2 = nb;
1236             for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
1237
1238 /*
1239                 Update and factorize the current diagonal block and test
1240                 for non-positive-definiteness.
1241
1242    Computing MIN
1243 */
1244                 i__3 = nb, i__4 = *n - j + 1;
1245                 jb = min(i__3,i__4);
1246                 i__3 = j - 1;
1247                 ssyrk_("Upper", "Transpose", &jb, &i__3, &c_b181, &a[j *
1248                         a_dim1 + 1], lda, &c_b164, &a[j + j * a_dim1], lda);
1249                 spotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info);
1250                 if (*info != 0) {
1251                     goto L30;
1252                 }
1253                 if (j + jb <= *n) {
1254
1255 /*                 Compute the current block row. */
1256
1257                     i__3 = *n - j - jb + 1;
1258                     i__4 = j - 1;
1259                     sgemm_("Transpose", "No transpose", &jb, &i__3, &i__4, &
1260                             c_b181, &a[j * a_dim1 + 1], lda, &a[(j + jb) *
1261                             a_dim1 + 1], lda, &c_b164, &a[j + (j + jb) *
1262                             a_dim1], lda);
1263                     i__3 = *n - j - jb + 1;
1264                     strsm_("Left", "Upper", "Transpose", "Non-unit", &jb, &
1265                             i__3, &c_b164, &a[j + j * a_dim1], lda, &a[j + (j
1266                             + jb) * a_dim1], lda);
1267                 }
1268 /* L10: */
1269             }
1270
1271         } else {
1272
1273 /*           Compute the Cholesky factorization A = L*L'. */
1274
1275             i__2 = *n;
1276             i__1 = nb;
1277             for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
1278
1279 /*
1280                 Update and factorize the current diagonal block and test
1281                 for non-positive-definiteness.
1282
1283    Computing MIN
1284 */
1285                 i__3 = nb, i__4 = *n - j + 1;
1286                 jb = min(i__3,i__4);
1287                 i__3 = j - 1;
1288                 ssyrk_("Lower", "No transpose", &jb, &i__3, &c_b181, &a[j +
1289                         a_dim1], lda, &c_b164, &a[j + j * a_dim1], lda);
1290                 spotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info);
1291                 if (*info != 0) {
1292                     goto L30;
1293                 }
1294                 if (j + jb <= *n) {
1295
1296 /*                 Compute the current block column. */
1297
1298                     i__3 = *n - j - jb + 1;
1299                     i__4 = j - 1;
1300                     sgemm_("No transpose", "Transpose", &i__3, &jb, &i__4, &
1301                             c_b181, &a[j + jb + a_dim1], lda, &a[j + a_dim1],
1302                             lda, &c_b164, &a[j + jb + j * a_dim1], lda);
1303                     i__3 = *n - j - jb + 1;
1304                     strsm_("Right", "Lower", "Transpose", "Non-unit", &i__3, &
1305                             jb, &c_b164, &a[j + j * a_dim1], lda, &a[j + jb +
1306                             j * a_dim1], lda);
1307                 }
1308 /* L20: */
1309             }
1310         }
1311     }
1312     goto L40;
1313
1314 L30:
1315     *info = *info + j - 1;
1316
1317 L40:
1318     return 0;
1319
1320 /*     End of SPOTRF */
1321
1322 } /* spotrf_ */
1323
1324 /* Subroutine */ int spotrs_(char *uplo, integer *n, integer *nrhs, real *a,
1325         integer *lda, real *b, integer *ldb, integer *info)
1326 {
1327     /* System generated locals */
1328     integer a_dim1, a_offset, b_dim1, b_offset, i__1;
1329
1330     /* Local variables */
1331     extern logical lsame_(char *, char *);
1332     static logical upper;
1333     extern /* Subroutine */ int strsm_(char *, char *, char *, char *,
1334             integer *, integer *, real *, real *, integer *, real *, integer *
1335             ), xerbla_(char *, integer *);
1336
1337
1338 /*
1339     -- LAPACK routine (version 3.0) --
1340        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1341        Courant Institute, Argonne National Lab, and Rice University
1342        March 31, 1993
1343
1344
1345     Purpose
1346     =======
1347
1348     SPOTRS solves a system of linear equations A*X = B with a symmetric
1349     positive definite matrix A using the Cholesky factorization
1350     A = U**T*U or A = L*L**T computed by SPOTRF.
1351
1352     Arguments
1353     =========
1354
1355     UPLO    (input) CHARACTER*1
1356             = 'U':  Upper triangle of A is stored;
1357             = 'L':  Lower triangle of A is stored.
1358
1359     N       (input) INTEGER
1360             The order of the matrix A.  N >= 0.
1361
1362     NRHS    (input) INTEGER
1363             The number of right hand sides, i.e., the number of columns
1364             of the matrix B.  NRHS >= 0.
1365
1366     A       (input) REAL array, dimension (LDA,N)
1367             The triangular factor U or L from the Cholesky factorization
1368             A = U**T*U or A = L*L**T, as computed by SPOTRF.
1369
1370     LDA     (input) INTEGER
1371             The leading dimension of the array A.  LDA >= max(1,N).
1372
1373     B       (input/output) REAL array, dimension (LDB,NRHS)
1374             On entry, the right hand side matrix B.
1375             On exit, the solution matrix X.
1376
1377     LDB     (input) INTEGER
1378             The leading dimension of the array B.  LDB >= max(1,N).
1379
1380     INFO    (output) INTEGER
1381             = 0:  successful exit
1382             < 0:  if INFO = -i, the i-th argument had an illegal value
1383
1384     =====================================================================
1385
1386
1387        Test the input parameters.
1388 */
1389
1390     /* Parameter adjustments */
1391     a_dim1 = *lda;
1392     a_offset = 1 + a_dim1;
1393     a -= a_offset;
1394     b_dim1 = *ldb;
1395     b_offset = 1 + b_dim1;
1396     b -= b_offset;
1397
1398     /* Function Body */
1399     *info = 0;
1400     upper = lsame_(uplo, "U");
1401     if (! upper && ! lsame_(uplo, "L")) {
1402         *info = -1;
1403     } else if (*n < 0) {
1404         *info = -2;
1405     } else if (*nrhs < 0) {
1406         *info = -3;
1407     } else if (*lda < max(1,*n)) {
1408         *info = -5;
1409     } else if (*ldb < max(1,*n)) {
1410         *info = -7;
1411     }
1412     if (*info != 0) {
1413         i__1 = -(*info);
1414         xerbla_("SPOTRS", &i__1);
1415         return 0;
1416     }
1417
1418 /*     Quick return if possible */
1419
1420     if (*n == 0 || *nrhs == 0) {
1421         return 0;
1422     }
1423
1424     if (upper) {
1425
1426 /*
1427           Solve A*X = B where A = U'*U.
1428
1429           Solve U'*X = B, overwriting B with X.
1430 */
1431
1432         strsm_("Left", "Upper", "Transpose", "Non-unit", n, nrhs, &c_b164, &a[
1433                 a_offset], lda, &b[b_offset], ldb);
1434
1435 /*        Solve U*X = B, overwriting B with X. */
1436
1437         strsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b164,
1438                 &a[a_offset], lda, &b[b_offset], ldb);
1439     } else {
1440
1441 /*
1442           Solve A*X = B where A = L*L'.
1443
1444           Solve L*X = B, overwriting B with X.
1445 */
1446
1447         strsm_("Left", "Lower", "No transpose", "Non-unit", n, nrhs, &c_b164,
1448                 &a[a_offset], lda, &b[b_offset], ldb);
1449
1450 /*        Solve L'*X = B, overwriting B with X. */
1451
1452         strsm_("Left", "Lower", "Transpose", "Non-unit", n, nrhs, &c_b164, &a[
1453                 a_offset], lda, &b[b_offset], ldb);
1454     }
1455
1456     return 0;
1457
1458 /*     End of SPOTRS */
1459
1460 } /* spotrs_ */
1461