LAPACK 3.3.1
Linear Algebra PACKage

dtrsm.f

Go to the documentation of this file.
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**T.
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**T.
00055 *
00056 *              TRANSA = 'C' or 'c'   op( A ) = A**T.
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**T )*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**T ).
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
 All Files Functions