cloudy trunk
Loading...
Searching...
No Matches
thirdparty_lapack.cpp
Go to the documentation of this file.
1/* These are wrappers for lapack linear algebra routines.
2 * There are two versions of the lapack routines - a fortran
3 * version that I converted to C with forc to use if nothing else is available
4 * (included in the Cloudy distribution),
5 * and an option to link into an external lapack library that may be optimized
6 * for your machine. By default the tralated C routines will be used.
7 * To use your machine's lapack library instead, define the macro
8 * LAPACK and link into your library. This is usually done with a command line
9 * switch "-DLAPACK" on the compiler command, and the linker option "-llapack"
10 */
11#include "cddefines.h"
12#include "thirdparty.h"
13/*lint -e725 expected pos indentation */
14/*lint -e801 goto is deprecated */
15
16#ifdef LAPACK
17/*********************The functions for FORTRAN version of the LAPACK calls *******************/
18/* dgetrf, dgetrs: lapack FORTRAN general full matrix solution using LU decomposition */
19
20extern "C" void dgetrf_(int32 *M, int32 *N, double *A, int32 *LDA, int32 *IPIV, int32 *INFO);
21extern "C" void dgetrs_(char *TRANS, int32 *N, int32 *NRHS, double *A, int32 *LDA, int32 *iPiv, double *B,
22 int32 *LDB, int32 *INFO, int32 translen);
23extern "C" void dgtsv_(int32 *n, int32 *nrhs, double *dl, double *d, double *du, double *b, int32 *ldb, int32 *info);
24
25#else
26
27/*********************The functions for C version of the LAPACK calls *******************/
28/*
29 * these are the routines that are, part of lapack, Some had their names slightly
30 * changed so as to not conflict with the special lapack that exists on our exemplar.
31 */
32
33/* DGETRF lapack, perform LU decomposition */
34STATIC void DGETRF(int32,int32,double*,int32,int32[],int32*);
35
36/* DGETRS lapack, solve linear system */
37STATIC void DGETRS(int32 TRANS,int32 N,int32 NRHS,double *A,int32 LDA,int32 IPIV[],double *B,int32 LDB,int32 *INFO);
38
39/* DGTSV lapack, solve tridiagonal system */
40/*static int32 DGTSV(int32 *n,int32 *nrhs,double *dl,double *d__,double *du,double *b,int32 *ldb,int32 *info);*/
41
42#endif
43
44
45/**********************************************************************************************/
46/* returns zero if successful termination */
47void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
48{
49 if( *info == 0 )
50 {
51 int32 M_loc, N_loc, lda_loc;
52
53 ASSERT( M < INT32_MAX && N < INT32_MAX && lda < INT32_MAX );
54
55 M_loc = (int32)M;
56 N_loc = (int32)N;
57 lda_loc = (int32)lda;
58
59# ifdef LAPACK
60 /* Calling the special version in library */
61 dgetrf_(&M_loc, &N_loc, A , &lda_loc, ipiv, info);
62# else
63 /* Calling the old slower one, included with cloudy */
64 DGETRF(M_loc, N_loc, A, lda_loc, ipiv, info);
65# endif
66 }
67}
68
69void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv,
70 double *B, long ldb, int32 *info)
71{
72 if( *info == 0 )
73 {
74 int32 N_loc, nrhs_loc, lda_loc, ldb_loc;
75
76 ASSERT( N < INT32_MAX && nrhs < INT32_MAX && lda < INT32_MAX && ldb < INT32_MAX );
77
78 N_loc = (int32)N;
79 nrhs_loc = (int32)nrhs;
80 lda_loc = (int32)lda;
81 ldb_loc = (int32)ldb;
82
83# ifdef LAPACK
84 /* Calling the special version in library */
85 dgetrs_(&trans, &N_loc, &nrhs_loc, A, &lda_loc, ipiv, B, &ldb_loc, info, sizeof(char));
86# else
87 /* Calling the old slower one, included with cloudy */
88 DGETRS(trans, N_loc, nrhs_loc, A, lda_loc, ipiv, B, ldb_loc, info);
89# endif
90 }
91}
92
93#if 0
94void dgtsv_wrapper(long N, long nrhs, double *dl, double *d__, double *du, double *b, long ldb, int32 *info)
95{
96 printf("Inside dgtsv\n");
98 if( *info == 0 )
99 {
100 int32 N_loc, nrhs_loc, ldb_loc;
101
102 ASSERT( N < INT32_MAX && nrhs < INT32_MAX && ldb < INT32_MAX );
103
104 N_loc = (int32)N;
105 nrhs_loc = (int32)nrhs;
106 ldb_loc = (int32)ldb;
107
108# ifdef LAPACK
109 /* Calling the special version in library */
110 dgtsv_(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
111# else
112 /* Calling the old slower one, included with cloudy */
113 /* DGTSV always returns zero, so it is safe to ignore the return value */
114 (void)DGTSV(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
115# endif
116 }
117}
118#endif
119
120
121#ifndef LAPACK
122
123#define ONE 1.0e0
124#define ZERO 0.0e0
125
126#ifdef AA
127# undef AA
128#endif
129#ifdef BB
130# undef BB
131#endif
132#ifdef CC
133# undef CC
134#endif
135
136#define AA(I_,J_) (*(A+(I_)*(LDA)+(J_)))
137#define BB(I_,J_) (*(B+(I_)*(LDB)+(J_)))
138#define CC(I_,J_) (*(C+(I_)*(LDC)+(J_)))
139
140/*
141 * these are the routines that are, part of lapack, Some had their names slightly
142 * changed so as to not conflict with the special lapack that exits on our exemplar.
143 */
144
145/* dgtsv, dgetrf, dgetrs: lapack general tridiagonal solution */
146/*int32 DGTSV(int32 *n, int32 *nrhs, double *dl,
147 double *d__, double *du, double *b, int32 *ldb, int32 *info);*/
148
149/* DGETRF lapack service routine */
150/*void DGETRF(int32,int32,double*,int32,int32[],int32*);*/
151
152/*DGETRS lapack matrix inversion routine */
153/*void DGETRS(int32 TRANS,
154 int32 N,
155 int32 NRHS,
156 double *A,
157 int32 LDA,
158 int32 IPIV[],
159 double *B,
160 int32 LDB,
161 int32 *INFO);
162*/
163/* DGEMM matrix inversion helper routine*/
164STATIC void DGEMM(int32 TRANSA,
165 int32 TRANSB,
166 int32 M,
167 int32 N,
168 int32 K,
169 double ALPHA,
170 double *A,
171 int32 LDA,
172 double *B,
173 int32 LDB,
174 double BETA,
175 double *C,
176 int32 LDC);
177
178/*LSAME LAPACK auxiliary routine */
179STATIC int32 LSAME(int32 CA,
180 int32 CB);
181
182/*IDAMAX lapack service routine */
183STATIC int32 IDAMAX(int32 n,
184 double dx[],
185 int32 incx);
186
187/*DTRSM lapack service routine */
188STATIC void DTRSM(int32 SIDE,
189 int32 UPLO,
190 int32 TRANSA,
191 int32 DIAG,
192 int32 M,
193 int32 N,
194 double ALPHA,
195 double *A,
196 int32 LDA,
197 double *B,
198 int32 LDB);
199
200/* ILAENV lapack helper routine */
201STATIC int32 ILAENV(int32 ISPEC,
202 const char *NAME,
203 /*char *OPTS, */
204 int32 N1,
205 int32 N2,
206 /*int32 N3, */
207 int32 N4);
208
209/*DSWAP lapack routine */
210STATIC void DSWAP(int32 n,
211 double dx[],
212 int32 incx,
213 double dy[],
214 int32 incy);
215
216/*DSCAL lapack routine */
217STATIC void DSCAL(int32 n,
218 double da,
219 double dx[],
220 int32 incx);
221
222/*DLASWP -- LAPACK auxiliary routine (version 2.0) --*/
223STATIC void DLASWP(int32 N,
224 double *A,
225 int32 LDA,
226 int32 K1,
227 int32 K2,
228 int32 IPIV[],
229 int32 INCX);
230
231/*DGETF2 lapack service routine */
232STATIC void DGETF2(int32 M,
233 int32 N,
234 double *A,
235 int32 LDA,
236 int32 IPIV[],
237 int32 *INFO);
238
239/*DGER service routine for matrix inversion */
240STATIC void DGER(int32 M,
241 int32 N,
242 double ALPHA,
243 double X[],
244 int32 INCX,
245 double Y[],
246 int32 INCY,
247 double *A,
248 int32 LDA);
249
250/*XERBLA -- LAPACK auxiliary routine (version 2.0) -- */
251STATIC void XERBLA(const char *SRNAME,
252 int32 INFO)
253{
254
255 DEBUG_ENTRY( "XERBLA()" );
256
257 /* -- LAPACK auxiliary routine (version 2.0) --
258 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
259 * Courant Institute, Argonne National Lab, and Rice University
260 * September 30, 1994
261 *
262 * .. Scalar Arguments .. */
263 /* ..
264 *
265 * Purpose
266 * =======
267 *
268 * XERBLA is an error handler for the LAPACK routines.
269 * It is called by an LAPACK routine if an input parameter has an
270 * invalid value. A message is printed and execution stops.
271 *
272 * Installers may consider modifying the STOP statement in order to
273 * call system-specific exception-handling facilities.
274 *
275 * Arguments
276 * =========
277 *
278 * SRNAME (input) CHARACTER*6
279 * The name of the routine which called XERBLA.
280 *
281 * INFO (input) INTEGER
282 * The position of the invalid parameter in the parameter list
283 * of the calling routine.
284 *
285 * =====================================================================
286 *
287 * .. Executable Statements ..
288 * */
289 fprintf( ioQQQ, " ** On entry to %6.6s parameter number %2ld had an illegal value\n",
290 SRNAME, (long)INFO );
291
293}
294
295
297 /* number of rows of the matrix */
298 int32 M,
299 /* number of columns of the matrix
300 * M=N for square matrix */
301 int32 N,
302 /* double precision matrix */
303 double *A,
304 /* LDA is right dimension of matrix */
305 int32 LDA,
306 /* following must dimensions the smaller of M or N */
307 int32 IPIV[],
308 /* following is zero for successful exit */
309 int32 *INFO)
310{
311
312 char chL1,
313 chL2,
314 chL3,
315 chL4;
316 int32 I,
317 IINFO,
318 I_,
319 J,
320 JB,
321 J_,
322 NB,
323
324 limit,
325 limit2;
326 /*double _d_l;*/
327
328 DEBUG_ENTRY( "DGETRF()" );
329
330 /* -- LAPACK routine (version 2.0) --
331 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
332 * Courant Institute, Argonne National Lab, and Rice University
333 * March 31, 1993
334 *
335 * Purpose
336 * =======
337 *
338 * DGETRF computes an LU factorization of a general M-by-N matrix A
339 * using partial pivoting with row interchanges.
340 *
341 * The factorization has the form
342 * A = P * L * U
343 * where P is a permutation matrix, L is lower triangular with unit
344 * diagonal elements (lower trapezoidal if m > n), and U is upper
345 * triangular (upper trapezoidal if m < n).
346 *
347 * This is the right-looking Level 3 BLAS version of the algorithm.
348 *
349 * Arguments
350 * =========
351 *
352 * M (input) INTEGER
353 * The number of rows of the matrix A. M >= 0.
354 *
355 * N (input) INTEGER
356 * The number of columns of the matrix A. N >= 0.
357 *
358 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
359 * On entry, the M-by-N matrix to be factored.
360 * On exit, the factors L and U from the factorization
361 * A = P*L*U; the unit diagonal elements of L are not stored.
362 *
363 * LDA (input) INTEGER
364 * The leading dimension of the array A. LDA >= MAX(1,M).
365 *
366 * IPIV (output) INTEGER array, dimension (MIN(M,N))
367 * The pivot indices; for 1 <= i <= MIN(M,N), row i of the
368 * matrix was interchanged with row IPIV(i).
369 *
370 * INFO (output) INTEGER
371 * = 0: successful exit
372 * < 0: if INFO = -i, the i-th argument had an illegal value
373 * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
374 * has been completed, but the factor U is exactly
375 * singular, and division by zero will occur if it is used
376 * to solve a system of equations.
377 *
378 * =====================================================================
379 *
380 * .. Parameters .. */
381 /* ..
382 * .. Local Scalars .. */
383 /* ..
384 * .. External Subroutines .. */
385 /* ..
386 * .. External Functions .. */
387 /* ..
388 * .. Intrinsic Functions .. */
389 /* ..
390 * .. Executable Statements ..
391 *
392 * Test the input parameters.
393 * */
394 *INFO = 0;
395 if( M < 0 )
396 {
397 *INFO = -1;
398 }
399 else if( N < 0 )
400 {
401 *INFO = -2;
402 }
403 else if( LDA < MAX2(1,M) )
404 {
405 *INFO = -4;
406 }
407 if( *INFO != 0 )
408 {
409 XERBLA("DGETRF",-*INFO);
410 /* XERBLA does not return */
411 }
412
413 /* Quick return if possible
414 * */
415 if( M == 0 || N == 0 )
416 {
417 return;
418 }
419
420 /* Determine the block size for this environment.
421 * */
422 /* >>chng 01 oct 22, remove two parameters since not used */
423 NB = ILAENV(1,"DGETRF",/*" ",*/M,N,/*-1,*/-1);
424 if( NB <= 1 || NB >= MIN2(M,N) )
425 {
426 /* Use unblocked code.
427 * */
428 DGETF2(M,N,A,LDA,IPIV,INFO);
429 }
430 else
431 {
432
433 /* Use blocked code.
434 * */
435 limit = MIN2(M,N);
436 /*for( J=1, _do0=DOCNT(J,limit,_do1 = NB); _do0 > 0; J += _do1, _do0-- )*/
437 /*do J = 1, limit , NB */
438 for( J=1; J<=limit; J += NB )
439 {
440 J_ = J - 1;
441 JB = MIN2(MIN2(M,N)-J+1,NB);
442
443 /* Factor diagonal and subdiagonal blocks and test for exact
444 * singularity.
445 * */
446 DGETF2(M-J+1,JB,&AA(J_,J_),LDA,&IPIV[J_],&IINFO);
447
448 /* Adjust INFO and the pivot indices.
449 * */
450 if( *INFO == 0 && IINFO > 0 )
451 *INFO = IINFO + J - 1;
452 limit2 = MIN2(M,J+JB-1);
453 for( I=J; I <= limit2; I++ )
454 {
455 I_ = I - 1;
456 IPIV[I_] += J - 1;
457 }
458
459 /* Apply interchanges to columns 1:J-1.
460 * */
461 DLASWP(J-1,A,LDA,J,J+JB-1,IPIV,1);
462
463 if( J + JB <= N )
464 {
465
466 /* Apply interchanges to columns J+JB:N.
467 * */
468 DLASWP(N-J-JB+1,&AA(J_+JB,0),LDA,J,J+JB-1,IPIV,1);
469
470 /* Compute block row of U.
471 * */
472 chL1 = 'L';
473 chL2 = 'L';
474 chL3 = 'N';
475 chL4 = 'U';
476 /* CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, */
477 DTRSM(chL1,chL2,chL3,chL4,JB,N-J-JB+1,ONE,&AA(J_,J_),
478 LDA,&AA(J_+JB,J_),LDA);
479 if( J + JB <= M )
480 {
481
482 /* Update trailing submatrix.
483 * */
484 chL1 = 'N';
485 chL2 = 'N';
486 /* CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1, */
487 DGEMM(chL1,chL2,M-J-JB+1,N-J-JB+1,JB,-ONE,&AA(J_,J_+JB),
488 LDA,&AA(J_+JB,J_),LDA,ONE,&AA(J_+JB,J_+JB),LDA);
489 }
490 }
491 }
492 }
493 return;
494
495 /* End of DGETRF
496 * */
497#undef A
498}
499
500/*DGETRS lapack matrix inversion routine */
501/*****************************************************************
502 *****************************************************************
503 *
504 * matrix inversion routine
505 *
506 * solves Ax=B A is an nXn matrix, C and B are nX1 matrices
507 * dim A(n,n), B(n,1) C overwrites B.
508 * integer ipiv(n)
509 * integer info see below:
510 * INFO (output) INTEGER
511 * = 0: successful exit
512 * < 0: if INFO = -i, the i-th argument had an illegal value
513 * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
514 * has been completed, but the factor U is exactly
515 * singular, and division by zero will occur if it is used
516 * to solve a system of equations.
517 *
518 *
519 *
520 *must call the routines in the following order:
521 * call dgetrf(n,n,A,n,ipiv,info)
522 * call dgetrs('N',n,1,A,n,ipiv,B,n,info)
523 *
524 *
525 *************************************************************************** */
526
528 /* 1 ch var saying what to do */
529 int32 TRANS,
530 /* order of the matrix */
531 int32 N,
532 /* number of right hand sides */
533 int32 NRHS,
534 /* double [N][LDA] */
535 double *A,
536 /* second dim of array */
537 int32 LDA,
538 /* helper vector, dimensioned N*/
539 int32 IPIV[],
540 /* on input the ri=hs vector, on output, the result */
541 double *B,
542 /* dimension of B, 1 if one vector */
543 int32 LDB,
544 /* = 0 if ok */
545 int32 *INFO)
546{
547/*#define A(I_,J_) (*(A+(I_)*(LDA)+(J_)))*/
548/*#define B(I_,J_) (*(B+(I_)*(LDB)+(J_)))*/
549 int32 NOTRAN;
550 char chL1,
551 chL2,
552 chL3,
553 chL4;
554
555 DEBUG_ENTRY( "DGETRS()" );
556
557 /* -- LAPACK routine (version 2.0) --
558 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
559 * Courant Institute, Argonne National Lab, and Rice University
560 * March 31, 1993
561 *
562 *
563 * Purpose
564 * =======
565 *
566 * DGETRS solves a system of linear equations
567 * A * X = B or A' * X = B
568 * with a general N-by-N matrix A using the LU factorization computed
569 * by DGETRF.
570 *
571 * Arguments
572 * =========
573 *
574 * TRANS (input) CHARACTER*1
575 * Specifies the form of the system of equations:
576 * = 'N': A * X = B (No transpose)
577 * = 'T': A'* X = B (Transpose)
578 * = 'C': A'* X = B (Conjugate transpose = Transpose)
579 *
580 * N (input) INTEGER
581 * The order of the matrix A. N >= 0.
582 *
583 * NRHS (input) INTEGER
584 * The number of right hand sides, i.e., the number of columns
585 * of the matrix B. NRHS >= 0.
586 *
587 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
588 * The factors L and U from the factorization A = P*L*U
589 * as computed by DGETRF.
590 *
591 * LDA (input) INTEGER
592 * The leading dimension of the array A. LDA >= MAX(1,N).
593 *
594 * IPIV (input) INTEGER array, dimension (N)
595 * The pivot indices from DGETRF; for 1<=i<=N, row i of the
596 * matrix was interchanged with row IPIV(i).
597 *
598 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
599 * On entry, the right hand side matrix B.
600 * On exit, the solution matrix X.
601 *
602 * LDB (input) INTEGER
603 * The leading dimension of the array B. LDB >= MAX(1,N).
604 *
605 * INFO (output) INTEGER
606 * = 0: successful exit
607 * < 0: if INFO = -i, the i-th argument had an illegal value
608 *
609 * =====================================================================
610 *
611 * .. Parameters .. */
612 /* ..
613 * .. Local Scalars .. */
614 /* ..
615 * .. External Functions .. */
616 /* ..
617 * .. External Subroutines .. */
618 /* ..
619 * .. Intrinsic Functions .. */
620 /* ..
621 * .. Executable Statements ..
622 *
623 * Test the input parameters.
624 * */
625 *INFO = 0;
626 NOTRAN = LSAME(TRANS,'N');
627 if( (!NOTRAN && !LSAME(TRANS,'T')) && !LSAME(TRANS,'C') )
628 {
629 *INFO = -1;
630 }
631 else if( N < 0 )
632 {
633 *INFO = -2;
634 }
635 else if( NRHS < 0 )
636 {
637 *INFO = -3;
638 }
639 else if( LDA < MAX2(1,N) )
640 {
641 *INFO = -5;
642 }
643 else if( LDB < MAX2(1,N) )
644 {
645 *INFO = -8;
646 }
647 if( *INFO != 0 )
648 {
649 XERBLA("DGETRS",-*INFO);
650 /* XERBLA does not return */
651 }
652
653 /* Quick return if possible
654 * */
655 if( N == 0 || NRHS == 0 )
656 {
657 return;
658 }
659
660 if( NOTRAN )
661 {
662
663 /* Solve A * X = B.
664 *
665 * Apply row interchanges to the right hand sides.
666 * */
667 DLASWP(NRHS,B,LDB,1,N,IPIV,1);
668
669 /* Solve L*X = B, overwriting B with X.
670 * */
671 chL1 = 'L';
672 chL2 = 'L';
673 chL3 = 'N';
674 chL4 = 'U';
675 /* CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, */
676 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
677
678 /* Solve U*X = B, overwriting B with X.
679 * */
680 chL1 = 'L';
681 chL2 = 'U';
682 chL3 = 'N';
683 chL4 = 'N';
684 /* CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, */
685 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
686 }
687 else
688 {
689
690 /* Solve A' * X = B.
691 *
692 * Solve U'*X = B, overwriting B with X.
693 * */
694 chL1 = 'L';
695 chL2 = 'U';
696 chL3 = 'T';
697 chL4 = 'N';
698 /* CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, */
699 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
700
701 /* Solve L'*X = B, overwriting B with X.
702 * */
703 chL1 = 'L';
704 chL2 = 'L';
705 chL3 = 'T';
706 chL4 = 'U';
707 /* CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE, */
708 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
709
710 /* Apply row interchanges to the solution vectors.
711 * */
712 DLASWP(NRHS,B,LDB,1,N,IPIV,-1);
713 }
714
715 return;
716
717 /* End of DGETRS
718 * */
719#undef B
720#undef A
721}
722
723/*LSAME LAPACK auxiliary routine */
724
725STATIC int32 LSAME(int32 CA,
726 int32 CB)
727{
728 /*
729 *
730 * -- LAPACK auxiliary routine (version 2.0) --
731 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
732 * Courant Institute, Argonne National Lab, and Rice University
733 * September 30, 1994
734 *
735 * .. Scalar Arguments ..
736 */
737 int32 LSAME_v;
738 int32 INTA,
739 INTB,
740 ZCODE;
741
742 DEBUG_ENTRY( "LSAME()" );
743 /* ..
744 *
745 * Purpose
746 * =======
747 *
748 * LSAME returns .true. if CA is the same letter as CB regardless of
749 * case.
750 *
751 * Arguments
752 * =========
753 *
754 * CA (input) CHARACTER*1
755 * CB (input) CHARACTER*1
756 * CA and CB specify the single characters to be compared.
757 *
758 * =====================================================================
759 *
760 * .. Intrinsic Functions .. */
761 /* ..
762 * .. Local Scalars .. */
763 /* ..
764 * .. Executable Statements ..
765 *
766 * Test if the characters are equal
767 * */
768 LSAME_v = CA == CB;
769 if( LSAME_v )
770 {
771 return LSAME_v;
772 }
773
774 /* Now test for equivalence if both characters are alphabetic.
775 * */
776 ZCODE = 'Z';
777
778 /* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
779 * machines, on which ICHAR returns a value with bit 8 set.
780 * ICHAR('A') on Prime machines returns 193 which is the same as
781 * ICHAR('A') on an EBCDIC machine.
782 * */
783 INTA = (CA);
784 INTB = (CB);
785
786 if( ZCODE == 90 || ZCODE == 122 )
787 {
788
789 /* ASCII is assumed - ZCODE is the ASCII code of either lower or
790 * upper case 'Z'.
791 * */
792 if( INTA >= 97 && INTA <= 122 )
793 INTA -= 32;
794 if( INTB >= 97 && INTB <= 122 )
795 INTB -= 32;
796
797 }
798 else if( ZCODE == 233 || ZCODE == 169 )
799 {
800
801 /* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
802 * upper case 'Z'.
803 * */
804 if( ((INTA >= 129 && INTA <= 137) || (INTA >= 145 && INTA <=
805 153)) || (INTA >= 162 && INTA <= 169) )
806 INTA += 64;
807 if( ((INTB >= 129 && INTB <= 137) || (INTB >= 145 && INTB <=
808 153)) || (INTB >= 162 && INTB <= 169) )
809 INTB += 64;
810
811 }
812 else if( ZCODE == 218 || ZCODE == 250 )
813 {
814
815 /* ASCII is assumed, on Prime machines - ZCODE is the ASCII code
816 * plus 128 of either lower or upper case 'Z'.
817 * */
818 if( INTA >= 225 && INTA <= 250 )
819 INTA -= 32;
820 if( INTB >= 225 && INTB <= 250 )
821 INTB -= 32;
822 }
823 LSAME_v = INTA == INTB;
824
825 /* RETURN
826 *
827 * End of LSAME
828 * */
829 return LSAME_v;
830}
831
832/*IDAMAX lapack service routine */
833
834STATIC int32 IDAMAX(int32 n,
835 double dx[],
836 int32 incx)
837{
838 /*
839 * finds the index of element having max. absolute value.
840 * jack dongarra, lapack, 3/11/78.
841 * modified 3/93 to return if incx .le. 0.
842 * modified 12/3/93, array(1) declarations changed to array(*)
843 *
844 */
845 int32 IDAMAX_v,
846 i,
847 ix;
848 double dmax;
849
850 DEBUG_ENTRY( "IDAMAX()" );
851
852 IDAMAX_v = 0;
853
854 if( n < 1 || incx <= 0 )
855 {
856 return IDAMAX_v;
857 }
858
859 IDAMAX_v = 1;
860
861 if( n == 1 )
862 {
863 return IDAMAX_v;
864 }
865
866 if( incx == 1 )
867 goto L_20;
868
869 /* code for increment not equal to 1
870 * */
871 ix = 1;
872 dmax = fabs(dx[0]);
873 ix += incx;
874 for( i=2; i <= n; i++ )
875 {
876 /* if(ABS(dx(ix)).le.dmax) go to 5 */
877 if( fabs(dx[ix-1]) > dmax )
878 {
879 IDAMAX_v = i;
880 dmax = fabs(dx[ix-1]);
881 }
882 ix += incx;
883 }
884 return IDAMAX_v;
885
886 /* code for increment equal to 1
887 * */
888L_20:
889 dmax = fabs(dx[0]);
890 for( i=1; i < n; i++ )
891 {
892 /* if(ABS(dx(i)).le.dmax) go to 30 */
893 if( fabs(dx[i]) > dmax )
894 {
895 IDAMAX_v = i+1;
896 dmax = fabs(dx[i]);
897 }
898 }
899 return IDAMAX_v;
900}
901
902/*DTRSM lapack service routine */
903STATIC void DTRSM(int32 SIDE,
904 int32 UPLO,
905 int32 TRANSA,
906 int32 DIAG,
907 int32 M,
908 int32 N,
909 double ALPHA,
910 double *A,
911 int32 LDA,
912 double *B,
913 int32 LDB)
914{
915 int32 LSIDE,
916 NOUNIT,
917 UPPER;
918 int32 I,
919 INFO,
920 I_,
921 J,
922 J_,
923 K,
924 K_,
925 NROWA;
926 double TEMP;
927
928 DEBUG_ENTRY( "DTRSM()" );
929 /* .. Scalar Arguments .. */
930 /* .. Array Arguments .. */
931 /* ..
932 *
933 * Purpose
934 * =======
935 *
936 * DTRSM solves one of the matrix equations
937 *
938 * op( A )*X = alpha*B, or X*op( A ) = alpha*B,
939 *
940 * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
941 * non-unit, upper or lower triangular matrix and op( A ) is one of
942 *
943 * op( A ) = A or op( A ) = A'.
944 *
945 * The matrix X is overwritten on B.
946 *
947 * Parameters
948 * ==========
949 *
950 * SIDE - CHARACTER*1.
951 * On entry, SIDE specifies whether op( A ) appears on the left
952 * or right of X as follows:
953 *
954 * SIDE = 'L' or 'l' op( A )*X = alpha*B.
955 *
956 * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
957 *
958 * Unchanged on exit.
959 *
960 * UPLO - CHARACTER*1.
961 * On entry, UPLO specifies whether the matrix A is an upper or
962 * lower triangular matrix as follows:
963 *
964 * UPLO = 'U' or 'u' A is an upper triangular matrix.
965 *
966 * UPLO = 'L' or 'l' A is a lower triangular matrix.
967 *
968 * Unchanged on exit.
969 *
970 * TRANSA - CHARACTER*1.
971 * On entry, TRANSA specifies the form of op( A ) to be used in
972 * the matrix multiplication as follows:
973 *
974 * TRANSA = 'N' or 'n' op( A ) = A.
975 *
976 * TRANSA = 'T' or 't' op( A ) = A'.
977 *
978 * TRANSA = 'C' or 'c' op( A ) = A'.
979 *
980 * Unchanged on exit.
981 *
982 * DIAG - CHARACTER*1.
983 * On entry, DIAG specifies whether or not A is unit triangular
984 * as follows:
985 *
986 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
987 *
988 * DIAG = 'N' or 'n' A is not assumed to be unit
989 * triangular.
990 *
991 * Unchanged on exit.
992 *
993 * M - INTEGER.
994 * On entry, M specifies the number of rows of B. M must be at
995 * least zero.
996 * Unchanged on exit.
997 *
998 * N - INTEGER.
999 * On entry, N specifies the number of columns of B. N must be
1000 * at least zero.
1001 * Unchanged on exit.
1002 *
1003 * ALPHA - DOUBLE PRECISION.
1004 * On entry, ALPHA specifies the scalar alpha. When alpha is
1005 * zero then A is not referenced and B need not be set before
1006 * entry.
1007 * Unchanged on exit.
1008 *
1009 * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
1010 * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
1011 * Before entry with UPLO = 'U' or 'u', the leading k by k
1012 * upper triangular part of the array A must contain the upper
1013 * triangular matrix and the strictly lower triangular part of
1014 * A is not referenced.
1015 * Before entry with UPLO = 'L' or 'l', the leading k by k
1016 * lower triangular part of the array A must contain the lower
1017 * triangular matrix and the strictly upper triangular part of
1018 * A is not referenced.
1019 * Note that when DIAG = 'U' or 'u', the diagonal elements of
1020 * A are not referenced either, but are assumed to be unity.
1021 * Unchanged on exit.
1022 *
1023 * LDA - INTEGER.
1024 * On entry, LDA specifies the first dimension of A as declared
1025 * in the calling (sub) program. When SIDE = 'L' or 'l' then
1026 * LDA must be at least MAX( 1, m ), when SIDE = 'R' or 'r'
1027 * then LDA must be at least MAX( 1, n ).
1028 * Unchanged on exit.
1029 *
1030 * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
1031 * Before entry, the leading m by n part of the array B must
1032 * contain the right-hand side matrix B, and on exit is
1033 * overwritten by the solution matrix X.
1034 *
1035 * LDB - INTEGER.
1036 * On entry, LDB specifies the first dimension of B as declared
1037 * in the calling (sub) program. LDB must be at least
1038 * MAX( 1, m ).
1039 * Unchanged on exit.
1040 *
1041 *
1042 * Level 3 Blas routine.
1043 *
1044 *
1045 * -- Written on 8-February-1989.
1046 * Jack Dongarra, Argonne National Laboratory.
1047 * Iain Duff, AERE Harwell.
1048 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
1049 * Sven Hammarling, Numerical Algorithms Group Ltd.
1050 *
1051 *
1052 * .. External Functions .. */
1053 /* .. External Subroutines .. */
1054 /* .. Intrinsic Functions .. */
1055 /* .. Local Scalars .. */
1056 /* .. Parameters .. */
1057 /* ..
1058 * .. Executable Statements ..
1059 *
1060 * Test the input parameters.
1061 * */
1062 LSIDE = LSAME(SIDE,'L');
1063 if( LSIDE )
1064 {
1065 NROWA = M;
1066 }
1067 else
1068 {
1069 NROWA = N;
1070 }
1071 NOUNIT = LSAME(DIAG,'N');
1072 UPPER = LSAME(UPLO,'U');
1073
1074 INFO = 0;
1075 if( (!LSIDE) && (!LSAME(SIDE,'R')) )
1076 {
1077 INFO = 1;
1078 }
1079 else if( (!UPPER) && (!LSAME(UPLO,'L')) )
1080 {
1081 INFO = 2;
1082 }
1083 else if( ((!LSAME(TRANSA,'N')) && (!LSAME(TRANSA,'T'))) && (!LSAME(TRANSA,
1084 'C')) )
1085 {
1086 INFO = 3;
1087 }
1088 else if( (!LSAME(DIAG,'U')) && (!LSAME(DIAG,'N')) )
1089 {
1090 INFO = 4;
1091 }
1092 else if( M < 0 )
1093 {
1094 INFO = 5;
1095 }
1096 else if( N < 0 )
1097 {
1098 INFO = 6;
1099 }
1100 else if( LDA < MAX2(1,NROWA) )
1101 {
1102 INFO = 9;
1103 }
1104 else if( LDB < MAX2(1,M) )
1105 {
1106 INFO = 11;
1107 }
1108 if( INFO != 0 )
1109 {
1110 XERBLA("DTRSM ",INFO);
1111 /* XERBLA does not return */
1112 }
1113
1114 /* Quick return if possible.
1115 * */
1116 if( N == 0 )
1117 {
1118 return;}
1119
1120 /* And when alpha.eq.zero.
1121 * */
1122 if( ALPHA == ZERO )
1123 {
1124 for( J=1; J <= N; J++ )
1125 {
1126 J_ = J - 1;
1127 for( I=1; I <= M; I++ )
1128 {
1129 I_ = I - 1;
1130 BB(J_,I_) = ZERO;
1131 }
1132 }
1133 return;
1134 }
1135
1136 /* Start the operations.
1137 * */
1138 if( LSIDE )
1139 {
1140 if( LSAME(TRANSA,'N') )
1141 {
1142
1143 /* Form B := alpha*inv( A )*B.
1144 * */
1145 if( UPPER )
1146 {
1147 for( J=1; J <= N; J++ )
1148 {
1149 J_ = J - 1;
1150 if( ALPHA != ONE )
1151 {
1152 for( I=1; I <= M; I++ )
1153 {
1154 I_ = I - 1;
1155 BB(J_,I_) *= ALPHA;
1156 }
1157 }
1158 for( K=M; K >= 1; K-- )
1159 {
1160 K_ = K - 1;
1161 if( BB(J_,K_) != ZERO )
1162 {
1163 if( NOUNIT )
1164 BB(J_,K_) /= AA(K_,K_);
1165 for( I=1; I <= (K - 1); I++ )
1166 {
1167 I_ = I - 1;
1168 BB(J_,I_) += -BB(J_,K_)*AA(K_,I_);
1169 }
1170 }
1171 }
1172 }
1173 }
1174 else
1175 {
1176 for( J=1; J <= N; J++ )
1177 {
1178 J_ = J - 1;
1179 if( ALPHA != ONE )
1180 {
1181 for( I=1; I <= M; I++ )
1182 {
1183 I_ = I - 1;
1184 BB(J_,I_) *= ALPHA;
1185 }
1186 }
1187 for( K=1; K <= M; K++ )
1188 {
1189 K_ = K - 1;
1190 if( BB(J_,K_) != ZERO )
1191 {
1192 if( NOUNIT )
1193 BB(J_,K_) /= AA(K_,K_);
1194 for( I=K + 1; I <= M; I++ )
1195 {
1196 I_ = I - 1;
1197 BB(J_,I_) += -BB(J_,K_)*AA(K_,I_);
1198 }
1199 }
1200 }
1201 }
1202 }
1203 }
1204 else
1205 {
1206
1207 /* Form B := alpha*inv( A' )*B.
1208 * */
1209 if( UPPER )
1210 {
1211 for( J=1; J <= N; J++ )
1212 {
1213 J_ = J - 1;
1214 for( I=1; I <= M; I++ )
1215 {
1216 I_ = I - 1;
1217 TEMP = ALPHA*BB(J_,I_);
1218 for( K=1; K <= (I - 1); K++ )
1219 {
1220 K_ = K - 1;
1221 TEMP += -AA(I_,K_)*BB(J_,K_);
1222 }
1223 if( NOUNIT )
1224 TEMP /= AA(I_,I_);
1225 BB(J_,I_) = TEMP;
1226 }
1227 }
1228 }
1229 else
1230 {
1231 for( J=1; J <= N; J++ )
1232 {
1233 J_ = J - 1;
1234 for( I=M; I >= 1; I-- )
1235 {
1236 I_ = I - 1;
1237 TEMP = ALPHA*BB(J_,I_);
1238 for( K=I + 1; K <= M; K++ )
1239 {
1240 K_ = K - 1;
1241 TEMP += -AA(I_,K_)*BB(J_,K_);
1242 }
1243 if( NOUNIT )
1244 TEMP /= AA(I_,I_);
1245 BB(J_,I_) = TEMP;
1246 }
1247 }
1248 }
1249 }
1250 }
1251 else
1252 {
1253 if( LSAME(TRANSA,'N') )
1254 {
1255
1256 /* Form B := alpha*B*inv( A ).
1257 * */
1258 if( UPPER )
1259 {
1260 for( J=1; J <= N; J++ )
1261 {
1262 J_ = J - 1;
1263 if( ALPHA != ONE )
1264 {
1265 for( I=1; I <= M; I++ )
1266 {
1267 I_ = I - 1;
1268 BB(J_,I_) *= ALPHA;
1269 }
1270 }
1271 for( K=1; K <= (J - 1); K++ )
1272 {
1273 K_ = K - 1;
1274 if( AA(J_,K_) != ZERO )
1275 {
1276 for( I=1; I <= M; I++ )
1277 {
1278 I_ = I - 1;
1279 BB(J_,I_) += -AA(J_,K_)*BB(K_,I_);
1280 }
1281 }
1282 }
1283 if( NOUNIT )
1284 {
1285 TEMP = ONE/AA(J_,J_);
1286 for( I=1; I <= M; I++ )
1287 {
1288 I_ = I - 1;
1289 BB(J_,I_) *= TEMP;
1290 }
1291 }
1292 }
1293 }
1294 else
1295 {
1296 for( J=N; J >= 1; J-- )
1297 {
1298 J_ = J - 1;
1299 if( ALPHA != ONE )
1300 {
1301 for( I=1; I <= M; I++ )
1302 {
1303 I_ = I - 1;
1304 BB(J_,I_) *= ALPHA;
1305 }
1306 }
1307 for( K=J + 1; K <= N; K++ )
1308 {
1309 K_ = K - 1;
1310 if( AA(J_,K_) != ZERO )
1311 {
1312 for( I=1; I <= M; I++ )
1313 {
1314 I_ = I - 1;
1315 BB(J_,I_) += -AA(J_,K_)*BB(K_,I_);
1316 }
1317 }
1318 }
1319 if( NOUNIT )
1320 {
1321 TEMP = ONE/AA(J_,J_);
1322 for( I=1; I <= M; I++ )
1323 {
1324 I_ = I - 1;
1325 BB(J_,I_) *= TEMP;
1326 }
1327 }
1328 }
1329 }
1330 }
1331 else
1332 {
1333
1334 /* Form B := alpha*B*inv( A' ).
1335 * */
1336 if( UPPER )
1337 {
1338 for( K=N; K >= 1; K-- )
1339 {
1340 K_ = K - 1;
1341 if( NOUNIT )
1342 {
1343 TEMP = ONE/AA(K_,K_);
1344 for( I=1; I <= M; I++ )
1345 {
1346 I_ = I - 1;
1347 BB(K_,I_) *= TEMP;
1348 }
1349 }
1350 for( J=1; J <= (K - 1); J++ )
1351 {
1352 J_ = J - 1;
1353 if( AA(K_,J_) != ZERO )
1354 {
1355 TEMP = AA(K_,J_);
1356 for( I=1; I <= M; I++ )
1357 {
1358 I_ = I - 1;
1359 BB(J_,I_) += -TEMP*BB(K_,I_);
1360 }
1361 }
1362 }
1363 if( ALPHA != ONE )
1364 {
1365 for( I=1; I <= M; I++ )
1366 {
1367 I_ = I - 1;
1368 BB(K_,I_) *= ALPHA;
1369 }
1370 }
1371 }
1372 }
1373 else
1374 {
1375 for( K=1; K <= N; K++ )
1376 {
1377 K_ = K - 1;
1378 if( NOUNIT )
1379 {
1380 TEMP = ONE/AA(K_,K_);
1381 for( I=1; I <= M; I++ )
1382 {
1383 I_ = I - 1;
1384 BB(K_,I_) *= TEMP;
1385 }
1386 }
1387 for( J=K + 1; J <= N; J++ )
1388 {
1389 J_ = J - 1;
1390 if( AA(K_,J_) != ZERO )
1391 {
1392 TEMP = AA(K_,J_);
1393 for( I=1; I <= M; I++ )
1394 {
1395 I_ = I - 1;
1396 BB(J_,I_) += -TEMP*BB(K_,I_);
1397 }
1398 }
1399 }
1400 if( ALPHA != ONE )
1401 {
1402 for( I=1; I <= M; I++ )
1403 {
1404 I_ = I - 1;
1405 BB(K_,I_) *= ALPHA;
1406 }
1407 }
1408 }
1409 }
1410 }
1411 }
1412
1413 return;
1414
1415 /* End of DTRSM .
1416 * */
1417#undef B
1418#undef A
1419}
1420
1421/* ILAENV lapack helper routine */
1422
1423/* >>chng 01 oct 22, remove two parameters since not used */
1424STATIC int32 ILAENV(int32 ISPEC,
1425 const char *NAME,
1426 /*char *OPTS, */
1427 int32 N1,
1428 int32 N2,
1429 /*int32 N3, */
1430 int32 N4)
1431{
1432 /*
1433 *
1434 * -- LAPACK auxiliary routine (version 2.0) --
1435 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1436 * Courant Institute, Argonne National Lab, and Rice University
1437 * September 30, 1994
1438 *
1439 * .. Scalar Arguments ..
1440 */
1441 char C2[3],
1442 C3[4],
1443 C4[3],
1444 SUBNAM[7];
1445 int32 CNAME,
1446 SNAME;
1447 char C1;
1448 int32 I,
1449 IC,
1450 ILAENV_v,
1451 IZ,
1452 NB,
1453 NBMIN,
1454 NX;
1455
1456 DEBUG_ENTRY( "ILAENV()" );
1457 /* ..
1458 *
1459 * Purpose
1460 * =======
1461 *
1462 * ILAENV is called from the LAPACK routines to choose problem-dependent
1463 * parameters for the local environment. See ISPEC for a description of
1464 * the parameters.
1465 *
1466 * This version provides a set of parameters which should give good,
1467 * but not optimal, performance on many of the currently available
1468 * computers. Users are encouraged to modify this subroutine to set
1469 * the tuning parameters for their particular machine using the option
1470 * and problem size information in the arguments.
1471 *
1472 * This routine will not function correctly if it is converted to all
1473 * lower case. Converting it to all upper case is allowed.
1474 *
1475 * Arguments
1476 * =========
1477 *
1478 * ISPEC (input) INTEGER
1479 * Specifies the parameter to be returned as the value of
1480 * ILAENV.
1481 * = 1: the optimal blocksize; if this value is 1, an unblocked
1482 * algorithm will give the best performance.
1483 * = 2: the minimum block size for which the block routine
1484 * should be used; if the usable block size is less than
1485 * this value, an unblocked routine should be used.
1486 * = 3: the crossover point (in a block routine, for N less
1487 * than this value, an unblocked routine should be used)
1488 * = 4: the number of shifts, used in the nonsymmetric
1489 * eigenvalue routines
1490 * = 5: the minimum column dimension for blocking to be used;
1491 * rectangular blocks must have dimension at least k by m,
1492 * where k is given by ILAENV(2,...) and m by ILAENV(5,...)
1493 * = 6: the crossover point for the SVD (when reducing an m by n
1494 * matrix to bidiagonal form, if MAX(m,n)/MIN(m,n) exceeds
1495 * this value, a QR factorization is used first to reduce
1496 * the matrix to a triangular form.)
1497 * = 7: the number of processors
1498 * = 8: the crossover point for the multishift QR and QZ methods
1499 * for nonsymmetric eigenvalue problems.
1500 *
1501 * NAME (input) CHARACTER*(*)
1502 * The name of the calling subroutine, in either upper case or
1503 * lower case.
1504 *
1505 * OPTS (input) CHARACTER*(*)
1506 * The character options to the subroutine NAME, concatenated
1507 * into a single character string. For example, UPLO = 'U',
1508 * TRANS = 'T', and DIAG = 'N' for a triangular routine would
1509 * be specified as OPTS = 'UTN'.
1510 *
1511 * N1 (input) INTEGER
1512 * N2 (input) INTEGER
1513 * N3 (input) INTEGER
1514 * N4 (input) INTEGER
1515 * Problem dimensions for the subroutine NAME; these may not all
1516 * be required.
1517 *
1518 * (ILAENV) (output) INTEGER
1519 * >= 0: the value of the parameter specified by ISPEC
1520 * < 0: if ILAENV = -k, the k-th argument had an illegal value.
1521 *
1522 * Further Details
1523 * ===============
1524 *
1525 * The following conventions have been used when calling ILAENV from the
1526 * LAPACK routines:
1527 * 1) OPTS is a concatenation of all of the character options to
1528 * subroutine NAME, in the same order that they appear in the
1529 * argument list for NAME, even if they are not used in determining
1530 * the value of the parameter specified by ISPEC.
1531 * 2) The problem dimensions N1, N2, N3, N4 are specified in the order
1532 * that they appear in the argument list for NAME. N1 is used
1533 * first, N2 second, and so on, and unused problem dimensions are
1534 * passed a value of -1.
1535 * 3) The parameter value returned by ILAENV is checked for validity in
1536 * the calling subroutine. For example, ILAENV is used to retrieve
1537 * the optimal blocksize for STRTRI as follows:
1538 *
1539 * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
1540 * IF( NB.LE.1 ) NB = MAX( 1, N )
1541 *
1542 * =====================================================================
1543 *
1544 * .. Local Scalars .. */
1545 /* ..
1546 * .. Intrinsic Functions .. */
1547 /* ..
1548 * .. Executable Statements ..
1549 * */
1550 switch( ISPEC )
1551 {
1552 case 1:
1553 {
1554 goto L_100;
1555 }
1556 case 2:
1557 {
1558 goto L_100;
1559 }
1560 case 3:
1561 {
1562 goto L_100;
1563 }
1564 case 4:
1565 {
1566 goto L_400;
1567 }
1568 case 5:
1569 {
1570 goto L_500;
1571 }
1572 case 6:
1573 {
1574 goto L_600;
1575 }
1576 case 7:
1577 {
1578 goto L_700;
1579 }
1580 case 8:
1581 {
1582 goto L_800;
1583 }
1584 /* this is impossible, added by gjf to stop lint from complaining */
1585 default:
1586 {
1587 /* Invalid value for ISPEC
1588 * */
1589 ILAENV_v = -1;
1590 return ILAENV_v;
1591 }
1592 }
1593
1594L_100:
1595
1596 /* Convert NAME to upper case if the first character is lower case.
1597 * */
1598 ILAENV_v = 1;
1599 strncpy( SUBNAM, NAME, 6 );
1600 SUBNAM[6] = '\0';
1601 IC = (SUBNAM[0]);
1602 IZ = 'Z';
1603 if( IZ == 90 || IZ == 122 )
1604 {
1605
1606 /* ASCII character set
1607 * */
1608 if( IC >= 97 && IC <= 122 )
1609 {
1610 SUBNAM[0] = (char)(IC - 32);
1611 for( I=2; I <= 6; I++ )
1612 {
1613 IC = (SUBNAM[I-1]);
1614 if( IC >= 97 && IC <= 122 )
1615 SUBNAM[I - 1] = (char)(IC - 32);
1616 }
1617 }
1618
1619 }
1620 else if( IZ == 233 || IZ == 169 )
1621 {
1622
1623 /* EBCDIC character set
1624 * */
1625 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <= 153)) ||
1626 (IC >= 162 && IC <= 169) )
1627 {
1628 SUBNAM[0] = (char)(IC + 64);
1629 for( I=2; I <= 6; I++ )
1630 {
1631 IC = (SUBNAM[I-1]);
1632 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <=
1633 153)) || (IC >= 162 && IC <= 169) )
1634 SUBNAM[I - 1] = (char)(IC + 64);
1635 }
1636 }
1637
1638 }
1639 else if( IZ == 218 || IZ == 250 )
1640 {
1641
1642 /* Prime machines: ASCII+128
1643 * */
1644 if( IC >= 225 && IC <= 250 )
1645 {
1646 SUBNAM[0] = (char)(IC - 32);
1647 for( I=2; I <= 6; I++ )
1648 {
1649 IC = (SUBNAM[I-1]);
1650 if( IC >= 225 && IC <= 250 )
1651 SUBNAM[I - 1] = (char)(IC - 32);
1652 }
1653 }
1654 }
1655
1656 C1 = SUBNAM[0];
1657 SNAME = C1 == 'S' || C1 == 'D';
1658 CNAME = C1 == 'C' || C1 == 'Z';
1659 if( !(CNAME || SNAME) )
1660 {
1661 return ILAENV_v;
1662 }
1663
1664# if 0
1665 strncpy( C2, SUBNAM+1, 2 );
1666 strncpy( C3, SUBNAM+3, 3 );
1667 strncpy( C4, C3+1, 2 );
1668# endif
1669
1670 /* >>chng 00 nov 08, from above to below, insure had run
1671 * off end of string*/
1672 strncpy( C2, SUBNAM+1, 2 );
1673 C2[2] = '\0';
1674 strncpy( C3, SUBNAM+3, 3 );
1675 C3[3] = '\0';
1676 strncpy( C4, C3+1, 2 );
1677 C4[2] = '\0';
1678
1679 switch( ISPEC )
1680 {
1681 case 1: goto L_110;
1682 case 2: goto L_200;
1683 case 3: goto L_300;
1684 }
1685
1686L_110:
1687
1688 /* ISPEC = 1: block size
1689 *
1690 * In these examples, separate code is provided for setting NB for
1691 * real and complex. We assume that NB will take the same value in
1692 * single or double precision.
1693 * */
1694 NB = 1;
1695
1696 if( strcmp(C2,"GE") == 0 )
1697 {
1698 if( strcmp(C3,"TRF") == 0 )
1699 {
1700 if( SNAME )
1701 {
1702 NB = 64;
1703 }
1704 else
1705 {
1706 NB = 64;
1707 }
1708 }
1709 else if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) ||
1710 strcmp(C3,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
1711 {
1712 if( SNAME )
1713 {
1714 NB = 32;
1715 }
1716 else
1717 {
1718 NB = 32;
1719 }
1720 }
1721 else if( strcmp(C3,"HRD") == 0 )
1722 {
1723 if( SNAME )
1724 {
1725 NB = 32;
1726 }
1727 else
1728 {
1729 NB = 32;
1730 }
1731 }
1732 else if( strcmp(C3,"BRD") == 0 )
1733 {
1734 if( SNAME )
1735 {
1736 NB = 32;
1737 }
1738 else
1739 {
1740 NB = 32;
1741 }
1742 }
1743 else if( strcmp(C3,"TRI") == 0 )
1744 {
1745 if( SNAME )
1746 {
1747 NB = 64;
1748 }
1749 else
1750 {
1751 NB = 64;
1752 }
1753 }
1754 }
1755 else if( strcmp(C2,"PO") == 0 )
1756 {
1757 if( strcmp(C3,"TRF") == 0 )
1758 {
1759 if( SNAME )
1760 {
1761 NB = 64;
1762 }
1763 else
1764 {
1765 NB = 64;
1766 }
1767 }
1768 }
1769 else if( strcmp(C2,"SY") == 0 )
1770 {
1771 if( strcmp(C3,"TRF") == 0 )
1772 {
1773 if( SNAME )
1774 {
1775 NB = 64;
1776 }
1777 else
1778 {
1779 NB = 64;
1780 }
1781 }
1782 else if( SNAME && strcmp(C3,"TRD") == 0 )
1783 {
1784 NB = 1;
1785 }
1786 else if( SNAME && strcmp(C3,"GST") == 0 )
1787 {
1788 NB = 64;
1789 }
1790 }
1791 else if( CNAME && strcmp(C2,"HE") == 0 )
1792 {
1793 if( strcmp(C3,"TRF") == 0 )
1794 {
1795 NB = 64;
1796 }
1797 else if( strcmp(C3,"TRD") == 0 )
1798 {
1799 NB = 1;
1800 }
1801 else if( strcmp(C3,"GST") == 0 )
1802 {
1803 NB = 64;
1804 }
1805 }
1806 else if( SNAME && strcmp(C2,"OR") == 0 )
1807 {
1808 if( C3[0] == 'G' )
1809 {
1810 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1811 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1812 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1813 0 )
1814 {
1815 NB = 32;
1816 }
1817 }
1818 else if( C3[0] == 'M' )
1819 {
1820 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1821 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1822 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1823 0 )
1824 {
1825 NB = 32;
1826 }
1827 }
1828 }
1829 else if( CNAME && strcmp(C2,"UN") == 0 )
1830 {
1831 if( C3[0] == 'G' )
1832 {
1833 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1834 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1835 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1836 0 )
1837 {
1838 NB = 32;
1839 }
1840 }
1841 else if( C3[0] == 'M' )
1842 {
1843 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1844 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1845 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1846 0 )
1847 {
1848 NB = 32;
1849 }
1850 }
1851 }
1852 else if( strcmp(C2,"GB") == 0 )
1853 {
1854 if( strcmp(C3,"TRF") == 0 )
1855 {
1856 if( SNAME )
1857 {
1858 if( N4 <= 64 )
1859 {
1860 NB = 1;
1861 }
1862 else
1863 {
1864 NB = 32;
1865 }
1866 }
1867 else
1868 {
1869 if( N4 <= 64 )
1870 {
1871 NB = 1;
1872 }
1873 else
1874 {
1875 NB = 32;
1876 }
1877 }
1878 }
1879 }
1880 else if( strcmp(C2,"PB") == 0 )
1881 {
1882 if( strcmp(C3,"TRF") == 0 )
1883 {
1884 if( SNAME )
1885 {
1886 if( N2 <= 64 )
1887 {
1888 NB = 1;
1889 }
1890 else
1891 {
1892 NB = 32;
1893 }
1894 }
1895 else
1896 {
1897 if( N2 <= 64 )
1898 {
1899 NB = 1;
1900 }
1901 else
1902 {
1903 NB = 32;
1904 }
1905 }
1906 }
1907 }
1908 else if( strcmp(C2,"TR") == 0 )
1909 {
1910 if( strcmp(C3,"TRI") == 0 )
1911 {
1912 if( SNAME )
1913 {
1914 NB = 64;
1915 }
1916 else
1917 {
1918 NB = 64;
1919 }
1920 }
1921 }
1922 else if( strcmp(C2,"LA") == 0 )
1923 {
1924 if( strcmp(C3,"UUM") == 0 )
1925 {
1926 if( SNAME )
1927 {
1928 NB = 64;
1929 }
1930 else
1931 {
1932 NB = 64;
1933 }
1934 }
1935 }
1936 else if( SNAME && strcmp(C2,"ST") == 0 )
1937 {
1938 if( strcmp(C3,"EBZ") == 0 )
1939 {
1940 NB = 1;
1941 }
1942 }
1943 ILAENV_v = NB;
1944 return ILAENV_v;
1945
1946L_200:
1947
1948 /* ISPEC = 2: minimum block size
1949 * */
1950 NBMIN = 2;
1951 if( strcmp(C2,"GE") == 0 )
1952 {
1953 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
1954 ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
1955 {
1956 if( SNAME )
1957 {
1958 NBMIN = 2;
1959 }
1960 else
1961 {
1962 NBMIN = 2;
1963 }
1964 }
1965 else if( strcmp(C3,"HRD") == 0 )
1966 {
1967 if( SNAME )
1968 {
1969 NBMIN = 2;
1970 }
1971 else
1972 {
1973 NBMIN = 2;
1974 }
1975 }
1976 else if( strcmp(C3,"BRD") == 0 )
1977 {
1978 if( SNAME )
1979 {
1980 NBMIN = 2;
1981 }
1982 else
1983 {
1984 NBMIN = 2;
1985 }
1986 }
1987 else if( strcmp(C3,"TRI") == 0 )
1988 {
1989 if( SNAME )
1990 {
1991 NBMIN = 2;
1992 }
1993 else
1994 {
1995 NBMIN = 2;
1996 }
1997 }
1998 }
1999 else if( strcmp(C2,"SY") == 0 )
2000 {
2001 if( strcmp(C3,"TRF") == 0 )
2002 {
2003 if( SNAME )
2004 {
2005 NBMIN = 8;
2006 }
2007 else
2008 {
2009 NBMIN = 8;
2010 }
2011 }
2012 else if( SNAME && strcmp(C3,"TRD") == 0 )
2013 {
2014 NBMIN = 2;
2015 }
2016 }
2017 else if( CNAME && strcmp(C2,"HE") == 0 )
2018 {
2019 if( strcmp(C3,"TRD") == 0 )
2020 {
2021 NBMIN = 2;
2022 }
2023 }
2024 else if( SNAME && strcmp(C2,"OR") == 0 )
2025 {
2026 if( C3[0] == 'G' )
2027 {
2028 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2029 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2030 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2031 0 )
2032 {
2033 NBMIN = 2;
2034 }
2035 }
2036 else if( C3[0] == 'M' )
2037 {
2038 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2039 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2040 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2041 0 )
2042 {
2043 NBMIN = 2;
2044 }
2045 }
2046 }
2047 else if( CNAME && strcmp(C2,"UN") == 0 )
2048 {
2049 if( C3[0] == 'G' )
2050 {
2051 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2052 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2053 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2054 0 )
2055 {
2056 NBMIN = 2;
2057 }
2058 }
2059 else if( C3[0] == 'M' )
2060 {
2061 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2062 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2063 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2064 0 )
2065 {
2066 NBMIN = 2;
2067 }
2068 }
2069 }
2070 ILAENV_v = NBMIN;
2071 return ILAENV_v;
2072
2073L_300:
2074
2075 /* ISPEC = 3: crossover point
2076 * */
2077 NX = 0;
2078 if( strcmp(C2,"GE") == 0 )
2079 {
2080 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
2081 ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
2082 {
2083 if( SNAME )
2084 {
2085 NX = 128;
2086 }
2087 else
2088 {
2089 NX = 128;
2090 }
2091 }
2092 else if( strcmp(C3,"HRD") == 0 )
2093 {
2094 if( SNAME )
2095 {
2096 NX = 128;
2097 }
2098 else
2099 {
2100 NX = 128;
2101 }
2102 }
2103 else if( strcmp(C3,"BRD") == 0 )
2104 {
2105 if( SNAME )
2106 {
2107 NX = 128;
2108 }
2109 else
2110 {
2111 NX = 128;
2112 }
2113 }
2114 }
2115 else if( strcmp(C2,"SY") == 0 )
2116 {
2117 if( SNAME && strcmp(C3,"TRD") == 0 )
2118 {
2119 NX = 1;
2120 }
2121 }
2122 else if( CNAME && strcmp(C2,"HE") == 0 )
2123 {
2124 if( strcmp(C3,"TRD") == 0 )
2125 {
2126 NX = 1;
2127 }
2128 }
2129 else if( SNAME && strcmp(C2,"OR") == 0 )
2130 {
2131 if( C3[0] == 'G' )
2132 {
2133 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2134 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2135 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2136 0 )
2137 {
2138 NX = 128;
2139 }
2140 }
2141 }
2142 else if( CNAME && strcmp(C2,"UN") == 0 )
2143 {
2144 if( C3[0] == 'G' )
2145 {
2146 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2147 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2148 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2149 0 )
2150 {
2151 NX = 128;
2152 }
2153 }
2154 }
2155 ILAENV_v = NX;
2156 return ILAENV_v;
2157
2158L_400:
2159
2160 /* ISPEC = 4: number of shifts (used by xHSEQR)
2161 * */
2162 ILAENV_v = 6;
2163 return ILAENV_v;
2164
2165L_500:
2166
2167 /* ISPEC = 5: minimum column dimension (not used)
2168 * */
2169 ILAENV_v = 2;
2170 return ILAENV_v;
2171
2172L_600:
2173
2174 /* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
2175 * */
2176 ILAENV_v = (int32)((realnum)(MIN2(N1,N2))*1.6e0);
2177 return ILAENV_v;
2178
2179L_700:
2180
2181 /* ISPEC = 7: number of processors (not used)
2182 * */
2183 ILAENV_v = 1;
2184 return ILAENV_v;
2185
2186L_800:
2187
2188 /* ISPEC = 8: crossover point for multishift (used by xHSEQR)
2189 * */
2190 ILAENV_v = 50;
2191 return ILAENV_v;
2192
2193 /* End of ILAENV
2194 * */
2195}
2196
2197/*DSWAP lapack routine */
2198
2199STATIC void DSWAP(int32 n,
2200 double dx[],
2201 int32 incx,
2202 double dy[],
2203 int32 incy)
2204{
2205 int32 i,
2206 ix,
2207 iy,
2208 m;
2209 double dtemp;
2210
2211 DEBUG_ENTRY( "DSWAP()" );
2212
2213 /* interchanges two vectors.
2214 * uses unrolled loops for increments equal one.
2215 * jack dongarra, lapack, 3/11/78.
2216 * modified 12/3/93, array(1) declarations changed to array(*)
2217 * */
2218
2219 if( n <= 0 )
2220 {
2221 return;}
2222 if( incx == 1 && incy == 1 )
2223 goto L_20;
2224
2225 /* code for unequal increments or equal increments not equal
2226 * to 1
2227 * */
2228 ix = 1;
2229 iy = 1;
2230
2231 if( incx < 0 )
2232 ix = (-n + 1)*incx + 1;
2233
2234 if( incy < 0 )
2235 iy = (-n + 1)*incy + 1;
2236
2237 for( i=0; i < n; i++ )
2238 {
2239 dtemp = dx[ix-1];
2240 dx[ix-1] = dy[iy-1];
2241 dy[iy-1] = dtemp;
2242 ix += incx;
2243 iy += incy;
2244 }
2245 return;
2246
2247 /* code for both increments equal to 1
2248 *
2249 *
2250 * clean-up loop
2251 * */
2252L_20:
2253 m = n%3;
2254 if( m == 0 )
2255 goto L_40;
2256
2257 for( i=0; i < m; i++ )
2258 {
2259 dtemp = dx[i];
2260 dx[i] = dy[i];
2261 dy[i] = dtemp;
2262 }
2263
2264 if( n < 3 )
2265 {
2266 return;
2267 }
2268
2269L_40:
2270 for( i=m; i < n; i += 3 )
2271 {
2272 dtemp = dx[i];
2273 dx[i] = dy[i];
2274 dy[i] = dtemp;
2275 dtemp = dx[i+1];
2276 dx[i+1] = dy[i+1];
2277 dy[i+1] = dtemp;
2278 dtemp = dx[i+2];
2279 dx[i+2] = dy[i+2];
2280 dy[i+2] = dtemp;
2281 }
2282 return;
2283}
2284
2285/*DSCAL lapack routine */
2286
2287STATIC void DSCAL(int32 n,
2288 double da,
2289 double dx[],
2290 int32 incx)
2291{
2292 int32 i,
2293 m,
2294 nincx;
2295
2296 DEBUG_ENTRY( "DSCAL()" );
2297
2298 /* scales a vector by a constant.
2299 * uses unrolled loops for increment equal to one.
2300 * jack dongarra, lapack, 3/11/78.
2301 * modified 3/93 to return if incx .le. 0.
2302 * modified 12/3/93, array(1) declarations changed to array(*)
2303 * */
2304
2305 if( n <= 0 || incx <= 0 )
2306 {
2307 return;}
2308 if( incx == 1 )
2309 goto L_20;
2310
2311 /* code for increment not equal to 1
2312 * */
2313 nincx = n*incx;
2314 /*for( i=1, _do0=DOCNT(i,nincx,_do1 = incx); _do0 > 0; i += _do1, _do0-- )*/
2315 for( i=0; i<nincx; i = i + incx)
2316 {
2317 dx[i] *= da;
2318 }
2319 return;
2320
2321 /* code for increment equal to 1
2322 *
2323 *
2324 * clean-up loop
2325 * */
2326L_20:
2327 m = n%5;
2328 if( m == 0 )
2329 goto L_40;
2330
2331 for( i=0; i < m; i++ )
2332 {
2333 dx[i] *= da;
2334 }
2335
2336 if( n < 5 )
2337 {
2338 return;
2339 }
2340
2341L_40:
2342
2343 for( i=m; i < n; i += 5 )
2344 {
2345 dx[i] *= da;
2346 dx[i+1] *= da;
2347 dx[i+2] *= da;
2348 dx[i+3] *= da;
2349 dx[i+4] *= da;
2350 }
2351 return;
2352}
2353
2354/*DLASWP -- LAPACK auxiliary routine (version 2.0) --*/
2355
2356STATIC void DLASWP(int32 N,
2357 double *A,
2358 int32 LDA,
2359 int32 K1,
2360 int32 K2,
2361 int32 IPIV[],
2362 int32 INCX)
2363{
2364 int32 I,
2365 IP,
2366 IX,
2367 I_;
2368
2369 DEBUG_ENTRY( "DLASWP()" );
2370
2371 /* -- LAPACK auxiliary routine (version 2.0) --
2372 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2373 * Courant Institute, Argonne National Lab, and Rice University
2374 * October 31, 1992
2375 *
2376 * .. Scalar Arguments .. */
2377 /* ..
2378 * .. Array Arguments .. */
2379 /* ..
2380 *
2381 * Purpose
2382 * =======
2383 *
2384 * DLASWP performs a series of row interchanges on the matrix A.
2385 * One row interchange is initiated for each of rows K1 through K2 of A.
2386 *
2387 * Arguments
2388 * =========
2389 *
2390 * N (input) INTEGER
2391 * The number of columns of the matrix A.
2392 *
2393 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
2394 * On entry, the matrix of column dimension N to which the row
2395 * interchanges will be applied.
2396 * On exit, the permuted matrix.
2397 *
2398 * LDA (input) INTEGER
2399 * The leading dimension of the array A.
2400 *
2401 * K1 (input) INTEGER
2402 * The first element of IPIV for which a row interchange will
2403 * be done.
2404 *
2405 * K2 (input) INTEGER
2406 * The last element of IPIV for which a row interchange will
2407 * be done.
2408 *
2409 * IPIV (input) INTEGER array, dimension (M*ABS(INCX))
2410 * The vector of pivot indices. Only the elements in positions
2411 * K1 through K2 of IPIV are accessed.
2412 * IPIV(K) = L implies rows K and L are to be interchanged.
2413 *
2414 * INCX (input) INTEGER
2415 * The increment between successive values of IPIV. If IPIV
2416 * is negative, the pivots are applied in reverse order.
2417 *
2418 * =====================================================================
2419 *
2420 * .. Local Scalars .. */
2421 /* ..
2422 * .. External Subroutines .. */
2423 /* ..
2424 * .. Executable Statements ..
2425 *
2426 * Interchange row I with row IPIV(I) for each of rows K1 through K2.
2427 * */
2428 if( INCX == 0 )
2429 {
2430 return;}
2431 if( INCX > 0 )
2432 {
2433 IX = K1;
2434 }
2435 else
2436 {
2437 IX = 1 + (1 - K2)*INCX;
2438 }
2439 if( INCX == 1 )
2440 {
2441 for( I=K1; I <= K2; I++ )
2442 {
2443 I_ = I - 1;
2444 IP = IPIV[I_];
2445 if( IP != I )
2446 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2447 }
2448 }
2449 else if( INCX > 1 )
2450 {
2451 for( I=K1; I <= K2; I++ )
2452 {
2453 I_ = I - 1;
2454 IP = IPIV[IX-1];
2455 if( IP != I )
2456 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2457 IX += INCX;
2458 }
2459 }
2460 else if( INCX < 0 )
2461 {
2462 for( I=K2; I >= K1; I-- )
2463 {
2464 I_ = I - 1;
2465 IP = IPIV[IX-1];
2466 if( IP != I )
2467 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2468 IX += INCX;
2469 }
2470 }
2471
2472 return;
2473
2474 /* End of DLASWP
2475 * */
2476#undef A
2477}
2478
2479/*DGETF2 lapack service routine */
2480
2481STATIC void DGETF2(int32 M,
2482 int32 N,
2483 double *A,
2484 int32 LDA,
2485 int32 IPIV[],
2486 int32 *INFO)
2487{
2488 int32 J,
2489 JP,
2490 J_,
2491 limit;
2492
2493 DEBUG_ENTRY( "DGETF2()" );
2494
2495 /* -- LAPACK routine (version 2.0) --
2496 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2497 * Courant Institute, Argonne National Lab, and Rice University
2498 * June 30, 1992
2499 *
2500 * .. Scalar Arguments .. */
2501 /* ..
2502 * .. Array Arguments .. */
2503 /* ..
2504 *
2505 * Purpose
2506 * =======
2507 *
2508 * DGETF2 computes an LU factorization of a general m-by-n matrix A
2509 * using partial pivoting with row interchanges.
2510 *
2511 * The factorization has the form
2512 * A = P * L * U
2513 * where P is a permutation matrix, L is lower triangular with unit
2514 * diagonal elements (lower trapezoidal if m > n), and U is upper
2515 * triangular (upper trapezoidal if m < n).
2516 *
2517 * This is the right-looking Level 2 BLAS version of the algorithm.
2518 *
2519 * Arguments
2520 * =========
2521 *
2522 * M (input) INTEGER
2523 * The number of rows of the matrix A. M >= 0.
2524 *
2525 * N (input) INTEGER
2526 * The number of columns of the matrix A. N >= 0.
2527 *
2528 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
2529 * On entry, the m by n matrix to be factored.
2530 * On exit, the factors L and U from the factorization
2531 * A = P*L*U; the unit diagonal elements of L are not stored.
2532 *
2533 * LDA (input) INTEGER
2534 * The leading dimension of the array A. LDA >= MAX(1,M).
2535 *
2536 * IPIV (output) INTEGER array, dimension (MIN(M,N))
2537 * The pivot indices; for 1 <= i <= MIN(M,N), row i of the
2538 * matrix was interchanged with row IPIV(i).
2539 *
2540 * INFO (output) INTEGER
2541 * = 0: successful exit
2542 * < 0: if INFO = -k, the k-th argument had an illegal value
2543 * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
2544 * has been completed, but the factor U is exactly
2545 * singular, and division by zero will occur if it is used
2546 * to solve a system of equations.
2547 *
2548 * =====================================================================
2549 *
2550 * .. Parameters .. */
2551 /* ..
2552 * .. Local Scalars .. */
2553 /* ..
2554 * .. External Functions .. */
2555 /* ..
2556 * .. External Subroutines .. */
2557 /* ..
2558 * .. Intrinsic Functions .. */
2559 /* ..
2560 * .. Executable Statements ..
2561 *
2562 * Test the input parameters.
2563 * */
2564 *INFO = 0;
2565 if( M < 0 )
2566 {
2567 *INFO = -1;
2568 }
2569 else if( N < 0 )
2570 {
2571 *INFO = -2;
2572 }
2573 else if( LDA < MAX2(1,M) )
2574 {
2575 *INFO = -4;
2576 }
2577 if( *INFO != 0 )
2578 {
2579 XERBLA("DGETF2",-*INFO);
2580 /* XERBLA does not return */
2581 }
2582
2583 /* Quick return if possible
2584 * */
2585 if( M == 0 || N == 0 )
2586 {
2587 return;}
2588
2589 limit = MIN2(M,N);
2590 for( J=1; J <= limit; J++ )
2591 {
2592 J_ = J - 1;
2593
2594 /* Find pivot and test for singularity.
2595 * */
2596 JP = J - 1 + IDAMAX(M-J+1,&AA(J_,J_),1);
2597 IPIV[J_] = JP;
2598 if( AA(J_,JP-1) != ZERO )
2599 {
2600 /* Apply the interchange to columns 1:N.
2601 * */
2602 if( JP != J )
2603 DSWAP(N,&AA(0,J_),LDA,&AA(0,JP-1),LDA);
2604
2605 /* Compute elements J+1:M of J-th column.
2606 * */
2607 if( J < M )
2608 DSCAL(M-J,ONE/AA(J_,J_),&AA(J_,J_+1),1);
2609 }
2610 else if( *INFO == 0 )
2611 {
2612 *INFO = J;
2613 }
2614
2615 if( J < MIN2(M,N) )
2616 {
2617 /* Update trailing submatrix.
2618 * */
2619 DGER(M-J,N-J,-ONE,&AA(J_,J_+1),1,&AA(J_+1,J_),LDA,&AA(J_+1,J_+1),
2620 LDA);
2621 }
2622 }
2623 return;
2624
2625 /* End of DGETF2
2626 * */
2627#undef A
2628}
2629
2630/*DGER service routine for matrix inversion */
2631
2632STATIC void DGER(int32 M,
2633 int32 N,
2634 double ALPHA,
2635 double X[],
2636 int32 INCX,
2637 double Y[],
2638 int32 INCY,
2639 double *A,
2640 int32 LDA)
2641{
2642 int32 I,
2643 INFO,
2644 IX,
2645 I_,
2646 J,
2647 JY,
2648 J_,
2649 KX;
2650 double TEMP;
2651
2652 DEBUG_ENTRY( "DGER()" );
2653 /* .. Scalar Arguments .. */
2654 /* .. Array Arguments .. */
2655 /* ..
2656 *
2657 * Purpose
2658 * =======
2659 *
2660 * DGER performs the rank 1 operation
2661 *
2662 * A := alpha*x*y' + A,
2663 *
2664 * where alpha is a scalar, x is an m element vector, y is an n element
2665 * vector and A is an m by n matrix.
2666 *
2667 * Parameters
2668 * ==========
2669 *
2670 * M - INTEGER.
2671 * On entry, M specifies the number of rows of the matrix A.
2672 * M must be at least zero.
2673 * Unchanged on exit.
2674 *
2675 * N - INTEGER.
2676 * On entry, N specifies the number of columns of the matrix A.
2677 * N must be at least zero.
2678 * Unchanged on exit.
2679 *
2680 * ALPHA - DOUBLE PRECISION.
2681 * On entry, ALPHA specifies the scalar alpha.
2682 * Unchanged on exit.
2683 *
2684 * X - DOUBLE PRECISION array of dimension at least
2685 * ( 1 + ( m - 1 )*ABS( INCX ) ).
2686 * Before entry, the incremented array X must contain the m
2687 * element vector x.
2688 * Unchanged on exit.
2689 *
2690 * INCX - INTEGER.
2691 * On entry, INCX specifies the increment for the elements of
2692 * X. INCX must not be zero.
2693 * Unchanged on exit.
2694 *
2695 * Y - DOUBLE PRECISION array of dimension at least
2696 * ( 1 + ( n - 1 )*ABS( INCY ) ).
2697 * Before entry, the incremented array Y must contain the n
2698 * element vector y.
2699 * Unchanged on exit.
2700 *
2701 * INCY - INTEGER.
2702 * On entry, INCY specifies the increment for the elements of
2703 * Y. INCY must not be zero.
2704 * Unchanged on exit.
2705 *
2706 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2707 * Before entry, the leading m by n part of the array A must
2708 * contain the matrix of coefficients. On exit, A is
2709 * overwritten by the updated matrix.
2710 *
2711 * LDA - INTEGER.
2712 * On entry, LDA specifies the first dimension of A as declared
2713 * in the calling (sub) program. LDA must be at least
2714 * MAX( 1, m ).
2715 * Unchanged on exit.
2716 *
2717 *
2718 * Level 2 Blas routine.
2719 *
2720 * -- Written on 22-October-1986.
2721 * Jack Dongarra, Argonne National Lab.
2722 * Jeremy Du Croz, Nag Central Office.
2723 * Sven Hammarling, Nag Central Office.
2724 * Richard Hanson, Sandia National Labs.
2725 *
2726 *
2727 * .. Parameters .. */
2728 /* .. Local Scalars .. */
2729 /* .. External Subroutines .. */
2730 /* .. Intrinsic Functions .. */
2731 /* ..
2732 * .. Executable Statements ..
2733 *
2734 * Test the input parameters.
2735 * */
2736 INFO = 0;
2737 if( M < 0 )
2738 {
2739 INFO = 1;
2740 }
2741 else if( N < 0 )
2742 {
2743 INFO = 2;
2744 }
2745 else if( INCX == 0 )
2746 {
2747 INFO = 5;
2748 }
2749 else if( INCY == 0 )
2750 {
2751 INFO = 7;
2752 }
2753 else if( LDA < MAX2(1,M) )
2754 {
2755 INFO = 9;
2756 }
2757 if( INFO != 0 )
2758 {
2759 XERBLA("DGER ",INFO);
2760 /* XERBLA does not return */
2761 }
2762
2763 /* Quick return if possible.
2764 * */
2765 if( ((M == 0) || (N == 0)) || (ALPHA == ZERO) )
2766 {
2767 return;}
2768
2769 /* Start the operations. In this version the elements of A are
2770 * accessed sequentially with one pass through A.
2771 * */
2772 if( INCY > 0 )
2773 {
2774 JY = 1;
2775 }
2776 else
2777 {
2778 JY = 1 - (N - 1)*INCY;
2779 }
2780 if( INCX == 1 )
2781 {
2782 for( J=1; J <= N; J++ )
2783 {
2784 J_ = J - 1;
2785 if( Y[JY-1] != ZERO )
2786 {
2787 TEMP = ALPHA*Y[JY-1];
2788 for( I=1; I <= M; I++ )
2789 {
2790 I_ = I - 1;
2791 AA(J_,I_) += X[I_]*TEMP;
2792 }
2793 }
2794 JY += INCY;
2795 }
2796 }
2797 else
2798 {
2799 if( INCX > 0 )
2800 {
2801 KX = 1;
2802 }
2803 else
2804 {
2805 KX = 1 - (M - 1)*INCX;
2806 }
2807 for( J=1; J <= N; J++ )
2808 {
2809 J_ = J - 1;
2810 if( Y[JY-1] != ZERO )
2811 {
2812 TEMP = ALPHA*Y[JY-1];
2813 IX = KX;
2814 for( I=1; I <= M; I++ )
2815 {
2816 I_ = I - 1;
2817 AA(J_,I_) += X[IX-1]*TEMP;
2818 IX += INCX;
2819 }
2820 }
2821 JY += INCY;
2822 }
2823 }
2824
2825 return;
2826
2827 /* End of DGER .
2828 * */
2829#undef A
2830}
2831
2832/* DGEMM matrix inversion helper routine*/
2833
2834STATIC void DGEMM(int32 TRANSA,
2835 int32 TRANSB,
2836 int32 M,
2837 int32 N,
2838 int32 K,
2839 double ALPHA,
2840 double *A,
2841 int32 LDA,
2842 double *B,
2843 int32 LDB,
2844 double BETA,
2845 double *C,
2846 int32 LDC)
2847{
2848 int32 NOTA,
2849 NOTB;
2850 int32 I,
2851 INFO,
2852 I_,
2853 J,
2854 J_,
2855 L,
2856 L_,
2857 NROWA,
2858 NROWB;
2859 double TEMP;
2860
2861 DEBUG_ENTRY( "DGEMM()" );
2862 /* .. Scalar Arguments .. */
2863 /* .. Array Arguments .. */
2864 /* ..
2865 *
2866 * Purpose
2867 * =======
2868 *
2869 * DGEMM performs one of the matrix-matrix operations
2870 *
2871 * C := alpha*op( A )*op( B ) + beta*C,
2872 *
2873 * where op( X ) is one of
2874 *
2875 * op( X ) = X or op( X ) = X',
2876 *
2877 * alpha and beta are scalars, and A, B and C are matrices, with op( A )
2878 * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
2879 *
2880 * Parameters
2881 * ==========
2882 *
2883 * TRANSA - CHARACTER*1.
2884 * On entry, TRANSA specifies the form of op( A ) to be used in
2885 * the matrix multiplication as follows:
2886 *
2887 * TRANSA = 'N' or 'n', op( A ) = A.
2888 *
2889 * TRANSA = 'T' or 't', op( A ) = A'.
2890 *
2891 * TRANSA = 'C' or 'c', op( A ) = A'.
2892 *
2893 * Unchanged on exit.
2894 *
2895 * TRANSB - CHARACTER*1.
2896 * On entry, TRANSB specifies the form of op( B ) to be used in
2897 * the matrix multiplication as follows:
2898 *
2899 * TRANSB = 'N' or 'n', op( B ) = B.
2900 *
2901 * TRANSB = 'T' or 't', op( B ) = B'.
2902 *
2903 * TRANSB = 'C' or 'c', op( B ) = B'.
2904 *
2905 * Unchanged on exit.
2906 *
2907 * M - INTEGER.
2908 * On entry, M specifies the number of rows of the matrix
2909 * op( A ) and of the matrix C. M must be at least zero.
2910 * Unchanged on exit.
2911 *
2912 * N - INTEGER.
2913 * On entry, N specifies the number of columns of the matrix
2914 * op( B ) and the number of columns of the matrix C. N must be
2915 * at least zero.
2916 * Unchanged on exit.
2917 *
2918 * K - INTEGER.
2919 * On entry, K specifies the number of columns of the matrix
2920 * op( A ) and the number of rows of the matrix op( B ). K must
2921 * be at least zero.
2922 * Unchanged on exit.
2923 *
2924 * ALPHA - DOUBLE PRECISION.
2925 * On entry, ALPHA specifies the scalar alpha.
2926 * Unchanged on exit.
2927 *
2928 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
2929 * k when TRANSA = 'N' or 'n', and is m otherwise.
2930 * Before entry with TRANSA = 'N' or 'n', the leading m by k
2931 * part of the array A must contain the matrix A, otherwise
2932 * the leading k by m part of the array A must contain the
2933 * matrix A.
2934 * Unchanged on exit.
2935 *
2936 * LDA - INTEGER.
2937 * On entry, LDA specifies the first dimension of A as declared
2938 * in the calling (sub) program. When TRANSA = 'N' or 'n' then
2939 * LDA must be at least MAX( 1, m ), otherwise LDA must be at
2940 * least MAX( 1, k ).
2941 * Unchanged on exit.
2942 *
2943 * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
2944 * n when TRANSB = 'N' or 'n', and is k otherwise.
2945 * Before entry with TRANSB = 'N' or 'n', the leading k by n
2946 * part of the array B must contain the matrix B, otherwise
2947 * the leading n by k part of the array B must contain the
2948 * matrix B.
2949 * Unchanged on exit.
2950 *
2951 * LDB - INTEGER.
2952 * On entry, LDB specifies the first dimension of B as declared
2953 * in the calling (sub) program. When TRANSB = 'N' or 'n' then
2954 * LDB must be at least MAX( 1, k ), otherwise LDB must be at
2955 * least MAX( 1, n ).
2956 * Unchanged on exit.
2957 *
2958 * BETA - DOUBLE PRECISION.
2959 * On entry, BETA specifies the scalar beta. When BETA is
2960 * supplied as zero then C need not be set on input.
2961 * Unchanged on exit.
2962 *
2963 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
2964 * Before entry, the leading m by n part of the array C must
2965 * contain the matrix C, except when beta is zero, in which
2966 * case C need not be set on entry.
2967 * On exit, the array C is overwritten by the m by n matrix
2968 * ( alpha*op( A )*op( B ) + beta*C ).
2969 *
2970 * LDC - INTEGER.
2971 * On entry, LDC specifies the first dimension of C as declared
2972 * in the calling (sub) program. LDC must be at least
2973 * MAX( 1, m ).
2974 * Unchanged on exit.
2975 *
2976 *
2977 * Level 3 Blas routine.
2978 *
2979 * -- Written on 8-February-1989.
2980 * Jack Dongarra, Argonne National Laboratory.
2981 * Iain Duff, AERE Harwell.
2982 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
2983 * Sven Hammarling, Numerical Algorithms Group Ltd.
2984 *
2985 *
2986 * .. External Functions .. */
2987 /* .. External Subroutines .. */
2988 /* .. Intrinsic Functions .. */
2989 /* .. Local Scalars .. */
2990 /* .. Parameters .. */
2991 /* ..
2992 * .. Executable Statements ..
2993 *
2994 * Set NOTA and NOTB as true if A and B respectively are not
2995 * transposed and set NROWA, NROWB as the number of rows
2996 * and columns of A and the number of rows of B respectively.
2997 * */
2998 NOTA = LSAME(TRANSA,'N');
2999 NOTB = LSAME(TRANSB,'N');
3000
3001 if( NOTA )
3002 {
3003 NROWA = M;
3004 }
3005 else
3006 {
3007 NROWA = K;
3008 }
3009
3010 if( NOTB )
3011 {
3012 NROWB = K;
3013 }
3014 else
3015 {
3016 NROWB = N;
3017 }
3018
3019 /* Test the input parameters.
3020 * */
3021 INFO = 0;
3022 if( ((!NOTA) && (!LSAME(TRANSA,'C'))) && (!LSAME(TRANSA,'T')) )
3023 {
3024 INFO = 1;
3025 }
3026 else if(
3027 ((!NOTB) && (!LSAME(TRANSB,'C'))) && (!LSAME(TRANSB,'T')) )
3028 {
3029 INFO = 2;
3030 }
3031
3032 else if( M < 0 )
3033 {
3034 INFO = 3;
3035 }
3036
3037 else if( N < 0 )
3038 {
3039 INFO = 4;
3040 }
3041
3042 else if( K < 0 )
3043 {
3044 INFO = 5;
3045 }
3046
3047 else if( LDA < MAX2(1,NROWA) )
3048 {
3049 INFO = 8;
3050 }
3051
3052 else if( LDB < MAX2(1,NROWB) )
3053 {
3054 INFO = 10;
3055 }
3056
3057 else if( LDC < MAX2(1,M) )
3058 {
3059 INFO = 13;
3060 }
3061
3062 if( INFO != 0 )
3063 {
3064 XERBLA("DGEMM ",INFO);
3065 /* XERBLA does not return */
3066 }
3067
3068 /* Quick return if possible.
3069 * */
3070 if( ((M == 0) || (N == 0)) || (((ALPHA == ZERO) || (K == 0)) &&
3071 (BETA == ONE)) )
3072 {
3073 return;
3074 }
3075
3076 /* And if alpha.eq.zero. */
3077 if( ALPHA == ZERO )
3078 {
3079 if( BETA == ZERO )
3080 {
3081 for( J=1; J <= N; J++ )
3082 {
3083 J_ = J - 1;
3084 for( I=1; I <= M; I++ )
3085 {
3086 I_ = I - 1;
3087 CC(J_,I_) = ZERO;
3088 }
3089 }
3090 }
3091
3092 else
3093 {
3094 for( J=1; J <= N; J++ )
3095 {
3096 J_ = J - 1;
3097 for( I=1; I <= M; I++ )
3098 {
3099 I_ = I - 1;
3100 CC(J_,I_) *= BETA;
3101 }
3102 }
3103 }
3104 return;
3105 }
3106
3107 /* Start the operations.*/
3108 if( NOTB )
3109 {
3110
3111 if( NOTA )
3112 {
3113
3114 /* Form C := alpha*A*B + beta*C.
3115 * */
3116 for( J=1; J <= N; J++ )
3117 {
3118 J_ = J - 1;
3119 if( BETA == ZERO )
3120 {
3121 for( I=1; I <= M; I++ )
3122 {
3123 I_ = I - 1;
3124 CC(J_,I_) = ZERO;
3125 }
3126 }
3127
3128 else if( BETA != ONE )
3129 {
3130 for( I=1; I <= M; I++ )
3131 {
3132 I_ = I - 1;
3133 CC(J_,I_) *= BETA;
3134 }
3135 }
3136
3137 for( L=1; L <= K; L++ )
3138 {
3139 L_ = L - 1;
3140 if( BB(J_,L_) != ZERO )
3141 {
3142 TEMP = ALPHA*BB(J_,L_);
3143 for( I=1; I <= M; I++ )
3144 {
3145 I_ = I - 1;
3146 CC(J_,I_) += TEMP*AA(L_,I_);
3147 }
3148 }
3149 }
3150 }
3151 }
3152 else
3153 {
3154
3155 /* Form C := alpha*A'*B + beta*C */
3156 for( J=1; J <= N; J++ )
3157 {
3158 J_ = J - 1;
3159 for( I=1; I <= M; I++ )
3160 {
3161 I_ = I - 1;
3162 TEMP = ZERO;
3163 for( L=1; L <= K; L++ )
3164 {
3165 L_ = L - 1;
3166 TEMP += AA(I_,L_)*BB(J_,L_);
3167 }
3168
3169 if( BETA == ZERO )
3170 {
3171 CC(J_,I_) = ALPHA*TEMP;
3172 }
3173 else
3174 {
3175 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
3176 }
3177 }
3178 }
3179 }
3180 }
3181 else
3182 {
3183 if( NOTA )
3184 {
3185
3186 /* Form C := alpha*A*B' + beta*C
3187 * */
3188 for( J=1; J <= N; J++ )
3189 {
3190 J_ = J - 1;
3191 if( BETA == ZERO )
3192 {
3193 for( I=1; I <= M; I++ )
3194 {
3195 I_ = I - 1;
3196 CC(J_,I_) = ZERO;
3197 }
3198 }
3199
3200 else if( BETA != ONE )
3201 {
3202 for( I=1; I <= M; I++ )
3203 {
3204 I_ = I - 1;
3205 CC(J_,I_) *= BETA;
3206 }
3207 }
3208
3209 for( L=1; L <= K; L++ )
3210 {
3211 L_ = L - 1;
3212 if( BB(L_,J_) != ZERO )
3213 {
3214 TEMP = ALPHA*BB(L_,J_);
3215 for( I=1; I <= M; I++ )
3216 {
3217 I_ = I - 1;
3218 CC(J_,I_) += TEMP*AA(L_,I_);
3219 }
3220 }
3221 }
3222 }
3223 }
3224
3225 else
3226 {
3227
3228 /* Form C := alpha*A'*B' + beta*C */
3229 for( J=1; J <= N; J++ )
3230 {
3231 J_ = J - 1;
3232
3233 for( I=1; I <= M; I++ )
3234 {
3235 I_ = I - 1;
3236 TEMP = ZERO;
3237
3238 for( L=1; L <= K; L++ )
3239 {
3240 L_ = L - 1;
3241 TEMP += AA(I_,L_)*BB(L_,J_);
3242 }
3243
3244 if( BETA == ZERO )
3245 {
3246 CC(J_,I_) = ALPHA*TEMP;
3247 }
3248 else
3249 {
3250 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
3251 }
3252
3253 }
3254 }
3255 }
3256 }
3257
3258 return;
3259
3260 /* End of DGEMM .*/
3261#undef C
3262#undef B
3263#undef A
3264}
3265#undef AA
3266#undef BB
3267#undef CC
3268
3269/* Subroutine */
3270#if 0
3271STATIC int32 DGTSV(int32 *n, int32 *nrhs, double *dl,
3272 double *d__, double *du, double *b, int32 *ldb, int32
3273 *info)
3274{
3275 /* System generated locals */
3276 int32 b_dim1, b_offset, i__1, i__2;
3277
3278 /* Local variables */
3279 double fact, temp;
3280 int32 i__, j;
3281
3282
3283#define b_ref(a_1,a_2) b[(a_2)*(b_dim1) + (a_1)]
3284
3285
3286/* -- LAPACK routine (version 3.0) --
3287 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3288 Courant Institute, Argonne National Lab, and Rice University
3289 October 31, 1999
3290
3291
3292 Purpose
3293 =======
3294
3295 DGTSV solves the equation
3296
3297 A*X = B,
3298
3299 where A is an n by n tridiagonal matrix, by Gaussian elimination with
3300 partial pivoting.
3301
3302 Note that the equation A'*X = B may be solved by interchanging the
3303 order of the arguments DU and DL.
3304
3305 Arguments
3306 =========
3307
3308 N (input) INTEGER
3309 The order of the matrix A. N >= 0.
3310
3311 NRHS (input) INTEGER
3312 The number of right hand sides, i.e., the number of columns
3313 of the matrix B. NRHS >= 0.
3314
3315 DL (input/output) DOUBLE PRECISION array, dimension (N-1)
3316 On entry, DL must contain the (n-1) sub-diagonal elements of
3317 A.
3318
3319 On exit, DL is overwritten by the (n-2) elements of the
3320 second super-diagonal of the upper triangular matrix U from
3321 the LU factorization of A, in DL(1), ..., DL(n-2).
3322
3323 D (input/output) DOUBLE PRECISION array, dimension (N)
3324 On entry, D must contain the diagonal elements of A.
3325
3326 On exit, D is overwritten by the n diagonal elements of U.
3327
3328 DU (input/output) DOUBLE PRECISION array, dimension (N-1)
3329 On entry, DU must contain the (n-1) super-diagonal elements
3330 of A.
3331
3332 On exit, DU is overwritten by the (n-1) elements of the first
3333 super-diagonal of U.
3334
3335 B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
3336 On entry, the N by NRHS matrix of right hand side matrix B.
3337 On exit, if INFO = 0, the N by NRHS solution matrix X.
3338
3339 LDB (input) INTEGER
3340 The leading dimension of the array B. LDB >= max(1,N).
3341
3342 INFO (output) INTEGER
3343 = 0: successful exit
3344 < 0: if INFO = -i, the i-th argument had an illegal value
3345 > 0: if INFO = i, U(i,i) is exactly zero, and the solution
3346 has not been computed. The factorization has not been
3347 completed unless i = N.
3348
3349 =====================================================================
3350
3351
3352 Parameter adjustments */
3353 --dl;
3354 --d__;
3355 --du;
3356 b_dim1 = *ldb;
3357 b_offset = 1 + b_dim1 * 1;
3358 b -= b_offset;
3359
3360 /* Function Body */
3361 *info = 0;
3362 if(*n < 0) {
3363 *info = -1;
3364 } else if(*nrhs < 0) {
3365 *info = -2;
3366 } else if(*ldb < *n && *ldb < 1) {
3367 *info = -7;
3368 }
3369 if(*info != 0) {
3370 i__1 = -(*info);
3371 XERBLA("DGTSV ", i__1);
3372 /* XERBLA does not return */
3373 }
3374
3375 if(*n == 0) {
3376 return 0;
3377 }
3378
3379 if(*nrhs == 1) {
3380 i__1 = *n - 2;
3381 for(i__ = 1; i__ <= i__1; ++i__) {
3382 if(fabs(d__[i__]) >= fabs(dl[i__])) {
3383
3384/* No row interchange required */
3385
3386 if(d__[i__] != 0.) {
3387 fact = dl[i__] / d__[i__];
3388 d__[i__ + 1] -= fact * du[i__];
3389 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3390 1);
3391 } else {
3392 *info = i__;
3393 return 0;
3394 }
3395 dl[i__] = 0.;
3396 } else {
3397
3398/* Interchange rows I and I+1 */
3399
3400 fact = d__[i__] / dl[i__];
3401 d__[i__] = dl[i__];
3402 temp = d__[i__ + 1];
3403 d__[i__ + 1] = du[i__] - fact * temp;
3404 dl[i__] = du[i__ + 1];
3405 du[i__ + 1] = -fact * dl[i__];
3406 du[i__] = temp;
3407 temp = b_ref(i__, 1);
3408 b_ref(i__, 1) = b_ref(i__ + 1, 1);
3409 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3410 }
3411/* L10: */
3412 }
3413 if(*n > 1) {
3414 i__ = *n - 1;
3415 if(fabs(d__[i__]) >= fabs(dl[i__])) {
3416 if(d__[i__] != 0.) {
3417 fact = dl[i__] / d__[i__];
3418 d__[i__ + 1] -= fact * du[i__];
3419 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3420 1);
3421 } else {
3422 *info = i__;
3423 return 0;
3424 }
3425 } else {
3426 fact = d__[i__] / dl[i__];
3427 d__[i__] = dl[i__];
3428 temp = d__[i__ + 1];
3429 d__[i__ + 1] = du[i__] - fact * temp;
3430 du[i__] = temp;
3431 temp = b_ref(i__, 1);
3432 b_ref(i__, 1) = b_ref(i__ + 1, 1);
3433 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3434 }
3435 }
3436 if(d__[*n] == 0.) {
3437 *info = *n;
3438 return 0;
3439 }
3440 } else {
3441 i__1 = *n - 2;
3442 for(i__ = 1; i__ <= i__1; ++i__) {
3443 if(fabs(d__[i__]) >= fabs(dl[i__])) {
3444
3445/* No row interchange required */
3446
3447 if(d__[i__] != 0.) {
3448 fact = dl[i__] / d__[i__];
3449 d__[i__ + 1] -= fact * du[i__];
3450 i__2 = *nrhs;
3451 for(j = 1; j <= i__2; ++j) {
3452 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3453 i__, j);
3454/* L20: */
3455 }
3456 } else {
3457 *info = i__;
3458 return 0;
3459 }
3460 dl[i__] = 0.;
3461 } else {
3462
3463/* Interchange rows I and I+1 */
3464
3465 fact = d__[i__] / dl[i__];
3466 d__[i__] = dl[i__];
3467 temp = d__[i__ + 1];
3468 d__[i__ + 1] = du[i__] - fact * temp;
3469 dl[i__] = du[i__ + 1];
3470 du[i__ + 1] = -fact * dl[i__];
3471 du[i__] = temp;
3472 i__2 = *nrhs;
3473 for(j = 1; j <= i__2; ++j) {
3474 temp = b_ref(i__, j);
3475 b_ref(i__, j) = b_ref(i__ + 1, j);
3476 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3477/* L30: */
3478 }
3479 }
3480/* L40: */
3481 }
3482 if(*n > 1) {
3483 i__ = *n - 1;
3484 if( fabs(d__[i__]) >= fabs(dl[i__]))
3485 {
3486 if(d__[i__] != 0.)
3487 {
3488 fact = dl[i__] / d__[i__];
3489 d__[i__ + 1] -= fact * du[i__];
3490 i__1 = *nrhs;
3491 for(j = 1; j <= i__1; ++j) {
3492 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3493 i__, j);
3494 /* L50: */
3495 }
3496 }
3497 else
3498 {
3499 *info = i__;
3500 return 0;
3501 }
3502 } else {
3503 fact = d__[i__] / dl[i__];
3504 d__[i__] = dl[i__];
3505 temp = d__[i__ + 1];
3506 d__[i__ + 1] = du[i__] - fact * temp;
3507 du[i__] = temp;
3508 i__1 = *nrhs;
3509 for(j = 1; j <= i__1; ++j) {
3510 temp = b_ref(i__, j);
3511 b_ref(i__, j) = b_ref(i__ + 1, j);
3512 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3513/* L60: */
3514 }
3515 }
3516 }
3517 if(d__[*n] == 0.) {
3518 *info = *n;
3519 return 0;
3520 }
3521 }
3522
3523/* Back solve with the matrix U from the factorization. */
3524
3525 if(*nrhs <= 2) {
3526 j = 1;
3527L70:
3528 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3529 if(*n > 1) {
3530 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, j))
3531 / d__[*n - 1];
3532 }
3533 for(i__ = *n - 2; i__ >= 1; --i__) {
3534 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) - dl[
3535 i__] * b_ref(i__ + 2, j)) / d__[i__];
3536/* L80: */
3537 }
3538 if(j < *nrhs) {
3539 ++j;
3540 goto L70;
3541 }
3542 } else {
3543 i__1 = *nrhs;
3544 for(j = 1; j <= i__1; ++j) {
3545 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3546 if(*n > 1) {
3547 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n,
3548 j)) / d__[*n - 1];
3549 }
3550 for(i__ = *n - 2; i__ >= 1; --i__) {
3551 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j)
3552 - dl[i__] * b_ref(i__ + 2, j)) / d__[i__];
3553/* L90: */
3554 }
3555/* L100: */
3556 }
3557 }
3558
3559 return 0;
3560
3561/* End of DGTSV */
3562
3563} /* dgtsv_ */
3564#endif
3565#undef b_ref
3566#endif
3567/*lint +e725 expected pos indentation */
3568/*lint +e801 goto is deprecated */
STATIC double da(double z, double temp, double eden)
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
#define INT32_MAX
Definition cpu.h:49
#define J_(A_)
Definition iso.h:23
#define L_(A_)
Definition iso.h:21
static const double C1
static const int M
static const int N
#define CC(I_, J_)
STATIC void DGETRS(int32 TRANS, int32 N, int32 NRHS, double *A, int32 LDA, int32 IPIV[], double *B, int32 LDB, int32 *INFO)
#define AA(I_, J_)
#define ONE
STATIC void DTRSM(int32 SIDE, int32 UPLO, int32 TRANSA, int32 DIAG, int32 M, int32 N, double ALPHA, double *A, int32 LDA, double *B, int32 LDB)
STATIC void XERBLA(const char *SRNAME, int32 INFO)
STATIC void DSCAL(int32 n, double da, double dx[], int32 incx)
STATIC int32 LSAME(int32 CA, int32 CB)
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
#define BB(I_, J_)
STATIC int32 ILAENV(int32 ISPEC, const char *NAME, int32 N1, int32 N2, int32 N4)
STATIC int32 IDAMAX(int32 n, double dx[], int32 incx)
STATIC void DGEMM(int32 TRANSA, int32 TRANSB, int32 M, int32 N, int32 K, double ALPHA, double *A, int32 LDA, double *B, int32 LDB, double BETA, double *C, int32 LDC)
STATIC void DLASWP(int32 N, double *A, int32 LDA, int32 K1, int32 K2, int32 IPIV[], int32 INCX)
STATIC void DGER(int32 M, int32 N, double ALPHA, double X[], int32 INCX, double Y[], int32 INCY, double *A, int32 LDA)
#define ZERO
STATIC void DGETRF(int32, int32, double *, int32, int32[], int32 *)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
STATIC void DGETF2(int32 M, int32 N, double *A, int32 LDA, int32 IPIV[], int32 *INFO)
STATIC void DSWAP(int32 n, double dx[], int32 incx, double dy[], int32 incy)