cloudy trunk
|
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 */