share/lapack/blas/fortran/dtrsm.f

   1       SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
   2      $                   B, LDB )
   3 *     .. Scalar Arguments ..
   4       CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
   5       INTEGER            M, N, LDA, LDB
   6       DOUBLE PRECISION   ALPHA
   7 *     .. Array Arguments ..
   8       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  DTRSM  solves one of the matrix equations
  15 *
  16 *     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
  17 *
  18 *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
  19 *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  20 *
  21 *     op( A ) = A   or   op( A ) = A'.
  22 *
  23 *  The matrix X is overwritten on B.
  24 *
  25 *  Parameters
  26 *  ==========
  27 *
  28 *  SIDE   - CHARACTER*1.
  29 *           On entry, SIDE specifies whether op( A ) appears on the left
  30 *           or right of X as follows:
  31 *
  32 *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
  33 *
  34 *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
  35 *
  36 *           Unchanged on exit.
  37 *
  38 *  UPLO   - CHARACTER*1.
  39 *           On entry, UPLO specifies whether the matrix A is an upper or
  40 *           lower triangular matrix as follows:
  41 *
  42 *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  43 *
  44 *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  45 *
  46 *           Unchanged on exit.
  47 *
  48 *  TRANSA - CHARACTER*1.
  49 *           On entry, TRANSA specifies the form of op( A ) to be used in
  50 *           the matrix multiplication as follows:
  51 *
  52 *              TRANSA = 'N' or 'n'   op( A ) = A.
  53 *
  54 *              TRANSA = 'T' or 't'   op( A ) = A'.
  55 *
  56 *              TRANSA = 'C' or 'c'   op( A ) = A'.
  57 *
  58 *           Unchanged on exit.
  59 *
  60 *  DIAG   - CHARACTER*1.
  61 *           On entry, DIAG specifies whether or not A is unit triangular
  62 *           as follows:
  63 *
  64 *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  65 *
  66 *              DIAG = 'N' or 'n'   A is not assumed to be unit
  67 *                                  triangular.
  68 *
  69 *           Unchanged on exit.
  70 *
  71 *  M      - INTEGER.
  72 *           On entry, M specifies the number of rows of B. M must be at
  73 *           least zero.
  74 *           Unchanged on exit.
  75 *
  76 *  N      - INTEGER.
  77 *           On entry, N specifies the number of columns of B.  N must be
  78 *           at least zero.
  79 *           Unchanged on exit.
  80 *
  81 *  ALPHA  - DOUBLE PRECISION.
  82 *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
  83 *           zero then  A is not referenced and  B need not be set before
  84 *           entry.
  85 *           Unchanged on exit.
  86 *
  87 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
  88 *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
  89 *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
  90 *           upper triangular part of the array  A must contain the upper
  91 *           triangular matrix  and the strictly lower triangular part of
  92 *           A is not referenced.
  93 *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
  94 *           lower triangular part of the array  A must contain the lower
  95 *           triangular matrix  and the strictly upper triangular part of
  96 *           A is not referenced.
  97 *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
  98 *           A  are not referenced either,  but are assumed to be  unity.
  99 *           Unchanged on exit.
 100 *
 101 *  LDA    - INTEGER.
 102 *           On entry, LDA specifies the first dimension of A as declared
 103 *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
 104 *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
 105 *           then LDA must be at least max( 1, n ).
 106 *           Unchanged on exit.
 107 *
 108 *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
 109 *           Before entry,  the leading  m by n part of the array  B must
 110 *           contain  the  right-hand  side  matrix  B,  and  on exit  is
 111 *           overwritten by the solution matrix  X.
 112 *
 113 *  LDB    - INTEGER.
 114 *           On entry, LDB specifies the first dimension of B as declared
 115 *           in  the  calling  (sub)  program.   LDB  must  be  at  least
 116 *           max( 1, m ).
 117 *           Unchanged on exit.
 118 *
 119 *
 120 *  Level 3 Blas routine.
 121 *
 122 *
 123 *  -- Written on 8-February-1989.
 124 *     Jack Dongarra, Argonne National Laboratory.
 125 *     Iain Duff, AERE Harwell.
 126 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 127 *     Sven Hammarling, Numerical Algorithms Group Ltd.
 128 *
 129 *
 130 *     .. External Functions ..
 131       LOGICAL            LSAME
 132       EXTERNAL           LSAME
 133 *     .. External Subroutines ..
 134       EXTERNAL           XERBLA
 135 *     .. Intrinsic Functions ..
 136       INTRINSIC          MAX
 137 *     .. Local Scalars ..
 138       LOGICAL            LSIDE, NOUNIT, UPPER
 139       INTEGER            I, INFO, J, K, NROWA
 140       DOUBLE PRECISION   TEMP
 141 *     .. Parameters ..
 142       DOUBLE PRECISION   ONE         , ZERO
 143       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 144 *     ..
 145 *     .. Executable Statements ..
 146 *
 147 *     Test the input parameters.
 148 *
 149       LSIDE  = LSAME( SIDE  , 'L' )
 150       IF( LSIDE )THEN
 151          NROWA = M
 152       ELSE
 153          NROWA = N
 154       END IF
 155       NOUNIT = LSAME( DIAG  , 'N' )
 156       UPPER  = LSAME( UPLO  , 'U' )
 157 *
 158       INFO   = 0
 159       IF(      ( .NOT.LSIDE                ).AND.
 160      $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
 161          INFO = 1
 162       ELSE IF( ( .NOT.UPPER                ).AND.
 163      $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
 164          INFO = 2
 165       ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
 166      $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
 167      $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
 168          INFO = 3
 169       ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
 170      $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
 171          INFO = 4
 172       ELSE IF( M  .LT.0               )THEN
 173          INFO = 5
 174       ELSE IF( N  .LT.0               )THEN
 175          INFO = 6
 176       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
 177          INFO = 9
 178       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
 179          INFO = 11
 180       END IF
 181       IF( INFO.NE.0 )THEN
 182          CALL XERBLA( 'DTRSM ', INFO )
 183          RETURN
 184       END IF
 185 *
 186 *     Quick return if possible.
 187 *
 188       IF( N.EQ.0 )
 189      $   RETURN
 190 *
 191 *     And when  alpha.eq.zero.
 192 *
 193       IF( ALPHA.EQ.ZERO )THEN
 194          DO 20, J = 1, N
 195             DO 10, I = 1, M
 196                B( I, J ) = ZERO
 197    10       CONTINUE
 198    20    CONTINUE
 199          RETURN
 200       END IF
 201 *
 202 *     Start the operations.
 203 *
 204       IF( LSIDE )THEN
 205          IF( LSAME( TRANSA, 'N' ) )THEN
 206 *
 207 *           Form  B := alpha*inv( A )*B.
 208 *
 209             IF( UPPER )THEN
 210                DO 60, J = 1, N
 211                   IF( ALPHA.NE.ONE )THEN
 212                      DO 30, I = 1, M
 213                         B( I, J ) = ALPHA*B( I, J )
 214    30                CONTINUE
 215                   END IF
 216                   DO 50, K = M, 1, -1
 217                      IF( B( K, J ).NE.ZERO )THEN
 218                         IF( NOUNIT )
 219      $                     B( K, J ) = B( K, J )/A( K, K )
 220                         DO 40, I = 1, K - 1
 221                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
 222    40                   CONTINUE
 223                      END IF
 224    50             CONTINUE
 225    60          CONTINUE
 226             ELSE
 227                DO 100, J = 1, N
 228                   IF( ALPHA.NE.ONE )THEN
 229                      DO 70, I = 1, M
 230                         B( I, J ) = ALPHA*B( I, J )
 231    70                CONTINUE
 232                   END IF
 233                   DO 90 K = 1, M
 234                      IF( B( K, J ).NE.ZERO )THEN
 235                         IF( NOUNIT )
 236      $                     B( K, J ) = B( K, J )/A( K, K )
 237                         DO 80, I = K + 1, M
 238                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
 239    80                   CONTINUE
 240                      END IF
 241    90             CONTINUE
 242   100          CONTINUE
 243             END IF
 244          ELSE
 245 *
 246 *           Form  B := alpha*inv( A' )*B.
 247 *
 248             IF( UPPER )THEN
 249                DO 130, J = 1, N
 250                   DO 120, I = 1, M
 251                      TEMP = ALPHA*B( I, J )
 252                      DO 110, K = 1, I - 1
 253                         TEMP = TEMP - A( K, I )*B( K, J )
 254   110                CONTINUE
 255                      IF( NOUNIT )
 256      $                  TEMP = TEMP/A( I, I )
 257                      B( I, J ) = TEMP
 258   120             CONTINUE
 259   130          CONTINUE
 260             ELSE
 261                DO 160, J = 1, N
 262                   DO 150, I = M, 1, -1
 263                      TEMP = ALPHA*B( I, J )
 264                      DO 140, K = I + 1, M
 265                         TEMP = TEMP - A( K, I )*B( K, J )
 266   140                CONTINUE
 267                      IF( NOUNIT )
 268      $                  TEMP = TEMP/A( I, I )
 269                      B( I, J ) = TEMP
 270   150             CONTINUE
 271   160          CONTINUE
 272             END IF
 273          END IF
 274       ELSE
 275          IF( LSAME( TRANSA, 'N' ) )THEN
 276 *
 277 *           Form  B := alpha*B*inv( A ).
 278 *
 279             IF( UPPER )THEN
 280                DO 210, J = 1, N
 281                   IF( ALPHA.NE.ONE )THEN
 282                      DO 170, I = 1, M
 283                         B( I, J ) = ALPHA*B( I, J )
 284   170                CONTINUE
 285                   END IF
 286                   DO 190, K = 1, J - 1
 287                      IF( A( K, J ).NE.ZERO )THEN
 288                         DO 180, I = 1, M
 289                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
 290   180                   CONTINUE
 291                      END IF
 292   190             CONTINUE
 293                   IF( NOUNIT )THEN
 294                      TEMP = ONE/A( J, J )
 295                      DO 200, I = 1, M
 296                         B( I, J ) = TEMP*B( I, J )
 297   200                CONTINUE
 298                   END IF
 299   210          CONTINUE
 300             ELSE
 301                DO 260, J = N, 1, -1
 302                   IF( ALPHA.NE.ONE )THEN
 303                      DO 220, I = 1, M
 304                         B( I, J ) = ALPHA*B( I, J )
 305   220                CONTINUE
 306                   END IF
 307                   DO 240, K = J + 1, N
 308                      IF( A( K, J ).NE.ZERO )THEN
 309                         DO 230, I = 1, M
 310                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
 311   230                   CONTINUE
 312                      END IF
 313   240             CONTINUE
 314                   IF( NOUNIT )THEN
 315                      TEMP = ONE/A( J, J )
 316                      DO 250, I = 1, M
 317                        B( I, J ) = TEMP*B( I, J )
 318   250                CONTINUE
 319                   END IF
 320   260          CONTINUE
 321             END IF
 322          ELSE
 323 *
 324 *           Form  B := alpha*B*inv( A' ).
 325 *
 326             IF( UPPER )THEN
 327                DO 310, K = N, 1, -1
 328                   IF( NOUNIT )THEN
 329                      TEMP = ONE/A( K, K )
 330                      DO 270, I = 1, M
 331                         B( I, K ) = TEMP*B( I, K )
 332   270                CONTINUE
 333                   END IF
 334                   DO 290, J = 1, K - 1
 335                      IF( A( J, K ).NE.ZERO )THEN
 336                         TEMP = A( J, K )
 337                         DO 280, I = 1, M
 338                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
 339   280                   CONTINUE
 340                      END IF
 341   290             CONTINUE
 342                   IF( ALPHA.NE.ONE )THEN
 343                      DO 300, I = 1, M
 344                         B( I, K ) = ALPHA*B( I, K )
 345   300                CONTINUE
 346                   END IF
 347   310          CONTINUE
 348             ELSE
 349                DO 360, K = 1, N
 350                   IF( NOUNIT )THEN
 351                      TEMP = ONE/A( K, K )
 352                      DO 320, I = 1, M
 353                         B( I, K ) = TEMP*B( I, K )
 354   320                CONTINUE
 355                   END IF
 356                   DO 340, J = K + 1, N
 357                      IF( A( J, K ).NE.ZERO )THEN
 358                         TEMP = A( J, K )
 359                         DO 330, I = 1, M
 360                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
 361   330                   CONTINUE
 362                      END IF
 363   340             CONTINUE
 364                   IF( ALPHA.NE.ONE )THEN
 365                      DO 350, I = 1, M
 366                         B( I, K ) = ALPHA*B( I, K )
 367   350                CONTINUE
 368                   END IF
 369   360          CONTINUE
 370             END IF
 371          END IF
 372       END IF
 373 *
 374       RETURN
 375 *
 376 *     End of DTRSM .
 377 *
 378       END