SphinxBase  5prealpha
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