|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) 00002 * .. Scalar Arguments .. 00003 DOUBLE PRECISION ALPHA 00004 INTEGER LDA,LDB,M,N 00005 CHARACTER DIAG,SIDE,TRANSA,UPLO 00006 * .. 00007 * .. Array Arguments .. 00008 DOUBLE PRECISION A(LDA,*),B(LDB,*) 00009 * .. 00010 * 00011 * Purpose 00012 * ======= 00013 * 00014 * DTRSM solves one of the matrix equations 00015 * 00016 * op( A )*X = alpha*B, or X*op( A ) = alpha*B, 00017 * 00018 * where alpha is a scalar, X and B are m by n matrices, A is a unit, or 00019 * non-unit, upper or lower triangular matrix and op( A ) is one of 00020 * 00021 * op( A ) = A or op( A ) = A'. 00022 * 00023 * The matrix X is overwritten on B. 00024 * 00025 * Arguments 00026 * ========== 00027 * 00028 * SIDE - CHARACTER*1. 00029 * On entry, SIDE specifies whether op( A ) appears on the left 00030 * or right of X as follows: 00031 * 00032 * SIDE = 'L' or 'l' op( A )*X = alpha*B. 00033 * 00034 * SIDE = 'R' or 'r' X*op( A ) = alpha*B. 00035 * 00036 * Unchanged on exit. 00037 * 00038 * UPLO - CHARACTER*1. 00039 * On entry, UPLO specifies whether the matrix A is an upper or 00040 * lower triangular matrix as follows: 00041 * 00042 * UPLO = 'U' or 'u' A is an upper triangular matrix. 00043 * 00044 * UPLO = 'L' or 'l' A is a lower triangular matrix. 00045 * 00046 * Unchanged on exit. 00047 * 00048 * TRANSA - CHARACTER*1. 00049 * On entry, TRANSA specifies the form of op( A ) to be used in 00050 * the matrix multiplication as follows: 00051 * 00052 * TRANSA = 'N' or 'n' op( A ) = A. 00053 * 00054 * TRANSA = 'T' or 't' op( A ) = A'. 00055 * 00056 * TRANSA = 'C' or 'c' op( A ) = A'. 00057 * 00058 * Unchanged on exit. 00059 * 00060 * DIAG - CHARACTER*1. 00061 * On entry, DIAG specifies whether or not A is unit triangular 00062 * as follows: 00063 * 00064 * DIAG = 'U' or 'u' A is assumed to be unit triangular. 00065 * 00066 * DIAG = 'N' or 'n' A is not assumed to be unit 00067 * triangular. 00068 * 00069 * Unchanged on exit. 00070 * 00071 * M - INTEGER. 00072 * On entry, M specifies the number of rows of B. M must be at 00073 * least zero. 00074 * Unchanged on exit. 00075 * 00076 * N - INTEGER. 00077 * On entry, N specifies the number of columns of B. N must be 00078 * at least zero. 00079 * Unchanged on exit. 00080 * 00081 * ALPHA - DOUBLE PRECISION. 00082 * On entry, ALPHA specifies the scalar alpha. When alpha is 00083 * zero then A is not referenced and B need not be set before 00084 * entry. 00085 * Unchanged on exit. 00086 * 00087 * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m 00088 * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. 00089 * Before entry with UPLO = 'U' or 'u', the leading k by k 00090 * upper triangular part of the array A must contain the upper 00091 * triangular matrix and the strictly lower triangular part of 00092 * A is not referenced. 00093 * Before entry with UPLO = 'L' or 'l', the leading k by k 00094 * lower triangular part of the array A must contain the lower 00095 * triangular matrix and the strictly upper triangular part of 00096 * A is not referenced. 00097 * Note that when DIAG = 'U' or 'u', the diagonal elements of 00098 * A are not referenced either, but are assumed to be unity. 00099 * Unchanged on exit. 00100 * 00101 * LDA - INTEGER. 00102 * On entry, LDA specifies the first dimension of A as declared 00103 * in the calling (sub) program. When SIDE = 'L' or 'l' then 00104 * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' 00105 * then LDA must be at least max( 1, n ). 00106 * Unchanged on exit. 00107 * 00108 * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). 00109 * Before entry, the leading m by n part of the array B must 00110 * contain the right-hand side matrix B, and on exit is 00111 * overwritten by the solution matrix X. 00112 * 00113 * LDB - INTEGER. 00114 * On entry, LDB specifies the first dimension of B as declared 00115 * in the calling (sub) program. LDB must be at least 00116 * max( 1, m ). 00117 * Unchanged on exit. 00118 * 00119 * Further Details 00120 * =============== 00121 * 00122 * Level 3 Blas routine. 00123 * 00124 * 00125 * -- Written on 8-February-1989. 00126 * Jack Dongarra, Argonne National Laboratory. 00127 * Iain Duff, AERE Harwell. 00128 * Jeremy Du Croz, Numerical Algorithms Group Ltd. 00129 * Sven Hammarling, Numerical Algorithms Group Ltd. 00130 * 00131 * ===================================================================== 00132 * 00133 * .. External Functions .. 00134 LOGICAL LSAME 00135 EXTERNAL LSAME 00136 * .. 00137 * .. External Subroutines .. 00138 EXTERNAL XERBLA 00139 * .. 00140 * .. Intrinsic Functions .. 00141 INTRINSIC MAX 00142 * .. 00143 * .. Local Scalars .. 00144 DOUBLE PRECISION TEMP 00145 INTEGER I,INFO,J,K,NROWA 00146 LOGICAL LSIDE,NOUNIT,UPPER 00147 * .. 00148 * .. Parameters .. 00149 DOUBLE PRECISION ONE,ZERO 00150 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) 00151 * .. 00152 * 00153 * Test the input parameters. 00154 * 00155 LSIDE = LSAME(SIDE,'L') 00156 IF (LSIDE) THEN 00157 NROWA = M 00158 ELSE 00159 NROWA = N 00160 END IF 00161 NOUNIT = LSAME(DIAG,'N') 00162 UPPER = LSAME(UPLO,'U') 00163 * 00164 INFO = 0 00165 IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN 00166 INFO = 1 00167 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN 00168 INFO = 2 00169 ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND. 00170 + (.NOT.LSAME(TRANSA,'T')) .AND. 00171 + (.NOT.LSAME(TRANSA,'C'))) THEN 00172 INFO = 3 00173 ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN 00174 INFO = 4 00175 ELSE IF (M.LT.0) THEN 00176 INFO = 5 00177 ELSE IF (N.LT.0) THEN 00178 INFO = 6 00179 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN 00180 INFO = 9 00181 ELSE IF (LDB.LT.MAX(1,M)) THEN 00182 INFO = 11 00183 END IF 00184 IF (INFO.NE.0) THEN 00185 CALL XERBLA('DTRSM ',INFO) 00186 RETURN 00187 END IF 00188 * 00189 * Quick return if possible. 00190 * 00191 IF (M.EQ.0 .OR. N.EQ.0) RETURN 00192 * 00193 * And when alpha.eq.zero. 00194 * 00195 IF (ALPHA.EQ.ZERO) THEN 00196 DO 20 J = 1,N 00197 DO 10 I = 1,M 00198 B(I,J) = ZERO 00199 10 CONTINUE 00200 20 CONTINUE 00201 RETURN 00202 END IF 00203 * 00204 * Start the operations. 00205 * 00206 IF (LSIDE) THEN 00207 IF (LSAME(TRANSA,'N')) THEN 00208 * 00209 * Form B := alpha*inv( A )*B. 00210 * 00211 IF (UPPER) THEN 00212 DO 60 J = 1,N 00213 IF (ALPHA.NE.ONE) THEN 00214 DO 30 I = 1,M 00215 B(I,J) = ALPHA*B(I,J) 00216 30 CONTINUE 00217 END IF 00218 DO 50 K = M,1,-1 00219 IF (B(K,J).NE.ZERO) THEN 00220 IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) 00221 DO 40 I = 1,K - 1 00222 B(I,J) = B(I,J) - B(K,J)*A(I,K) 00223 40 CONTINUE 00224 END IF 00225 50 CONTINUE 00226 60 CONTINUE 00227 ELSE 00228 DO 100 J = 1,N 00229 IF (ALPHA.NE.ONE) THEN 00230 DO 70 I = 1,M 00231 B(I,J) = ALPHA*B(I,J) 00232 70 CONTINUE 00233 END IF 00234 DO 90 K = 1,M 00235 IF (B(K,J).NE.ZERO) THEN 00236 IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) 00237 DO 80 I = K + 1,M 00238 B(I,J) = B(I,J) - B(K,J)*A(I,K) 00239 80 CONTINUE 00240 END IF 00241 90 CONTINUE 00242 100 CONTINUE 00243 END IF 00244 ELSE 00245 * 00246 * Form B := alpha*inv( A' )*B. 00247 * 00248 IF (UPPER) THEN 00249 DO 130 J = 1,N 00250 DO 120 I = 1,M 00251 TEMP = ALPHA*B(I,J) 00252 DO 110 K = 1,I - 1 00253 TEMP = TEMP - A(K,I)*B(K,J) 00254 110 CONTINUE 00255 IF (NOUNIT) TEMP = TEMP/A(I,I) 00256 B(I,J) = TEMP 00257 120 CONTINUE 00258 130 CONTINUE 00259 ELSE 00260 DO 160 J = 1,N 00261 DO 150 I = M,1,-1 00262 TEMP = ALPHA*B(I,J) 00263 DO 140 K = I + 1,M 00264 TEMP = TEMP - A(K,I)*B(K,J) 00265 140 CONTINUE 00266 IF (NOUNIT) TEMP = TEMP/A(I,I) 00267 B(I,J) = TEMP 00268 150 CONTINUE 00269 160 CONTINUE 00270 END IF 00271 END IF 00272 ELSE 00273 IF (LSAME(TRANSA,'N')) THEN 00274 * 00275 * Form B := alpha*B*inv( A ). 00276 * 00277 IF (UPPER) THEN 00278 DO 210 J = 1,N 00279 IF (ALPHA.NE.ONE) THEN 00280 DO 170 I = 1,M 00281 B(I,J) = ALPHA*B(I,J) 00282 170 CONTINUE 00283 END IF 00284 DO 190 K = 1,J - 1 00285 IF (A(K,J).NE.ZERO) THEN 00286 DO 180 I = 1,M 00287 B(I,J) = B(I,J) - A(K,J)*B(I,K) 00288 180 CONTINUE 00289 END IF 00290 190 CONTINUE 00291 IF (NOUNIT) THEN 00292 TEMP = ONE/A(J,J) 00293 DO 200 I = 1,M 00294 B(I,J) = TEMP*B(I,J) 00295 200 CONTINUE 00296 END IF 00297 210 CONTINUE 00298 ELSE 00299 DO 260 J = N,1,-1 00300 IF (ALPHA.NE.ONE) THEN 00301 DO 220 I = 1,M 00302 B(I,J) = ALPHA*B(I,J) 00303 220 CONTINUE 00304 END IF 00305 DO 240 K = J + 1,N 00306 IF (A(K,J).NE.ZERO) THEN 00307 DO 230 I = 1,M 00308 B(I,J) = B(I,J) - A(K,J)*B(I,K) 00309 230 CONTINUE 00310 END IF 00311 240 CONTINUE 00312 IF (NOUNIT) THEN 00313 TEMP = ONE/A(J,J) 00314 DO 250 I = 1,M 00315 B(I,J) = TEMP*B(I,J) 00316 250 CONTINUE 00317 END IF 00318 260 CONTINUE 00319 END IF 00320 ELSE 00321 * 00322 * Form B := alpha*B*inv( A' ). 00323 * 00324 IF (UPPER) THEN 00325 DO 310 K = N,1,-1 00326 IF (NOUNIT) THEN 00327 TEMP = ONE/A(K,K) 00328 DO 270 I = 1,M 00329 B(I,K) = TEMP*B(I,K) 00330 270 CONTINUE 00331 END IF 00332 DO 290 J = 1,K - 1 00333 IF (A(J,K).NE.ZERO) THEN 00334 TEMP = A(J,K) 00335 DO 280 I = 1,M 00336 B(I,J) = B(I,J) - TEMP*B(I,K) 00337 280 CONTINUE 00338 END IF 00339 290 CONTINUE 00340 IF (ALPHA.NE.ONE) THEN 00341 DO 300 I = 1,M 00342 B(I,K) = ALPHA*B(I,K) 00343 300 CONTINUE 00344 END IF 00345 310 CONTINUE 00346 ELSE 00347 DO 360 K = 1,N 00348 IF (NOUNIT) THEN 00349 TEMP = ONE/A(K,K) 00350 DO 320 I = 1,M 00351 B(I,K) = TEMP*B(I,K) 00352 320 CONTINUE 00353 END IF 00354 DO 340 J = K + 1,N 00355 IF (A(J,K).NE.ZERO) THEN 00356 TEMP = A(J,K) 00357 DO 330 I = 1,M 00358 B(I,J) = B(I,J) - TEMP*B(I,K) 00359 330 CONTINUE 00360 END IF 00361 340 CONTINUE 00362 IF (ALPHA.NE.ONE) THEN 00363 DO 350 I = 1,M 00364 B(I,K) = ALPHA*B(I,K) 00365 350 CONTINUE 00366 END IF 00367 360 CONTINUE 00368 END IF 00369 END IF 00370 END IF 00371 * 00372 RETURN 00373 * 00374 * End of DTRSM . 00375 * 00376 END
1.7.3 
