cloudy: /home66/gary/public_html/cloudy/c08_branch/source/thirdparty

00001 /* These are wrappers for lapack linear algebra routines.
00002  * There are two versions of the lapack routines - a fortran
00003  * version that I converted to C with forc to use if nothing else is available
00004  * (included in the Cloudy distribution),
00005  * and an option to link into an external lapack library that may be optimized
00006  * for your machine.  By default the tralated C routines will be used.
00007  * To use your machine's lapack library instead, define the macro
00008  * LAPACK and link into your library.  This is usually done with a command line
00009  * switch "-DLAPACK" on the compiler command, and the linker option "-llapack"
00010  */
00011 #include "cddefines.h"
00012 #include "thirdparty.h"
00013 /*lint -e725 expected pos indentation */
00014 /*lint -e801 goto is deprecated */
00015 
00016 #ifdef LAPACK
00017 /*********************The functions for FORTRAN version of the LAPACK calls *******************/
00018 /* dgetrf, dgetrs: lapack FORTRAN general full matrix solution using LU decomposition */
00019 
00020 extern "C" void dgetrf_(int32 *M, int32 *N, double *A, int32 *LDA, int32 *IPIV, int32 *INFO);
00021 extern "C" void dgetrs_(char *TRANS, int32 *N, int32 *NRHS, double *A, int32 *LDA, int32 *iPiv, double *B,
00022              int32 *LDB, int32 *INFO, int32 translen);
00023 extern "C" void dgtsv_(int32 *n, int32 *nrhs, double *dl, double *d, double *du, double *b, int32 *ldb, int32 *info);
00024 
00025 #else
00026 
00027 /*********************The functions for C version of the LAPACK calls *******************/
00028 /*
00029  * these are the routines that are, part of lapack, Some had their names slightly
00030  * changed so as to not conflict with the special lapack that exists on our exemplar.
00031  */
00032 
00033 /* DGETRF lapack, perform LU decomposition */
00034 STATIC void DGETRF(int32,int32,double*,int32,int32[],int32*);
00035 
00036 /* DGETRS lapack, solve linear system */
00037 STATIC void DGETRS(int32 TRANS,int32 N,int32 NRHS,double *A,int32 LDA,int32 IPIV[],double *B,int32 LDB,int32 *INFO);
00038 
00039 /* DGTSV lapack, solve tridiagonal system */
00040 /*static int32 DGTSV(int32 *n,int32 *nrhs,double *dl,double *d__,double *du,double *b,int32 *ldb,int32 *info);*/
00041 
00042 #endif
00043 
00044 
00045 /**********************************************************************************************/
00046 /* returns zero if successful termination */
00047 void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
00048 {
00049         if( *info == 0 )
00050         {
00051                 int32 M_loc, N_loc, lda_loc;
00052 
00053                 ASSERT( M < INT32_MAX && N < INT32_MAX && lda < INT32_MAX );
00054 
00055                 M_loc = (int32)M;
00056                 N_loc = (int32)N;
00057                 lda_loc = (int32)lda;
00058 
00059 #               ifdef  LAPACK
00060                 /* Calling the special version in library */
00061                 dgetrf_(&M_loc, &N_loc, A , &lda_loc, ipiv, info);
00062 #               else
00063                 /* Calling the old slower one, included with cloudy */
00064                 DGETRF(M_loc, N_loc, A, lda_loc, ipiv, info);
00065 #               endif
00066         }
00067 }
00068 
00069 void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, 
00070                    double *B, long ldb, int32 *info)
00071 {
00072         if( *info == 0 )
00073         {
00074                 int32 N_loc, nrhs_loc, lda_loc, ldb_loc;
00075 
00076                 ASSERT( N < INT32_MAX && nrhs < INT32_MAX && lda < INT32_MAX && ldb < INT32_MAX );
00077 
00078                 N_loc = (int32)N;
00079                 nrhs_loc = (int32)nrhs;
00080                 lda_loc = (int32)lda;
00081                 ldb_loc = (int32)ldb;
00082 
00083 #               ifdef LAPACK
00084                 /* Calling the special version in library */
00085                 dgetrs_(&trans, &N_loc, &nrhs_loc, A, &lda_loc, ipiv, B, &ldb_loc, info, sizeof(char));
00086 #               else
00087                 /* Calling the old slower one, included with cloudy */
00088                 DGETRS(trans, N_loc, nrhs_loc, A, lda_loc, ipiv, B, ldb_loc, info);
00089 #               endif
00090         }
00091 }
00092 
00093 #if 0
00094 void dgtsv_wrapper(long N, long nrhs, double *dl, double *d__, double *du, double *b, long ldb, int32 *info)
00095 {
00096         printf("Inside dgtsv\n");
00097         cdEXIT( EXIT_FAILURE );
00098         if( *info == 0 )
00099         {
00100                 int32 N_loc, nrhs_loc, ldb_loc;
00101 
00102                 ASSERT( N < INT32_MAX && nrhs < INT32_MAX && ldb < INT32_MAX );
00103 
00104                 N_loc = (int32)N;
00105                 nrhs_loc = (int32)nrhs;
00106                 ldb_loc = (int32)ldb;
00107 
00108 #               ifdef LAPACK
00109                 /* Calling the special version in library */
00110                 dgtsv_(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
00111 #               else
00112                 /* Calling the old slower one, included with cloudy */
00113                 /* DGTSV always returns zero, so it is safe to ignore the return value */
00114                 (void)DGTSV(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
00115 #               endif
00116         }
00117 }
00118 #endif
00119 
00120 
00121 #ifndef LAPACK
00122 
00123 #define ONE     1.0e0
00124 #define ZERO    0.0e0
00125 
00126 #ifdef AA
00127 #       undef AA
00128 #endif
00129 #ifdef BB
00130 #       undef BB
00131 #endif
00132 #ifdef CC
00133 #       undef CC
00134 #endif
00135 
00136 #define AA(I_,J_)       (*(A+(I_)*(LDA)+(J_)))
00137 #define BB(I_,J_)       (*(B+(I_)*(LDB)+(J_)))
00138 #define CC(I_,J_)       (*(C+(I_)*(LDC)+(J_)))
00139 
00140 /* 
00141  * these are the routines that are, part of lapack, Some had their names slightly
00142  * changed so as to not conflict with the special lapack that exits on our exemplar.
00143  */
00144 
00145 /* dgtsv, dgetrf, dgetrs: lapack general tridiagonal solution */
00146 /*int32 DGTSV(int32 *n, int32 *nrhs, double *dl, 
00147         double *d__, double *du, double *b, int32 *ldb, int32 *info);*/
00148 
00149 /* DGETRF lapack service routine */
00150 /*void DGETRF(int32,int32,double*,int32,int32[],int32*);*/
00151 
00152 /*DGETRS lapack matrix inversion routine */
00153 /*void DGETRS(int32 TRANS, 
00154           int32 N, 
00155           int32 NRHS, 
00156           double *A, 
00157           int32 LDA, 
00158           int32 IPIV[], 
00159           double *B, 
00160           int32 LDB, 
00161           int32 *INFO);
00162 */
00163 /* DGEMM matrix inversion helper routine*/
00164 STATIC void DGEMM(int32 TRANSA, 
00165           int32 TRANSB, 
00166           int32 M, 
00167           int32 N, 
00168           int32 K, 
00169           double ALPHA, 
00170           double *A, 
00171           int32 LDA, 
00172           double *B, 
00173           int32 LDB, 
00174           double BETA, 
00175           double *C, 
00176           int32 LDC);
00177 
00178 /*LSAME LAPACK auxiliary routine  */
00179 STATIC int32 LSAME(int32 CA, 
00180           int32 CB);
00181 
00182 /*IDAMAX lapack service routine */
00183 STATIC int32 IDAMAX(int32 n, 
00184           double dx[], 
00185           int32 incx);
00186 
00187 /*DTRSM lapack service routine */
00188 STATIC void DTRSM(int32 SIDE, 
00189           int32 UPLO, 
00190           int32 TRANSA, 
00191           int32 DIAG, 
00192           int32 M, 
00193           int32 N, 
00194           double ALPHA, 
00195           double *A, 
00196           int32 LDA, 
00197           double *B, 
00198           int32 LDB);
00199 
00200 /* ILAENV lapack helper routine */
00201 STATIC int32 ILAENV(int32 ISPEC, 
00202           const char *NAME, 
00203           /*char *OPTS, */
00204           int32 N1, 
00205           int32 N2, 
00206           /*int32 N3, */
00207           int32 N4);
00208 
00209 /*DSWAP lapack routine */
00210 STATIC void DSWAP(int32 n, 
00211           double dx[], 
00212           int32 incx, 
00213           double dy[], 
00214           int32 incy);
00215 
00216 /*DSCAL lapack routine */
00217 STATIC void DSCAL(int32 n, 
00218           double da, 
00219           double dx[], 
00220           int32 incx);
00221 
00222 /*DLASWP  -- LAPACK auxiliary routine (version 2.0) --*/
00223 STATIC void DLASWP(int32 N, 
00224           double *A, 
00225           int32 LDA, 
00226           int32 K1, 
00227           int32 K2, 
00228           int32 IPIV[], 
00229           int32 INCX);
00230 
00231 /*DGETF2 lapack service routine */
00232 STATIC void DGETF2(int32 M, 
00233           int32 N, 
00234           double *A, 
00235           int32 LDA, 
00236           int32 IPIV[], 
00237           int32 *INFO);
00238 
00239 /*DGER service routine for matrix inversion */
00240 STATIC void DGER(int32 M, 
00241           int32 N, 
00242           double ALPHA, 
00243           double X[], 
00244           int32 INCX, 
00245           double Y[], 
00246           int32 INCY, 
00247           double *A, 
00248           int32 LDA);
00249 
00250 /*XERBLA  -- LAPACK auxiliary routine (version 2.0) -- */
00251 STATIC void XERBLA(const char *SRNAME, 
00252                    int32 INFO)
00253 {
00254 
00255         DEBUG_ENTRY( "XERBLA()" );
00256 
00257         /*  -- LAPACK auxiliary routine (version 2.0) --
00258          * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00259          * Courant Institute, Argonne National Lab, and Rice University
00260          * September 30, 1994
00261          *
00262          * .. Scalar Arguments .. */
00263         /* ..
00264          *
00265          *  Purpose
00266          *  =======
00267          *
00268          *  XERBLA  is an error handler for the LAPACK routines.
00269          *  It is called by an LAPACK routine if an input parameter has an
00270          *  invalid value.  A message is printed and execution stops.
00271          *
00272          *  Installers may consider modifying the STOP statement in order to
00273          *  call system-specific exception-handling facilities.
00274          *
00275          *  Arguments
00276          *  =========
00277          *
00278          *  SRNAME  (input) CHARACTER*6
00279          *  The name of the routine which called XERBLA.
00280          *
00281          *  INFO    (input) INTEGER
00282          *  The position of the invalid parameter in the parameter list
00283          *  of the calling routine.
00284          *
00285          * =====================================================================
00286          *
00287          * .. Executable Statements ..
00288          * */
00289         fprintf( ioQQQ, " ** On entry to %6.6s parameter number %2ld had an illegal value\n",
00290                  SRNAME, (long)INFO );
00291 
00292         cdEXIT(EXIT_FAILURE);
00293 }
00294 
00295 
00296 STATIC void DGETRF(
00297                 /* number of rows of the matrix */
00298           int32 M, 
00299                 /* number of columns of the matrix
00300                  * M=N for square matrix */
00301           int32 N, 
00302           /* double precision matrix */
00303           double *A, 
00304           /* LDA is right dimension of matrix */
00305           int32 LDA, 
00306           /* following must dimensions the smaller of M or N */
00307           int32 IPIV[], 
00308           /* following is zero for successful exit */
00309           int32 *INFO)
00310 {
00311 
00312         char chL1, 
00313           chL2, 
00314           chL3, 
00315           chL4;
00316         int32 I, 
00317           IINFO, 
00318           I_, 
00319           J, 
00320           JB, 
00321           J_, 
00322           NB, 
00323 
00324           limit, 
00325           limit2;
00326         /*double _d_l;*/
00327 
00328         DEBUG_ENTRY( "DGETRF()" );
00329 
00330         /*  -- LAPACK routine (version 2.0) --
00331          * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00332          * Courant Institute, Argonne National Lab, and Rice University
00333          * March 31, 1993
00334          *
00335          *  Purpose
00336          *  =======
00337          *
00338          *  DGETRF computes an LU factorization of a general M-by-N matrix A
00339          *  using partial pivoting with row interchanges.
00340          *
00341          *  The factorization has the form
00342          * A = P * L * U
00343          *  where P is a permutation matrix, L is lower triangular with unit
00344          *  diagonal elements (lower trapezoidal if m > n), and U is upper
00345          *  triangular (upper trapezoidal if m < n).
00346          *
00347          *  This is the right-looking Level 3 BLAS version of the algorithm.
00348          *
00349          *  Arguments
00350          *  =========
00351          *
00352          *  M       (input) INTEGER
00353          *  The number of rows of the matrix A.  M >= 0.
00354          *
00355          *  N       (input) INTEGER
00356          *  The number of columns of the matrix A.  N >= 0.
00357          *
00358          *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00359          *  On entry, the M-by-N matrix to be factored.
00360          *  On exit, the factors L and U from the factorization
00361          *  A = P*L*U; the unit diagonal elements of L are not stored.
00362          *
00363          *  LDA     (input) INTEGER
00364          *  The leading dimension of the array A.  LDA >= MAX(1,M).
00365          *
00366          *  IPIV    (output) INTEGER array, dimension (MIN(M,N))
00367          *  The pivot indices; for 1 <= i <= MIN(M,N), row i of the
00368          *  matrix was interchanged with row IPIV(i).
00369          *
00370          *  INFO    (output) INTEGER
00371          *  = 0:  successful exit
00372          *  < 0:  if INFO = -i, the i-th argument had an illegal value
00373          *  > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
00374          *      has been completed, but the factor U is exactly
00375          *      singular, and division by zero will occur if it is used
00376          *      to solve a system of equations.
00377          *
00378          *  =====================================================================
00379          *
00380          * .. Parameters .. */
00381         /* ..
00382          * .. Local Scalars .. */
00383         /* ..
00384          * .. External Subroutines .. */
00385         /* ..
00386          * .. External Functions .. */
00387         /* ..
00388          * .. Intrinsic Functions .. */
00389         /* ..
00390          * .. Executable Statements ..
00391          *
00392          * Test the input parameters.
00393          * */
00394         *INFO = 0;
00395         if( M < 0 )
00396         {
00397                 *INFO = -1;
00398         }
00399         else if( N < 0 )
00400         {
00401                 *INFO = -2;
00402         }
00403         else if( LDA < MAX2(1,M) )
00404         {
00405                 *INFO = -4;
00406         }
00407         if( *INFO != 0 )
00408         {
00409                 XERBLA("DGETRF",-*INFO);
00410                 /* XERBLA does not return */
00411         }
00412 
00413         /* Quick return if possible
00414          * */
00415         if( M == 0 || N == 0 )
00416         { 
00417                 return;
00418         }
00419 
00420         /* Determine the block size for this environment.
00421          * */
00422                 /* >>chng 01 oct 22, remove two parameters since not used */
00423         NB = ILAENV(1,"DGETRF",/*" ",*/M,N,/*-1,*/-1);
00424         if( NB <= 1 || NB >= MIN2(M,N) )
00425         {
00426                 /*  Use unblocked code.
00427                  * */
00428                 DGETF2(M,N,A,LDA,IPIV,INFO);
00429         }
00430         else
00431         {
00432 
00433                 /*  Use blocked code.
00434                  * */
00435                 limit = MIN2(M,N);
00436                 /*for( J=1, _do0=DOCNT(J,limit,_do1 = NB); _do0 > 0; J += _do1, _do0-- )*/
00437                 /*do J = 1, limit , NB */
00438                 for( J=1; J<=limit; J += NB )
00439                 {
00440                         J_ = J - 1;
00441                         JB = MIN2(MIN2(M,N)-J+1,NB);
00442 
00443                         /* Factor diagonal and subdiagonal blocks and test for exact
00444                          * singularity.
00445                          * */
00446                         DGETF2(M-J+1,JB,&AA(J_,J_),LDA,&IPIV[J_],&IINFO);
00447 
00448                         /* Adjust INFO and the pivot indices.
00449                          * */
00450                         if( *INFO == 0 && IINFO > 0 )
00451                                 *INFO = IINFO + J - 1;
00452                         limit2 = MIN2(M,J+JB-1);
00453                         for( I=J; I <= limit2; I++ )
00454                         {
00455                                 I_ = I - 1;
00456                                 IPIV[I_] += J - 1;
00457                         }
00458 
00459                         /* Apply interchanges to columns 1:J-1.
00460                          * */
00461                         DLASWP(J-1,A,LDA,J,J+JB-1,IPIV,1);
00462 
00463                         if( J + JB <= N )
00464                         {
00465 
00466                                 /*    Apply interchanges to columns J+JB:N.
00467                                  * */
00468                                 DLASWP(N-J-JB+1,&AA(J_+JB,0),LDA,J,J+JB-1,IPIV,1);
00469 
00470                                 /*    Compute block row of U.
00471                                  * */
00472                                 chL1 = 'L';
00473                                 chL2 = 'L';
00474                                 chL3 = 'N';
00475                                 chL4 = 'U';
00476                                 /*    CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, */
00477                                 DTRSM(chL1,chL2,chL3,chL4,JB,N-J-JB+1,ONE,&AA(J_,J_),
00478                                   LDA,&AA(J_+JB,J_),LDA);
00479                                 if( J + JB <= M )
00480                                 {
00481 
00482                                         /*       Update trailing submatrix.
00483                                          * */
00484                                         chL1 = 'N';
00485                                         chL2 = 'N';
00486                                         /*       CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1, */
00487                                         DGEMM(chL1,chL2,M-J-JB+1,N-J-JB+1,JB,-ONE,&AA(J_,J_+JB),
00488                                           LDA,&AA(J_+JB,J_),LDA,ONE,&AA(J_+JB,J_+JB),LDA);
00489                                 }
00490                         }
00491                 }
00492         }
00493         return;
00494 
00495         /* End of DGETRF
00496          * */
00497 #undef  A
00498 }
00499 
00500 /*DGETRS lapack matrix inversion routine */
00501 /*****************************************************************
00502  *****************************************************************
00503  *
00504  * matrix inversion routine
00505  *
00506  * solves Ax=B  A is an nXn matrix, C and B are nX1 matrices
00507  * dim A(n,n), B(n,1)  C overwrites B.
00508  * integer ipiv(n)
00509  * integer info  see below: 
00510  *  INFO    (output) INTEGER
00511  *  = 0:  successful exit
00512  *  < 0:  if INFO = -i, the i-th argument had an illegal value
00513  *  > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
00514  *      has been completed, but the factor U is exactly
00515  *      singular, and division by zero will occur if it is used
00516  *      to solve a system of equations.
00517  *
00518  *
00519  *
00520  *must call the routines in the following order:
00521  * call dgetrf(n,n,A,n,ipiv,info)
00522  * call dgetrs('N',n,1,A,n,ipiv,B,n,info)
00523  *
00524  *
00525  *************************************************************************** */
00526 
00527 STATIC void DGETRS(
00528           /* 1 ch var saying what to do */
00529           int32 TRANS, 
00530           /* order of the matrix */
00531           int32 N, 
00532           /* number of right hand sides */
00533           int32 NRHS, 
00534           /* double [N][LDA] */
00535           double *A, 
00536           /* second dim of array */
00537           int32 LDA, 
00538           /* helper vector, dimensioned N*/
00539           int32 IPIV[], 
00540           /* on input the ri=hs vector, on output, the result */
00541           double *B, 
00542           /* dimension of B, 1 if one vector */
00543           int32 LDB, 
00544           /* = 0 if ok */
00545           int32 *INFO)
00546 {
00547 /*#define A(I_,J_)      (*(A+(I_)*(LDA)+(J_)))*/
00548 /*#define B(I_,J_)      (*(B+(I_)*(LDB)+(J_)))*/
00549         int32 NOTRAN;
00550         char chL1, 
00551           chL2, 
00552           chL3, 
00553           chL4;
00554 
00555         DEBUG_ENTRY( "DGETRS()" );
00556 
00557         /*  -- LAPACK routine (version 2.0) --
00558          * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00559          * Courant Institute, Argonne National Lab, and Rice University
00560          * March 31, 1993
00561          *
00562          *
00563          *  Purpose
00564          *  =======
00565          *
00566          *  DGETRS solves a system of linear equations
00567          * A * X = B  or  A' * X = B
00568          *  with a general N-by-N matrix A using the LU factorization computed
00569          *  by DGETRF.
00570          *
00571          *  Arguments
00572          *  =========
00573          *
00574          *  TRANS   (input) CHARACTER*1
00575          *  Specifies the form of the system of equations:
00576          *  = 'N':  A * X = B  (No transpose)
00577          *  = 'T':  A'* X = B  (Transpose)
00578          *  = 'C':  A'* X = B  (Conjugate transpose = Transpose)
00579          *
00580          *  N       (input) INTEGER
00581          *  The order of the matrix A.  N >= 0.
00582          *
00583          *  NRHS    (input) INTEGER
00584          *  The number of right hand sides, i.e., the number of columns
00585          *  of the matrix B.  NRHS >= 0.
00586          *
00587          *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00588          *  The factors L and U from the factorization A = P*L*U
00589          *  as computed by DGETRF.
00590          *
00591          *  LDA     (input) INTEGER
00592          *  The leading dimension of the array A.  LDA >= MAX(1,N).
00593          *
00594          *  IPIV    (input) INTEGER array, dimension (N)
00595          *  The pivot indices from DGETRF; for 1<=i<=N, row i of the
00596          *  matrix was interchanged with row IPIV(i).
00597          *
00598          *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00599          *  On entry, the right hand side matrix B.
00600          *  On exit, the solution matrix X.
00601          *
00602          *  LDB     (input) INTEGER
00603          *  The leading dimension of the array B.  LDB >= MAX(1,N).
00604          *
00605          *  INFO    (output) INTEGER
00606          *  = 0:  successful exit
00607          *  < 0:  if INFO = -i, the i-th argument had an illegal value
00608          *
00609          *  =====================================================================
00610          *
00611          * .. Parameters .. */
00612         /* ..
00613          * .. Local Scalars .. */
00614         /* ..
00615          * .. External Functions .. */
00616         /* ..
00617          * .. External Subroutines .. */
00618         /* ..
00619          * .. Intrinsic Functions .. */
00620         /* ..
00621          * .. Executable Statements ..
00622          *
00623          * Test the input parameters.
00624          * */
00625         *INFO = 0;
00626         NOTRAN = LSAME(TRANS,'N');
00627         if( (!NOTRAN && !LSAME(TRANS,'T')) && !LSAME(TRANS,'C') )
00628         {
00629                 *INFO = -1;
00630         }
00631         else if( N < 0 )
00632         {
00633                 *INFO = -2;
00634         }
00635         else if( NRHS < 0 )
00636         {
00637                 *INFO = -3;
00638         }
00639         else if( LDA < MAX2(1,N) )
00640         {
00641                 *INFO = -5;
00642         }
00643         else if( LDB < MAX2(1,N) )
00644         {
00645                 *INFO = -8;
00646         }
00647         if( *INFO != 0 )
00648         {
00649                 XERBLA("DGETRS",-*INFO);
00650                 /* XERBLA does not return */
00651         }
00652 
00653         /* Quick return if possible
00654          * */
00655         if( N == 0 || NRHS == 0 )
00656         { 
00657                 return;
00658         }
00659 
00660         if( NOTRAN )
00661         {
00662 
00663                 /*  Solve A * X = B.
00664                  *
00665                  *  Apply row interchanges to the right hand sides.
00666                  * */
00667                 DLASWP(NRHS,B,LDB,1,N,IPIV,1);
00668 
00669                 /*  Solve L*X = B, overwriting B with X.
00670                  * */
00671                 chL1 = 'L';
00672                 chL2 = 'L';
00673                 chL3 = 'N';
00674                 chL4 = 'U';
00675                 /*  CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, */
00676                 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00677 
00678                 /*  Solve U*X = B, overwriting B with X.
00679                  * */
00680                 chL1 = 'L';
00681                 chL2 = 'U';
00682                 chL3 = 'N';
00683                 chL4 = 'N';
00684                 /*  CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, */
00685                 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00686         }
00687         else
00688         {
00689 
00690                 /*  Solve A' * X = B.
00691                  *
00692                  *  Solve U'*X = B, overwriting B with X.
00693                  * */
00694                 chL1 = 'L';
00695                 chL2 = 'U';
00696                 chL3 = 'T';
00697                 chL4 = 'N';
00698                 /*  CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, */
00699                 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00700 
00701                 /*  Solve L'*X = B, overwriting B with X.
00702                  * */
00703                 chL1 = 'L';
00704                 chL2 = 'L';
00705                 chL3 = 'T';
00706                 chL4 = 'U';
00707                 /*  CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE, */
00708                 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00709 
00710                 /*  Apply row interchanges to the solution vectors.
00711                  * */
00712                 DLASWP(NRHS,B,LDB,1,N,IPIV,-1);
00713         }
00714 
00715         return;
00716 
00717         /* End of DGETRS
00718          * */
00719 #undef  B
00720 #undef  A
00721 }
00722 
00723 /*LSAME LAPACK auxiliary routine  */
00724 
00725 STATIC int32 LSAME(int32 CA, 
00726           int32 CB)
00727 {
00728         /*
00729          *
00730          *  -- LAPACK auxiliary routine (version 2.0) --
00731          * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00732          * Courant Institute, Argonne National Lab, and Rice University
00733          * September 30, 1994
00734          *
00735          * .. Scalar Arguments ..
00736          */
00737         int32 LSAME_v;
00738         int32 INTA, 
00739           INTB, 
00740           ZCODE;
00741 
00742         DEBUG_ENTRY( "LSAME()" );
00743         /* ..
00744          *
00745          *  Purpose
00746          *  =======
00747          *
00748          *  LSAME returns .true. if CA is the same letter as CB regardless of
00749          *  case.
00750          *
00751          *  Arguments
00752          *  =========
00753          *
00754          *  CA      (input) CHARACTER*1
00755          *  CB      (input) CHARACTER*1
00756          *  CA and CB specify the single characters to be compared.
00757          *
00758          * =====================================================================
00759          *
00760          * .. Intrinsic Functions .. */
00761         /* ..
00762          * .. Local Scalars .. */
00763         /* ..
00764          * .. Executable Statements ..
00765          *
00766          * Test if the characters are equal
00767          * */
00768         LSAME_v = CA == CB;
00769         if( LSAME_v )
00770         { 
00771                 return LSAME_v;
00772         }
00773 
00774         /* Now test for equivalence if both characters are alphabetic.
00775          * */
00776         ZCODE = 'Z';
00777 
00778         /* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
00779          * machines, on which ICHAR returns a value with bit 8 set.
00780          * ICHAR('A') on Prime machines returns 193 which is the same as
00781          * ICHAR('A') on an EBCDIC machine.
00782          * */
00783         INTA = (CA);
00784         INTB = (CB);
00785 
00786         if( ZCODE == 90 || ZCODE == 122 )
00787         {
00788 
00789                 /*  ASCII is assumed - ZCODE is the ASCII code of either lower or
00790                  *  upper case 'Z'.
00791                  * */
00792                 if( INTA >= 97 && INTA <= 122 )
00793                         INTA -= 32;
00794                 if( INTB >= 97 && INTB <= 122 )
00795                         INTB -= 32;
00796 
00797         }
00798         else if( ZCODE == 233 || ZCODE == 169 )
00799         {
00800 
00801                 /*  EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
00802                  *  upper case 'Z'.
00803                  * */
00804                 if( ((INTA >= 129 && INTA <= 137) || (INTA >= 145 && INTA <= 
00805                   153)) || (INTA >= 162 && INTA <= 169) )
00806                         INTA += 64;
00807                 if( ((INTB >= 129 && INTB <= 137) || (INTB >= 145 && INTB <= 
00808                   153)) || (INTB >= 162 && INTB <= 169) )
00809                         INTB += 64;
00810 
00811         }
00812         else if( ZCODE == 218 || ZCODE == 250 )
00813         {
00814 
00815                 /*  ASCII is assumed, on Prime machines - ZCODE is the ASCII code
00816                  *  plus 128 of either lower or upper case 'Z'.
00817                  * */
00818                 if( INTA >= 225 && INTA <= 250 )
00819                         INTA -= 32;
00820                 if( INTB >= 225 && INTB <= 250 )
00821                         INTB -= 32;
00822         }
00823         LSAME_v = INTA == INTB;
00824 
00825         /* RETURN
00826          *
00827          * End of LSAME
00828          * */
00829         return LSAME_v;
00830 }
00831 
00832 /*IDAMAX lapack service routine */
00833 
00834 STATIC int32 IDAMAX(int32 n, 
00835           double dx[], 
00836           int32 incx)
00837 {
00838         /*
00839          * finds the index of element having max. absolute value.
00840          * jack dongarra, lapack, 3/11/78.
00841          * modified 3/93 to return if incx .le. 0.
00842          * modified 12/3/93, array(1) declarations changed to array(*)
00843          *
00844          */
00845         int32 IDAMAX_v, 
00846           i, 
00847           ix;
00848         double dmax;
00849 
00850         DEBUG_ENTRY( "IDAMAX()" );
00851 
00852         IDAMAX_v = 0;
00853 
00854         if( n < 1 || incx <= 0 )
00855         { 
00856                 return IDAMAX_v;
00857         }
00858 
00859         IDAMAX_v = 1;
00860 
00861         if( n == 1 )
00862         { 
00863                 return IDAMAX_v;
00864         }
00865 
00866         if( incx == 1 )
00867                 goto L_20;
00868 
00869         /*  code for increment not equal to 1
00870          * */
00871         ix = 1;
00872         dmax = fabs(dx[0]);
00873         ix += incx;
00874         for( i=2; i <= n; i++ )
00875         {
00876                 /*  if(ABS(dx(ix)).le.dmax) go to 5 */
00877                 if( fabs(dx[ix-1]) > dmax )
00878                 {
00879                         IDAMAX_v = i;
00880                         dmax = fabs(dx[ix-1]);
00881                 }
00882                 ix += incx;
00883         }
00884         return IDAMAX_v;
00885 
00886         /*  code for increment equal to 1
00887          * */
00888 L_20:
00889         dmax = fabs(dx[0]);
00890         for( i=1; i < n; i++ )
00891         {
00892                 /*  if(ABS(dx(i)).le.dmax) go to 30 */
00893                 if( fabs(dx[i]) > dmax )
00894                 {
00895                         IDAMAX_v = i+1;
00896                         dmax = fabs(dx[i]);
00897                 }
00898         }
00899         return IDAMAX_v;
00900 }
00901 
00902 /*DTRSM lapack service routine */
00903 STATIC void DTRSM(int32 SIDE, 
00904           int32 UPLO, 
00905           int32 TRANSA, 
00906           int32 DIAG, 
00907           int32 M, 
00908           int32 N, 
00909           double ALPHA, 
00910           double *A, 
00911           int32 LDA, 
00912           double *B, 
00913           int32 LDB)
00914 {
00915         int32 LSIDE, 
00916           NOUNIT, 
00917           UPPER;
00918         int32 I, 
00919           INFO, 
00920           I_, 
00921           J, 
00922           J_, 
00923           K, 
00924           K_, 
00925           NROWA;
00926         double TEMP;
00927 
00928         DEBUG_ENTRY( "DTRSM()" );
00929         /* .. Scalar Arguments .. */
00930         /* .. Array Arguments .. */
00931         /* ..
00932          *
00933          *  Purpose
00934          *  =======
00935          *
00936          *  DTRSM  solves one of the matrix equations
00937          *
00938          * op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
00939          *
00940          *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
00941          *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
00942          *
00943          * op( A ) = A   or   op( A ) = A'.
00944          *
00945          *  The matrix X is overwritten on B.
00946          *
00947          *  Parameters
00948          *  ==========
00949          *
00950          *  SIDE   - CHARACTER*1.
00951          * On entry, SIDE specifies whether op( A ) appears on the left
00952          * or right of X as follows:
00953          *
00954          *    SIDE = 'L' or 'l'   op( A )*X = alpha*B.
00955          *
00956          *    SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
00957          *
00958          * Unchanged on exit.
00959          *
00960          *  UPLO   - CHARACTER*1.
00961          * On entry, UPLO specifies whether the matrix A is an upper or
00962          * lower triangular matrix as follows:
00963          *
00964          *    UPLO = 'U' or 'u'   A is an upper triangular matrix.
00965          *
00966          *    UPLO = 'L' or 'l'   A is a lower triangular matrix.
00967          *
00968          * Unchanged on exit.
00969          *
00970          *  TRANSA - CHARACTER*1.
00971          * On entry, TRANSA specifies the form of op( A ) to be used in
00972          * the matrix multiplication as follows:
00973          *
00974          *    TRANSA = 'N' or 'n'   op( A ) = A.
00975          *
00976          *    TRANSA = 'T' or 't'   op( A ) = A'.
00977          *
00978          *    TRANSA = 'C' or 'c'   op( A ) = A'.
00979          *
00980          * Unchanged on exit.
00981          *
00982          *  DIAG   - CHARACTER*1.
00983          * On entry, DIAG specifies whether or not A is unit triangular
00984          * as follows:
00985          *
00986          *    DIAG = 'U' or 'u'   A is assumed to be unit triangular.
00987          *
00988          *    DIAG = 'N' or 'n'   A is not assumed to be unit
00989          *                        triangular.
00990          *
00991          * Unchanged on exit.
00992          *
00993          *  M      - INTEGER.
00994          * On entry, M specifies the number of rows of B. M must be at
00995          * least zero.
00996          * Unchanged on exit.
00997          *
00998          *  N      - INTEGER.
00999          * On entry, N specifies the number of columns of B.  N must be
01000          * at least zero.
01001          * Unchanged on exit.
01002          *
01003          *  ALPHA  - DOUBLE PRECISION.
01004          * On entry,  ALPHA specifies the scalar  alpha. When  alpha is
01005          * zero then  A is not referenced and  B need not be set before
01006          * entry.
01007          * Unchanged on exit.
01008          *
01009          *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
01010          * when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
01011          * Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
01012          * upper triangular part of the array  A must contain the upper
01013          * triangular matrix  and the strictly lower triangular part of
01014          * A is not referenced.
01015          * Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
01016          * lower triangular part of the array  A must contain the lower
01017          * triangular matrix  and the strictly upper triangular part of
01018          * A is not referenced.
01019          * Note that when  DIAG = 'U' or 'u',  the diagonal elements of
01020          * A  are not referenced either,  but are assumed to be  unity.
01021          * Unchanged on exit.
01022          *
01023          *  LDA    - INTEGER.
01024          * On entry, LDA specifies the first dimension of A as declared
01025          * in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
01026          * LDA  must be at least  MAX( 1, m ),  when  SIDE = 'R' or 'r'
01027          * then LDA must be at least MAX( 1, n ).
01028          * Unchanged on exit.
01029          *
01030          *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
01031          * Before entry,  the leading  m by n part of the array  B must
01032          * contain  the  right-hand  side  matrix  B,  and  on exit  is
01033          * overwritten by the solution matrix  X.
01034          *
01035          *  LDB    - INTEGER.
01036          * On entry, LDB specifies the first dimension of B as declared
01037          * in  the  calling  (sub)  program.   LDB  must  be  at  least
01038          * MAX( 1, m ).
01039          * Unchanged on exit.
01040          *
01041          *
01042          *  Level 3 Blas routine.
01043          *
01044          *
01045          *  -- Written on 8-February-1989.
01046          * Jack Dongarra, Argonne National Laboratory.
01047          * Iain Duff, AERE Harwell.
01048          * Jeremy Du Croz, Numerical Algorithms Group Ltd.
01049          * Sven Hammarling, Numerical Algorithms Group Ltd.
01050          *
01051          *
01052          * .. External Functions .. */
01053         /* .. External Subroutines .. */
01054         /* .. Intrinsic Functions .. */
01055         /* .. Local Scalars .. */
01056         /* .. Parameters .. */
01057         /* ..
01058          * .. Executable Statements ..
01059          *
01060          * Test the input parameters.
01061          * */
01062         LSIDE = LSAME(SIDE,'L');
01063         if( LSIDE )
01064         {
01065                 NROWA = M;
01066         }
01067         else
01068         {
01069                 NROWA = N;
01070         }
01071         NOUNIT = LSAME(DIAG,'N');
01072         UPPER = LSAME(UPLO,'U');
01073 
01074         INFO = 0;
01075         if( (!LSIDE) && (!LSAME(SIDE,'R')) )
01076         {
01077                 INFO = 1;
01078         }
01079         else if( (!UPPER) && (!LSAME(UPLO,'L')) )
01080         {
01081                 INFO = 2;
01082         }
01083         else if( ((!LSAME(TRANSA,'N')) && (!LSAME(TRANSA,'T'))) && (!LSAME(TRANSA,
01084           'C')) )
01085         {
01086                 INFO = 3;
01087         }
01088         else if( (!LSAME(DIAG,'U')) && (!LSAME(DIAG,'N')) )
01089         {
01090                 INFO = 4;
01091         }
01092         else if( M < 0 )
01093         {
01094                 INFO = 5;
01095         }
01096         else if( N < 0 )
01097         {
01098                 INFO = 6;
01099         }
01100         else if( LDA < MAX2(1,NROWA) )
01101         {
01102                 INFO = 9;
01103         }
01104         else if( LDB < MAX2(1,M) )
01105         {
01106                 INFO = 11;
01107         }
01108         if( INFO != 0 )
01109         {
01110                 XERBLA("DTRSM ",INFO);
01111                 /* XERBLA does not return */
01112         }
01113 
01114         /* Quick return if possible.
01115          * */
01116         if( N == 0 )
01117                 { 
01118                 return;}
01119 
01120         /* And when  alpha.eq.zero.
01121          * */
01122         if( ALPHA == ZERO )
01123         {
01124                 for( J=1; J <= N; J++ )
01125                 {
01126                         J_ = J - 1;
01127                         for( I=1; I <= M; I++ )
01128                         {
01129                                 I_ = I - 1;
01130                                 BB(J_,I_) = ZERO;
01131                         }
01132                 }
01133                 return;
01134         }
01135 
01136         /* Start the operations.
01137          * */
01138         if( LSIDE )
01139         {
01140                 if( LSAME(TRANSA,'N') )
01141                 {
01142 
01143                         /* Form  B := alpha*inv( A )*B.
01144                          * */
01145                         if( UPPER )
01146                         {
01147                                 for( J=1; J <= N; J++ )
01148                                 {
01149                                         J_ = J - 1;
01150                                         if( ALPHA != ONE )
01151                                         {
01152                                                 for( I=1; I <= M; I++ )
01153                                                 {
01154                                                         I_ = I - 1;
01155                                                         BB(J_,I_) *= ALPHA;
01156                                                 }
01157                                         }
01158                                         for( K=M; K >= 1; K-- )
01159                                         {
01160                                                 K_ = K - 1;
01161                                                 if( BB(J_,K_) != ZERO )
01162                                                 {
01163                                                         if( NOUNIT )
01164                                                                 BB(J_,K_) /= AA(K_,K_);
01165                                                         for( I=1; I <= (K - 1); I++ )
01166                                                         {
01167                                                                 I_ = I - 1;
01168                                                                 BB(J_,I_) += -BB(J_,K_)*AA(K_,I_);
01169                                                         }
01170                                                 }
01171                                         }
01172                                 }
01173                         }
01174                         else
01175                         {
01176                                 for( J=1; J <= N; J++ )
01177                                 {
01178                                         J_ = J - 1;
01179                                         if( ALPHA != ONE )
01180                                         {
01181                                                 for( I=1; I <= M; I++ )
01182                                                 {
01183                                                         I_ = I - 1;
01184                                                         BB(J_,I_) *= ALPHA;
01185                                                 }
01186                                         }
01187                                         for( K=1; K <= M; K++ )
01188                                         {
01189                                                 K_ = K - 1;
01190                                                 if( BB(J_,K_) != ZERO )
01191                                                 {
01192                                                         if( NOUNIT )
01193                                                                 BB(J_,K_) /= AA(K_,K_);
01194                                                         for( I=K + 1; I <= M; I++ )
01195                                                         {
01196                                                                 I_ = I - 1;
01197                                                                 BB(J_,I_) += -BB(J_,K_)*AA(K_,I_);
01198                                                         }
01199                                                 }
01200                                         }
01201                                 }
01202                         }
01203                 }
01204                 else
01205                 {
01206 
01207                         /* Form  B := alpha*inv( A' )*B.
01208                          * */
01209                         if( UPPER )
01210                         {
01211                                 for( J=1; J <= N; J++ )
01212                                 {
01213                                         J_ = J - 1;
01214                                         for( I=1; I <= M; I++ )
01215                                         {
01216                                                 I_ = I - 1;
01217                                                 TEMP = ALPHA*BB(J_,I_);
01218                                                 for( K=1; K <= (I - 1); K++ )
01219                                                 {
01220                                                         K_ = K - 1;
01221                                                         TEMP += -AA(I_,K_)*BB(J_,K_);
01222                                                 }
01223                                                 if( NOUNIT )
01224                                                         TEMP /= AA(I_,I_);
01225                                                 BB(J_,I_) = TEMP;
01226                                         }
01227                                 }
01228                         }
01229                         else
01230                         {
01231                                 for( J=1; J <= N; J++ )
01232                                 {
01233                                         J_ = J - 1;
01234                                         for( I=M; I >= 1; I-- )
01235                                         {
01236                                                 I_ = I - 1;
01237                                                 TEMP = ALPHA*BB(J_,I_);
01238                                                 for( K=I + 1; K <= M; K++ )
01239                                                 {
01240                                                         K_ = K - 1;
01241                                                         TEMP += -AA(I_,K_)*BB(J_,K_);
01242                                                 }
01243                                                 if( NOUNIT )
01244                                                         TEMP /= AA(I_,I_);
01245                                                 BB(J_,I_) = TEMP;
01246                                         }
01247                                 }
01248                         }
01249                 }
01250         }
01251         else
01252         {
01253                 if( LSAME(TRANSA,'N') )
01254                 {
01255 
01256                         /* Form  B := alpha*B*inv( A ).
01257                          * */
01258                         if( UPPER )
01259                         {
01260                                 for( J=1; J <= N; J++ )
01261                                 {
01262                                         J_ = J - 1;
01263                                         if( ALPHA != ONE )
01264                                         {
01265                                                 for( I=1; I <= M; I++ )
01266                                                 {
01267                                                         I_ = I - 1;
01268                                                         BB(J_,I_) *= ALPHA;
01269                                                 }
01270                                         }
01271                                         for( K=1; K <= (J - 1); K++ )
01272                                         {
01273                                                 K_ = K - 1;
01274                                                 if( AA(J_,K_) != ZERO )
01275                                                 {
01276                                                         for( I=1; I <= M; I++ )
01277                                                         {
01278                                                                 I_ = I - 1;
01279                                                                 BB(J_,I_) += -AA(J_,K_)*BB(K_,I_);
01280                                                         }
01281                                                 }
01282                                         }
01283                                         if( NOUNIT )
01284                                         {
01285                                                 TEMP = ONE/AA(J_,J_);
01286                                                 for( I=1; I <= M; I++ )
01287                                                 {
01288                                                         I_ = I - 1;
01289                                                         BB(J_,I_) *= TEMP;
01290                                                 }
01291                                         }
01292                                 }
01293                         }
01294                         else
01295                         {
01296                                 for( J=N; J >= 1; J-- )
01297                                 {
01298                                         J_ = J - 1;
01299                                         if( ALPHA != ONE )
01300                                         {
01301                                                 for( I=1; I <= M; I++ )
01302                                                 {
01303                                                         I_ = I - 1;
01304                                                         BB(J_,I_) *= ALPHA;
01305                                                 }
01306                                         }
01307                                         for( K=J + 1; K <= N; K++ )
01308                                         {
01309                                                 K_ = K - 1;
01310                                                 if( AA(J_,K_) != ZERO )
01311                                                 {
01312                                                         for( I=1; I <= M; I++ )
01313                                                         {
01314                                                                 I_ = I - 1;
01315                                                                 BB(J_,I_) += -AA(J_,K_)*BB(K_,I_);
01316                                                         }
01317                                                 }
01318                                         }
01319                                         if( NOUNIT )
01320                                         {
01321                                                 TEMP = ONE/AA(J_,J_);
01322                                                 for( I=1; I <= M; I++ )
01323                                                 {
01324                                                         I_ = I - 1;
01325                                                         BB(J_,I_) *= TEMP;
01326                                                 }
01327                                         }
01328                                 }
01329                         }
01330                 }
01331                 else
01332                 {
01333 
01334                         /* Form  B := alpha*B*inv( A' ).
01335                          * */
01336                         if( UPPER )
01337                         {
01338                                 for( K=N; K >= 1; K-- )
01339                                 {
01340                                         K_ = K - 1;
01341                                         if( NOUNIT )
01342                                         {
01343                                                 TEMP = ONE/AA(K_,K_);
01344                                                 for( I=1; I <= M; I++ )
01345                                                 {
01346                                                         I_ = I - 1;
01347                                                         BB(K_,I_) *= TEMP;
01348                                                 }
01349                                         }
01350                                         for( J=1; J <= (K - 1); J++ )
01351                                         {
01352                                                 J_ = J - 1;
01353                                                 if( AA(K_,J_) != ZERO )
01354                                                 {
01355                                                         TEMP = AA(K_,J_);
01356                                                         for( I=1; I <= M; I++ )
01357                                                         {
01358                                                                 I_ = I - 1;
01359                                                                 BB(J_,I_) += -TEMP*BB(K_,I_);
01360                                                         }
01361                                                 }
01362                                         }
01363                                         if( ALPHA != ONE )
01364                                         {
01365                                                 for( I=1; I <= M; I++ )
01366                                                 {
01367                                                         I_ = I - 1;
01368                                                         BB(K_,I_) *= ALPHA;
01369                                                 }
01370                                         }
01371                                 }
01372                         }
01373                         else
01374                         {
01375                                 for( K=1; K <= N; K++ )
01376                                 {
01377                                         K_ = K - 1;
01378                                         if( NOUNIT )
01379                                         {
01380                                                 TEMP = ONE/AA(K_,K_);
01381                                                 for( I=1; I <= M; I++ )
01382                                                 {
01383                                                         I_ = I - 1;
01384                                                         BB(K_,I_) *= TEMP;
01385                                                 }
01386                                         }
01387                                         for( J=K + 1; J <= N; J++ )
01388                                         {
01389                                                 J_ = J - 1;
01390                                                 if( AA(K_,J_) != ZERO )
01391                                                 {
01392                                                         TEMP = AA(K_,J_);
01393                                                         for( I=1; I <= M; I++ )
01394                                                         {
01395                                                                 I_ = I - 1;
01396                                                                 BB(J_,I_) += -TEMP*BB(K_,I_);
01397                                                         }
01398                                                 }
01399                                         }
01400                                         if( ALPHA != ONE )
01401                                         {
01402                                                 for( I=1; I <= M; I++ )
01403                                                 {
01404                                                         I_ = I - 1;
01405                                                         BB(K_,I_) *= ALPHA;
01406                                                 }
01407                                         }
01408                                 }
01409                         }
01410                 }
01411         }
01412 
01413         return;
01414 
01415         /* End of DTRSM .
01416          * */
01417 #undef  B
01418 #undef  A
01419 }
01420 
01421 /* ILAENV lapack helper routine */
01422 
01423 /* >>chng 01 oct 22, remove two parameters since not used */
01424 STATIC int32 ILAENV(int32 ISPEC, 
01425           const char *NAME, 
01426           /*char *OPTS, */
01427           int32 N1, 
01428           int32 N2, 
01429           /*int32 N3, */
01430           int32 N4)
01431 {
01432         /*
01433          *
01434          *  -- LAPACK auxiliary routine (version 2.0) --
01435          * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
01436          * Courant Institute, Argonne National Lab, and Rice University
01437          * September 30, 1994
01438          *
01439          * .. Scalar Arguments ..
01440          */
01441         char C2[3], 
01442           C3[4], 
01443           C4[3], 
01444           SUBNAM[7];
01445         int32 CNAME, 
01446           SNAME;
01447         char C1;
01448         int32 I, 
01449           IC, 
01450           ILAENV_v, 
01451           IZ, 
01452           NB, 
01453           NBMIN, 
01454           NX;
01455 
01456         DEBUG_ENTRY( "ILAENV()" );
01457         /* ..
01458          *
01459          *  Purpose
01460          *  =======
01461          *
01462          *  ILAENV is called from the LAPACK routines to choose problem-dependent
01463          *  parameters for the local environment.  See ISPEC for a description of
01464          *  the parameters.
01465          *
01466          *  This version provides a set of parameters which should give good,
01467          *  but not optimal, performance on many of the currently available
01468          *  computers.  Users are encouraged to modify this subroutine to set
01469          *  the tuning parameters for their particular machine using the option
01470          *  and problem size information in the arguments.
01471          *
01472          *  This routine will not function correctly if it is converted to all
01473          *  lower case.  Converting it to all upper case is allowed.
01474          *
01475          *  Arguments
01476          *  =========
01477          *
01478          *  ISPEC   (input) INTEGER
01479          *  Specifies the parameter to be returned as the value of
01480          *  ILAENV.
01481          *  = 1: the optimal blocksize; if this value is 1, an unblocked
01482          *     algorithm will give the best performance.
01483          *  = 2: the minimum block size for which the block routine
01484          *     should be used; if the usable block size is less than
01485          *     this value, an unblocked routine should be used.
01486          *  = 3: the crossover point (in a block routine, for N less
01487          *     than this value, an unblocked routine should be used)
01488          *  = 4: the number of shifts, used in the nonsymmetric
01489          *     eigenvalue routines
01490          *  = 5: the minimum column dimension for blocking to be used;
01491          *     rectangular blocks must have dimension at least k by m,
01492          *     where k is given by ILAENV(2,...) and m by ILAENV(5,...)
01493          *  = 6: the crossover point for the SVD (when reducing an m by n
01494          *     matrix to bidiagonal form, if MAX(m,n)/MIN(m,n) exceeds
01495          *     this value, a QR factorization is used first to reduce
01496          *     the matrix to a triangular form.)
01497          *  = 7: the number of processors
01498          *  = 8: the crossover point for the multishift QR and QZ methods
01499          *     for nonsymmetric eigenvalue problems.
01500          *
01501          *  NAME    (input) CHARACTER*(*)
01502          *  The name of the calling subroutine, in either upper case or
01503          *  lower case.
01504          *
01505          *  OPTS    (input) CHARACTER*(*)
01506          *  The character options to the subroutine NAME, concatenated
01507          *  into a single character string.  For example, UPLO = 'U',
01508          *  TRANS = 'T', and DIAG = 'N' for a triangular routine would
01509          *  be specified as OPTS = 'UTN'.
01510          *
01511          *  N1      (input) INTEGER
01512          *  N2      (input) INTEGER
01513          *  N3      (input) INTEGER
01514          *  N4      (input) INTEGER
01515          *  Problem dimensions for the subroutine NAME; these may not all
01516          *  be required.
01517          *
01518          * (ILAENV) (output) INTEGER
01519          *  >= 0: the value of the parameter specified by ISPEC
01520          *  < 0:  if ILAENV = -k, the k-th argument had an illegal value.
01521          *
01522          *  Further Details
01523          *  ===============
01524          *
01525          *  The following conventions have been used when calling ILAENV from the
01526          *  LAPACK routines:
01527          *  1)  OPTS is a concatenation of all of the character options to
01528          *  subroutine NAME, in the same order that they appear in the
01529          *  argument list for NAME, even if they are not used in determining
01530          *  the value of the parameter specified by ISPEC.
01531          *  2)  The problem dimensions N1, N2, N3, N4 are specified in the order
01532          *  that they appear in the argument list for NAME.  N1 is used
01533          *  first, N2 second, and so on, and unused problem dimensions are
01534          *  passed a value of -1.
01535          *  3)  The parameter value returned by ILAENV is checked for validity in
01536          *  the calling subroutine.  For example, ILAENV is used to retrieve
01537          *  the optimal blocksize for STRTRI as follows:
01538          *
01539          *  NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
01540          *  IF( NB.LE.1 ) NB = MAX( 1, N )
01541          *
01542          *  =====================================================================
01543          *
01544          * .. Local Scalars .. */
01545         /* ..
01546          * .. Intrinsic Functions .. */
01547         /* ..
01548          * .. Executable Statements ..
01549          * */
01550         switch( ISPEC )
01551         {
01552                 case 1: 
01553                         {
01554                                 goto L_100;
01555                         }
01556                 case 2: 
01557                         {
01558                                 goto L_100;
01559                         }
01560                 case 3: 
01561                         {
01562                                 goto L_100;
01563                         }
01564                 case 4: 
01565                         {
01566                                 goto L_400;
01567                         }
01568                 case 5: 
01569                         {
01570                                 goto L_500;
01571                         }
01572                 case 6: 
01573                         {
01574                                 goto L_600;
01575                         }
01576                 case 7: 
01577                         {
01578                                 goto L_700;
01579                         }
01580                 case 8: 
01581                         {
01582                                 goto L_800;
01583                         }
01584                 /* this is impossible, added by gjf to stop lint from complaining */
01585                 default: 
01586                 {
01587                         /* Invalid value for ISPEC
01588                          * */
01589                         ILAENV_v = -1;
01590                         return ILAENV_v;
01591                 }
01592         }
01593 
01594 L_100:
01595 
01596         /* Convert NAME to upper case if the first character is lower case.
01597          * */
01598         ILAENV_v = 1;
01599         strncpy( SUBNAM, NAME, 6 );
01600         IC = (SUBNAM[0]);
01601         IZ = 'Z';
01602         if( IZ == 90 || IZ == 122 )
01603         {
01604 
01605                 /*  ASCII character set
01606                  * */
01607                 if( IC >= 97 && IC <= 122 )
01608                 {
01609                         SUBNAM[0] = (char)(IC - 32);
01610                         for( I=2; I <= 6; I++ )
01611                         {
01612                                 IC = (SUBNAM[I-1]);
01613                                 if( IC >= 97 && IC <= 122 )
01614                                         SUBNAM[I - 1] = (char)(IC - 32);
01615                         }
01616                 }
01617 
01618         }
01619         else if( IZ == 233 || IZ == 169 )
01620         {
01621 
01622                 /*  EBCDIC character set
01623                  * */
01624                 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <= 153)) || 
01625                   (IC >= 162 && IC <= 169) )
01626                 {
01627                         SUBNAM[0] = (char)(IC + 64);
01628                         for( I=2; I <= 6; I++ )
01629                         {
01630                                 IC = (SUBNAM[I-1]);
01631                                 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <= 
01632                                   153)) || (IC >= 162 && IC <= 169) )
01633                                         SUBNAM[I - 1] = (char)(IC + 64);
01634                         }
01635                 }
01636 
01637         }
01638         else if( IZ == 218 || IZ == 250 )
01639         {
01640 
01641                 /*  Prime machines:  ASCII+128
01642                  * */
01643                 if( IC >= 225 && IC <= 250 )
01644                 {
01645                         SUBNAM[0] = (char)(IC - 32);
01646                         for( I=2; I <= 6; I++ )
01647                         {
01648                                 IC = (SUBNAM[I-1]);
01649                                 if( IC >= 225 && IC <= 250 )
01650                                         SUBNAM[I - 1] = (char)(IC - 32);
01651                         }
01652                 }
01653         }
01654 
01655         C1 = SUBNAM[0];
01656         SNAME = C1 == 'S' || C1 == 'D';
01657         CNAME = C1 == 'C' || C1 == 'Z';
01658         if( !(CNAME || SNAME) )
01659         { 
01660                 return ILAENV_v;
01661         }
01662 
01663 #       if 0
01664         strncpy( C2, SUBNAM+1, 2 );
01665         strncpy( C3, SUBNAM+3, 3 );
01666         strncpy( C4, C3+1, 2 );
01667 #       endif
01668 
01669         /* >>chng 00 nov 08, from above to below, insure had run
01670          * off end of string*/
01671         strncpy( C2, SUBNAM+1, 2 );
01672         C2[2] = '\0'; 
01673         strncpy( C3, SUBNAM+3, 3 );
01674         C3[3] = '\0'; 
01675         strncpy( C4, C3+1, 2 );
01676         C4[2] = '\0'; 
01677 
01678         switch( ISPEC )
01679         {
01680                 case 1: goto L_110;
01681                 case 2: goto L_200;
01682                 case 3: goto L_300;
01683         }
01684 
01685 L_110:
01686 
01687         /* ISPEC = 1:  block size
01688          *
01689          * In these examples, separate code is provided for setting NB for
01690          * real and complex.  We assume that NB will take the same value in
01691          * single or double precision.
01692          * */
01693         NB = 1;
01694 
01695         if( strcmp(C2,"GE") == 0 )
01696         {
01697                 if( strcmp(C3,"TRF") == 0 )
01698                 {
01699                         if( SNAME )
01700                         {
01701                                 NB = 64;
01702                         }
01703                         else
01704                         {
01705                                 NB = 64;
01706                         }
01707                 }
01708                 else if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || 
01709                   strcmp(C3,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
01710                 {
01711                         if( SNAME )
01712                         {
01713                                 NB = 32;
01714                         }
01715                         else
01716                         {
01717                                 NB = 32;
01718                         }
01719                 }
01720                 else if( strcmp(C3,"HRD") == 0 )
01721                 {
01722                         if( SNAME )
01723                         {
01724                                 NB = 32;
01725                         }
01726                         else
01727                         {
01728                                 NB = 32;
01729                         }
01730                 }
01731                 else if( strcmp(C3,"BRD") == 0 )
01732                 {
01733                         if( SNAME )
01734                         {
01735                                 NB = 32;
01736                         }
01737                         else
01738                         {
01739                                 NB = 32;
01740                         }
01741                 }
01742                 else if( strcmp(C3,"TRI") == 0 )
01743                 {
01744                         if( SNAME )
01745                         {
01746                                 NB = 64;
01747                         }
01748                         else
01749                         {
01750                                 NB = 64;
01751                         }
01752                 }
01753         }
01754         else if( strcmp(C2,"PO") == 0 )
01755         {
01756                 if( strcmp(C3,"TRF") == 0 )
01757                 {
01758                         if( SNAME )
01759                         {
01760                                 NB = 64;
01761                         }
01762                         else
01763                         {
01764                                 NB = 64;
01765                         }
01766                 }
01767         }
01768         else if( strcmp(C2,"SY") == 0 )
01769         {
01770                 if( strcmp(C3,"TRF") == 0 )
01771                 {
01772                         if( SNAME )
01773                         {
01774                                 NB = 64;
01775                         }
01776                         else
01777                         {
01778                                 NB = 64;
01779                         }
01780                 }
01781                 else if( SNAME && strcmp(C3,"TRD") == 0 )
01782                 {
01783                         NB = 1;
01784                 }
01785                 else if( SNAME && strcmp(C3,"GST") == 0 )
01786                 {
01787                         NB = 64;
01788                 }
01789         }
01790         else if( CNAME && strcmp(C2,"HE") == 0 )
01791         {
01792                 if( strcmp(C3,"TRF") == 0 )
01793                 {
01794                         NB = 64;
01795                 }
01796                 else if( strcmp(C3,"TRD") == 0 )
01797                 {
01798                         NB = 1;
01799                 }
01800                 else if( strcmp(C3,"GST") == 0 )
01801                 {
01802                         NB = 64;
01803                 }
01804         }
01805         else if( SNAME && strcmp(C2,"OR") == 0 )
01806         {
01807                 if( C3[0] == 'G' )
01808                 {
01809                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
01810                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01811                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
01812                           0 )
01813                         {
01814                                 NB = 32;
01815                         }
01816                 }
01817                 else if( C3[0] == 'M' )
01818                 {
01819                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
01820                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01821                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
01822                           0 )
01823                         {
01824                                 NB = 32;
01825                         }
01826                 }
01827         }
01828         else if( CNAME && strcmp(C2,"UN") == 0 )
01829         {
01830                 if( C3[0] == 'G' )
01831                 {
01832                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
01833                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01834                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
01835                           0 )
01836                         {
01837                                 NB = 32;
01838                         }
01839                 }
01840                 else if( C3[0] == 'M' )
01841                 {
01842                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
01843                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01844                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
01845                           0 )
01846                         {
01847                                 NB = 32;
01848                         }
01849                 }
01850         }
01851         else if( strcmp(C2,"GB") == 0 )
01852         {
01853                 if( strcmp(C3,"TRF") == 0 )
01854                 {
01855                         if( SNAME )
01856                         {
01857                                 if( N4 <= 64 )
01858                                 {
01859                                         NB = 1;
01860                                 }
01861                                 else
01862                                 {
01863                                         NB = 32;
01864                                 }
01865                         }
01866                         else
01867                         {
01868                                 if( N4 <= 64 )
01869                                 {
01870                                         NB = 1;
01871                                 }
01872                                 else
01873                                 {
01874                                         NB = 32;
01875                                 }
01876                         }
01877                 }
01878         }
01879         else if( strcmp(C2,"PB") == 0 )
01880         {
01881                 if( strcmp(C3,"TRF") == 0 )
01882                 {
01883                         if( SNAME )
01884                         {
01885                                 if( N2 <= 64 )
01886                                 {
01887                                         NB = 1;
01888                                 }
01889                                 else
01890                                 {
01891                                         NB = 32;
01892                                 }
01893                         }
01894                         else
01895                         {
01896                                 if( N2 <= 64 )
01897                                 {
01898                                         NB = 1;
01899                                 }
01900                                 else
01901                                 {
01902                                         NB = 32;
01903                                 }
01904                         }
01905                 }
01906         }
01907         else if( strcmp(C2,"TR") == 0 )
01908         {
01909                 if( strcmp(C3,"TRI") == 0 )
01910                 {
01911                         if( SNAME )
01912                         {
01913                                 NB = 64;
01914                         }
01915                         else
01916                         {
01917                                 NB = 64;
01918                         }
01919                 }
01920         }
01921         else if( strcmp(C2,"LA") == 0 )
01922         {
01923                 if( strcmp(C3,"UUM") == 0 )
01924                 {
01925                         if( SNAME )
01926                         {
01927                                 NB = 64;
01928                         }
01929                         else
01930                         {
01931                                 NB = 64;
01932                         }
01933                 }
01934         }
01935         else if( SNAME && strcmp(C2,"ST") == 0 )
01936         {
01937                 if( strcmp(C3,"EBZ") == 0 )
01938                 {
01939                         NB = 1;
01940                 }
01941         }
01942         ILAENV_v = NB;
01943         return ILAENV_v;
01944 
01945 L_200:
01946 
01947         /* ISPEC = 2:  minimum block size
01948          * */
01949         NBMIN = 2;
01950         if( strcmp(C2,"GE") == 0 )
01951         {
01952                 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
01953                   ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
01954                 {
01955                         if( SNAME )
01956                         {
01957                                 NBMIN = 2;
01958                         }
01959                         else
01960                         {
01961                                 NBMIN = 2;
01962                         }
01963                 }
01964                 else if( strcmp(C3,"HRD") == 0 )
01965                 {
01966                         if( SNAME )
01967                         {
01968                                 NBMIN = 2;
01969                         }
01970                         else
01971                         {
01972                                 NBMIN = 2;
01973                         }
01974                 }
01975                 else if( strcmp(C3,"BRD") == 0 )
01976                 {
01977                         if( SNAME )
01978                         {
01979                                 NBMIN = 2;
01980                         }
01981                         else
01982                         {
01983                                 NBMIN = 2;
01984                         }
01985                 }
01986                 else if( strcmp(C3,"TRI") == 0 )
01987                 {
01988                         if( SNAME )
01989                         {
01990                                 NBMIN = 2;
01991                         }
01992                         else
01993                         {
01994                                 NBMIN = 2;
01995                         }
01996                 }
01997         }
01998         else if( strcmp(C2,"SY") == 0 )
01999         {
02000                 if( strcmp(C3,"TRF") == 0 )
02001                 {
02002                         if( SNAME )
02003                         {
02004                                 NBMIN = 8;
02005                         }
02006                         else
02007                         {
02008                                 NBMIN = 8;
02009                         }
02010                 }
02011                 else if( SNAME && strcmp(C3,"TRD") == 0 )
02012                 {
02013                         NBMIN = 2;
02014                 }
02015         }
02016         else if( CNAME && strcmp(C2,"HE") == 0 )
02017         {
02018                 if( strcmp(C3,"TRD") == 0 )
02019                 {
02020                         NBMIN = 2;
02021                 }
02022         }
02023         else if( SNAME && strcmp(C2,"OR") == 0 )
02024         {
02025                 if( C3[0] == 'G' )
02026                 {
02027                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
02028                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02029                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
02030                           0 )
02031                         {
02032                                 NBMIN = 2;
02033                         }
02034                 }
02035                 else if( C3[0] == 'M' )
02036                 {
02037                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
02038                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02039                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
02040                           0 )
02041                         {
02042                                 NBMIN = 2;
02043                         }
02044                 }
02045         }
02046         else if( CNAME && strcmp(C2,"UN") == 0 )
02047         {
02048                 if( C3[0] == 'G' )
02049                 {
02050                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
02051                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02052                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
02053                           0 )
02054                         {
02055                                 NBMIN = 2;
02056                         }
02057                 }
02058                 else if( C3[0] == 'M' )
02059                 {
02060                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
02061                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02062                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
02063                           0 )
02064                         {
02065                                 NBMIN = 2;
02066                         }
02067                 }
02068         }
02069         ILAENV_v = NBMIN;
02070         return ILAENV_v;
02071 
02072 L_300:
02073 
02074         /* ISPEC = 3:  crossover point
02075          * */
02076         NX = 0;
02077         if( strcmp(C2,"GE") == 0 )
02078         {
02079                 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
02080                   ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
02081                 {
02082                         if( SNAME )
02083                         {
02084                                 NX = 128;
02085                         }
02086                         else
02087                         {
02088                                 NX = 128;
02089                         }
02090                 }
02091                 else if( strcmp(C3,"HRD") == 0 )
02092                 {
02093                         if( SNAME )
02094                         {
02095                                 NX = 128;
02096                         }
02097                         else
02098                         {
02099                                 NX = 128;
02100                         }
02101                 }
02102                 else if( strcmp(C3,"BRD") == 0 )
02103                 {
02104                         if( SNAME )
02105                         {
02106                                 NX = 128;
02107                         }
02108                         else
02109                         {
02110                                 NX = 128;
02111                         }
02112                 }
02113         }
02114         else if( strcmp(C2,"SY") == 0 )
02115         {
02116                 if( SNAME && strcmp(C3,"TRD") == 0 )
02117                 {
02118                         NX = 1;
02119                 }
02120         }
02121         else if( CNAME && strcmp(C2,"HE") == 0 )
02122         {
02123                 if( strcmp(C3,"TRD") == 0 )
02124                 {
02125                         NX = 1;
02126                 }
02127         }
02128         else if( SNAME && strcmp(C2,"OR") == 0 )
02129         {
02130                 if( C3[0] == 'G' )
02131                 {
02132                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
02133                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02134                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
02135                           0 )
02136                         {
02137                                 NX = 128;
02138                         }
02139                 }
02140         }
02141         else if( CNAME && strcmp(C2,"UN") == 0 )
02142         {
02143                 if( C3[0] == 'G' )
02144                 {
02145                         if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) || 
02146                           strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02147                           ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") == 
02148                           0 )
02149                         {
02150                                 NX = 128;
02151                         }
02152                 }
02153         }
02154         ILAENV_v = NX;
02155         return ILAENV_v;
02156 
02157 L_400:
02158 
02159         /* ISPEC = 4:  number of shifts (used by xHSEQR)
02160          * */
02161         ILAENV_v = 6;
02162         return ILAENV_v;
02163 
02164 L_500:
02165 
02166         /* ISPEC = 5:  minimum column dimension (not used)
02167          * */
02168         ILAENV_v = 2;
02169         return ILAENV_v;
02170 
02171 L_600:
02172 
02173         /* ISPEC = 6:  crossover point for SVD (used by xGELSS and xGESVD)
02174          * */
02175         ILAENV_v = (int32)((realnum)(MIN2(N1,N2))*1.6e0);
02176         return ILAENV_v;
02177 
02178 L_700:
02179 
02180         /* ISPEC = 7:  number of processors (not used)
02181          * */
02182         ILAENV_v = 1;
02183         return ILAENV_v;
02184 
02185 L_800:
02186 
02187         /* ISPEC = 8:  crossover point for multishift (used by xHSEQR)
02188          * */
02189         ILAENV_v = 50;
02190         return ILAENV_v;
02191 
02192         /* End of ILAENV
02193          * */
02194 }
02195 
02196 /*DSWAP lapack routine */
02197 
02198 STATIC void DSWAP(int32 n, 
02199           double dx[], 
02200           int32 incx, 
02201           double dy[], 
02202           int32 incy)
02203 {
02204         int32 i, 
02205           ix, 
02206           iy, 
02207           m;
02208         double dtemp;
02209 
02210         DEBUG_ENTRY( "DSWAP()" );
02211 
02212         /* interchanges two vectors.
02213          * uses unrolled loops for increments equal one.
02214          * jack dongarra, lapack, 3/11/78.
02215          * modified 12/3/93, array(1) declarations changed to array(*)
02216          * */
02217 
02218         if( n <= 0 )
02219                 { 
02220                 return;}
02221         if( incx == 1 && incy == 1 )
02222                 goto L_20;
02223 
02224         /* code for unequal increments or equal increments not equal
02225          * to 1
02226          * */
02227         ix = 1;
02228         iy = 1;
02229 
02230         if( incx < 0 )
02231                 ix = (-n + 1)*incx + 1;
02232 
02233         if( incy < 0 )
02234                 iy = (-n + 1)*incy + 1;
02235 
02236         for( i=0; i < n; i++ )
02237         {
02238                 dtemp = dx[ix-1];
02239                 dx[ix-1] = dy[iy-1];
02240                 dy[iy-1] = dtemp;
02241                 ix += incx;
02242                 iy += incy;
02243         }
02244         return;
02245 
02246         /* code for both increments equal to 1
02247          *
02248          *
02249          * clean-up loop
02250          * */
02251 L_20:
02252         m = n%3;
02253         if( m == 0 )
02254                 goto L_40;
02255 
02256         for( i=0; i < m; i++ )
02257         {
02258                 dtemp = dx[i];
02259                 dx[i] = dy[i];
02260                 dy[i] = dtemp;
02261         }
02262 
02263         if( n < 3 )
02264         { 
02265                 return;
02266         }
02267 
02268 L_40:
02269         for( i=m; i < n; i += 3 )
02270         {
02271                 dtemp = dx[i];
02272                 dx[i] = dy[i];
02273                 dy[i] = dtemp;
02274                 dtemp = dx[i+1];
02275                 dx[i+1] = dy[i+1];
02276                 dy[i+1] = dtemp;
02277                 dtemp = dx[i+2];
02278                 dx[i+2] = dy[i+2];
02279                 dy[i+2] = dtemp;
02280         }
02281         return;
02282 }
02283 
02284 /*DSCAL lapack routine */
02285 
02286 STATIC void DSCAL(int32 n, 
02287           double da, 
02288           double dx[], 
02289           int32 incx)
02290 {
02291         int32 i, 
02292           m, 
02293           nincx;
02294 
02295         DEBUG_ENTRY( "DSCAL()" );
02296 
02297         /* scales a vector by a constant.
02298          * uses unrolled loops for increment equal to one.
02299          * jack dongarra, lapack, 3/11/78.
02300          * modified 3/93 to return if incx .le. 0.
02301          * modified 12/3/93, array(1) declarations changed to array(*)
02302          * */
02303 
02304         if( n <= 0 || incx <= 0 )
02305                 { 
02306                 return;}
02307         if( incx == 1 )
02308                 goto L_20;
02309 
02310         /*  code for increment not equal to 1
02311          * */
02312         nincx = n*incx;
02313         /*for( i=1, _do0=DOCNT(i,nincx,_do1 = incx); _do0 > 0; i += _do1, _do0-- )*/
02314         for( i=0; i<nincx; i = i + incx)
02315         {
02316                 dx[i] *= da;
02317         }
02318         return;
02319 
02320         /*  code for increment equal to 1
02321          *
02322          *
02323          *  clean-up loop
02324          * */
02325 L_20:
02326         m = n%5;
02327         if( m == 0 )
02328                 goto L_40;
02329 
02330         for( i=0; i < m; i++ )
02331         {
02332                 dx[i] *= da;
02333         }
02334 
02335         if( n < 5 )
02336         { 
02337                 return;
02338         }
02339 
02340 L_40:
02341 
02342         for( i=m; i < n; i += 5 )
02343         {
02344                 dx[i] *= da;
02345                 dx[i+1] *= da;
02346                 dx[i+2] *= da;
02347                 dx[i+3] *= da;
02348                 dx[i+4] *= da;
02349         }
02350         return;
02351 }
02352 
02353 /*DLASWP  -- LAPACK auxiliary routine (version 2.0) --*/
02354 
02355 STATIC void DLASWP(int32 N, 
02356           double *A, 
02357           int32 LDA, 
02358           int32 K1, 
02359           int32 K2, 
02360           int32 IPIV[], 
02361           int32 INCX)
02362 {
02363         int32 I, 
02364           IP, 
02365           IX, 
02366           I_;
02367 
02368         DEBUG_ENTRY( "DLASWP()" );
02369 
02370         /*  -- LAPACK auxiliary routine (version 2.0) --
02371          * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
02372          * Courant Institute, Argonne National Lab, and Rice University
02373          * October 31, 1992
02374          *
02375          * .. Scalar Arguments .. */
02376         /* ..
02377          * .. Array Arguments .. */
02378         /* ..
02379          *
02380          *  Purpose
02381          *  =======
02382          *
02383          *  DLASWP performs a series of row interchanges on the matrix A.
02384          *  One row interchange is initiated for each of rows K1 through K2 of A.
02385          *
02386          *  Arguments
02387          *  =========
02388          *
02389          *  N       (input) INTEGER
02390          *  The number of columns of the matrix A.
02391          *
02392          *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
02393          *  On entry, the matrix of column dimension N to which the row
02394          *  interchanges will be applied.
02395          *  On exit, the permuted matrix.
02396          *
02397          *  LDA     (input) INTEGER
02398          *  The leading dimension of the array A.
02399          *
02400          *  K1      (input) INTEGER
02401          *  The first element of IPIV for which a row interchange will
02402          *  be done.
02403          *
02404          *  K2      (input) INTEGER
02405          *  The last element of IPIV for which a row interchange will
02406          *  be done.
02407          *
02408          *  IPIV    (input) INTEGER array, dimension (M*ABS(INCX))
02409          *  The vector of pivot indices.  Only the elements in positions
02410          *  K1 through K2 of IPIV are accessed.
02411          *  IPIV(K) = L implies rows K and L are to be interchanged.
02412          *
02413          *  INCX    (input) INTEGER
02414          *  The increment between successive values of IPIV.  If IPIV
02415          *  is negative, the pivots are applied in reverse order.
02416          *
02417          * =====================================================================
02418          *
02419          * .. Local Scalars .. */
02420         /* ..
02421          * .. External Subroutines .. */
02422         /* ..
02423          * .. Executable Statements ..
02424          *
02425          * Interchange row I with row IPIV(I) for each of rows K1 through K2.
02426          * */
02427         if( INCX == 0 )
02428                 { 
02429                 return;}
02430         if( INCX > 0 )
02431         {
02432                 IX = K1;
02433         }
02434         else
02435         {
02436                 IX = 1 + (1 - K2)*INCX;
02437         }
02438         if( INCX == 1 )
02439         {
02440                 for( I=K1; I <= K2; I++ )
02441                 {
02442                         I_ = I - 1;
02443                         IP = IPIV[I_];
02444                         if( IP != I )
02445                                 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02446                 }
02447         }
02448         else if( INCX > 1 )
02449         {
02450                 for( I=K1; I <= K2; I++ )
02451                 {
02452                         I_ = I - 1;
02453                         IP = IPIV[IX-1];
02454                         if( IP != I )
02455                                 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02456                         IX += INCX;
02457                 }
02458         }
02459         else if( INCX < 0 )
02460         {
02461                 for( I=K2; I >= K1; I-- )
02462                 {
02463                         I_ = I - 1;
02464                         IP = IPIV[IX-1];
02465                         if( IP != I )
02466                                 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02467                         IX += INCX;
02468                 }
02469         }
02470 
02471         return;
02472 
02473         /* End of DLASWP
02474          * */
02475 #undef  A
02476 }
02477 
02478 /*DGETF2 lapack service routine */
02479 
02480 STATIC void DGETF2(int32 M, 
02481           int32 N, 
02482           double *A, 
02483           int32 LDA, 
02484           int32 IPIV[], 
02485           int32 *INFO)
02486 {
02487         int32 J, 
02488           JP, 
02489           J_, 
02490           limit;
02491 
02492         DEBUG_ENTRY( "DGETF2()" );
02493 
02494         /*  -- LAPACK routine (version 2.0) --
02495          * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
02496          * Courant Institute, Argonne National Lab, and Rice University
02497          * June 30, 1992
02498          *
02499          * .. Scalar Arguments .. */
02500         /* ..
02501          * .. Array Arguments .. */
02502         /* ..
02503          *
02504          *  Purpose
02505          *  =======
02506          *
02507          *  DGETF2 computes an LU factorization of a general m-by-n matrix A
02508          *  using partial pivoting with row interchanges.
02509          *
02510          *  The factorization has the form
02511          * A = P * L * U
02512          *  where P is a permutation matrix, L is lower triangular with unit
02513          *  diagonal elements (lower trapezoidal if m > n), and U is upper
02514          *  triangular (upper trapezoidal if m < n).
02515          *
02516          *  This is the right-looking Level 2 BLAS version of the algorithm.
02517          *
02518          *  Arguments
02519          *  =========
02520          *
02521          *  M       (input) INTEGER
02522          *  The number of rows of the matrix A.  M >= 0.
02523          *
02524          *  N       (input) INTEGER
02525          *  The number of columns of the matrix A.  N >= 0.
02526          *
02527          *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
02528          *  On entry, the m by n matrix to be factored.
02529          *  On exit, the factors L and U from the factorization
02530          *  A = P*L*U; the unit diagonal elements of L are not stored.
02531          *
02532          *  LDA     (input) INTEGER
02533          *  The leading dimension of the array A.  LDA >= MAX(1,M).
02534          *
02535          *  IPIV    (output) INTEGER array, dimension (MIN(M,N))
02536          *  The pivot indices; for 1 <= i <= MIN(M,N), row i of the
02537          *  matrix was interchanged with row IPIV(i).
02538          *
02539          *  INFO    (output) INTEGER
02540          *  = 0: successful exit
02541          *  < 0: if INFO = -k, the k-th argument had an illegal value
02542          *  > 0: if INFO = k, U(k,k) is exactly zero. The factorization
02543          *     has been completed, but the factor U is exactly
02544          *     singular, and division by zero will occur if it is used
02545          *     to solve a system of equations.
02546          *
02547          *  =====================================================================
02548          *
02549          * .. Parameters .. */
02550         /* ..
02551          * .. Local Scalars .. */
02552         /* ..
02553          * .. External Functions .. */
02554         /* ..
02555          * .. External Subroutines .. */
02556         /* ..
02557          * .. Intrinsic Functions .. */
02558         /* ..
02559          * .. Executable Statements ..
02560          *
02561          * Test the input parameters.
02562          * */
02563         *INFO = 0;
02564         if( M < 0 )
02565         {
02566                 *INFO = -1;
02567         }
02568         else if( N < 0 )
02569         {
02570                 *INFO = -2;
02571         }
02572         else if( LDA < MAX2(1,M) )
02573         {
02574                 *INFO = -4;
02575         }
02576         if( *INFO != 0 )
02577         {
02578                 XERBLA("DGETF2",-*INFO);
02579                 /* XERBLA does not return */
02580         }
02581 
02582         /* Quick return if possible
02583          * */
02584         if( M == 0 || N == 0 )
02585                 { 
02586                 return;}
02587 
02588         limit = MIN2(M,N);
02589         for( J=1; J <= limit; J++ )
02590         {
02591                 J_ = J - 1;
02592 
02593                 /*  Find pivot and test for singularity.
02594                  * */
02595                 JP = J - 1 + IDAMAX(M-J+1,&AA(J_,J_),1);
02596                 IPIV[J_] = JP;
02597                 if( AA(J_,JP-1) != ZERO )
02598                 {
02599                         /* Apply the interchange to columns 1:N.
02600                          * */
02601                         if( JP != J )
02602                                 DSWAP(N,&AA(0,J_),LDA,&AA(0,JP-1),LDA);
02603 
02604                         /* Compute elements J+1:M of J-th column.
02605                          * */
02606                         if( J < M )
02607                                 DSCAL(M-J,ONE/AA(J_,J_),&AA(J_,J_+1),1);
02608                 }
02609                 else if( *INFO == 0 )
02610                 {
02611                         *INFO = J;
02612                 }
02613 
02614                 if( J < MIN2(M,N) )
02615                 {
02616                         /* Update trailing submatrix.
02617                          * */
02618                         DGER(M-J,N-J,-ONE,&AA(J_,J_+1),1,&AA(J_+1,J_),LDA,&AA(J_+1,J_+1),
02619                           LDA);
02620                 }
02621         }
02622         return;
02623 
02624         /* End of DGETF2
02625          * */
02626 #undef  A
02627 }
02628 
02629 /*DGER service routine for matrix inversion */
02630 
02631 STATIC void DGER(int32 M, 
02632           int32 N, 
02633           double ALPHA, 
02634           double X[], 
02635           int32 INCX, 
02636           double Y[], 
02637           int32 INCY, 
02638           double *A, 
02639           int32 LDA)
02640 {
02641         int32 I, 
02642           INFO, 
02643           IX, 
02644           I_, 
02645           J, 
02646           JY, 
02647           J_, 
02648           KX;
02649         double TEMP;
02650 
02651         DEBUG_ENTRY( "DGER()" );
02652         /* .. Scalar Arguments .. */
02653         /* .. Array Arguments .. */
02654         /* ..
02655          *
02656          *  Purpose
02657          *  =======
02658          *
02659          *  DGER   performs the rank 1 operation
02660          *
02661          * A := alpha*x*y' + A,
02662          *
02663          *  where alpha is a scalar, x is an m element vector, y is an n element
02664          *  vector and A is an m by n matrix.
02665          *
02666          *  Parameters
02667          *  ==========
02668          *
02669          *  M      - INTEGER.
02670          * On entry, M specifies the number of rows of the matrix A.
02671          * M must be at least zero.
02672          * Unchanged on exit.
02673          *
02674          *  N      - INTEGER.
02675          * On entry, N specifies the number of columns of the matrix A.
02676          * N must be at least zero.
02677          * Unchanged on exit.
02678          *
02679          *  ALPHA  - DOUBLE PRECISION.
02680          * On entry, ALPHA specifies the scalar alpha.
02681          * Unchanged on exit.
02682          *
02683          *  X      - DOUBLE PRECISION array of dimension at least
02684          * ( 1 + ( m - 1 )*ABS( INCX ) ).
02685          * Before entry, the incremented array X must contain the m
02686          * element vector x.
02687          * Unchanged on exit.
02688          *
02689          *  INCX   - INTEGER.
02690          * On entry, INCX specifies the increment for the elements of
02691          * X. INCX must not be zero.
02692          * Unchanged on exit.
02693          *
02694          *  Y      - DOUBLE PRECISION array of dimension at least
02695          * ( 1 + ( n - 1 )*ABS( INCY ) ).
02696          * Before entry, the incremented array Y must contain the n
02697          * element vector y.
02698          * Unchanged on exit.
02699          *
02700          *  INCY   - INTEGER.
02701          * On entry, INCY specifies the increment for the elements of
02702          * Y. INCY must not be zero.
02703          * Unchanged on exit.
02704          *
02705          *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
02706          * Before entry, the leading m by n part of the array A must
02707          * contain the matrix of coefficients. On exit, A is
02708          * overwritten by the updated matrix.
02709          *
02710          *  LDA    - INTEGER.
02711          * On entry, LDA specifies the first dimension of A as declared
02712          * in the calling (sub) program. LDA must be at least
02713          * MAX( 1, m ).
02714          * Unchanged on exit.
02715          *
02716          *
02717          *  Level 2 Blas routine.
02718          *
02719          *  -- Written on 22-October-1986.
02720          * Jack Dongarra, Argonne National Lab.
02721          * Jeremy Du Croz, Nag Central Office.
02722          * Sven Hammarling, Nag Central Office.
02723          * Richard Hanson, Sandia National Labs.
02724          *
02725          *
02726          * .. Parameters .. */
02727         /* .. Local Scalars .. */
02728         /* .. External Subroutines .. */
02729         /* .. Intrinsic Functions .. */
02730         /* ..
02731          * .. Executable Statements ..
02732          *
02733          * Test the input parameters.
02734          * */
02735         INFO = 0;
02736         if( M < 0 )
02737         {
02738                 INFO = 1;
02739         }
02740         else if( N < 0 )
02741         {
02742                 INFO = 2;
02743         }
02744         else if( INCX == 0 )
02745         {
02746                 INFO = 5;
02747         }
02748         else if( INCY == 0 )
02749         {
02750                 INFO = 7;
02751         }
02752         else if( LDA < MAX2(1,M) )
02753         {
02754                 INFO = 9;
02755         }
02756         if( INFO != 0 )
02757         {
02758                 XERBLA("DGER  ",INFO);
02759                 /* XERBLA does not return */
02760         }
02761 
02762         /* Quick return if possible.
02763          * */
02764         if( ((M == 0) || (N == 0)) || (ALPHA == ZERO) )
02765                 { 
02766                 return;}
02767 
02768         /* Start the operations. In this version the elements of A are
02769          * accessed sequentially with one pass through A.
02770          * */
02771         if( INCY > 0 )
02772         {
02773                 JY = 1;
02774         }
02775         else
02776         {
02777                 JY = 1 - (N - 1)*INCY;
02778         }
02779         if( INCX == 1 )
02780         {
02781                 for( J=1; J <= N; J++ )
02782                 {
02783                         J_ = J - 1;
02784                         if( Y[JY-1] != ZERO )
02785                         {
02786                                 TEMP = ALPHA*Y[JY-1];
02787                                 for( I=1; I <= M; I++ )
02788                                 {
02789                                         I_ = I - 1;
02790                                         AA(J_,I_) += X[I_]*TEMP;
02791                                 }
02792                         }
02793                         JY += INCY;
02794                 }
02795         }
02796         else
02797         {
02798                 if( INCX > 0 )
02799                 {
02800                         KX = 1;
02801                 }
02802                 else
02803                 {
02804                         KX = 1 - (M - 1)*INCX;
02805                 }
02806                 for( J=1; J <= N; J++ )
02807                 {
02808                         J_ = J - 1;
02809                         if( Y[JY-1] != ZERO )
02810                         {
02811                                 TEMP = ALPHA*Y[JY-1];
02812                                 IX = KX;
02813                                 for( I=1; I <= M; I++ )
02814                                 {
02815                                         I_ = I - 1;
02816                                         AA(J_,I_) += X[IX-1]*TEMP;
02817                                         IX += INCX;
02818                                 }
02819                         }
02820                         JY += INCY;
02821                 }
02822         }
02823 
02824         return;
02825 
02826         /* End of DGER  .
02827          * */
02828 #undef  A
02829 }
02830 
02831 /* DGEMM matrix inversion helper routine*/
02832 
02833 STATIC void DGEMM(int32 TRANSA, 
02834           int32 TRANSB, 
02835           int32 M, 
02836           int32 N, 
02837           int32 K, 
02838           double ALPHA, 
02839           double *A, 
02840           int32 LDA, 
02841           double *B, 
02842           int32 LDB, 
02843           double BETA, 
02844           double *C, 
02845           int32 LDC)
02846 {
02847         int32 NOTA, 
02848           NOTB;
02849         int32 I, 
02850           INFO, 
02851           I_, 
02852           J, 
02853           J_, 
02854           L, 
02855           L_, 
02856           NROWA, 
02857           NROWB;
02858         double TEMP;
02859 
02860         DEBUG_ENTRY( "DGEMM()" );
02861         /* .. Scalar Arguments .. */
02862         /* .. Array Arguments .. */
02863         /* ..
02864          *
02865          *  Purpose
02866          *  =======
02867          *
02868          *  DGEMM  performs one of the matrix-matrix operations
02869          *
02870          * C := alpha*op( A )*op( B ) + beta*C,
02871          *
02872          *  where  op( X ) is one of
02873          *
02874          * op( X ) = X   or   op( X ) = X',
02875          *
02876          *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
02877          *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
02878          *
02879          *  Parameters
02880          *  ==========
02881          *
02882          *  TRANSA - CHARACTER*1.
02883          * On entry, TRANSA specifies the form of op( A ) to be used in
02884          * the matrix multiplication as follows:
02885          *
02886          *    TRANSA = 'N' or 'n',  op( A ) = A.
02887          *
02888          *    TRANSA = 'T' or 't',  op( A ) = A'.
02889          *
02890          *    TRANSA = 'C' or 'c',  op( A ) = A'.
02891          *
02892          * Unchanged on exit.
02893          *
02894          *  TRANSB - CHARACTER*1.
02895          * On entry, TRANSB specifies the form of op( B ) to be used in
02896          * the matrix multiplication as follows:
02897          *
02898          *    TRANSB = 'N' or 'n',  op( B ) = B.
02899          *
02900          *    TRANSB = 'T' or 't',  op( B ) = B'.
02901          *
02902          *    TRANSB = 'C' or 'c',  op( B ) = B'.
02903          *
02904          * Unchanged on exit.
02905          *
02906          *  M      - INTEGER.
02907          * On entry,  M  specifies  the number  of rows  of the  matrix
02908          * op( A )  and of the  matrix  C.  M  must  be at least  zero.
02909          * Unchanged on exit.
02910          *
02911          *  N      - INTEGER.
02912          * On entry,  N  specifies the number  of columns of the matrix
02913          * op( B ) and the number of columns of the matrix C. N must be
02914          * at least zero.
02915          * Unchanged on exit.
02916          *
02917          *  K      - INTEGER.
02918          * On entry,  K  specifies  the number of columns of the matrix
02919          * op( A ) and the number of rows of the matrix op( B ). K must
02920          * be at least  zero.
02921          * Unchanged on exit.
02922          *
02923          *  ALPHA  - DOUBLE PRECISION.
02924          * On entry, ALPHA specifies the scalar alpha.
02925          * Unchanged on exit.
02926          *
02927          *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
02928          * k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
02929          * Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
02930          * part of the array  A  must contain the matrix  A,  otherwise
02931          * the leading  k by m  part of the array  A  must contain  the
02932          * matrix A.
02933          * Unchanged on exit.
02934          *
02935          *  LDA    - INTEGER.
02936          * On entry, LDA specifies the first dimension of A as declared
02937          * in the calling (sub) program. When  TRANSA = 'N' or 'n' then
02938          * LDA must be at least  MAX( 1, m ), otherwise  LDA must be at
02939          * least  MAX( 1, k ).
02940          * Unchanged on exit.
02941          *
02942          *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
02943          * n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
02944          * Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
02945          * part of the array  B  must contain the matrix  B,  otherwise
02946          * the leading  n by k  part of the array  B  must contain  the
02947          * matrix B.
02948          * Unchanged on exit.
02949          *
02950          *  LDB    - INTEGER.
02951          * On entry, LDB specifies the first dimension of B as declared
02952          * in the calling (sub) program. When  TRANSB = 'N' or 'n' then
02953          * LDB must be at least  MAX( 1, k ), otherwise  LDB must be at
02954          * least  MAX( 1, n ).
02955          * Unchanged on exit.
02956          *
02957          *  BETA   - DOUBLE PRECISION.
02958          * On entry,  BETA  specifies the scalar  beta.  When  BETA  is
02959          * supplied as zero then C need not be set on input.
02960          * Unchanged on exit.
02961          *
02962          *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
02963          * Before entry, the leading  m by n  part of the array  C must
02964          * contain the matrix  C,  except when  beta  is zero, in which
02965          * case C need not be set on entry.
02966          * On exit, the array  C  is overwritten by the  m by n  matrix
02967          * ( alpha*op( A )*op( B ) + beta*C ).
02968          *
02969          *  LDC    - INTEGER.
02970          * On entry, LDC specifies the first dimension of C as declared
02971          * in  the  calling  (sub)  program.   LDC  must  be  at  least
02972          * MAX( 1, m ).
02973          * Unchanged on exit.
02974          *
02975          *
02976          *  Level 3 Blas routine.
02977          *
02978          *  -- Written on 8-February-1989.
02979          * Jack Dongarra, Argonne National Laboratory.
02980          * Iain Duff, AERE Harwell.
02981          * Jeremy Du Croz, Numerical Algorithms Group Ltd.
02982          * Sven Hammarling, Numerical Algorithms Group Ltd.
02983          *
02984          *
02985          * .. External Functions .. */
02986         /* .. External Subroutines .. */
02987         /* .. Intrinsic Functions .. */
02988         /* .. Local Scalars .. */
02989         /* .. Parameters .. */
02990         /* ..
02991          * .. Executable Statements ..
02992          *
02993          * Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
02994          * transposed and set  NROWA, NROWB  as the number of rows
02995          * and  columns of  A  and the  number of  rows  of  B  respectively.
02996          * */
02997         NOTA = LSAME(TRANSA,'N');
02998         NOTB = LSAME(TRANSB,'N');
02999 
03000         if( NOTA )
03001         {
03002                 NROWA = M;
03003         }
03004         else
03005         {
03006                 NROWA = K;
03007         }
03008 
03009         if( NOTB )
03010         {
03011                 NROWB = K;
03012         }
03013         else
03014         {
03015                 NROWB = N;
03016         }
03017 
03018         /* Test the input parameters.
03019          * */
03020         INFO = 0;
03021         if( ((!NOTA) && (!LSAME(TRANSA,'C'))) && (!LSAME(TRANSA,'T')) )
03022         {
03023                 INFO = 1;
03024         }
03025         else if( 
03026                 ((!NOTB) && (!LSAME(TRANSB,'C'))) && (!LSAME(TRANSB,'T')) )
03027         {
03028                 INFO = 2;
03029         }
03030 
03031         else if( M < 0 )
03032         {
03033                 INFO = 3;
03034         }
03035 
03036         else if( N < 0 )
03037         {
03038                 INFO = 4;
03039         }
03040 
03041         else if( K < 0 )
03042         {
03043                 INFO = 5;
03044         }
03045 
03046         else if( LDA < MAX2(1,NROWA) )
03047         {
03048                 INFO = 8;
03049         }
03050 
03051         else if( LDB < MAX2(1,NROWB) )
03052         {
03053                 INFO = 10;
03054         }
03055 
03056         else if( LDC < MAX2(1,M) )
03057         {
03058                 INFO = 13;
03059         }
03060 
03061         if( INFO != 0 )
03062         {
03063                 XERBLA("DGEMM ",INFO);
03064                 /* XERBLA does not return */
03065         }
03066 
03067         /* Quick return if possible.
03068          * */
03069         if( ((M == 0) || (N == 0)) || (((ALPHA == ZERO) || (K == 0)) && 
03070           (BETA == ONE)) )
03071         { 
03072                 return;
03073         }
03074 
03075         /* And if  alpha.eq.zero. */
03076         if( ALPHA == ZERO )
03077         {
03078                 if( BETA == ZERO )
03079                 {
03080                         for( J=1; J <= N; J++ )
03081                         {
03082                                 J_ = J - 1;
03083                                 for( I=1; I <= M; I++ )
03084                                 {
03085                                         I_ = I - 1;
03086                                         CC(J_,I_) = ZERO;
03087                                 }
03088                         }
03089                 }
03090 
03091                 else
03092                 {
03093                         for( J=1; J <= N; J++ )
03094                         {
03095                                 J_ = J - 1;
03096                                 for( I=1; I <= M; I++ )
03097                                 {
03098                                         I_ = I - 1;
03099                                         CC(J_,I_) *= BETA;
03100                                 }
03101                         }
03102                 }
03103                 return;
03104         }
03105 
03106         /* Start the operations.*/
03107         if( NOTB )
03108         {
03109 
03110                 if( NOTA )
03111                 {
03112 
03113                         /* Form  C := alpha*A*B + beta*C.
03114                          * */
03115                         for( J=1; J <= N; J++ )
03116                         {
03117                                 J_ = J - 1;
03118                                 if( BETA == ZERO )
03119                                 {
03120                                         for( I=1; I <= M; I++ )
03121                                         {
03122                                                 I_ = I - 1;
03123                                                 CC(J_,I_) = ZERO;
03124                                         }
03125                                 }
03126 
03127                                 else if( BETA != ONE )
03128                                 {
03129                                         for( I=1; I <= M; I++ )
03130                                         {
03131                                                 I_ = I - 1;
03132                                                 CC(J_,I_) *= BETA;
03133                                         }
03134                                 }
03135 
03136                                 for( L=1; L <= K; L++ )
03137                                 {
03138                                         L_ = L - 1;
03139                                         if( BB(J_,L_) != ZERO )
03140                                         {
03141                                                 TEMP = ALPHA*BB(J_,L_);
03142                                                 for( I=1; I <= M; I++ )
03143                                                 {
03144                                                         I_ = I - 1;
03145                                                         CC(J_,I_) += TEMP*AA(L_,I_);
03146                                                 }
03147                                         }
03148                                 }
03149                         }
03150                 }
03151                 else
03152                 {
03153 
03154                         /* Form  C := alpha*A'*B + beta*C */
03155                         for( J=1; J <= N; J++ )
03156                         {
03157                                 J_ = J - 1;
03158                                 for( I=1; I <= M; I++ )
03159                                 {
03160                                         I_ = I - 1;
03161                                         TEMP = ZERO;
03162                                         for( L=1; L <= K; L++ )
03163                                         {
03164                                                 L_ = L - 1;
03165                                                 TEMP += AA(I_,L_)*BB(J_,L_);
03166                                         }
03167 
03168                                         if( BETA == ZERO )
03169                                         {
03170                                                 CC(J_,I_) = ALPHA*TEMP;
03171                                         }
03172                                         else
03173                                         {
03174                                                 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
03175                                         }
03176                                 }
03177                         }
03178                 }
03179         }
03180         else
03181         {
03182                 if( NOTA )
03183                 {
03184 
03185                         /* Form  C := alpha*A*B' + beta*C
03186                          * */
03187                         for( J=1; J <= N; J++ )
03188                         {
03189                                 J_ = J - 1;
03190                                 if( BETA == ZERO )
03191                                 {
03192                                         for( I=1; I <= M; I++ )
03193                                         {
03194                                                 I_ = I - 1;
03195                                                 CC(J_,I_) = ZERO;
03196                                         }
03197                                 }
03198 
03199                                 else if( BETA != ONE )
03200                                 {
03201                                         for( I=1; I <= M; I++ )
03202                                         {
03203                                                 I_ = I - 1;
03204                                                 CC(J_,I_) *= BETA;
03205                                         }
03206                                 }
03207 
03208                                 for( L=1; L <= K; L++ )
03209                                 {
03210                                         L_ = L - 1;
03211                                         if( BB(L_,J_) != ZERO )
03212                                         {
03213                                                 TEMP = ALPHA*BB(L_,J_);
03214                                                 for( I=1; I <= M; I++ )
03215                                                 {
03216                                                         I_ = I - 1;
03217                                                         CC(J_,I_) += TEMP*AA(L_,I_);
03218                                                 }
03219                                         }
03220                                 }
03221                         }
03222                 }
03223 
03224                 else
03225                 {
03226 
03227                         /* Form  C := alpha*A'*B' + beta*C */
03228                         for( J=1; J <= N; J++ )
03229                         {
03230                                 J_ = J - 1;
03231 
03232                                 for( I=1; I <= M; I++ )
03233                                 {
03234                                         I_ = I - 1;
03235                                         TEMP = ZERO;
03236 
03237                                         for( L=1; L <= K; L++ )
03238                                         {
03239                                                 L_ = L - 1;
03240                                                 TEMP += AA(I_,L_)*BB(L_,J_);
03241                                         }
03242 
03243                                         if( BETA == ZERO )
03244                                         {
03245                                                 CC(J_,I_) = ALPHA*TEMP;
03246                                         }
03247                                         else
03248                                         {
03249                                                 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
03250                                         }
03251 
03252                                 }
03253                         }
03254                 }
03255         }
03256 
03257         return;
03258 
03259         /* End of DGEMM .*/
03260 #undef  C
03261 #undef  B
03262 #undef  A
03263 }
03264 #undef AA
03265 #undef BB
03266 #undef CC
03267 
03268 /* Subroutine */ 
03269 #if 0
03270 STATIC int32 DGTSV(int32 *n, int32 *nrhs, double *dl, 
03271         double *d__, double *du, double *b, int32 *ldb, int32 
03272         *info)
03273 {
03274     /* System generated locals */
03275     int32 b_dim1, b_offset, i__1, i__2;
03276 
03277     /* Local variables */
03278     double fact, temp;
03279     int32 i__, j;
03280 
03281 
03282 #define b_ref(a_1,a_2) b[(a_2)*(b_dim1) + (a_1)]
03283 
03284 
03285 /*  -- LAPACK routine (version 3.0) --   
03286        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
03287        Courant Institute, Argonne National Lab, and Rice University   
03288        October 31, 1999   
03289 
03290 
03291     Purpose   
03292     =======   
03293 
03294     DGTSV  solves the equation   
03295 
03296        A*X = B,   
03297 
03298     where A is an n by n tridiagonal matrix, by Gaussian elimination with   
03299     partial pivoting.   
03300 
03301     Note that the equation  A'*X = B  may be solved by interchanging the   
03302     order of the arguments DU and DL.   
03303 
03304     Arguments   
03305     =========   
03306 
03307     N       (input) INTEGER   
03308             The order of the matrix A.  N >= 0.   
03309 
03310     NRHS    (input) INTEGER   
03311             The number of right hand sides, i.e., the number of columns   
03312             of the matrix B.  NRHS >= 0.   
03313 
03314     DL      (input/output) DOUBLE PRECISION array, dimension (N-1)   
03315             On entry, DL must contain the (n-1) sub-diagonal elements of   
03316             A.   
03317 
03318             On exit, DL is overwritten by the (n-2) elements of the   
03319             second super-diagonal of the upper triangular matrix U from   
03320             the LU factorization of A, in DL(1), ..., DL(n-2).   
03321 
03322     D       (input/output) DOUBLE PRECISION array, dimension (N)   
03323             On entry, D must contain the diagonal elements of A.   
03324 
03325             On exit, D is overwritten by the n diagonal elements of U.   
03326 
03327     DU      (input/output) DOUBLE PRECISION array, dimension (N-1)   
03328             On entry, DU must contain the (n-1) super-diagonal elements   
03329             of A.   
03330 
03331             On exit, DU is overwritten by the (n-1) elements of the first   
03332             super-diagonal of U.   
03333 
03334     B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)   
03335             On entry, the N by NRHS matrix of right hand side matrix B.   
03336             On exit, if INFO = 0, the N by NRHS solution matrix X.   
03337 
03338     LDB     (input) INTEGER   
03339             The leading dimension of the array B.  LDB >= max(1,N).   
03340 
03341     INFO    (output) INTEGER   
03342             = 0: successful exit   
03343             < 0: if INFO = -i, the i-th argument had an illegal value   
03344             > 0: if INFO = i, U(i,i) is exactly zero, and the solution   
03345                  has not been computed.  The factorization has not been   
03346                  completed unless i = N.   
03347 
03348     =====================================================================   
03349 
03350 
03351        Parameter adjustments */
03352     --dl;
03353     --d__;
03354     --du;
03355     b_dim1 = *ldb;
03356     b_offset = 1 + b_dim1 * 1;
03357     b -= b_offset;
03358 
03359     /* Function Body */
03360     *info = 0;
03361     if(*n < 0) {
03362         *info = -1;
03363     } else if(*nrhs < 0) {
03364         *info = -2;
03365     } else if(*ldb < *n && *ldb < 1) {
03366         *info = -7;
03367     }
03368     if(*info != 0) {
03369         i__1 = -(*info);
03370         XERBLA("DGTSV ", i__1);
03371         /* XERBLA does not return */
03372     }
03373 
03374     if(*n == 0) {
03375         return 0;
03376     }
03377 
03378     if(*nrhs == 1) {
03379         i__1 = *n - 2;
03380         for(i__ = 1; i__ <= i__1; ++i__) {
03381             if(fabs(d__[i__]) >= fabs(dl[i__])) {
03382 
03383 /*              No row interchange required */
03384 
03385                 if(d__[i__] != 0.) {
03386                     fact = dl[i__] / d__[i__];
03387                     d__[i__ + 1] -= fact * du[i__];
03388                     b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__, 
03389                             1);
03390                 } else {
03391                     *info = i__;
03392                     return 0;
03393                 }
03394                 dl[i__] = 0.;
03395             } else {
03396 
03397 /*              Interchange rows I and I+1 */
03398 
03399                 fact = d__[i__] / dl[i__];
03400                 d__[i__] = dl[i__];
03401                 temp = d__[i__ + 1];
03402                 d__[i__ + 1] = du[i__] - fact * temp;
03403                 dl[i__] = du[i__ + 1];
03404                 du[i__ + 1] = -fact * dl[i__];
03405                 du[i__] = temp;
03406                 temp = b_ref(i__, 1);
03407                 b_ref(i__, 1) = b_ref(i__ + 1, 1);
03408                 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
03409             }
03410 /* L10: */
03411         }
03412         if(*n > 1) {
03413             i__ = *n - 1;
03414             if(fabs(d__[i__]) >= fabs(dl[i__])) {
03415                 if(d__[i__] != 0.) {
03416                     fact = dl[i__] / d__[i__];
03417                     d__[i__ + 1] -= fact * du[i__];
03418                     b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__, 
03419                             1);
03420                 } else {
03421                     *info = i__;
03422                     return 0;
03423                 }
03424             } else {
03425                 fact = d__[i__] / dl[i__];
03426                 d__[i__] = dl[i__];
03427                 temp = d__[i__ + 1];
03428                 d__[i__ + 1] = du[i__] - fact * temp;
03429                 du[i__] = temp;
03430                 temp = b_ref(i__, 1);
03431                 b_ref(i__, 1) = b_ref(i__ + 1, 1);
03432                 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
03433             }
03434         }
03435         if(d__[*n] == 0.) {
03436             *info = *n;
03437             return 0;
03438         }
03439     } else {
03440         i__1 = *n - 2;
03441         for(i__ = 1; i__ <= i__1; ++i__) {
03442             if(fabs(d__[i__]) >= fabs(dl[i__])) {
03443 
03444 /*              No row interchange required */
03445 
03446                 if(d__[i__] != 0.) {
03447                     fact = dl[i__] / d__[i__];
03448                     d__[i__ + 1] -= fact * du[i__];
03449                     i__2 = *nrhs;
03450                     for(j = 1; j <= i__2; ++j) {
03451                         b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
03452                                 i__, j);
03453 /* L20: */
03454                     }
03455                 } else {
03456                     *info = i__;
03457                     return 0;
03458                 }
03459                 dl[i__] = 0.;
03460             } else {
03461 
03462 /*              Interchange rows I and I+1 */
03463 
03464                 fact = d__[i__] / dl[i__];
03465                 d__[i__] = dl[i__];
03466                 temp = d__[i__ + 1];
03467                 d__[i__ + 1] = du[i__] - fact * temp;
03468                 dl[i__] = du[i__ + 1];
03469                 du[i__ + 1] = -fact * dl[i__];
03470                 du[i__] = temp;
03471                 i__2 = *nrhs;
03472                 for(j = 1; j <= i__2; ++j) {
03473                     temp = b_ref(i__, j);
03474                     b_ref(i__, j) = b_ref(i__ + 1, j);
03475                     b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
03476 /* L30: */
03477                 }
03478             }
03479 /* L40: */
03480         }
03481         if(*n > 1) {
03482             i__ = *n - 1;
03483             if( fabs(d__[i__]) >= fabs(dl[i__])) 
03484                 {
03485                         if(d__[i__] != 0.) 
03486                         {
03487                                 fact = dl[i__] / d__[i__];
03488                                 d__[i__ + 1] -= fact * du[i__];
03489                                 i__1 = *nrhs;
03490                                 for(j = 1; j <= i__1; ++j) {
03491                                 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
03492                                         i__, j);
03493         /* L50: */
03494                                 }
03495                         } 
03496                         else 
03497                         {
03498                                 *info = i__;
03499                                 return 0;
03500                         }
03501             } else {
03502                 fact = d__[i__] / dl[i__];
03503                 d__[i__] = dl[i__];
03504                 temp = d__[i__ + 1];
03505                 d__[i__ + 1] = du[i__] - fact * temp;
03506                 du[i__] = temp;
03507                 i__1 = *nrhs;
03508                 for(j = 1; j <= i__1; ++j) {
03509                     temp = b_ref(i__, j);
03510                     b_ref(i__, j) = b_ref(i__ + 1, j);
03511                     b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
03512 /* L60: */
03513                 }
03514             }
03515         }
03516         if(d__[*n] == 0.) {
03517             *info = *n;
03518             return 0;
03519         }
03520     }
03521 
03522 /*     Back solve with the matrix U from the factorization. */
03523 
03524     if(*nrhs <= 2) {
03525         j = 1;
03526 L70:
03527         b_ref(*n, j) = b_ref(*n, j) / d__[*n];
03528         if(*n > 1) {
03529             b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, j)) 
03530                     / d__[*n - 1];
03531         }
03532         for(i__ = *n - 2; i__ >= 1; --i__) {
03533             b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) - dl[
03534                     i__] * b_ref(i__ + 2, j)) / d__[i__];
03535 /* L80: */
03536         }
03537         if(j < *nrhs) {
03538             ++j;
03539             goto L70;
03540         }
03541     } else {
03542         i__1 = *nrhs;
03543         for(j = 1; j <= i__1; ++j) {
03544             b_ref(*n, j) = b_ref(*n, j) / d__[*n];
03545             if(*n > 1) {
03546                 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, 
03547                         j)) / d__[*n - 1];
03548             }
03549             for(i__ = *n - 2; i__ >= 1; --i__) {
03550                 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) 
03551                         - dl[i__] * b_ref(i__ + 2, j)) / d__[i__];
03552 /* L90: */
03553             }
03554 /* L100: */
03555         }
03556     }
03557 
03558     return 0;
03559 
03560 /*     End of DGTSV */
03561 
03562 } /* dgtsv_ */
03563 #endif
03564 #undef b_ref
03565 #endif
03566 /*lint +e725 expected pos indentation */
03567 /*lint +e801 goto is deprecated */
/home66/gary/public_html/cloudy/c08_branch/source/thirdparty_lapack.cpp
