1       SUBROUTINE CTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
  2 *     .. Scalar Arguments ..
  3       COMPLEX ALPHA
  4       INTEGER LDA,LDB,M,N
  5       CHARACTER DIAG,SIDE,TRANSA,UPLO
  6 *     ..
  7 *     .. Array Arguments ..
  8       COMPLEX A(LDA,*),B(LDB,*)
  9 *     ..
 10 *
 11 *  Purpose
 12 *  =======
 13 *
 14 *  CTRSM  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**T   or   op( A ) = A**H.
 22 *
 23 *  The matrix X is overwritten on B.
 24 *
 25 *  Arguments
 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**T.
 55 *
 56 *              TRANSA = 'C' or 'c'   op( A ) = A**H.
 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  - COMPLEX         .
 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      - COMPLEX          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      - COMPLEX          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 *  Further Details
120 *  ===============
121 *
122 *  Level 3 Blas routine.
123 *
124 *  -- Written on 8-February-1989.
125 *     Jack Dongarra, Argonne National Laboratory.
126 *     Iain Duff, AERE Harwell.
127 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
128 *     Sven Hammarling, Numerical Algorithms Group Ltd.
129 *
130 *  =====================================================================
131 *
132 *     .. External Functions ..
133       LOGICAL LSAME
134       EXTERNAL LSAME
135 *     ..
136 *     .. External Subroutines ..
137       EXTERNAL XERBLA
138 *     ..
139 *     .. Intrinsic Functions ..
140       INTRINSIC CONJG,MAX
141 *     ..
142 *     .. Local Scalars ..
143       COMPLEX TEMP
144       INTEGER I,INFO,J,K,NROWA
145       LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER
146 *     ..
147 *     .. Parameters ..
148       COMPLEX ONE
149       PARAMETER (ONE= (1.0E+0,0.0E+0))
150       COMPLEX ZERO
151       PARAMETER (ZERO= (0.0E+0,0.0E+0))
152 *     ..
153 *
154 *     Test the input parameters.
155 *
156       LSIDE = LSAME(SIDE,'L')
157       IF (LSIDE) THEN
158           NROWA = M
159       ELSE
160           NROWA = N
161       END IF
162       NOCONJ = LSAME(TRANSA,'T')
163       NOUNIT = LSAME(DIAG,'N')
164       UPPER = LSAME(UPLO,'U')
165 *
166       INFO = 0
167       IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
168           INFO = 1
169       ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
170           INFO = 2
171       ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
172      +         (.NOT.LSAME(TRANSA,'T')) .AND.
173      +         (.NOT.LSAME(TRANSA,'C'))) THEN
174           INFO = 3
175       ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
176           INFO = 4
177       ELSE IF (M.LT.0THEN
178           INFO = 5
179       ELSE IF (N.LT.0THEN
180           INFO = 6
181       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
182           INFO = 9
183       ELSE IF (LDB.LT.MAX(1,M)) THEN
184           INFO = 11
185       END IF
186       IF (INFO.NE.0THEN
187           CALL XERBLA('CTRSM ',INFO)
188           RETURN
189       END IF
190 *
191 *     Quick return if possible.
192 *
193       IF (M.EQ.0 .OR. N.EQ.0RETURN
194 *
195 *     And when  alpha.eq.zero.
196 *
197       IF (ALPHA.EQ.ZERO) THEN
198           DO 20 J = 1,N
199               DO 10 I = 1,M
200                   B(I,J) = ZERO
201    10         CONTINUE
202    20     CONTINUE
203           RETURN
204       END IF
205 *
206 *     Start the operations.
207 *
208       IF (LSIDE) THEN
209           IF (LSAME(TRANSA,'N')) THEN
210 *
211 *           Form  B := alpha*inv( A )*B.
212 *
213               IF (UPPER) THEN
214                   DO 60 J = 1,N
215                       IF (ALPHA.NE.ONE) THEN
216                           DO 30 I = 1,M
217                               B(I,J) = ALPHA*B(I,J)
218    30                     CONTINUE
219                       END IF
220                       DO 50 K = M,1,-1
221                           IF (B(K,J).NE.ZERO) THEN
222                               IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
223                               DO 40 I = 1,K - 1
224                                   B(I,J) = B(I,J) - B(K,J)*A(I,K)
225    40                         CONTINUE
226                           END IF
227    50                 CONTINUE
228    60             CONTINUE
229               ELSE
230                   DO 100 J = 1,N
231                       IF (ALPHA.NE.ONE) THEN
232                           DO 70 I = 1,M
233                               B(I,J) = ALPHA*B(I,J)
234    70                     CONTINUE
235                       END IF
236                       DO 90 K = 1,M
237                           IF (B(K,J).NE.ZERO) THEN
238                               IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
239                               DO 80 I = K + 1,M
240                                   B(I,J) = B(I,J) - B(K,J)*A(I,K)
241    80                         CONTINUE
242                           END IF
243    90                 CONTINUE
244   100             CONTINUE
245               END IF
246           ELSE
247 *
248 *           Form  B := alpha*inv( A**T )*B
249 *           or    B := alpha*inv( A**H )*B.
250 *
251               IF (UPPER) THEN
252                   DO 140 J = 1,N
253                       DO 130 I = 1,M
254                           TEMP = ALPHA*B(I,J)
255                           IF (NOCONJ) THEN
256                               DO 110 K = 1,I - 1
257                                   TEMP = TEMP - A(K,I)*B(K,J)
258   110                         CONTINUE
259                               IF (NOUNIT) TEMP = TEMP/A(I,I)
260                           ELSE
261                               DO 120 K = 1,I - 1
262                                   TEMP = TEMP - CONJG(A(K,I))*B(K,J)
263   120                         CONTINUE
264                               IF (NOUNIT) TEMP = TEMP/CONJG(A(I,I))
265                           END IF
266                           B(I,J) = TEMP
267   130                 CONTINUE
268   140             CONTINUE
269               ELSE
270                   DO 180 J = 1,N
271                       DO 170 I = M,1,-1
272                           TEMP = ALPHA*B(I,J)
273                           IF (NOCONJ) THEN
274                               DO 150 K = I + 1,M
275                                   TEMP = TEMP - A(K,I)*B(K,J)
276   150                         CONTINUE
277                               IF (NOUNIT) TEMP = TEMP/A(I,I)
278                           ELSE
279                               DO 160 K = I + 1,M
280                                   TEMP = TEMP - CONJG(A(K,I))*B(K,J)
281   160                         CONTINUE
282                               IF (NOUNIT) TEMP = TEMP/CONJG(A(I,I))
283                           END IF
284                           B(I,J) = TEMP
285   170                 CONTINUE
286   180             CONTINUE
287               END IF
288           END IF
289       ELSE
290           IF (LSAME(TRANSA,'N')) THEN
291 *
292 *           Form  B := alpha*B*inv( A ).
293 *
294               IF (UPPER) THEN
295                   DO 230 J = 1,N
296                       IF (ALPHA.NE.ONE) THEN
297                           DO 190 I = 1,M
298                               B(I,J) = ALPHA*B(I,J)
299   190                     CONTINUE
300                       END IF
301                       DO 210 K = 1,J - 1
302                           IF (A(K,J).NE.ZERO) THEN
303                               DO 200 I = 1,M
304                                   B(I,J) = B(I,J) - A(K,J)*B(I,K)
305   200                         CONTINUE
306                           END IF
307   210                 CONTINUE
308                       IF (NOUNIT) THEN
309                           TEMP = ONE/A(J,J)
310                           DO 220 I = 1,M
311                               B(I,J) = TEMP*B(I,J)
312   220                     CONTINUE
313                       END IF
314   230             CONTINUE
315               ELSE
316                   DO 280 J = N,1,-1
317                       IF (ALPHA.NE.ONE) THEN
318                           DO 240 I = 1,M
319                               B(I,J) = ALPHA*B(I,J)
320   240                     CONTINUE
321                       END IF
322                       DO 260 K = J + 1,N
323                           IF (A(K,J).NE.ZERO) THEN
324                               DO 250 I = 1,M
325                                   B(I,J) = B(I,J) - A(K,J)*B(I,K)
326   250                         CONTINUE
327                           END IF
328   260                 CONTINUE
329                       IF (NOUNIT) THEN
330                           TEMP = ONE/A(J,J)
331                           DO 270 I = 1,M
332                               B(I,J) = TEMP*B(I,J)
333   270                     CONTINUE
334                       END IF
335   280             CONTINUE
336               END IF
337           ELSE
338 *
339 *           Form  B := alpha*B*inv( A**T )
340 *           or    B := alpha*B*inv( A**H ).
341 *
342               IF (UPPER) THEN
343                   DO 330 K = N,1,-1
344                       IF (NOUNIT) THEN
345                           IF (NOCONJ) THEN
346                               TEMP = ONE/A(K,K)
347                           ELSE
348                               TEMP = ONE/CONJG(A(K,K))
349                           END IF
350                           DO 290 I = 1,M
351                               B(I,K) = TEMP*B(I,K)
352   290                     CONTINUE
353                       END IF
354                       DO 310 J = 1,K - 1
355                           IF (A(J,K).NE.ZERO) THEN
356                               IF (NOCONJ) THEN
357                                   TEMP = A(J,K)
358                               ELSE
359                                   TEMP = CONJG(A(J,K))
360                               END IF
361                               DO 300 I = 1,M
362                                   B(I,J) = B(I,J) - TEMP*B(I,K)
363   300                         CONTINUE
364                           END IF
365   310                 CONTINUE
366                       IF (ALPHA.NE.ONE) THEN
367                           DO 320 I = 1,M
368                               B(I,K) = ALPHA*B(I,K)
369   320                     CONTINUE
370                       END IF
371   330             CONTINUE
372               ELSE
373                   DO 380 K = 1,N
374                       IF (NOUNIT) THEN
375                           IF (NOCONJ) THEN
376                               TEMP = ONE/A(K,K)
377                           ELSE
378                               TEMP = ONE/CONJG(A(K,K))
379                           END IF
380                           DO 340 I = 1,M
381                               B(I,K) = TEMP*B(I,K)
382   340                     CONTINUE
383                       END IF
384                       DO 360 J = K + 1,N
385                           IF (A(J,K).NE.ZERO) THEN
386                               IF (NOCONJ) THEN
387                                   TEMP = A(J,K)
388                               ELSE
389                                   TEMP = CONJG(A(J,K))
390                               END IF
391                               DO 350 I = 1,M
392                                   B(I,J) = B(I,J) - TEMP*B(I,K)
393   350                         CONTINUE
394                           END IF
395   360                 CONTINUE
396                       IF (ALPHA.NE.ONE) THEN
397                           DO 370 I = 1,M
398                               B(I,K) = ALPHA*B(I,K)
399   370                     CONTINUE
400                       END IF
401   380             CONTINUE
402               END IF
403           END IF
404       END IF
405 *
406       RETURN
407 *
408 *     End of CTRSM .
409 *
410       END