BLAS/SRC/ctrsm.f

   1 *> \brief \b CTRSM
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE CTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
  12 *
  13 *       .. Scalar Arguments ..
  14 *       COMPLEX ALPHA
  15 *       INTEGER LDA,LDB,M,N
  16 *       CHARACTER DIAG,SIDE,TRANSA,UPLO
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       COMPLEX A(LDA,*),B(LDB,*)
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> CTRSM  solves one of the matrix equations
  29 *>
  30 *>    op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
  31 *>
  32 *> where alpha is a scalar, X and B are m by n matrices, A is a unit, or
  33 *> non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  34 *>
  35 *>    op( A ) = A   or   op( A ) = A**T   or   op( A ) = A**H.
  36 *>
  37 *> The matrix X is overwritten on B.
  38 *> \endverbatim
  39 *
  40 *  Arguments:
  41 *  ==========
  42 *
  43 *> \param[in] SIDE
  44 *> \verbatim
  45 *>          SIDE is CHARACTER*1
  46 *>           On entry, SIDE specifies whether op( A ) appears on the left
  47 *>           or right of X as follows:
  48 *>
  49 *>              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
  50 *>
  51 *>              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
  52 *> \endverbatim
  53 *>
  54 *> \param[in] UPLO
  55 *> \verbatim
  56 *>          UPLO is CHARACTER*1
  57 *>           On entry, UPLO specifies whether the matrix A is an upper or
  58 *>           lower triangular matrix as follows:
  59 *>
  60 *>              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  61 *>
  62 *>              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] TRANSA
  66 *> \verbatim
  67 *>          TRANSA is CHARACTER*1
  68 *>           On entry, TRANSA specifies the form of op( A ) to be used in
  69 *>           the matrix multiplication as follows:
  70 *>
  71 *>              TRANSA = 'N' or 'n'   op( A ) = A.
  72 *>
  73 *>              TRANSA = 'T' or 't'   op( A ) = A**T.
  74 *>
  75 *>              TRANSA = 'C' or 'c'   op( A ) = A**H.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] DIAG
  79 *> \verbatim
  80 *>          DIAG is CHARACTER*1
  81 *>           On entry, DIAG specifies whether or not A is unit triangular
  82 *>           as follows:
  83 *>
  84 *>              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  85 *>
  86 *>              DIAG = 'N' or 'n'   A is not assumed to be unit
  87 *>                                  triangular.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] M
  91 *> \verbatim
  92 *>          M is INTEGER
  93 *>           On entry, M specifies the number of rows of B. M must be at
  94 *>           least zero.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] N
  98 *> \verbatim
  99 *>          N is INTEGER
 100 *>           On entry, N specifies the number of columns of B.  N must be
 101 *>           at least zero.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] ALPHA
 105 *> \verbatim
 106 *>          ALPHA is COMPLEX
 107 *>           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
 108 *>           zero then  A is not referenced and  B need not be set before
 109 *>           entry.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] A
 113 *> \verbatim
 114 *>          A is COMPLEX array of DIMENSION ( LDA, k ),
 115 *>           where k is m when SIDE = 'L' or 'l'
 116 *>             and k is n when SIDE = 'R' or 'r'.
 117 *>           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
 118 *>           upper triangular part of the array  A must contain the upper
 119 *>           triangular matrix  and the strictly lower triangular part of
 120 *>           A is not referenced.
 121 *>           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
 122 *>           lower triangular part of the array  A must contain the lower
 123 *>           triangular matrix  and the strictly upper triangular part of
 124 *>           A is not referenced.
 125 *>           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
 126 *>           A  are not referenced either,  but are assumed to be  unity.
 127 *> \endverbatim
 128 *>
 129 *> \param[in] LDA
 130 *> \verbatim
 131 *>          LDA is INTEGER
 132 *>           On entry, LDA specifies the first dimension of A as declared
 133 *>           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
 134 *>           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
 135 *>           then LDA must be at least max( 1, n ).
 136 *> \endverbatim
 137 *>
 138 *> \param[in,out] B
 139 *> \verbatim
 140 *>          B is COMPLEX array of DIMENSION ( LDB, n ).
 141 *>           Before entry,  the leading  m by n part of the array  B must
 142 *>           contain  the  right-hand  side  matrix  B,  and  on exit  is
 143 *>           overwritten by the solution matrix  X.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] LDB
 147 *> \verbatim
 148 *>          LDB is INTEGER
 149 *>           On entry, LDB specifies the first dimension of B as declared
 150 *>           in  the  calling  (sub)  program.   LDB  must  be  at  least
 151 *>           max( 1, m ).
 152 *> \endverbatim
 153 *
 154 *  Authors:
 155 *  ========
 156 *
 157 *> \author Univ. of Tennessee
 158 *> \author Univ. of California Berkeley
 159 *> \author Univ. of Colorado Denver
 160 *> \author NAG Ltd.
 161 *
 162 *> \date November 2011
 163 *
 164 *> \ingroup complex_blas_level3
 165 *
 166 *> \par Further Details:
 167 *  =====================
 168 *>
 169 *> \verbatim
 170 *>
 171 *>  Level 3 Blas routine.
 172 *>
 173 *>  -- Written on 8-February-1989.
 174 *>     Jack Dongarra, Argonne National Laboratory.
 175 *>     Iain Duff, AERE Harwell.
 176 *>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 177 *>     Sven Hammarling, Numerical Algorithms Group Ltd.
 178 *> \endverbatim
 179 *>
 180 *  =====================================================================
 181       SUBROUTINE CTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
 182 *
 183 *  -- Reference BLAS level3 routine (version 3.4.0) --
 184 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 185 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 186 *     November 2011
 187 *
 188 *     .. Scalar Arguments ..
 189       COMPLEX ALPHA
 190       INTEGER LDA,LDB,M,N
 191       CHARACTER DIAG,SIDE,TRANSA,UPLO
 192 *     ..
 193 *     .. Array Arguments ..
 194       COMPLEX A(LDA,*),B(LDB,*)
 195 *     ..
 196 *
 197 *  =====================================================================
 198 *
 199 *     .. External Functions ..
 200       LOGICAL LSAME
 201       EXTERNAL LSAME
 202 *     ..
 203 *     .. External Subroutines ..
 204       EXTERNAL XERBLA
 205 *     ..
 206 *     .. Intrinsic Functions ..
 207       INTRINSIC CONJG,MAX
 208 *     ..
 209 *     .. Local Scalars ..
 210       COMPLEX TEMP
 211       INTEGER I,INFO,J,K,NROWA
 212       LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER
 213 *     ..
 214 *     .. Parameters ..
 215       COMPLEX ONE
 216       PARAMETER (ONE= (1.0E+0,0.0E+0))
 217       COMPLEX ZERO
 218       PARAMETER (ZERO= (0.0E+0,0.0E+0))
 219 *     ..
 220 *
 221 *     Test the input parameters.
 222 *
 223       LSIDE = LSAME(SIDE,'L')
 224       IF (LSIDE) THEN
 225           NROWA = M
 226       ELSE
 227           NROWA = N
 228       END IF
 229       NOCONJ = LSAME(TRANSA,'T')
 230       NOUNIT = LSAME(DIAG,'N')
 231       UPPER = LSAME(UPLO,'U')
 232 *
 233       INFO = 0
 234       IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
 235           INFO = 1
 236       ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
 237           INFO = 2
 238       ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
 239      +         (.NOT.LSAME(TRANSA,'T')) .AND.
 240      +         (.NOT.LSAME(TRANSA,'C'))) THEN
 241           INFO = 3
 242       ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
 243           INFO = 4
 244       ELSE IF (M.LT.0) THEN
 245           INFO = 5
 246       ELSE IF (N.LT.0) THEN
 247           INFO = 6
 248       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
 249           INFO = 9
 250       ELSE IF (LDB.LT.MAX(1,M)) THEN
 251           INFO = 11
 252       END IF
 253       IF (INFO.NE.0) THEN
 254           CALL XERBLA('CTRSM ',INFO)
 255           RETURN
 256       END IF
 257 *
 258 *     Quick return if possible.
 259 *
 260       IF (M.EQ.0 .OR. N.EQ.0) RETURN
 261 *
 262 *     And when  alpha.eq.zero.
 263 *
 264       IF (ALPHA.EQ.ZERO) THEN
 265           DO 20 J = 1,N
 266               DO 10 I = 1,M
 267                   B(I,J) = ZERO
 268    10         CONTINUE
 269    20     CONTINUE
 270           RETURN
 271       END IF
 272 *
 273 *     Start the operations.
 274 *
 275       IF (LSIDE) THEN
 276           IF (LSAME(TRANSA,'N')) THEN
 277 *
 278 *           Form  B := alpha*inv( A )*B.
 279 *
 280               IF (UPPER) THEN
 281                   DO 60 J = 1,N
 282                       IF (ALPHA.NE.ONE) THEN
 283                           DO 30 I = 1,M
 284                               B(I,J) = ALPHA*B(I,J)
 285    30                     CONTINUE
 286                       END IF
 287                       DO 50 K = M,1,-1
 288                           IF (B(K,J).NE.ZERO) THEN
 289                               IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
 290                               DO 40 I = 1,K - 1
 291                                   B(I,J) = B(I,J) - B(K,J)*A(I,K)
 292    40                         CONTINUE
 293                           END IF
 294    50                 CONTINUE
 295    60             CONTINUE
 296               ELSE
 297                   DO 100 J = 1,N
 298                       IF (ALPHA.NE.ONE) THEN
 299                           DO 70 I = 1,M
 300                               B(I,J) = ALPHA*B(I,J)
 301    70                     CONTINUE
 302                       END IF
 303                       DO 90 K = 1,M
 304                           IF (B(K,J).NE.ZERO) THEN
 305                               IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
 306                               DO 80 I = K + 1,M
 307                                   B(I,J) = B(I,J) - B(K,J)*A(I,K)
 308    80                         CONTINUE
 309                           END IF
 310    90                 CONTINUE
 311   100             CONTINUE
 312               END IF
 313           ELSE
 314 *
 315 *           Form  B := alpha*inv( A**T )*B
 316 *           or    B := alpha*inv( A**H )*B.
 317 *
 318               IF (UPPER) THEN
 319                   DO 140 J = 1,N
 320                       DO 130 I = 1,M
 321                           TEMP = ALPHA*B(I,J)
 322                           IF (NOCONJ) THEN
 323                               DO 110 K = 1,I - 1
 324                                   TEMP = TEMP - A(K,I)*B(K,J)
 325   110                         CONTINUE
 326                               IF (NOUNIT) TEMP = TEMP/A(I,I)
 327                           ELSE
 328                               DO 120 K = 1,I - 1
 329                                   TEMP = TEMP - CONJG(A(K,I))*B(K,J)
 330   120                         CONTINUE
 331                               IF (NOUNIT) TEMP = TEMP/CONJG(A(I,I))
 332                           END IF
 333                           B(I,J) = TEMP
 334   130                 CONTINUE
 335   140             CONTINUE
 336               ELSE
 337                   DO 180 J = 1,N
 338                       DO 170 I = M,1,-1
 339                           TEMP = ALPHA*B(I,J)
 340                           IF (NOCONJ) THEN
 341                               DO 150 K = I + 1,M
 342                                   TEMP = TEMP - A(K,I)*B(K,J)
 343   150                         CONTINUE
 344                               IF (NOUNIT) TEMP = TEMP/A(I,I)
 345                           ELSE
 346                               DO 160 K = I + 1,M
 347                                   TEMP = TEMP - CONJG(A(K,I))*B(K,J)
 348   160                         CONTINUE
 349                               IF (NOUNIT) TEMP = TEMP/CONJG(A(I,I))
 350                           END IF
 351                           B(I,J) = TEMP
 352   170                 CONTINUE
 353   180             CONTINUE
 354               END IF
 355           END IF
 356       ELSE
 357           IF (LSAME(TRANSA,'N')) THEN
 358 *
 359 *           Form  B := alpha*B*inv( A ).
 360 *
 361               IF (UPPER) THEN
 362                   DO 230 J = 1,N
 363                       IF (ALPHA.NE.ONE) THEN
 364                           DO 190 I = 1,M
 365                               B(I,J) = ALPHA*B(I,J)
 366   190                     CONTINUE
 367                       END IF
 368                       DO 210 K = 1,J - 1
 369                           IF (A(K,J).NE.ZERO) THEN
 370                               DO 200 I = 1,M
 371                                   B(I,J) = B(I,J) - A(K,J)*B(I,K)
 372   200                         CONTINUE
 373                           END IF
 374   210                 CONTINUE
 375                       IF (NOUNIT) THEN
 376                           TEMP = ONE/A(J,J)
 377                           DO 220 I = 1,M
 378                               B(I,J) = TEMP*B(I,J)
 379   220                     CONTINUE
 380                       END IF
 381   230             CONTINUE
 382               ELSE
 383                   DO 280 J = N,1,-1
 384                       IF (ALPHA.NE.ONE) THEN
 385                           DO 240 I = 1,M
 386                               B(I,J) = ALPHA*B(I,J)
 387   240                     CONTINUE
 388                       END IF
 389                       DO 260 K = J + 1,N
 390                           IF (A(K,J).NE.ZERO) THEN
 391                               DO 250 I = 1,M
 392                                   B(I,J) = B(I,J) - A(K,J)*B(I,K)
 393   250                         CONTINUE
 394                           END IF
 395   260                 CONTINUE
 396                       IF (NOUNIT) THEN
 397                           TEMP = ONE/A(J,J)
 398                           DO 270 I = 1,M
 399                               B(I,J) = TEMP*B(I,J)
 400   270                     CONTINUE
 401                       END IF
 402   280             CONTINUE
 403               END IF
 404           ELSE
 405 *
 406 *           Form  B := alpha*B*inv( A**T )
 407 *           or    B := alpha*B*inv( A**H ).
 408 *
 409               IF (UPPER) THEN
 410                   DO 330 K = N,1,-1
 411                       IF (NOUNIT) THEN
 412                           IF (NOCONJ) THEN
 413                               TEMP = ONE/A(K,K)
 414                           ELSE
 415                               TEMP = ONE/CONJG(A(K,K))
 416                           END IF
 417                           DO 290 I = 1,M
 418                               B(I,K) = TEMP*B(I,K)
 419   290                     CONTINUE
 420                       END IF
 421                       DO 310 J = 1,K - 1
 422                           IF (A(J,K).NE.ZERO) THEN
 423                               IF (NOCONJ) THEN
 424                                   TEMP = A(J,K)
 425                               ELSE
 426                                   TEMP = CONJG(A(J,K))
 427                               END IF
 428                               DO 300 I = 1,M
 429                                   B(I,J) = B(I,J) - TEMP*B(I,K)
 430   300                         CONTINUE
 431                           END IF
 432   310                 CONTINUE
 433                       IF (ALPHA.NE.ONE) THEN
 434                           DO 320 I = 1,M
 435                               B(I,K) = ALPHA*B(I,K)
 436   320                     CONTINUE
 437                       END IF
 438   330             CONTINUE
 439               ELSE
 440                   DO 380 K = 1,N
 441                       IF (NOUNIT) THEN
 442                           IF (NOCONJ) THEN
 443                               TEMP = ONE/A(K,K)
 444                           ELSE
 445                               TEMP = ONE/CONJG(A(K,K))
 446                           END IF
 447                           DO 340 I = 1,M
 448                               B(I,K) = TEMP*B(I,K)
 449   340                     CONTINUE
 450                       END IF
 451                       DO 360 J = K + 1,N
 452                           IF (A(J,K).NE.ZERO) THEN
 453                               IF (NOCONJ) THEN
 454                                   TEMP = A(J,K)
 455                               ELSE
 456                                   TEMP = CONJG(A(J,K))
 457                               END IF
 458                               DO 350 I = 1,M
 459                                   B(I,J) = B(I,J) - TEMP*B(I,K)
 460   350                         CONTINUE
 461                           END IF
 462   360                 CONTINUE
 463                       IF (ALPHA.NE.ONE) THEN
 464                           DO 370 I = 1,M
 465                               B(I,K) = ALPHA*B(I,K)
 466   370                     CONTINUE
 467                       END IF
 468   380             CONTINUE
 469               END IF
 470           END IF
 471       END IF
 472 *
 473       RETURN
 474 *
 475 *     End of CTRSM .
 476 *
 477       END