cloudy  trunk
 All Data Structures 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 
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 
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  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 */
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  IC = (SUBNAM[0]);
1601  IZ = 'Z';
1602  if( IZ == 90 || IZ == 122 )
1603  {
1604 
1605  /* ASCII character set
1606  * */
1607  if( IC >= 97 && IC <= 122 )
1608  {
1609  SUBNAM[0] = (char)(IC - 32);
1610  for( I=2; I <= 6; I++ )
1611  {
1612  IC = (SUBNAM[I-1]);
1613  if( IC >= 97 && IC <= 122 )
1614  SUBNAM[I - 1] = (char)(IC - 32);
1615  }
1616  }
1617 
1618  }
1619  else if( IZ == 233 || IZ == 169 )
1620  {
1621 
1622  /* EBCDIC character set
1623  * */
1624  if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <= 153)) ||
1625  (IC >= 162 && IC <= 169) )
1626  {
1627  SUBNAM[0] = (char)(IC + 64);
1628  for( I=2; I <= 6; I++ )
1629  {
1630  IC = (SUBNAM[I-1]);
1631  if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <=
1632  153)) || (IC >= 162 && IC <= 169) )
1633  SUBNAM[I - 1] = (char)(IC + 64);
1634  }
1635  }
1636 
1637  }
1638  else if( IZ == 218 || IZ == 250 )
1639  {
1640 
1641  /* Prime machines: ASCII+128
1642  * */
1643  if( IC >= 225 && IC <= 250 )
1644  {
1645  SUBNAM[0] = (char)(IC - 32);
1646  for( I=2; I <= 6; I++ )
1647  {
1648  IC = (SUBNAM[I-1]);
1649  if( IC >= 225 && IC <= 250 )
1650  SUBNAM[I - 1] = (char)(IC - 32);
1651  }
1652  }
1653  }
1654 
1655  C1 = SUBNAM[0];
1656  SNAME = C1 == 'S' || C1 == 'D';
1657  CNAME = C1 == 'C' || C1 == 'Z';
1658  if( !(CNAME || SNAME) )
1659  {
1660  return ILAENV_v;
1661  }
1662 
1663 # if 0
1664  strncpy( C2, SUBNAM+1, 2 );
1665  strncpy( C3, SUBNAM+3, 3 );
1666  strncpy( C4, C3+1, 2 );
1667 # endif
1668 
1669  /* >>chng 00 nov 08, from above to below, insure had run
1670  * off end of string*/
1671  strncpy( C2, SUBNAM+1, 2 );
1672  C2[2] = '\0';
1673  strncpy( C3, SUBNAM+3, 3 );
1674  C3[3] = '\0';
1675  strncpy( C4, C3+1, 2 );
1676  C4[2] = '\0';
1677 
1678  switch( ISPEC )
1679  {
1680  case 1: goto L_110;
1681  case 2: goto L_200;
1682  case 3: goto L_300;
1683  }
1684 
1685 L_110:
1686 
1687  /* ISPEC = 1: block size
1688  *
1689  * In these examples, separate code is provided for setting NB for
1690  * real and complex. We assume that NB will take the same value in
1691  * single or double precision.
1692  * */
1693  NB = 1;
1694 
1695  if( strcmp(C2,"GE") == 0 )
1696  {
1697  if( strcmp(C3,"TRF") == 0 )
1698  {
1699  if( SNAME )
1700  {
1701  NB = 64;
1702  }
1703  else
1704  {
1705  NB = 64;
1706  }
1707  }
1708  else if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) ||
1709  strcmp(C3,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
1710  {
1711  if( SNAME )
1712  {
1713  NB = 32;
1714  }
1715  else
1716  {
1717  NB = 32;
1718  }
1719  }
1720  else if( strcmp(C3,"HRD") == 0 )
1721  {
1722  if( SNAME )
1723  {
1724  NB = 32;
1725  }
1726  else
1727  {
1728  NB = 32;
1729  }
1730  }
1731  else if( strcmp(C3,"BRD") == 0 )
1732  {
1733  if( SNAME )
1734  {
1735  NB = 32;
1736  }
1737  else
1738  {
1739  NB = 32;
1740  }
1741  }
1742  else if( strcmp(C3,"TRI") == 0 )
1743  {
1744  if( SNAME )
1745  {
1746  NB = 64;
1747  }
1748  else
1749  {
1750  NB = 64;
1751  }
1752  }
1753  }
1754  else if( strcmp(C2,"PO") == 0 )
1755  {
1756  if( strcmp(C3,"TRF") == 0 )
1757  {
1758  if( SNAME )
1759  {
1760  NB = 64;
1761  }
1762  else
1763  {
1764  NB = 64;
1765  }
1766  }
1767  }
1768  else if( strcmp(C2,"SY") == 0 )
1769  {
1770  if( strcmp(C3,"TRF") == 0 )
1771  {
1772  if( SNAME )
1773  {
1774  NB = 64;
1775  }
1776  else
1777  {
1778  NB = 64;
1779  }
1780  }
1781  else if( SNAME && strcmp(C3,"TRD") == 0 )
1782  {
1783  NB = 1;
1784  }
1785  else if( SNAME && strcmp(C3,"GST") == 0 )
1786  {
1787  NB = 64;
1788  }
1789  }
1790  else if( CNAME && strcmp(C2,"HE") == 0 )
1791  {
1792  if( strcmp(C3,"TRF") == 0 )
1793  {
1794  NB = 64;
1795  }
1796  else if( strcmp(C3,"TRD") == 0 )
1797  {
1798  NB = 1;
1799  }
1800  else if( strcmp(C3,"GST") == 0 )
1801  {
1802  NB = 64;
1803  }
1804  }
1805  else if( SNAME && strcmp(C2,"OR") == 0 )
1806  {
1807  if( C3[0] == 'G' )
1808  {
1809  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1810  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1811  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1812  0 )
1813  {
1814  NB = 32;
1815  }
1816  }
1817  else if( C3[0] == 'M' )
1818  {
1819  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1820  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1821  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1822  0 )
1823  {
1824  NB = 32;
1825  }
1826  }
1827  }
1828  else if( CNAME && strcmp(C2,"UN") == 0 )
1829  {
1830  if( C3[0] == 'G' )
1831  {
1832  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1833  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1834  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1835  0 )
1836  {
1837  NB = 32;
1838  }
1839  }
1840  else if( C3[0] == 'M' )
1841  {
1842  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
1843  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
1844  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
1845  0 )
1846  {
1847  NB = 32;
1848  }
1849  }
1850  }
1851  else if( strcmp(C2,"GB") == 0 )
1852  {
1853  if( strcmp(C3,"TRF") == 0 )
1854  {
1855  if( SNAME )
1856  {
1857  if( N4 <= 64 )
1858  {
1859  NB = 1;
1860  }
1861  else
1862  {
1863  NB = 32;
1864  }
1865  }
1866  else
1867  {
1868  if( N4 <= 64 )
1869  {
1870  NB = 1;
1871  }
1872  else
1873  {
1874  NB = 32;
1875  }
1876  }
1877  }
1878  }
1879  else if( strcmp(C2,"PB") == 0 )
1880  {
1881  if( strcmp(C3,"TRF") == 0 )
1882  {
1883  if( SNAME )
1884  {
1885  if( N2 <= 64 )
1886  {
1887  NB = 1;
1888  }
1889  else
1890  {
1891  NB = 32;
1892  }
1893  }
1894  else
1895  {
1896  if( N2 <= 64 )
1897  {
1898  NB = 1;
1899  }
1900  else
1901  {
1902  NB = 32;
1903  }
1904  }
1905  }
1906  }
1907  else if( strcmp(C2,"TR") == 0 )
1908  {
1909  if( strcmp(C3,"TRI") == 0 )
1910  {
1911  if( SNAME )
1912  {
1913  NB = 64;
1914  }
1915  else
1916  {
1917  NB = 64;
1918  }
1919  }
1920  }
1921  else if( strcmp(C2,"LA") == 0 )
1922  {
1923  if( strcmp(C3,"UUM") == 0 )
1924  {
1925  if( SNAME )
1926  {
1927  NB = 64;
1928  }
1929  else
1930  {
1931  NB = 64;
1932  }
1933  }
1934  }
1935  else if( SNAME && strcmp(C2,"ST") == 0 )
1936  {
1937  if( strcmp(C3,"EBZ") == 0 )
1938  {
1939  NB = 1;
1940  }
1941  }
1942  ILAENV_v = NB;
1943  return ILAENV_v;
1944 
1945 L_200:
1946 
1947  /* ISPEC = 2: minimum block size
1948  * */
1949  NBMIN = 2;
1950  if( strcmp(C2,"GE") == 0 )
1951  {
1952  if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
1953  ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
1954  {
1955  if( SNAME )
1956  {
1957  NBMIN = 2;
1958  }
1959  else
1960  {
1961  NBMIN = 2;
1962  }
1963  }
1964  else if( strcmp(C3,"HRD") == 0 )
1965  {
1966  if( SNAME )
1967  {
1968  NBMIN = 2;
1969  }
1970  else
1971  {
1972  NBMIN = 2;
1973  }
1974  }
1975  else if( strcmp(C3,"BRD") == 0 )
1976  {
1977  if( SNAME )
1978  {
1979  NBMIN = 2;
1980  }
1981  else
1982  {
1983  NBMIN = 2;
1984  }
1985  }
1986  else if( strcmp(C3,"TRI") == 0 )
1987  {
1988  if( SNAME )
1989  {
1990  NBMIN = 2;
1991  }
1992  else
1993  {
1994  NBMIN = 2;
1995  }
1996  }
1997  }
1998  else if( strcmp(C2,"SY") == 0 )
1999  {
2000  if( strcmp(C3,"TRF") == 0 )
2001  {
2002  if( SNAME )
2003  {
2004  NBMIN = 8;
2005  }
2006  else
2007  {
2008  NBMIN = 8;
2009  }
2010  }
2011  else if( SNAME && strcmp(C3,"TRD") == 0 )
2012  {
2013  NBMIN = 2;
2014  }
2015  }
2016  else if( CNAME && strcmp(C2,"HE") == 0 )
2017  {
2018  if( strcmp(C3,"TRD") == 0 )
2019  {
2020  NBMIN = 2;
2021  }
2022  }
2023  else if( SNAME && strcmp(C2,"OR") == 0 )
2024  {
2025  if( C3[0] == 'G' )
2026  {
2027  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2028  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2029  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2030  0 )
2031  {
2032  NBMIN = 2;
2033  }
2034  }
2035  else if( C3[0] == 'M' )
2036  {
2037  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2038  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2039  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2040  0 )
2041  {
2042  NBMIN = 2;
2043  }
2044  }
2045  }
2046  else if( CNAME && strcmp(C2,"UN") == 0 )
2047  {
2048  if( C3[0] == 'G' )
2049  {
2050  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2051  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2052  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2053  0 )
2054  {
2055  NBMIN = 2;
2056  }
2057  }
2058  else if( C3[0] == 'M' )
2059  {
2060  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2061  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2062  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2063  0 )
2064  {
2065  NBMIN = 2;
2066  }
2067  }
2068  }
2069  ILAENV_v = NBMIN;
2070  return ILAENV_v;
2071 
2072 L_300:
2073 
2074  /* ISPEC = 3: crossover point
2075  * */
2076  NX = 0;
2077  if( strcmp(C2,"GE") == 0 )
2078  {
2079  if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
2080  ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
2081  {
2082  if( SNAME )
2083  {
2084  NX = 128;
2085  }
2086  else
2087  {
2088  NX = 128;
2089  }
2090  }
2091  else if( strcmp(C3,"HRD") == 0 )
2092  {
2093  if( SNAME )
2094  {
2095  NX = 128;
2096  }
2097  else
2098  {
2099  NX = 128;
2100  }
2101  }
2102  else if( strcmp(C3,"BRD") == 0 )
2103  {
2104  if( SNAME )
2105  {
2106  NX = 128;
2107  }
2108  else
2109  {
2110  NX = 128;
2111  }
2112  }
2113  }
2114  else if( strcmp(C2,"SY") == 0 )
2115  {
2116  if( SNAME && strcmp(C3,"TRD") == 0 )
2117  {
2118  NX = 1;
2119  }
2120  }
2121  else if( CNAME && strcmp(C2,"HE") == 0 )
2122  {
2123  if( strcmp(C3,"TRD") == 0 )
2124  {
2125  NX = 1;
2126  }
2127  }
2128  else if( SNAME && strcmp(C2,"OR") == 0 )
2129  {
2130  if( C3[0] == 'G' )
2131  {
2132  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2133  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2134  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2135  0 )
2136  {
2137  NX = 128;
2138  }
2139  }
2140  }
2141  else if( CNAME && strcmp(C2,"UN") == 0 )
2142  {
2143  if( C3[0] == 'G' )
2144  {
2145  if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
2146  strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
2147  ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
2148  0 )
2149  {
2150  NX = 128;
2151  }
2152  }
2153  }
2154  ILAENV_v = NX;
2155  return ILAENV_v;
2156 
2157 L_400:
2158 
2159  /* ISPEC = 4: number of shifts (used by xHSEQR)
2160  * */
2161  ILAENV_v = 6;
2162  return ILAENV_v;
2163 
2164 L_500:
2165 
2166  /* ISPEC = 5: minimum column dimension (not used)
2167  * */
2168  ILAENV_v = 2;
2169  return ILAENV_v;
2170 
2171 L_600:
2172 
2173  /* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
2174  * */
2175  ILAENV_v = (int32)((realnum)(MIN2(N1,N2))*1.6e0);
2176  return ILAENV_v;
2177 
2178 L_700:
2179 
2180  /* ISPEC = 7: number of processors (not used)
2181  * */
2182  ILAENV_v = 1;
2183  return ILAENV_v;
2184 
2185 L_800:
2186 
2187  /* ISPEC = 8: crossover point for multishift (used by xHSEQR)
2188  * */
2189  ILAENV_v = 50;
2190  return ILAENV_v;
2191 
2192  /* End of ILAENV
2193  * */
2194 }
2195 
2196 /*DSWAP lapack routine */
2197 
2198 STATIC void DSWAP(int32 n,
2199  double dx[],
2200  int32 incx,
2201  double dy[],
2202  int32 incy)
2203 {
2204  int32 i,
2205  ix,
2206  iy,
2207  m;
2208  double dtemp;
2209 
2210  DEBUG_ENTRY( "DSWAP()" );
2211 
2212  /* interchanges two vectors.
2213  * uses unrolled loops for increments equal one.
2214  * jack dongarra, lapack, 3/11/78.
2215  * modified 12/3/93, array(1) declarations changed to array(*)
2216  * */
2217 
2218  if( n <= 0 )
2219  {
2220  return;}
2221  if( incx == 1 && incy == 1 )
2222  goto L_20;
2223 
2224  /* code for unequal increments or equal increments not equal
2225  * to 1
2226  * */
2227  ix = 1;
2228  iy = 1;
2229 
2230  if( incx < 0 )
2231  ix = (-n + 1)*incx + 1;
2232 
2233  if( incy < 0 )
2234  iy = (-n + 1)*incy + 1;
2235 
2236  for( i=0; i < n; i++ )
2237  {
2238  dtemp = dx[ix-1];
2239  dx[ix-1] = dy[iy-1];
2240  dy[iy-1] = dtemp;
2241  ix += incx;
2242  iy += incy;
2243  }
2244  return;
2245 
2246  /* code for both increments equal to 1
2247  *
2248  *
2249  * clean-up loop
2250  * */
2251 L_20:
2252  m = n%3;
2253  if( m == 0 )
2254  goto L_40;
2255 
2256  for( i=0; i < m; i++ )
2257  {
2258  dtemp = dx[i];
2259  dx[i] = dy[i];
2260  dy[i] = dtemp;
2261  }
2262 
2263  if( n < 3 )
2264  {
2265  return;
2266  }
2267 
2268 L_40:
2269  for( i=m; i < n; i += 3 )
2270  {
2271  dtemp = dx[i];
2272  dx[i] = dy[i];
2273  dy[i] = dtemp;
2274  dtemp = dx[i+1];
2275  dx[i+1] = dy[i+1];
2276  dy[i+1] = dtemp;
2277  dtemp = dx[i+2];
2278  dx[i+2] = dy[i+2];
2279  dy[i+2] = dtemp;
2280  }
2281  return;
2282 }
2283 
2284 /*DSCAL lapack routine */
2285 
2286 STATIC void DSCAL(int32 n,
2287  double da,
2288  double dx[],
2289  int32 incx)
2290 {
2291  int32 i,
2292  m,
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  m = n%5;
2327  if( m == 0 )
2328  goto L_40;
2329 
2330  for( i=0; i < m; i++ )
2331  {
2332  dx[i] *= da;
2333  }
2334 
2335  if( n < 5 )
2336  {
2337  return;
2338  }
2339 
2340 L_40:
2341 
2342  for( i=m; i < n; i += 5 )
2343  {
2344  dx[i] *= da;
2345  dx[i+1] *= da;
2346  dx[i+2] *= da;
2347  dx[i+3] *= da;
2348  dx[i+4] *= da;
2349  }
2350  return;
2351 }
2352 
2353 /*DLASWP -- LAPACK auxiliary routine (version 2.0) --*/
2354 
2355 STATIC void DLASWP(int32 N,
2356  double *A,
2357  int32 LDA,
2358  int32 K1,
2359  int32 K2,
2360  int32 IPIV[],
2361  int32 INCX)
2362 {
2363  int32 I,
2364  IP,
2365  IX,
2366  I_;
2367 
2368  DEBUG_ENTRY( "DLASWP()" );
2369 
2370  /* -- LAPACK auxiliary routine (version 2.0) --
2371  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2372  * Courant Institute, Argonne National Lab, and Rice University
2373  * October 31, 1992
2374  *
2375  * .. Scalar Arguments .. */
2376  /* ..
2377  * .. Array Arguments .. */
2378  /* ..
2379  *
2380  * Purpose
2381  * =======
2382  *
2383  * DLASWP performs a series of row interchanges on the matrix A.
2384  * One row interchange is initiated for each of rows K1 through K2 of A.
2385  *
2386  * Arguments
2387  * =========
2388  *
2389  * N (input) INTEGER
2390  * The number of columns of the matrix A.
2391  *
2392  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
2393  * On entry, the matrix of column dimension N to which the row
2394  * interchanges will be applied.
2395  * On exit, the permuted matrix.
2396  *
2397  * LDA (input) INTEGER
2398  * The leading dimension of the array A.
2399  *
2400  * K1 (input) INTEGER
2401  * The first element of IPIV for which a row interchange will
2402  * be done.
2403  *
2404  * K2 (input) INTEGER
2405  * The last element of IPIV for which a row interchange will
2406  * be done.
2407  *
2408  * IPIV (input) INTEGER array, dimension (M*ABS(INCX))
2409  * The vector of pivot indices. Only the elements in positions
2410  * K1 through K2 of IPIV are accessed.
2411  * IPIV(K) = L implies rows K and L are to be interchanged.
2412  *
2413  * INCX (input) INTEGER
2414  * The increment between successive values of IPIV. If IPIV
2415  * is negative, the pivots are applied in reverse order.
2416  *
2417  * =====================================================================
2418  *
2419  * .. Local Scalars .. */
2420  /* ..
2421  * .. External Subroutines .. */
2422  /* ..
2423  * .. Executable Statements ..
2424  *
2425  * Interchange row I with row IPIV(I) for each of rows K1 through K2.
2426  * */
2427  if( INCX == 0 )
2428  {
2429  return;}
2430  if( INCX > 0 )
2431  {
2432  IX = K1;
2433  }
2434  else
2435  {
2436  IX = 1 + (1 - K2)*INCX;
2437  }
2438  if( INCX == 1 )
2439  {
2440  for( I=K1; I <= K2; I++ )
2441  {
2442  I_ = I - 1;
2443  IP = IPIV[I_];
2444  if( IP != I )
2445  DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2446  }
2447  }
2448  else if( INCX > 1 )
2449  {
2450  for( I=K1; I <= K2; I++ )
2451  {
2452  I_ = I - 1;
2453  IP = IPIV[IX-1];
2454  if( IP != I )
2455  DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2456  IX += INCX;
2457  }
2458  }
2459  else if( INCX < 0 )
2460  {
2461  for( I=K2; I >= K1; I-- )
2462  {
2463  I_ = I - 1;
2464  IP = IPIV[IX-1];
2465  if( IP != I )
2466  DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
2467  IX += INCX;
2468  }
2469  }
2470 
2471  return;
2472 
2473  /* End of DLASWP
2474  * */
2475 #undef A
2476 }
2477 
2478 /*DGETF2 lapack service routine */
2479 
2480 STATIC void DGETF2(int32 M,
2481  int32 N,
2482  double *A,
2483  int32 LDA,
2484  int32 IPIV[],
2485  int32 *INFO)
2486 {
2487  int32 J,
2488  JP,
2489  J_,
2490  limit;
2491 
2492  DEBUG_ENTRY( "DGETF2()" );
2493 
2494  /* -- LAPACK routine (version 2.0) --
2495  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2496  * Courant Institute, Argonne National Lab, and Rice University
2497  * June 30, 1992
2498  *
2499  * .. Scalar Arguments .. */
2500  /* ..
2501  * .. Array Arguments .. */
2502  /* ..
2503  *
2504  * Purpose
2505  * =======
2506  *
2507  * DGETF2 computes an LU factorization of a general m-by-n matrix A
2508  * using partial pivoting with row interchanges.
2509  *
2510  * The factorization has the form
2511  * A = P * L * U
2512  * where P is a permutation matrix, L is lower triangular with unit
2513  * diagonal elements (lower trapezoidal if m > n), and U is upper
2514  * triangular (upper trapezoidal if m < n).
2515  *
2516  * This is the right-looking Level 2 BLAS version of the algorithm.
2517  *
2518  * Arguments
2519  * =========
2520  *
2521  * M (input) INTEGER
2522  * The number of rows of the matrix A. M >= 0.
2523  *
2524  * N (input) INTEGER
2525  * The number of columns of the matrix A. N >= 0.
2526  *
2527  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
2528  * On entry, the m by n matrix to be factored.
2529  * On exit, the factors L and U from the factorization
2530  * A = P*L*U; the unit diagonal elements of L are not stored.
2531  *
2532  * LDA (input) INTEGER
2533  * The leading dimension of the array A. LDA >= MAX(1,M).
2534  *
2535  * IPIV (output) INTEGER array, dimension (MIN(M,N))
2536  * The pivot indices; for 1 <= i <= MIN(M,N), row i of the
2537  * matrix was interchanged with row IPIV(i).
2538  *
2539  * INFO (output) INTEGER
2540  * = 0: successful exit
2541  * < 0: if INFO = -k, the k-th argument had an illegal value
2542  * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
2543  * has been completed, but the factor U is exactly
2544  * singular, and division by zero will occur if it is used
2545  * to solve a system of equations.
2546  *
2547  * =====================================================================
2548  *
2549  * .. Parameters .. */
2550  /* ..
2551  * .. Local Scalars .. */
2552  /* ..
2553  * .. External Functions .. */
2554  /* ..
2555  * .. External Subroutines .. */
2556  /* ..
2557  * .. Intrinsic Functions .. */
2558  /* ..
2559  * .. Executable Statements ..
2560  *
2561  * Test the input parameters.
2562  * */
2563  *INFO = 0;
2564  if( M < 0 )
2565  {
2566  *INFO = -1;
2567  }
2568  else if( N < 0 )
2569  {
2570  *INFO = -2;
2571  }
2572  else if( LDA < MAX2(1,M) )
2573  {
2574  *INFO = -4;
2575  }
2576  if( *INFO != 0 )
2577  {
2578  XERBLA("DGETF2",-*INFO);
2579  /* XERBLA does not return */
2580  }
2581 
2582  /* Quick return if possible
2583  * */
2584  if( M == 0 || N == 0 )
2585  {
2586  return;}
2587 
2588  limit = MIN2(M,N);
2589  for( J=1; J <= limit; J++ )
2590  {
2591  J_ = J - 1;
2592 
2593  /* Find pivot and test for singularity.
2594  * */
2595  JP = J - 1 + IDAMAX(M-J+1,&AA(J_,J_),1);
2596  IPIV[J_] = JP;
2597  if( AA(J_,JP-1) != ZERO )
2598  {
2599  /* Apply the interchange to columns 1:N.
2600  * */
2601  if( JP != J )
2602  DSWAP(N,&AA(0,J_),LDA,&AA(0,JP-1),LDA);
2603 
2604  /* Compute elements J+1:M of J-th column.
2605  * */
2606  if( J < M )
2607  DSCAL(M-J,ONE/AA(J_,J_),&AA(J_,J_+1),1);
2608  }
2609  else if( *INFO == 0 )
2610  {
2611  *INFO = J;
2612  }
2613 
2614  if( J < MIN2(M,N) )
2615  {
2616  /* Update trailing submatrix.
2617  * */
2618  DGER(M-J,N-J,-ONE,&AA(J_,J_+1),1,&AA(J_+1,J_),LDA,&AA(J_+1,J_+1),
2619  LDA);
2620  }
2621  }
2622  return;
2623 
2624  /* End of DGETF2
2625  * */
2626 #undef A
2627 }
2628 
2629 /*DGER service routine for matrix inversion */
2630 
2631 STATIC void DGER(int32 M,
2632  int32 N,
2633  double ALPHA,
2634  double X[],
2635  int32 INCX,
2636  double Y[],
2637  int32 INCY,
2638  double *A,
2639  int32 LDA)
2640 {
2641  int32 I,
2642  INFO,
2643  IX,
2644  I_,
2645  J,
2646  JY,
2647  J_,
2648  KX;
2649  double TEMP;
2650 
2651  DEBUG_ENTRY( "DGER()" );
2652  /* .. Scalar Arguments .. */
2653  /* .. Array Arguments .. */
2654  /* ..
2655  *
2656  * Purpose
2657  * =======
2658  *
2659  * DGER performs the rank 1 operation
2660  *
2661  * A := alpha*x*y' + A,
2662  *
2663  * where alpha is a scalar, x is an m element vector, y is an n element
2664  * vector and A is an m by n matrix.
2665  *
2666  * Parameters
2667  * ==========
2668  *
2669  * M - INTEGER.
2670  * On entry, M specifies the number of rows of the matrix A.
2671  * M must be at least zero.
2672  * Unchanged on exit.
2673  *
2674  * N - INTEGER.
2675  * On entry, N specifies the number of columns of the matrix A.
2676  * N must be at least zero.
2677  * Unchanged on exit.
2678  *
2679  * ALPHA - DOUBLE PRECISION.
2680  * On entry, ALPHA specifies the scalar alpha.
2681  * Unchanged on exit.
2682  *
2683  * X - DOUBLE PRECISION array of dimension at least
2684  * ( 1 + ( m - 1 )*ABS( INCX ) ).
2685  * Before entry, the incremented array X must contain the m
2686  * element vector x.
2687  * Unchanged on exit.
2688  *
2689  * INCX - INTEGER.
2690  * On entry, INCX specifies the increment for the elements of
2691  * X. INCX must not be zero.
2692  * Unchanged on exit.
2693  *
2694  * Y - DOUBLE PRECISION array of dimension at least
2695  * ( 1 + ( n - 1 )*ABS( INCY ) ).
2696  * Before entry, the incremented array Y must contain the n
2697  * element vector y.
2698  * Unchanged on exit.
2699  *
2700  * INCY - INTEGER.
2701  * On entry, INCY specifies the increment for the elements of
2702  * Y. INCY must not be zero.
2703  * Unchanged on exit.
2704  *
2705  * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2706  * Before entry, the leading m by n part of the array A must
2707  * contain the matrix of coefficients. On exit, A is
2708  * overwritten by the updated matrix.
2709  *
2710  * LDA - INTEGER.
2711  * On entry, LDA specifies the first dimension of A as declared
2712  * in the calling (sub) program. LDA must be at least
2713  * MAX( 1, m ).
2714  * Unchanged on exit.
2715  *
2716  *
2717  * Level 2 Blas routine.
2718  *
2719  * -- Written on 22-October-1986.
2720  * Jack Dongarra, Argonne National Lab.
2721  * Jeremy Du Croz, Nag Central Office.
2722  * Sven Hammarling, Nag Central Office.
2723  * Richard Hanson, Sandia National Labs.
2724  *
2725  *
2726  * .. Parameters .. */
2727  /* .. Local Scalars .. */
2728  /* .. External Subroutines .. */
2729  /* .. Intrinsic Functions .. */
2730  /* ..
2731  * .. Executable Statements ..
2732  *
2733  * Test the input parameters.
2734  * */
2735  INFO = 0;
2736  if( M < 0 )
2737  {
2738  INFO = 1;
2739  }
2740  else if( N < 0 )
2741  {
2742  INFO = 2;
2743  }
2744  else if( INCX == 0 )
2745  {
2746  INFO = 5;
2747  }
2748  else if( INCY == 0 )
2749  {
2750  INFO = 7;
2751  }
2752  else if( LDA < MAX2(1,M) )
2753  {
2754  INFO = 9;
2755  }
2756  if( INFO != 0 )
2757  {
2758  XERBLA("DGER ",INFO);
2759  /* XERBLA does not return */
2760  }
2761 
2762  /* Quick return if possible.
2763  * */
2764  if( ((M == 0) || (N == 0)) || (ALPHA == ZERO) )
2765  {
2766  return;}
2767 
2768  /* Start the operations. In this version the elements of A are
2769  * accessed sequentially with one pass through A.
2770  * */
2771  if( INCY > 0 )
2772  {
2773  JY = 1;
2774  }
2775  else
2776  {
2777  JY = 1 - (N - 1)*INCY;
2778  }
2779  if( INCX == 1 )
2780  {
2781  for( J=1; J <= N; J++ )
2782  {
2783  J_ = J - 1;
2784  if( Y[JY-1] != ZERO )
2785  {
2786  TEMP = ALPHA*Y[JY-1];
2787  for( I=1; I <= M; I++ )
2788  {
2789  I_ = I - 1;
2790  AA(J_,I_) += X[I_]*TEMP;
2791  }
2792  }
2793  JY += INCY;
2794  }
2795  }
2796  else
2797  {
2798  if( INCX > 0 )
2799  {
2800  KX = 1;
2801  }
2802  else
2803  {
2804  KX = 1 - (M - 1)*INCX;
2805  }
2806  for( J=1; J <= N; J++ )
2807  {
2808  J_ = J - 1;
2809  if( Y[JY-1] != ZERO )
2810  {
2811  TEMP = ALPHA*Y[JY-1];
2812  IX = KX;
2813  for( I=1; I <= M; I++ )
2814  {
2815  I_ = I - 1;
2816  AA(J_,I_) += X[IX-1]*TEMP;
2817  IX += INCX;
2818  }
2819  }
2820  JY += INCY;
2821  }
2822  }
2823 
2824  return;
2825 
2826  /* End of DGER .
2827  * */
2828 #undef A
2829 }
2830 
2831 /* DGEMM matrix inversion helper routine*/
2832 
2833 STATIC void DGEMM(int32 TRANSA,
2834  int32 TRANSB,
2835  int32 M,
2836  int32 N,
2837  int32 K,
2838  double ALPHA,
2839  double *A,
2840  int32 LDA,
2841  double *B,
2842  int32 LDB,
2843  double BETA,
2844  double *C,
2845  int32 LDC)
2846 {
2847  int32 NOTA,
2848  NOTB;
2849  int32 I,
2850  INFO,
2851  I_,
2852  J,
2853  J_,
2854  L,
2855  L_,
2856  NROWA,
2857  NROWB;
2858  double TEMP;
2859 
2860  DEBUG_ENTRY( "DGEMM()" );
2861  /* .. Scalar Arguments .. */
2862  /* .. Array Arguments .. */
2863  /* ..
2864  *
2865  * Purpose
2866  * =======
2867  *
2868  * DGEMM performs one of the matrix-matrix operations
2869  *
2870  * C := alpha*op( A )*op( B ) + beta*C,
2871  *
2872  * where op( X ) is one of
2873  *
2874  * op( X ) = X or op( X ) = X',
2875  *
2876  * alpha and beta are scalars, and A, B and C are matrices, with op( A )
2877  * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
2878  *
2879  * Parameters
2880  * ==========
2881  *
2882  * TRANSA - CHARACTER*1.
2883  * On entry, TRANSA specifies the form of op( A ) to be used in
2884  * the matrix multiplication as follows:
2885  *
2886  * TRANSA = 'N' or 'n', op( A ) = A.
2887  *
2888  * TRANSA = 'T' or 't', op( A ) = A'.
2889  *
2890  * TRANSA = 'C' or 'c', op( A ) = A'.
2891  *
2892  * Unchanged on exit.
2893  *
2894  * TRANSB - CHARACTER*1.
2895  * On entry, TRANSB specifies the form of op( B ) to be used in
2896  * the matrix multiplication as follows:
2897  *
2898  * TRANSB = 'N' or 'n', op( B ) = B.
2899  *
2900  * TRANSB = 'T' or 't', op( B ) = B'.
2901  *
2902  * TRANSB = 'C' or 'c', op( B ) = B'.
2903  *
2904  * Unchanged on exit.
2905  *
2906  * M - INTEGER.
2907  * On entry, M specifies the number of rows of the matrix
2908  * op( A ) and of the matrix C. M must be at least zero.
2909  * Unchanged on exit.
2910  *
2911  * N - INTEGER.
2912  * On entry, N specifies the number of columns of the matrix
2913  * op( B ) and the number of columns of the matrix C. N must be
2914  * at least zero.
2915  * Unchanged on exit.
2916  *
2917  * K - INTEGER.
2918  * On entry, K specifies the number of columns of the matrix
2919  * op( A ) and the number of rows of the matrix op( B ). K must
2920  * be at least zero.
2921  * Unchanged on exit.
2922  *
2923  * ALPHA - DOUBLE PRECISION.
2924  * On entry, ALPHA specifies the scalar alpha.
2925  * Unchanged on exit.
2926  *
2927  * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
2928  * k when TRANSA = 'N' or 'n', and is m otherwise.
2929  * Before entry with TRANSA = 'N' or 'n', the leading m by k
2930  * part of the array A must contain the matrix A, otherwise
2931  * the leading k by m part of the array A must contain the
2932  * matrix A.
2933  * Unchanged on exit.
2934  *
2935  * LDA - INTEGER.
2936  * On entry, LDA specifies the first dimension of A as declared
2937  * in the calling (sub) program. When TRANSA = 'N' or 'n' then
2938  * LDA must be at least MAX( 1, m ), otherwise LDA must be at
2939  * least MAX( 1, k ).
2940  * Unchanged on exit.
2941  *
2942  * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
2943  * n when TRANSB = 'N' or 'n', and is k otherwise.
2944  * Before entry with TRANSB = 'N' or 'n', the leading k by n
2945  * part of the array B must contain the matrix B, otherwise
2946  * the leading n by k part of the array B must contain the
2947  * matrix B.
2948  * Unchanged on exit.
2949  *
2950  * LDB - INTEGER.
2951  * On entry, LDB specifies the first dimension of B as declared
2952  * in the calling (sub) program. When TRANSB = 'N' or 'n' then
2953  * LDB must be at least MAX( 1, k ), otherwise LDB must be at
2954  * least MAX( 1, n ).
2955  * Unchanged on exit.
2956  *
2957  * BETA - DOUBLE PRECISION.
2958  * On entry, BETA specifies the scalar beta. When BETA is
2959  * supplied as zero then C need not be set on input.
2960  * Unchanged on exit.
2961  *
2962  * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
2963  * Before entry, the leading m by n part of the array C must
2964  * contain the matrix C, except when beta is zero, in which
2965  * case C need not be set on entry.
2966  * On exit, the array C is overwritten by the m by n matrix
2967  * ( alpha*op( A )*op( B ) + beta*C ).
2968  *
2969  * LDC - INTEGER.
2970  * On entry, LDC specifies the first dimension of C as declared
2971  * in the calling (sub) program. LDC must be at least
2972  * MAX( 1, m ).
2973  * Unchanged on exit.
2974  *
2975  *
2976  * Level 3 Blas routine.
2977  *
2978  * -- Written on 8-February-1989.
2979  * Jack Dongarra, Argonne National Laboratory.
2980  * Iain Duff, AERE Harwell.
2981  * Jeremy Du Croz, Numerical Algorithms Group Ltd.
2982  * Sven Hammarling, Numerical Algorithms Group Ltd.
2983  *
2984  *
2985  * .. External Functions .. */
2986  /* .. External Subroutines .. */
2987  /* .. Intrinsic Functions .. */
2988  /* .. Local Scalars .. */
2989  /* .. Parameters .. */
2990  /* ..
2991  * .. Executable Statements ..
2992  *
2993  * Set NOTA and NOTB as true if A and B respectively are not
2994  * transposed and set NROWA, NROWB as the number of rows
2995  * and columns of A and the number of rows of B respectively.
2996  * */
2997  NOTA = LSAME(TRANSA,'N');
2998  NOTB = LSAME(TRANSB,'N');
2999 
3000  if( NOTA )
3001  {
3002  NROWA = M;
3003  }
3004  else
3005  {
3006  NROWA = K;
3007  }
3008 
3009  if( NOTB )
3010  {
3011  NROWB = K;
3012  }
3013  else
3014  {
3015  NROWB = N;
3016  }
3017 
3018  /* Test the input parameters.
3019  * */
3020  INFO = 0;
3021  if( ((!NOTA) && (!LSAME(TRANSA,'C'))) && (!LSAME(TRANSA,'T')) )
3022  {
3023  INFO = 1;
3024  }
3025  else if(
3026  ((!NOTB) && (!LSAME(TRANSB,'C'))) && (!LSAME(TRANSB,'T')) )
3027  {
3028  INFO = 2;
3029  }
3030 
3031  else if( M < 0 )
3032  {
3033  INFO = 3;
3034  }
3035 
3036  else if( N < 0 )
3037  {
3038  INFO = 4;
3039  }
3040 
3041  else if( K < 0 )
3042  {
3043  INFO = 5;
3044  }
3045 
3046  else if( LDA < MAX2(1,NROWA) )
3047  {
3048  INFO = 8;
3049  }
3050 
3051  else if( LDB < MAX2(1,NROWB) )
3052  {
3053  INFO = 10;
3054  }
3055 
3056  else if( LDC < MAX2(1,M) )
3057  {
3058  INFO = 13;
3059  }
3060 
3061  if( INFO != 0 )
3062  {
3063  XERBLA("DGEMM ",INFO);
3064  /* XERBLA does not return */
3065  }
3066 
3067  /* Quick return if possible.
3068  * */
3069  if( ((M == 0) || (N == 0)) || (((ALPHA == ZERO) || (K == 0)) &&
3070  (BETA == ONE)) )
3071  {
3072  return;
3073  }
3074 
3075  /* And if alpha.eq.zero. */
3076  if( ALPHA == ZERO )
3077  {
3078  if( BETA == ZERO )
3079  {
3080  for( J=1; J <= N; J++ )
3081  {
3082  J_ = J - 1;
3083  for( I=1; I <= M; I++ )
3084  {
3085  I_ = I - 1;
3086  CC(J_,I_) = ZERO;
3087  }
3088  }
3089  }
3090 
3091  else
3092  {
3093  for( J=1; J <= N; J++ )
3094  {
3095  J_ = J - 1;
3096  for( I=1; I <= M; I++ )
3097  {
3098  I_ = I - 1;
3099  CC(J_,I_) *= BETA;
3100  }
3101  }
3102  }
3103  return;
3104  }
3105 
3106  /* Start the operations.*/
3107  if( NOTB )
3108  {
3109 
3110  if( NOTA )
3111  {
3112 
3113  /* Form C := alpha*A*B + beta*C.
3114  * */
3115  for( J=1; J <= N; J++ )
3116  {
3117  J_ = J - 1;
3118  if( BETA == ZERO )
3119  {
3120  for( I=1; I <= M; I++ )
3121  {
3122  I_ = I - 1;
3123  CC(J_,I_) = ZERO;
3124  }
3125  }
3126 
3127  else if( BETA != ONE )
3128  {
3129  for( I=1; I <= M; I++ )
3130  {
3131  I_ = I - 1;
3132  CC(J_,I_) *= BETA;
3133  }
3134  }
3135 
3136  for( L=1; L <= K; L++ )
3137  {
3138  L_ = L - 1;
3139  if( BB(J_,L_) != ZERO )
3140  {
3141  TEMP = ALPHA*BB(J_,L_);
3142  for( I=1; I <= M; I++ )
3143  {
3144  I_ = I - 1;
3145  CC(J_,I_) += TEMP*AA(L_,I_);
3146  }
3147  }
3148  }
3149  }
3150  }
3151  else
3152  {
3153 
3154  /* Form C := alpha*A'*B + beta*C */
3155  for( J=1; J <= N; J++ )
3156  {
3157  J_ = J - 1;
3158  for( I=1; I <= M; I++ )
3159  {
3160  I_ = I - 1;
3161  TEMP = ZERO;
3162  for( L=1; L <= K; L++ )
3163  {
3164  L_ = L - 1;
3165  TEMP += AA(I_,L_)*BB(J_,L_);
3166  }
3167 
3168  if( BETA == ZERO )
3169  {
3170  CC(J_,I_) = ALPHA*TEMP;
3171  }
3172  else
3173  {
3174  CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
3175  }
3176  }
3177  }
3178  }
3179  }
3180  else
3181  {
3182  if( NOTA )
3183  {
3184 
3185  /* Form C := alpha*A*B' + beta*C
3186  * */
3187  for( J=1; J <= N; J++ )
3188  {
3189  J_ = J - 1;
3190  if( BETA == ZERO )
3191  {
3192  for( I=1; I <= M; I++ )
3193  {
3194  I_ = I - 1;
3195  CC(J_,I_) = ZERO;
3196  }
3197  }
3198 
3199  else if( BETA != ONE )
3200  {
3201  for( I=1; I <= M; I++ )
3202  {
3203  I_ = I - 1;
3204  CC(J_,I_) *= BETA;
3205  }
3206  }
3207 
3208  for( L=1; L <= K; L++ )
3209  {
3210  L_ = L - 1;
3211  if( BB(L_,J_) != ZERO )
3212  {
3213  TEMP = ALPHA*BB(L_,J_);
3214  for( I=1; I <= M; I++ )
3215  {
3216  I_ = I - 1;
3217  CC(J_,I_) += TEMP*AA(L_,I_);
3218  }
3219  }
3220  }
3221  }
3222  }
3223 
3224  else
3225  {
3226 
3227  /* Form C := alpha*A'*B' + beta*C */
3228  for( J=1; J <= N; J++ )
3229  {
3230  J_ = J - 1;
3231 
3232  for( I=1; I <= M; I++ )
3233  {
3234  I_ = I - 1;
3235  TEMP = ZERO;
3236 
3237  for( L=1; L <= K; L++ )
3238  {
3239  L_ = L - 1;
3240  TEMP += AA(I_,L_)*BB(L_,J_);
3241  }
3242 
3243  if( BETA == ZERO )
3244  {
3245  CC(J_,I_) = ALPHA*TEMP;
3246  }
3247  else
3248  {
3249  CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
3250  }
3251 
3252  }
3253  }
3254  }
3255  }
3256 
3257  return;
3258 
3259  /* End of DGEMM .*/
3260 #undef C
3261 #undef B
3262 #undef A
3263 }
3264 #undef AA
3265 #undef BB
3266 #undef CC
3267 
3268 /* Subroutine */
3269 #if 0
3270 STATIC int32 DGTSV(int32 *n, int32 *nrhs, double *dl,
3271  double *d__, double *du, double *b, int32 *ldb, int32
3272  *info)
3273 {
3274  /* System generated locals */
3275  int32 b_dim1, b_offset, i__1, i__2;
3276 
3277  /* Local variables */
3278  double fact, temp;
3279  int32 i__, j;
3280 
3281 
3282 #define b_ref(a_1,a_2) b[(a_2)*(b_dim1) + (a_1)]
3283 
3284 
3285 /* -- LAPACK routine (version 3.0) --
3286  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3287  Courant Institute, Argonne National Lab, and Rice University
3288  October 31, 1999
3289 
3290 
3291  Purpose
3292  =======
3293 
3294  DGTSV solves the equation
3295 
3296  A*X = B,
3297 
3298  where A is an n by n tridiagonal matrix, by Gaussian elimination with
3299  partial pivoting.
3300 
3301  Note that the equation A'*X = B may be solved by interchanging the
3302  order of the arguments DU and DL.
3303 
3304  Arguments
3305  =========
3306 
3307  N (input) INTEGER
3308  The order of the matrix A. N >= 0.
3309 
3310  NRHS (input) INTEGER
3311  The number of right hand sides, i.e., the number of columns
3312  of the matrix B. NRHS >= 0.
3313 
3314  DL (input/output) DOUBLE PRECISION array, dimension (N-1)
3315  On entry, DL must contain the (n-1) sub-diagonal elements of
3316  A.
3317 
3318  On exit, DL is overwritten by the (n-2) elements of the
3319  second super-diagonal of the upper triangular matrix U from
3320  the LU factorization of A, in DL(1), ..., DL(n-2).
3321 
3322  D (input/output) DOUBLE PRECISION array, dimension (N)
3323  On entry, D must contain the diagonal elements of A.
3324 
3325  On exit, D is overwritten by the n diagonal elements of U.
3326 
3327  DU (input/output) DOUBLE PRECISION array, dimension (N-1)
3328  On entry, DU must contain the (n-1) super-diagonal elements
3329  of A.
3330 
3331  On exit, DU is overwritten by the (n-1) elements of the first
3332  super-diagonal of U.
3333 
3334  B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
3335  On entry, the N by NRHS matrix of right hand side matrix B.
3336  On exit, if INFO = 0, the N by NRHS solution matrix X.
3337 
3338  LDB (input) INTEGER
3339  The leading dimension of the array B. LDB >= max(1,N).
3340 
3341  INFO (output) INTEGER
3342  = 0: successful exit
3343  < 0: if INFO = -i, the i-th argument had an illegal value
3344  > 0: if INFO = i, U(i,i) is exactly zero, and the solution
3345  has not been computed. The factorization has not been
3346  completed unless i = N.
3347 
3348  =====================================================================
3349 
3350 
3351  Parameter adjustments */
3352  --dl;
3353  --d__;
3354  --du;
3355  b_dim1 = *ldb;
3356  b_offset = 1 + b_dim1 * 1;
3357  b -= b_offset;
3358 
3359  /* Function Body */
3360  *info = 0;
3361  if(*n < 0) {
3362  *info = -1;
3363  } else if(*nrhs < 0) {
3364  *info = -2;
3365  } else if(*ldb < *n && *ldb < 1) {
3366  *info = -7;
3367  }
3368  if(*info != 0) {
3369  i__1 = -(*info);
3370  XERBLA("DGTSV ", i__1);
3371  /* XERBLA does not return */
3372  }
3373 
3374  if(*n == 0) {
3375  return 0;
3376  }
3377 
3378  if(*nrhs == 1) {
3379  i__1 = *n - 2;
3380  for(i__ = 1; i__ <= i__1; ++i__) {
3381  if(fabs(d__[i__]) >= fabs(dl[i__])) {
3382 
3383 /* No row interchange required */
3384 
3385  if(d__[i__] != 0.) {
3386  fact = dl[i__] / d__[i__];
3387  d__[i__ + 1] -= fact * du[i__];
3388  b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3389  1);
3390  } else {
3391  *info = i__;
3392  return 0;
3393  }
3394  dl[i__] = 0.;
3395  } else {
3396 
3397 /* Interchange rows I and I+1 */
3398 
3399  fact = d__[i__] / dl[i__];
3400  d__[i__] = dl[i__];
3401  temp = d__[i__ + 1];
3402  d__[i__ + 1] = du[i__] - fact * temp;
3403  dl[i__] = du[i__ + 1];
3404  du[i__ + 1] = -fact * dl[i__];
3405  du[i__] = temp;
3406  temp = b_ref(i__, 1);
3407  b_ref(i__, 1) = b_ref(i__ + 1, 1);
3408  b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3409  }
3410 /* L10: */
3411  }
3412  if(*n > 1) {
3413  i__ = *n - 1;
3414  if(fabs(d__[i__]) >= fabs(dl[i__])) {
3415  if(d__[i__] != 0.) {
3416  fact = dl[i__] / d__[i__];
3417  d__[i__ + 1] -= fact * du[i__];
3418  b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3419  1);
3420  } else {
3421  *info = i__;
3422  return 0;
3423  }
3424  } else {
3425  fact = d__[i__] / dl[i__];
3426  d__[i__] = dl[i__];
3427  temp = d__[i__ + 1];
3428  d__[i__ + 1] = du[i__] - fact * temp;
3429  du[i__] = temp;
3430  temp = b_ref(i__, 1);
3431  b_ref(i__, 1) = b_ref(i__ + 1, 1);
3432  b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3433  }
3434  }
3435  if(d__[*n] == 0.) {
3436  *info = *n;
3437  return 0;
3438  }
3439  } else {
3440  i__1 = *n - 2;
3441  for(i__ = 1; i__ <= i__1; ++i__) {
3442  if(fabs(d__[i__]) >= fabs(dl[i__])) {
3443 
3444 /* No row interchange required */
3445 
3446  if(d__[i__] != 0.) {
3447  fact = dl[i__] / d__[i__];
3448  d__[i__ + 1] -= fact * du[i__];
3449  i__2 = *nrhs;
3450  for(j = 1; j <= i__2; ++j) {
3451  b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3452  i__, j);
3453 /* L20: */
3454  }
3455  } else {
3456  *info = i__;
3457  return 0;
3458  }
3459  dl[i__] = 0.;
3460  } else {
3461 
3462 /* Interchange rows I and I+1 */
3463 
3464  fact = d__[i__] / dl[i__];
3465  d__[i__] = dl[i__];
3466  temp = d__[i__ + 1];
3467  d__[i__ + 1] = du[i__] - fact * temp;
3468  dl[i__] = du[i__ + 1];
3469  du[i__ + 1] = -fact * dl[i__];
3470  du[i__] = temp;
3471  i__2 = *nrhs;
3472  for(j = 1; j <= i__2; ++j) {
3473  temp = b_ref(i__, j);
3474  b_ref(i__, j) = b_ref(i__ + 1, j);
3475  b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3476 /* L30: */
3477  }
3478  }
3479 /* L40: */
3480  }
3481  if(*n > 1) {
3482  i__ = *n - 1;
3483  if( fabs(d__[i__]) >= fabs(dl[i__]))
3484  {
3485  if(d__[i__] != 0.)
3486  {
3487  fact = dl[i__] / d__[i__];
3488  d__[i__ + 1] -= fact * du[i__];
3489  i__1 = *nrhs;
3490  for(j = 1; j <= i__1; ++j) {
3491  b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3492  i__, j);
3493  /* L50: */
3494  }
3495  }
3496  else
3497  {
3498  *info = i__;
3499  return 0;
3500  }
3501  } else {
3502  fact = d__[i__] / dl[i__];
3503  d__[i__] = dl[i__];
3504  temp = d__[i__ + 1];
3505  d__[i__ + 1] = du[i__] - fact * temp;
3506  du[i__] = temp;
3507  i__1 = *nrhs;
3508  for(j = 1; j <= i__1; ++j) {
3509  temp = b_ref(i__, j);
3510  b_ref(i__, j) = b_ref(i__ + 1, j);
3511  b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3512 /* L60: */
3513  }
3514  }
3515  }
3516  if(d__[*n] == 0.) {
3517  *info = *n;
3518  return 0;
3519  }
3520  }
3521 
3522 /* Back solve with the matrix U from the factorization. */
3523 
3524  if(*nrhs <= 2) {
3525  j = 1;
3526 L70:
3527  b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3528  if(*n > 1) {
3529  b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, j))
3530  / d__[*n - 1];
3531  }
3532  for(i__ = *n - 2; i__ >= 1; --i__) {
3533  b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) - dl[
3534  i__] * b_ref(i__ + 2, j)) / d__[i__];
3535 /* L80: */
3536  }
3537  if(j < *nrhs) {
3538  ++j;
3539  goto L70;
3540  }
3541  } else {
3542  i__1 = *nrhs;
3543  for(j = 1; j <= i__1; ++j) {
3544  b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3545  if(*n > 1) {
3546  b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n,
3547  j)) / d__[*n - 1];
3548  }
3549  for(i__ = *n - 2; i__ >= 1; --i__) {
3550  b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j)
3551  - dl[i__] * b_ref(i__ + 2, j)) / d__[i__];
3552 /* L90: */
3553  }
3554 /* L100: */
3555  }
3556  }
3557 
3558  return 0;
3559 
3560 /* End of DGTSV */
3561 
3562 } /* dgtsv_ */
3563 #endif
3564 #undef b_ref
3565 #endif
3566 /*lint +e725 expected pos indentation */
3567 /*lint +e801 goto is deprecated */

Generated for cloudy by doxygen 1.8.1.2