LAPACK: BLAS/SRC/ctrsm.f Source File

Go to the documentation of this file.
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