cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
20 extern "C" void dgetrf_(int32 *M, int32 *N, double *A, int32 *LDA, int32 *IPIV, int32 *INFO);
21 extern "C" void dgetrs_(char *TRANS, int32 *N, int32 *NRHS, double *A, int32 *LDA, int32 *iPiv, double *B,
22  int32 *LDB, int32 *INFO, int32 translen);
23 extern "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 */
34 STATIC void DGETRF(int32,int32,double*,int32,int32[],int32*);
35 
36 /* DGETRS lapack, solve linear system */
37 STATIC 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 */
47 void 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 
69 void 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
94 void dgtsv_wrapper(long N, long nrhs, double *dl, double *d__, double *du, double *b, long ldb, int32 *info)
95 {
96  printf("Inside dgtsv\n");
97  cdEXIT( EXIT_FAILURE );
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*/
164 STATIC 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 */
179 STATIC int32 LSAME(int32 CA,
180  int32 CB);
181 
182 /*IDAMAX lapack service routine */
183 STATIC int32 IDAMAX(int32 n,
184  double dx[],
185  int32 incx);
186 
187 /*DTRSM lapack service routine */
188 STATIC 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 */
201 STATIC 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 */
210 STATIC void DSWAP(int32 n,
211  double dx[],
212  int32 incx,
213  double dy[],
214  int32 incy);
215 
216 /*DSCAL lapack routine */
217 STATIC void DSCAL(int32 n,
218  double da,
219  double dx[],
220  int32 incx);
221 
222 /*DLASWP -- LAPACK auxiliary routine (version 2.0) --*/
223 STATIC 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 */
232 STATIC 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 */
240 STATIC 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) -- */
251 STATIC 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 
292  cdEXIT(EXIT_FAILURE);
293 }
294 
295 
296 STATIC void DGETRF(
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 
527 STATIC void DGETRS(
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 
725 STATIC 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 
834 STATIC 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  * */
888 L_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 */
903 STATIC 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  TEMP = -BB(J_,K_);
1166  for( I=0; I < K-1; I++ )
1167  {
1168  BB(J_,I) += TEMP*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  TEMP = -BB(J_,K_);
1195  for( I=K; I < M; I++ )
1196  {
1197  BB(J_,I) += TEMP*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 */
1424 STATIC 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 
1594 L_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 
1686 L_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 
1946 L_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 
2073 L_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 
2158 L_400:
2159 
2160  /* ISPEC = 4: number of shifts (used by xHSEQR)
2161  * */
2162  ILAENV_v = 6;
2163  return ILAENV_v;
2164 
2165 L_500:
2166 
2167  /* ISPEC = 5: minimum column dimension (not used)
2168  * */
2169  ILAENV_v = 2;
2170  return ILAENV_v;
2171 
2172 L_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 
2179 L_700:
2180 
2181  /* ISPEC = 7: number of processors (not used)
2182  * */
2183  ILAENV_v = 1;
2184  return ILAENV_v;
2185 
2186 L_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 
2199 STATIC 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  * */
2252 L_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 
2269 L_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 
2287 STATIC void DSCAL(int32 n,
2288  double da,
2289  double dx[],
2290  int32 incx)
2291 {
2292  int32 i,
2293  nincx;
2294 
2295  DEBUG_ENTRY( "DSCAL()" );
2296 
2297  /* scales a vector by a constant.
2298  * uses unrolled loops for increment equal to one.
2299  * jack dongarra, lapack, 3/11/78.
2300  * modified 3/93 to return if incx .le. 0.
2301  * modified 12/3/93, array(1) declarations changed to array(*)
2302  * */
2303 
2304  if( n <= 0 || incx <= 0 )
2305  {
2306  return;}
2307  if( incx == 1 )
2308  goto L_20;
2309 
2310  /* code for increment not equal to 1
2311  * */
2312  nincx = n*incx;
2313  /*for( i=1, _do0=DOCNT(i,nincx,_do1 = incx); _do0 > 0; i += _do1, _do0-- )*/
2314  for( i=0; i<nincx; i = i + incx)
2315  {
2316  dx[i] *= da;
2317  }
2318  return;
2319 
2320  /* code for increment equal to 1
2321  *
2322  *
2323  * clean-up loop
2324  * */
2325 L_20:
2326  for( i=0; i < n; i++ )
2327  {
2328  dx[i] *= da;
2329  }
2330  return;
2331 }
2332 
2333 /*DLASWP -- LAPACK auxiliary routine (version 2.0) --*/
2334 
2335 STATIC void DLASWP(int32 N,
2336  double *A,
2337  int32 LDA,
2338  int32 K1,
2339  int32 K2,
2340  int32 IPIV[],
2341  int32 INCX)
2342 {
2343  int32 I,
2344  IP,
2345  IX,
2346  I_;
2347 
2348  DEBUG_ENTRY( "DLASWP()" );
2349 
2350  /* -- LAPACK auxiliary routine (version 2.0) --
2351  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2352  * Courant Institute, Argonne National Lab, and Rice University
2353  * October 31, 1992
2354  *
2355  * .. Scalar Arguments .. */
2356  /* ..
2357  * .. Array Arguments .. */
2358  /* ..
2359  *
2360  * Purpose
2361  * =======
2362  *
2363  * DLASWP performs a series of row interchanges on the matrix A.
2364  * One row interchange is initiated for each of rows K1 through K2 of A.
2365  *
2366  * Arguments
2367  * =========
2368  *
2369  * N (input) INTEGER
2370  * The number of columns of the matrix A.
2371  *
2372  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
2373  * On entry, the matrix of column dimension N to which the row
2374  * interchanges will be applied.
2375  * On exit, the permuted matrix.
2376  *
2377  * LDA (input) INTEGER
2378  * The leading dimension of the array A.
2379  *
2380  * K1 (input) INTEGER
2381  * The first element of IPIV for which a row interchange will
2382  * be done.
2383  *
2384  * K2 (input) INTEGER
2385  * The last element of IPIV for which a row interchange will
2386  * be done.
2387  *
2388  * IPIV (input) INTEGER array, dimension (M*ABS(INCX))
2389  * The vector of pivot indices. Only the elements in positions
2390  * K1 through K2 of IPIV are accessed.
2391  * IPIV(K) = L implies rows K and L are to be interchanged.
2392  *
2393  * INCX (input) INTEGER
2394  * The increment between successive values of IPIV. If IPIV
2395  * is negative, the pivots are applied in reverse order.
2396  *
2397  * =====================================================================
2398  *
2399  * .. Local Scalars .. */
2400  /* ..
2401  * .. External Subroutines .. */
2402  /* ..
2403  * .. Executable Statements ..
2404  *
2405  * Interchange row I with row IPIV(I) for each of rows K1 through K2.
2406  * */
2407  if( INCX == 0 )
2408  {
2409  return;}
2410  if( INCX > 0 )
2411  {
2412  IX = K1;
2413  }
2414  else
2415  {
2416  IX = 1 + (1 - K2)*INCX;
2417  }
2418  if( INCX == 1 )
2419  {
2420  for( I=K1; I <= K2; I++ )
2421  {
2422  I_ = I - 1;
2423  IP = IPIV[I_];
2424  if( IP != I )
2425  DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2426  }
2427  }
2428  else if( INCX > 1 )
2429  {
2430  for( I=K1; I <= K2; I++ )
2431  {
2432  I_ = I - 1;
2433  IP = IPIV[IX-1];
2434  if( IP != I )
2435  DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2436  IX += INCX;
2437  }
2438  }
2439  else if( INCX < 0 )
2440  {
2441  for( I=K2; I >= K1; I-- )
2442  {
2443  I_ = I - 1;
2444  IP = IPIV[IX-1];
2445  if( IP != I )
2446  DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2447  IX += INCX;
2448  }
2449  }
2450 
2451  return;
2452 
2453  /* End of DLASWP
2454  * */
2455 #undef A
2456 }
2457 
2458 /*DGETF2 lapack service routine */
2459 
2460 STATIC void DGETF2(int32 M,
2461  int32 N,
2462  double *A,
2463  int32 LDA,
2464  int32 IPIV[],
2465  int32 *INFO)
2466 {
2467  int32 J,
2468  JP,
2469  J_,
2470  limit;
2471 
2472  DEBUG_ENTRY( "DGETF2()" );
2473 
2474  /* -- LAPACK routine (version 2.0) --
2475  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2476  * Courant Institute, Argonne National Lab, and Rice University
2477  * June 30, 1992
2478  *
2479  * .. Scalar Arguments .. */
2480  /* ..
2481  * .. Array Arguments .. */
2482  /* ..
2483  *
2484  * Purpose
2485  * =======
2486  *
2487  * DGETF2 computes an LU factorization of a general m-by-n matrix A
2488  * using partial pivoting with row interchanges.
2489  *
2490  * The factorization has the form
2491  * A = P * L * U
2492  * where P is a permutation matrix, L is lower triangular with unit
2493  * diagonal elements (lower trapezoidal if m > n), and U is upper
2494  * triangular (upper trapezoidal if m < n).
2495  *
2496  * This is the right-looking Level 2 BLAS version of the algorithm.
2497  *
2498  * Arguments
2499  * =========
2500  *
2501  * M (input) INTEGER
2502  * The number of rows of the matrix A. M >= 0.
2503  *
2504  * N (input) INTEGER
2505  * The number of columns of the matrix A. N >= 0.
2506  *
2507  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
2508  * On entry, the m by n matrix to be factored.
2509  * On exit, the factors L and U from the factorization
2510  * A = P*L*U; the unit diagonal elements of L are not stored.
2511  *
2512  * LDA (input) INTEGER
2513  * The leading dimension of the array A. LDA >= MAX(1,M).
2514  *
2515  * IPIV (output) INTEGER array, dimension (MIN(M,N))
2516  * The pivot indices; for 1 <= i <= MIN(M,N), row i of the
2517  * matrix was interchanged with row IPIV(i).
2518  *
2519  * INFO (output) INTEGER
2520  * = 0: successful exit
2521  * < 0: if INFO = -k, the k-th argument had an illegal value
2522  * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
2523  * has been completed, but the factor U is exactly
2524  * singular, and division by zero will occur if it is used
2525  * to solve a system of equations.
2526  *
2527  * =====================================================================
2528  *
2529  * .. Parameters .. */
2530  /* ..
2531  * .. Local Scalars .. */
2532  /* ..
2533  * .. External Functions .. */
2534  /* ..
2535  * .. External Subroutines .. */
2536  /* ..
2537  * .. Intrinsic Functions .. */
2538  /* ..
2539  * .. Executable Statements ..
2540  *
2541  * Test the input parameters.
2542  * */
2543  *INFO = 0;
2544  if( M < 0 )
2545  {
2546  *INFO = -1;
2547  }
2548  else if( N < 0 )
2549  {
2550  *INFO = -2;
2551  }
2552  else if( LDA < MAX2(1,M) )
2553  {
2554  *INFO = -4;
2555  }
2556  if( *INFO != 0 )
2557  {
2558  XERBLA("DGETF2",-*INFO);
2559  /* XERBLA does not return */
2560  }
2561 
2562  /* Quick return if possible
2563  * */
2564  if( M == 0 || N == 0 )
2565  {
2566  return;}
2567 
2568  limit = MIN2(M,N);
2569  for( J=1; J <= limit; J++ )
2570  {
2571  J_ = J - 1;
2572 
2573  /* Find pivot and test for singularity.
2574  * */
2575  JP = J - 1 + IDAMAX(M-J+1,&AA(J_,J_),1);
2576  IPIV[J_] = JP;
2577  if( AA(J_,JP-1) != ZERO )
2578  {
2579  /* Apply the interchange to columns 1:N.
2580  * */
2581  if( JP != J )
2582  DSWAP(N,&AA(0,J_),LDA,&AA(0,JP-1),LDA);
2583 
2584  /* Compute elements J+1:M of J-th column.
2585  * */
2586  if( J < M )
2587  DSCAL(M-J,ONE/AA(J_,J_),&AA(J_,J_+1),1);
2588  }
2589  else if( *INFO == 0 )
2590  {
2591  *INFO = J;
2592  }
2593 
2594  if( J < MIN2(M,N) )
2595  {
2596  /* Update trailing submatrix.
2597  * */
2598  DGER(M-J,N-J,-ONE,&AA(J_,J_+1),1,&AA(J_+1,J_),LDA,&AA(J_+1,J_+1),
2599  LDA);
2600  }
2601  }
2602  return;
2603 
2604  /* End of DGETF2
2605  * */
2606 #undef A
2607 }
2608 
2609 /*DGER service routine for matrix inversion */
2610 
2611 STATIC void DGER(int32 M,
2612  int32 N,
2613  double ALPHA,
2614  double X[],
2615  int32 INCX,
2616  double Y[],
2617  int32 INCY,
2618  double *A,
2619  int32 LDA)
2620 {
2621  int32 I,
2622  INFO,
2623  IX,
2624  I_,
2625  J,
2626  JY,
2627  J_,
2628  KX;
2629  double TEMP;
2630 
2631  DEBUG_ENTRY( "DGER()" );
2632  /* .. Scalar Arguments .. */
2633  /* .. Array Arguments .. */
2634  /* ..
2635  *
2636  * Purpose
2637  * =======
2638  *
2639  * DGER performs the rank 1 operation
2640  *
2641  * A := alpha*x*y' + A,
2642  *
2643  * where alpha is a scalar, x is an m element vector, y is an n element
2644  * vector and A is an m by n matrix.
2645  *
2646  * Parameters
2647  * ==========
2648  *
2649  * M - INTEGER.
2650  * On entry, M specifies the number of rows of the matrix A.
2651  * M must be at least zero.
2652  * Unchanged on exit.
2653  *
2654  * N - INTEGER.
2655  * On entry, N specifies the number of columns of the matrix A.
2656  * N must be at least zero.
2657  * Unchanged on exit.
2658  *
2659  * ALPHA - DOUBLE PRECISION.
2660  * On entry, ALPHA specifies the scalar alpha.
2661  * Unchanged on exit.
2662  *
2663  * X - DOUBLE PRECISION array of dimension at least
2664  * ( 1 + ( m - 1 )*ABS( INCX ) ).
2665  * Before entry, the incremented array X must contain the m
2666  * element vector x.
2667  * Unchanged on exit.
2668  *
2669  * INCX - INTEGER.
2670  * On entry, INCX specifies the increment for the elements of
2671  * X. INCX must not be zero.
2672  * Unchanged on exit.
2673  *
2674  * Y - DOUBLE PRECISION array of dimension at least
2675  * ( 1 + ( n - 1 )*ABS( INCY ) ).
2676  * Before entry, the incremented array Y must contain the n
2677  * element vector y.
2678  * Unchanged on exit.
2679  *
2680  * INCY - INTEGER.
2681  * On entry, INCY specifies the increment for the elements of
2682  * Y. INCY must not be zero.
2683  * Unchanged on exit.
2684  *
2685  * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2686  * Before entry, the leading m by n part of the array A must
2687  * contain the matrix of coefficients. On exit, A is
2688  * overwritten by the updated matrix.
2689  *
2690  * LDA - INTEGER.
2691  * On entry, LDA specifies the first dimension of A as declared
2692  * in the calling (sub) program. LDA must be at least
2693  * MAX( 1, m ).
2694  * Unchanged on exit.
2695  *
2696  *
2697  * Level 2 Blas routine.
2698  *
2699  * -- Written on 22-October-1986.
2700  * Jack Dongarra, Argonne National Lab.
2701  * Jeremy Du Croz, Nag Central Office.
2702  * Sven Hammarling, Nag Central Office.
2703  * Richard Hanson, Sandia National Labs.
2704  *
2705  *
2706  * .. Parameters .. */
2707  /* .. Local Scalars .. */
2708  /* .. External Subroutines .. */
2709  /* .. Intrinsic Functions .. */
2710  /* ..
2711  * .. Executable Statements ..
2712  *
2713  * Test the input parameters.
2714  * */
2715  INFO = 0;
2716  if( M < 0 )
2717  {
2718  INFO = 1;
2719  }
2720  else if( N < 0 )
2721  {
2722  INFO = 2;
2723  }
2724  else if( INCX == 0 )
2725  {
2726  INFO = 5;
2727  }
2728  else if( INCY == 0 )
2729  {
2730  INFO = 7;
2731  }
2732  else if( LDA < MAX2(1,M) )
2733  {
2734  INFO = 9;
2735  }
2736  if( INFO != 0 )
2737  {
2738  XERBLA("DGER ",INFO);
2739  /* XERBLA does not return */
2740  }
2741 
2742  /* Quick return if possible.
2743  * */
2744  if( ((M == 0) || (N == 0)) || (ALPHA == ZERO) )
2745  {
2746  return;}
2747 
2748  /* Start the operations. In this version the elements of A are
2749  * accessed sequentially with one pass through A.
2750  * */
2751  if( INCY > 0 )
2752  {
2753  JY = 1;
2754  }
2755  else
2756  {
2757  JY = 1 - (N - 1)*INCY;
2758  }
2759  if( INCX == 1 )
2760  {
2761  for( J=1; J <= N; J++ )
2762  {
2763  J_ = J - 1;
2764  if( Y[JY-1] != ZERO )
2765  {
2766  TEMP = ALPHA*Y[JY-1];
2767  for( I=0; I < M; I++ )
2768  {
2769  AA(J_,I) += X[I]*TEMP;
2770  }
2771  }
2772  JY += INCY;
2773  }
2774  }
2775  else
2776  {
2777  if( INCX > 0 )
2778  {
2779  KX = 1;
2780  }
2781  else
2782  {
2783  KX = 1 - (M - 1)*INCX;
2784  }
2785  for( J=1; J <= N; J++ )
2786  {
2787  J_ = J - 1;
2788  if( Y[JY-1] != ZERO )
2789  {
2790  TEMP = ALPHA*Y[JY-1];
2791  IX = KX;
2792  for( I=1; I <= M; I++ )
2793  {
2794  I_ = I - 1;
2795  AA(J_,I_) += X[IX-1]*TEMP;
2796  IX += INCX;
2797  }
2798  }
2799  JY += INCY;
2800  }
2801  }
2802 
2803  return;
2804 
2805  /* End of DGER .
2806  * */
2807 #undef A
2808 }
2809 
2810 /* DGEMM matrix inversion helper routine*/
2811 
2812 STATIC void DGEMM(int32 TRANSA,
2813  int32 TRANSB,
2814  int32 M,
2815  int32 N,
2816  int32 K,
2817  double ALPHA,
2818  double *A,
2819  int32 LDA,
2820  double *B,
2821  int32 LDB,
2822  double BETA,
2823  double *C,
2824  int32 LDC)
2825 {
2826  int32 NOTA,
2827  NOTB;
2828  int32 I,
2829  INFO,
2830  J,
2831  L,
2832  NROWA,
2833  NROWB;
2834  double TEMP;
2835 
2836  DEBUG_ENTRY( "DGEMM()" );
2837  /* .. Scalar Arguments .. */
2838  /* .. Array Arguments .. */
2839  /* ..
2840  *
2841  * Purpose
2842  * =======
2843  *
2844  * DGEMM performs one of the matrix-matrix operations
2845  *
2846  * C := alpha*op( A )*op( B ) + beta*C,
2847  *
2848  * where op( X ) is one of
2849  *
2850  * op( X ) = X or op( X ) = X',
2851  *
2852  * alpha and beta are scalars, and A, B and C are matrices, with op( A )
2853  * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
2854  *
2855  * Parameters
2856  * ==========
2857  *
2858  * TRANSA - CHARACTER*1.
2859  * On entry, TRANSA specifies the form of op( A ) to be used in
2860  * the matrix multiplication as follows:
2861  *
2862  * TRANSA = 'N' or 'n', op( A ) = A.
2863  *
2864  * TRANSA = 'T' or 't', op( A ) = A'.
2865  *
2866  * TRANSA = 'C' or 'c', op( A ) = A'.
2867  *
2868  * Unchanged on exit.
2869  *
2870  * TRANSB - CHARACTER*1.
2871  * On entry, TRANSB specifies the form of op( B ) to be used in
2872  * the matrix multiplication as follows:
2873  *
2874  * TRANSB = 'N' or 'n', op( B ) = B.
2875  *
2876  * TRANSB = 'T' or 't', op( B ) = B'.
2877  *
2878  * TRANSB = 'C' or 'c', op( B ) = B'.
2879  *
2880  * Unchanged on exit.
2881  *
2882  * M - INTEGER.
2883  * On entry, M specifies the number of rows of the matrix
2884  * op( A ) and of the matrix C. M must be at least zero.
2885  * Unchanged on exit.
2886  *
2887  * N - INTEGER.
2888  * On entry, N specifies the number of columns of the matrix
2889  * op( B ) and the number of columns of the matrix C. N must be
2890  * at least zero.
2891  * Unchanged on exit.
2892  *
2893  * K - INTEGER.
2894  * On entry, K specifies the number of columns of the matrix
2895  * op( A ) and the number of rows of the matrix op( B ). K must
2896  * be at least zero.
2897  * Unchanged on exit.
2898  *
2899  * ALPHA - DOUBLE PRECISION.
2900  * On entry, ALPHA specifies the scalar alpha.
2901  * Unchanged on exit.
2902  *
2903  * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
2904  * k when TRANSA = 'N' or 'n', and is m otherwise.
2905  * Before entry with TRANSA = 'N' or 'n', the leading m by k
2906  * part of the array A must contain the matrix A, otherwise
2907  * the leading k by m part of the array A must contain the
2908  * matrix A.
2909  * Unchanged on exit.
2910  *
2911  * LDA - INTEGER.
2912  * On entry, LDA specifies the first dimension of A as declared
2913  * in the calling (sub) program. When TRANSA = 'N' or 'n' then
2914  * LDA must be at least MAX( 1, m ), otherwise LDA must be at
2915  * least MAX( 1, k ).
2916  * Unchanged on exit.
2917  *
2918  * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
2919  * n when TRANSB = 'N' or 'n', and is k otherwise.
2920  * Before entry with TRANSB = 'N' or 'n', the leading k by n
2921  * part of the array B must contain the matrix B, otherwise
2922  * the leading n by k part of the array B must contain the
2923  * matrix B.
2924  * Unchanged on exit.
2925  *
2926  * LDB - INTEGER.
2927  * On entry, LDB specifies the first dimension of B as declared
2928  * in the calling (sub) program. When TRANSB = 'N' or 'n' then
2929  * LDB must be at least MAX( 1, k ), otherwise LDB must be at
2930  * least MAX( 1, n ).
2931  * Unchanged on exit.
2932  *
2933  * BETA - DOUBLE PRECISION.
2934  * On entry, BETA specifies the scalar beta. When BETA is
2935  * supplied as zero then C need not be set on input.
2936  * Unchanged on exit.
2937  *
2938  * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
2939  * Before entry, the leading m by n part of the array C must
2940  * contain the matrix C, except when beta is zero, in which
2941  * case C need not be set on entry.
2942  * On exit, the array C is overwritten by the m by n matrix
2943  * ( alpha*op( A )*op( B ) + beta*C ).
2944  *
2945  * LDC - INTEGER.
2946  * On entry, LDC specifies the first dimension of C as declared
2947  * in the calling (sub) program. LDC must be at least
2948  * MAX( 1, m ).
2949  * Unchanged on exit.
2950  *
2951  *
2952  * Level 3 Blas routine.
2953  *
2954  * -- Written on 8-February-1989.
2955  * Jack Dongarra, Argonne National Laboratory.
2956  * Iain Duff, AERE Harwell.
2957  * Jeremy Du Croz, Numerical Algorithms Group Ltd.
2958  * Sven Hammarling, Numerical Algorithms Group Ltd.
2959  *
2960  *
2961  * .. External Functions .. */
2962  /* .. External Subroutines .. */
2963  /* .. Intrinsic Functions .. */
2964  /* .. Local Scalars .. */
2965  /* .. Parameters .. */
2966  /* ..
2967  * .. Executable Statements ..
2968  *
2969  * Set NOTA and NOTB as true if A and B respectively are not
2970  * transposed and set NROWA, NROWB as the number of rows
2971  * and columns of A and the number of rows of B respectively.
2972  * */
2973  NOTA = LSAME(TRANSA,'N');
2974  NOTB = LSAME(TRANSB,'N');
2975 
2976  if( NOTA )
2977  {
2978  NROWA = M;
2979  }
2980  else
2981  {
2982  NROWA = K;
2983  }
2984 
2985  if( NOTB )
2986  {
2987  NROWB = K;
2988  }
2989  else
2990  {
2991  NROWB = N;
2992  }
2993 
2994  /* Test the input parameters.
2995  * */
2996  INFO = 0;
2997  if( ((!NOTA) && (!LSAME(TRANSA,'C'))) && (!LSAME(TRANSA,'T')) )
2998  {
2999  INFO = 1;
3000  }
3001  else if(
3002  ((!NOTB) && (!LSAME(TRANSB,'C'))) && (!LSAME(TRANSB,'T')) )
3003  {
3004  INFO = 2;
3005  }
3006 
3007  else if( M < 0 )
3008  {
3009  INFO = 3;
3010  }
3011 
3012  else if( N < 0 )
3013  {
3014  INFO = 4;
3015  }
3016 
3017  else if( K < 0 )
3018  {
3019  INFO = 5;
3020  }
3021 
3022  else if( LDA < MAX2(1,NROWA) )
3023  {
3024  INFO = 8;
3025  }
3026 
3027  else if( LDB < MAX2(1,NROWB) )
3028  {
3029  INFO = 10;
3030  }
3031 
3032  else if( LDC < MAX2(1,M) )
3033  {
3034  INFO = 13;
3035  }
3036 
3037  if( INFO != 0 )
3038  {
3039  XERBLA("DGEMM ",INFO);
3040  /* XERBLA does not return */
3041  }
3042 
3043  /* Quick return if possible.
3044  * */
3045  if( ((M == 0) || (N == 0)) || (((ALPHA == ZERO) || (K == 0)) &&
3046  (BETA == ONE)) )
3047  {
3048  return;
3049  }
3050 
3051  /* And if alpha.eq.zero. */
3052  if( ALPHA == ZERO )
3053  {
3054  if( BETA == ZERO )
3055  {
3056  for( J=0; J < N; J++ )
3057  {
3058  for( I=0; I < M; I++ )
3059  {
3060  CC(J,I) = ZERO;
3061  }
3062  }
3063  }
3064 
3065  else
3066  {
3067  for( J=0; J < N; J++ )
3068  {
3069  for( I=0; I < M; I++ )
3070  {
3071  CC(J,I) *= BETA;
3072  }
3073  }
3074  }
3075  return;
3076  }
3077 
3078  /* Start the operations.*/
3079  if( NOTB )
3080  {
3081 
3082  if( NOTA )
3083  {
3084 
3085  /* Form C := alpha*A*B + beta*C.
3086  * */
3087  for( J=0; J < N; J++ )
3088  {
3089  if( BETA == ZERO )
3090  {
3091  for( I=0; I < M; I++ )
3092  {
3093  CC(J,I) = ZERO;
3094  }
3095  }
3096 
3097  else if( BETA != ONE )
3098  {
3099  for( I=0; I < M; I++ )
3100  {
3101  CC(J,I) *= BETA;
3102  }
3103  }
3104 
3105  for( L=0; L < K; L++ )
3106  {
3107  if( BB(J,L) != ZERO )
3108  {
3109  TEMP = ALPHA*BB(J,L);
3110  for( I=0; I < M; I++ )
3111  {
3112  CC(J,I) += TEMP*AA(L,I);
3113  }
3114  }
3115  }
3116  }
3117  }
3118  else
3119  {
3120 
3121  /* Form C := alpha*A'*B + beta*C */
3122  for( J=0; J < N; J++ )
3123  {
3124  for( I=0; I < M; I++ )
3125  {
3126  TEMP = ZERO;
3127  for( L=0; L < K; L++ )
3128  {
3129  TEMP += AA(I,L)*BB(J,L);
3130  }
3131 
3132  if( BETA == ZERO )
3133  {
3134  CC(J,I) = ALPHA*TEMP;
3135  }
3136  else
3137  {
3138  CC(J,I) = ALPHA*TEMP + BETA*CC(J,I);
3139  }
3140  }
3141  }
3142  }
3143  }
3144  else
3145  {
3146  if( NOTA )
3147  {
3148 
3149  /* Form C := alpha*A*B' + beta*C
3150  * */
3151  for( J=0; J < N; J++ )
3152  {
3153  if( BETA == ZERO )
3154  {
3155  for( I=0; I < M; I++ )
3156  {
3157  CC(J,I) = ZERO;
3158  }
3159  }
3160 
3161  else if( BETA != ONE )
3162  {
3163  for( I=0; I < M; I++ )
3164  {
3165  CC(J,I) *= BETA;
3166  }
3167  }
3168 
3169  for( L=0; L < K; L++ )
3170  {
3171  if( BB(L,J) != ZERO )
3172  {
3173  TEMP = ALPHA*BB(L,J);
3174  for( I=0; I < M; I++ )
3175  {
3176  CC(J,I) += TEMP*AA(L,I);
3177  }
3178  }
3179  }
3180  }
3181  }
3182 
3183  else
3184  {
3185 
3186  /* Form C := alpha*A'*B' + beta*C */
3187  for( J=0; J < N; J++ )
3188  {
3189 
3190  for( I=0; I < M; I++ )
3191  {
3192  TEMP = ZERO;
3193 
3194  for( L=0; L < K; L++ )
3195  {
3196  TEMP += AA(I,L)*BB(L,J);
3197  }
3198 
3199  if( BETA == ZERO )
3200  {
3201  CC(J,I) = ALPHA*TEMP;
3202  }
3203  else
3204  {
3205  CC(J,I) = ALPHA*TEMP + BETA*CC(J,I);
3206  }
3207 
3208  }
3209  }
3210  }
3211  }
3212 
3213  return;
3214 
3215  /* End of DGEMM .*/
3216 #undef C
3217 #undef B
3218 #undef A
3219 }
3220 #undef AA
3221 #undef BB
3222 #undef CC
3223 
3224 /* Subroutine */
3225 #if 0
3226 STATIC int32 DGTSV(int32 *n, int32 *nrhs, double *dl,
3227  double *d__, double *du, double *b, int32 *ldb, int32
3228  *info)
3229 {
3230  /* System generated locals */
3231  int32 b_dim1, b_offset, i__1, i__2;
3232 
3233  /* Local variables */
3234  double fact, temp;
3235  int32 i__, j;
3236 
3237 
3238 #define b_ref(a_1,a_2) b[(a_2)*(b_dim1) + (a_1)]
3239 
3240 
3241 /* -- LAPACK routine (version 3.0) --
3242  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3243  Courant Institute, Argonne National Lab, and Rice University
3244  October 31, 1999
3245 
3246 
3247  Purpose
3248  =======
3249 
3250  DGTSV solves the equation
3251 
3252  A*X = B,
3253 
3254  where A is an n by n tridiagonal matrix, by Gaussian elimination with
3255  partial pivoting.
3256 
3257  Note that the equation A'*X = B may be solved by interchanging the
3258  order of the arguments DU and DL.
3259 
3260  Arguments
3261  =========
3262 
3263  N (input) INTEGER
3264  The order of the matrix A. N >= 0.
3265 
3266  NRHS (input) INTEGER
3267  The number of right hand sides, i.e., the number of columns
3268  of the matrix B. NRHS >= 0.
3269 
3270  DL (input/output) DOUBLE PRECISION array, dimension (N-1)
3271  On entry, DL must contain the (n-1) sub-diagonal elements of
3272  A.
3273 
3274  On exit, DL is overwritten by the (n-2) elements of the
3275  second super-diagonal of the upper triangular matrix U from
3276  the LU factorization of A, in DL(1), ..., DL(n-2).
3277 
3278  D (input/output) DOUBLE PRECISION array, dimension (N)
3279  On entry, D must contain the diagonal elements of A.
3280 
3281  On exit, D is overwritten by the n diagonal elements of U.
3282 
3283  DU (input/output) DOUBLE PRECISION array, dimension (N-1)
3284  On entry, DU must contain the (n-1) super-diagonal elements
3285  of A.
3286 
3287  On exit, DU is overwritten by the (n-1) elements of the first
3288  super-diagonal of U.
3289 
3290  B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
3291  On entry, the N by NRHS matrix of right hand side matrix B.
3292  On exit, if INFO = 0, the N by NRHS solution matrix X.
3293 
3294  LDB (input) INTEGER
3295  The leading dimension of the array B. LDB >= max(1,N).
3296 
3297  INFO (output) INTEGER
3298  = 0: successful exit
3299  < 0: if INFO = -i, the i-th argument had an illegal value
3300  > 0: if INFO = i, U(i,i) is exactly zero, and the solution
3301  has not been computed. The factorization has not been
3302  completed unless i = N.
3303 
3304  =====================================================================
3305 
3306 
3307  Parameter adjustments */
3308  --dl;
3309  --d__;
3310  --du;
3311  b_dim1 = *ldb;
3312  b_offset = 1 + b_dim1 * 1;
3313  b -= b_offset;
3314 
3315  /* Function Body */
3316  *info = 0;
3317  if(*n < 0) {
3318  *info = -1;
3319  } else if(*nrhs < 0) {
3320  *info = -2;
3321  } else if(*ldb < *n && *ldb < 1) {
3322  *info = -7;
3323  }
3324  if(*info != 0) {
3325  i__1 = -(*info);
3326  XERBLA("DGTSV ", i__1);
3327  /* XERBLA does not return */
3328  }
3329 
3330  if(*n == 0) {
3331  return 0;
3332  }
3333 
3334  if(*nrhs == 1) {
3335  i__1 = *n - 2;
3336  for(i__ = 1; i__ <= i__1; ++i__) {
3337  if(fabs(d__[i__]) >= fabs(dl[i__])) {
3338 
3339 /* No row interchange required */
3340 
3341  if(d__[i__] != 0.) {
3342  fact = dl[i__] / d__[i__];
3343  d__[i__ + 1] -= fact * du[i__];
3344  b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3345  1);
3346  } else {
3347  *info = i__;
3348  return 0;
3349  }
3350  dl[i__] = 0.;
3351  } else {
3352 
3353 /* Interchange rows I and I+1 */
3354 
3355  fact = d__[i__] / dl[i__];
3356  d__[i__] = dl[i__];
3357  temp = d__[i__ + 1];
3358  d__[i__ + 1] = du[i__] - fact * temp;
3359  dl[i__] = du[i__ + 1];
3360  du[i__ + 1] = -fact * dl[i__];
3361  du[i__] = temp;
3362  temp = b_ref(i__, 1);
3363  b_ref(i__, 1) = b_ref(i__ + 1, 1);
3364  b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3365  }
3366 /* L10: */
3367  }
3368  if(*n > 1) {
3369  i__ = *n - 1;
3370  if(fabs(d__[i__]) >= fabs(dl[i__])) {
3371  if(d__[i__] != 0.) {
3372  fact = dl[i__] / d__[i__];
3373  d__[i__ + 1] -= fact * du[i__];
3374  b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3375  1);
3376  } else {
3377  *info = i__;
3378  return 0;
3379  }
3380  } else {
3381  fact = d__[i__] / dl[i__];
3382  d__[i__] = dl[i__];
3383  temp = d__[i__ + 1];
3384  d__[i__ + 1] = du[i__] - fact * temp;
3385  du[i__] = temp;
3386  temp = b_ref(i__, 1);
3387  b_ref(i__, 1) = b_ref(i__ + 1, 1);
3388  b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3389  }
3390  }
3391  if(d__[*n] == 0.) {
3392  *info = *n;
3393  return 0;
3394  }
3395  } else {
3396  i__1 = *n - 2;
3397  for(i__ = 1; i__ <= i__1; ++i__) {
3398  if(fabs(d__[i__]) >= fabs(dl[i__])) {
3399 
3400 /* No row interchange required */
3401 
3402  if(d__[i__] != 0.) {
3403  fact = dl[i__] / d__[i__];
3404  d__[i__ + 1] -= fact * du[i__];
3405  i__2 = *nrhs;
3406  for(j = 1; j <= i__2; ++j) {
3407  b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3408  i__, j);
3409 /* L20: */
3410  }
3411  } else {
3412  *info = i__;
3413  return 0;
3414  }
3415  dl[i__] = 0.;
3416  } else {
3417 
3418 /* Interchange rows I and I+1 */
3419 
3420  fact = d__[i__] / dl[i__];
3421  d__[i__] = dl[i__];
3422  temp = d__[i__ + 1];
3423  d__[i__ + 1] = du[i__] - fact * temp;
3424  dl[i__] = du[i__ + 1];
3425  du[i__ + 1] = -fact * dl[i__];
3426  du[i__] = temp;
3427  i__2 = *nrhs;
3428  for(j = 1; j <= i__2; ++j) {
3429  temp = b_ref(i__, j);
3430  b_ref(i__, j) = b_ref(i__ + 1, j);
3431  b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3432 /* L30: */
3433  }
3434  }
3435 /* L40: */
3436  }
3437  if(*n > 1) {
3438  i__ = *n - 1;
3439  if( fabs(d__[i__]) >= fabs(dl[i__]))
3440  {
3441  if(d__[i__] != 0.)
3442  {
3443  fact = dl[i__] / d__[i__];
3444  d__[i__ + 1] -= fact * du[i__];
3445  i__1 = *nrhs;
3446  for(j = 1; j <= i__1; ++j) {
3447  b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3448  i__, j);
3449  /* L50: */
3450  }
3451  }
3452  else
3453  {
3454  *info = i__;
3455  return 0;
3456  }
3457  } else {
3458  fact = d__[i__] / dl[i__];
3459  d__[i__] = dl[i__];
3460  temp = d__[i__ + 1];
3461  d__[i__ + 1] = du[i__] - fact * temp;
3462  du[i__] = temp;
3463  i__1 = *nrhs;
3464  for(j = 1; j <= i__1; ++j) {
3465  temp = b_ref(i__, j);
3466  b_ref(i__, j) = b_ref(i__ + 1, j);
3467  b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3468 /* L60: */
3469  }
3470  }
3471  }
3472  if(d__[*n] == 0.) {
3473  *info = *n;
3474  return 0;
3475  }
3476  }
3477 
3478 /* Back solve with the matrix U from the factorization. */
3479 
3480  if(*nrhs <= 2) {
3481  j = 1;
3482 L70:
3483  b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3484  if(*n > 1) {
3485  b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, j))
3486  / d__[*n - 1];
3487  }
3488  for(i__ = *n - 2; i__ >= 1; --i__) {
3489  b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) - dl[
3490  i__] * b_ref(i__ + 2, j)) / d__[i__];
3491 /* L80: */
3492  }
3493  if(j < *nrhs) {
3494  ++j;
3495  goto L70;
3496  }
3497  } else {
3498  i__1 = *nrhs;
3499  for(j = 1; j <= i__1; ++j) {
3500  b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3501  if(*n > 1) {
3502  b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n,
3503  j)) / d__[*n - 1];
3504  }
3505  for(i__ = *n - 2; i__ >= 1; --i__) {
3506  b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j)
3507  - dl[i__] * b_ref(i__ + 2, j)) / d__[i__];
3508 /* L90: */
3509  }
3510 /* L100: */
3511  }
3512  }
3513 
3514  return 0;
3515 
3516 /* End of DGTSV */
3517 
3518 } /* dgtsv_ */
3519 #endif
3520 #undef b_ref
3521 #endif
3522 /*lint +e725 expected pos indentation */
3523 /*lint +e801 goto is deprecated */
ZERO
#define ZERO
Definition: thirdparty_lapack.cpp:124
DLASWP
STATIC void DLASWP(int32 N, double *A, int32 LDA, int32 K1, int32 K2, int32 IPIV[], int32 INCX)
Definition: thirdparty_lapack.cpp:2335
DGETRF
STATIC void DGETRF(int32, int32, double *, int32, int32[], int32 *)
Definition: thirdparty_lapack.cpp:296
getrs_wrapper
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
Definition: thirdparty_lapack.cpp:69
ILAENV
STATIC int32 ILAENV(int32 ISPEC, const char *NAME, int32 N1, int32 N2, int32 N4)
Definition: thirdparty_lapack.cpp:1424
J_
#define J_(A_)
Definition: iso.h:25
DSWAP
STATIC void DSWAP(int32 n, double dx[], int32 incx, double dy[], int32 incy)
Definition: thirdparty_lapack.cpp:2199
CC
#define CC(I_, J_)
Definition: thirdparty_lapack.cpp:138
N
static const int N
Definition: thirdparty.cpp:3954
INT32_MAX
#define INT32_MAX
Definition: cpu.h:49
IDAMAX
STATIC int32 IDAMAX(int32 n, double dx[], int32 incx)
Definition: thirdparty_lapack.cpp:834
C1
static const double C1
Definition: thirdparty.cpp:1674
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
MIN2
#define MIN2(a, b)
Definition: cddefines.h:803
ONE
#define ONE
Definition: thirdparty_lapack.cpp:123
getrf_wrapper
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
Definition: thirdparty_lapack.cpp:47
STATIC
#define STATIC
Definition: cddefines.h:118
M
static const int M
Definition: thirdparty.cpp:3955
realnum
float realnum
Definition: cddefines.h:124
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
DGER
STATIC void DGER(int32 M, int32 N, double ALPHA, double X[], int32 INCX, double Y[], int32 INCY, double *A, int32 LDA)
Definition: thirdparty_lapack.cpp:2611
da
STATIC double da(double z, double temp, double eden)
Definition: atmdat_3body.cpp:111
DGEMM
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)
Definition: thirdparty_lapack.cpp:2812
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
DGETF2
STATIC void DGETF2(int32 M, int32 N, double *A, int32 LDA, int32 IPIV[], int32 *INFO)
Definition: thirdparty_lapack.cpp:2460
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
LSAME
STATIC int32 LSAME(int32 CA, int32 CB)
Definition: thirdparty_lapack.cpp:725
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
AA
#define AA(I_, J_)
Definition: thirdparty_lapack.cpp:136
DGETRS
STATIC void DGETRS(int32 TRANS, int32 N, int32 NRHS, double *A, int32 LDA, int32 IPIV[], double *B, int32 LDB, int32 *INFO)
Definition: thirdparty_lapack.cpp:527
XERBLA
STATIC void XERBLA(const char *SRNAME, int32 INFO)
Definition: thirdparty_lapack.cpp:251
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
DTRSM
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)
Definition: thirdparty_lapack.cpp:903
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
DSCAL
STATIC void DSCAL(int32 n, double da, double dx[], int32 incx)
Definition: thirdparty_lapack.cpp:2287
BB
#define BB(I_, J_)
Definition: thirdparty_lapack.cpp:137
Generated by   doxygen 1.8.5