|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) 00002 * .. Scalar Arguments .. 00003 COMPLEX ALPHA 00004 INTEGER LDA,LDB,M,N 00005 CHARACTER DIAG,SIDE,TRANSA,UPLO 00006 * .. 00007 * .. Array Arguments .. 00008 COMPLEX A(LDA,*),B(LDB,*) 00009 * .. 00010 * 00011 * Purpose 00012 * ======= 00013 * 00014 * CTRSM 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' or op( A ) = conjg( 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 ) = conjg( 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 - COMPLEX . 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 - COMPLEX 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 - COMPLEX 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 * -- Written on 8-February-1989. 00125 * Jack Dongarra, Argonne National Laboratory. 00126 * Iain Duff, AERE Harwell. 00127 * Jeremy Du Croz, Numerical Algorithms Group Ltd. 00128 * Sven Hammarling, Numerical Algorithms Group Ltd. 00129 * 00130 * ===================================================================== 00131 * 00132 * .. External Functions .. 00133 LOGICAL LSAME 00134 EXTERNAL LSAME 00135 * .. 00136 * .. External Subroutines .. 00137 EXTERNAL XERBLA 00138 * .. 00139 * .. Intrinsic Functions .. 00140 INTRINSIC CONJG,MAX 00141 * .. 00142 * .. Local Scalars .. 00143 COMPLEX TEMP 00144 INTEGER I,INFO,J,K,NROWA 00145 LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER 00146 * .. 00147 * .. Parameters .. 00148 COMPLEX ONE 00149 PARAMETER (ONE= (1.0E+0,0.0E+0)) 00150 COMPLEX ZERO 00151 PARAMETER (ZERO= (0.0E+0,0.0E+0)) 00152 * .. 00153 * 00154 * Test the input parameters. 00155 * 00156 LSIDE = LSAME(SIDE,'L') 00157 IF (LSIDE) THEN 00158 NROWA = M 00159 ELSE 00160 NROWA = N 00161 END IF 00162 NOCONJ = LSAME(TRANSA,'T') 00163 NOUNIT = LSAME(DIAG,'N') 00164 UPPER = LSAME(UPLO,'U') 00165 * 00166 INFO = 0 00167 IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN 00168 INFO = 1 00169 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN 00170 INFO = 2 00171 ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND. 00172 + (.NOT.LSAME(TRANSA,'T')) .AND. 00173 + (.NOT.LSAME(TRANSA,'C'))) THEN 00174 INFO = 3 00175 ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN 00176 INFO = 4 00177 ELSE IF (M.LT.0) THEN 00178 INFO = 5 00179 ELSE IF (N.LT.0) THEN 00180 INFO = 6 00181 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN 00182 INFO = 9 00183 ELSE IF (LDB.LT.MAX(1,M)) THEN 00184 INFO = 11 00185 END IF 00186 IF (INFO.NE.0) THEN 00187 CALL XERBLA('CTRSM ',INFO) 00188 RETURN 00189 END IF 00190 * 00191 * Quick return if possible. 00192 * 00193 IF (M.EQ.0 .OR. N.EQ.0) RETURN 00194 * 00195 * And when alpha.eq.zero. 00196 * 00197 IF (ALPHA.EQ.ZERO) THEN 00198 DO 20 J = 1,N 00199 DO 10 I = 1,M 00200 B(I,J) = ZERO 00201 10 CONTINUE 00202 20 CONTINUE 00203 RETURN 00204 END IF 00205 * 00206 * Start the operations. 00207 * 00208 IF (LSIDE) THEN 00209 IF (LSAME(TRANSA,'N')) THEN 00210 * 00211 * Form B := alpha*inv( A )*B. 00212 * 00213 IF (UPPER) THEN 00214 DO 60 J = 1,N 00215 IF (ALPHA.NE.ONE) THEN 00216 DO 30 I = 1,M 00217 B(I,J) = ALPHA*B(I,J) 00218 30 CONTINUE 00219 END IF 00220 DO 50 K = M,1,-1 00221 IF (B(K,J).NE.ZERO) THEN 00222 IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) 00223 DO 40 I = 1,K - 1 00224 B(I,J) = B(I,J) - B(K,J)*A(I,K) 00225 40 CONTINUE 00226 END IF 00227 50 CONTINUE 00228 60 CONTINUE 00229 ELSE 00230 DO 100 J = 1,N 00231 IF (ALPHA.NE.ONE) THEN 00232 DO 70 I = 1,M 00233 B(I,J) = ALPHA*B(I,J) 00234 70 CONTINUE 00235 END IF 00236 DO 90 K = 1,M 00237 IF (B(K,J).NE.ZERO) THEN 00238 IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) 00239 DO 80 I = K + 1,M 00240 B(I,J) = B(I,J) - B(K,J)*A(I,K) 00241 80 CONTINUE 00242 END IF 00243 90 CONTINUE 00244 100 CONTINUE 00245 END IF 00246 ELSE 00247 * 00248 * Form B := alpha*inv( A' )*B 00249 * or B := alpha*inv( conjg( A' ) )*B. 00250 * 00251 IF (UPPER) THEN 00252 DO 140 J = 1,N 00253 DO 130 I = 1,M 00254 TEMP = ALPHA*B(I,J) 00255 IF (NOCONJ) THEN 00256 DO 110 K = 1,I - 1 00257 TEMP = TEMP - A(K,I)*B(K,J) 00258 110 CONTINUE 00259 IF (NOUNIT) TEMP = TEMP/A(I,I) 00260 ELSE 00261 DO 120 K = 1,I - 1 00262 TEMP = TEMP - CONJG(A(K,I))*B(K,J) 00263 120 CONTINUE 00264 IF (NOUNIT) TEMP = TEMP/CONJG(A(I,I)) 00265 END IF 00266 B(I,J) = TEMP 00267 130 CONTINUE 00268 140 CONTINUE 00269 ELSE 00270 DO 180 J = 1,N 00271 DO 170 I = M,1,-1 00272 TEMP = ALPHA*B(I,J) 00273 IF (NOCONJ) THEN 00274 DO 150 K = I + 1,M 00275 TEMP = TEMP - A(K,I)*B(K,J) 00276 150 CONTINUE 00277 IF (NOUNIT) TEMP = TEMP/A(I,I) 00278 ELSE 00279 DO 160 K = I + 1,M 00280 TEMP = TEMP - CONJG(A(K,I))*B(K,J) 00281 160 CONTINUE 00282 IF (NOUNIT) TEMP = TEMP/CONJG(A(I,I)) 00283 END IF 00284 B(I,J) = TEMP 00285 170 CONTINUE 00286 180 CONTINUE 00287 END IF 00288 END IF 00289 ELSE 00290 IF (LSAME(TRANSA,'N')) THEN 00291 * 00292 * Form B := alpha*B*inv( A ). 00293 * 00294 IF (UPPER) THEN 00295 DO 230 J = 1,N 00296 IF (ALPHA.NE.ONE) THEN 00297 DO 190 I = 1,M 00298 B(I,J) = ALPHA*B(I,J) 00299 190 CONTINUE 00300 END IF 00301 DO 210 K = 1,J - 1 00302 IF (A(K,J).NE.ZERO) THEN 00303 DO 200 I = 1,M 00304 B(I,J) = B(I,J) - A(K,J)*B(I,K) 00305 200 CONTINUE 00306 END IF 00307 210 CONTINUE 00308 IF (NOUNIT) THEN 00309 TEMP = ONE/A(J,J) 00310 DO 220 I = 1,M 00311 B(I,J) = TEMP*B(I,J) 00312 220 CONTINUE 00313 END IF 00314 230 CONTINUE 00315 ELSE 00316 DO 280 J = N,1,-1 00317 IF (ALPHA.NE.ONE) THEN 00318 DO 240 I = 1,M 00319 B(I,J) = ALPHA*B(I,J) 00320 240 CONTINUE 00321 END IF 00322 DO 260 K = J + 1,N 00323 IF (A(K,J).NE.ZERO) THEN 00324 DO 250 I = 1,M 00325 B(I,J) = B(I,J) - A(K,J)*B(I,K) 00326 250 CONTINUE 00327 END IF 00328 260 CONTINUE 00329 IF (NOUNIT) THEN 00330 TEMP = ONE/A(J,J) 00331 DO 270 I = 1,M 00332 B(I,J) = TEMP*B(I,J) 00333 270 CONTINUE 00334 END IF 00335 280 CONTINUE 00336 END IF 00337 ELSE 00338 * 00339 * Form B := alpha*B*inv( A' ) 00340 * or B := alpha*B*inv( conjg( A' ) ). 00341 * 00342 IF (UPPER) THEN 00343 DO 330 K = N,1,-1 00344 IF (NOUNIT) THEN 00345 IF (NOCONJ) THEN 00346 TEMP = ONE/A(K,K) 00347 ELSE 00348 TEMP = ONE/CONJG(A(K,K)) 00349 END IF 00350 DO 290 I = 1,M 00351 B(I,K) = TEMP*B(I,K) 00352 290 CONTINUE 00353 END IF 00354 DO 310 J = 1,K - 1 00355 IF (A(J,K).NE.ZERO) THEN 00356 IF (NOCONJ) THEN 00357 TEMP = A(J,K) 00358 ELSE 00359 TEMP = CONJG(A(J,K)) 00360 END IF 00361 DO 300 I = 1,M 00362 B(I,J) = B(I,J) - TEMP*B(I,K) 00363 300 CONTINUE 00364 END IF 00365 310 CONTINUE 00366 IF (ALPHA.NE.ONE) THEN 00367 DO 320 I = 1,M 00368 B(I,K) = ALPHA*B(I,K) 00369 320 CONTINUE 00370 END IF 00371 330 CONTINUE 00372 ELSE 00373 DO 380 K = 1,N 00374 IF (NOUNIT) THEN 00375 IF (NOCONJ) THEN 00376 TEMP = ONE/A(K,K) 00377 ELSE 00378 TEMP = ONE/CONJG(A(K,K)) 00379 END IF 00380 DO 340 I = 1,M 00381 B(I,K) = TEMP*B(I,K) 00382 340 CONTINUE 00383 END IF 00384 DO 360 J = K + 1,N 00385 IF (A(J,K).NE.ZERO) THEN 00386 IF (NOCONJ) THEN 00387 TEMP = A(J,K) 00388 ELSE 00389 TEMP = CONJG(A(J,K)) 00390 END IF 00391 DO 350 I = 1,M 00392 B(I,J) = B(I,J) - TEMP*B(I,K) 00393 350 CONTINUE 00394 END IF 00395 360 CONTINUE 00396 IF (ALPHA.NE.ONE) THEN 00397 DO 370 I = 1,M 00398 B(I,K) = ALPHA*B(I,K) 00399 370 CONTINUE 00400 END IF 00401 380 CONTINUE 00402 END IF 00403 END IF 00404 END IF 00405 * 00406 RETURN 00407 * 00408 * End of CTRSM . 00409 * 00410 END
1.7.3 
