SphinxBase 5prealpha
slapack_lite.c
1/*
2NOTE: 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
10extern 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
18extern doublereal slapy2_(real *, real *);
19
20
21
22/* Table of constant values */
23
24static integer c__0 = 0;
25static real c_b163 = 0.f;
26static real c_b164 = 1.f;
27static integer c__1 = 1;
28static real c_b181 = -1.f;
29static integer c_n1 = -1;
30
31integer 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
180integer 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
319L100:
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
391L110:
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
562L200:
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
654L300:
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
716L400:
717
718/* ISPEC = 4: number of shifts (used by xHSEQR) */
719
720 ret_val = 6;
721 return ret_val;
722
723L500:
724
725/* ISPEC = 5: minimum column dimension (not used) */
726
727 ret_val = 2;
728 return ret_val;
729
730L600:
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
737L700:
738
739/* ISPEC = 7: number of processors (not used) */
740
741 ret_val = 1;
742 return ret_val;
743
744L800:
745
746/* ISPEC = 8: crossover point for multishift (used by xHSEQR) */
747
748 ret_val = 50;
749 return ret_val;
750
751L900:
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
762L1000:
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
775L1100:
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
1099L30:
1100 *info = j;
1101
1102L40:
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
1314L30:
1315 *info = *info + j - 1;
1316
1317L40:
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