001:       SUBROUTINE CTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
002: *     .. Scalar Arguments ..
003:       COMPLEX ALPHA
004:       INTEGER LDA,LDB,M,N
005:       CHARACTER DIAG,SIDE,TRANSA,UPLO
006: *     ..
007: *     .. Array Arguments ..
008:       COMPLEX A(LDA,*),B(LDB,*)
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  CTRSM  solves one of the matrix equations
015: *
016: *     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
017: *
018: *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
019: *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
020: *
021: *     op( A ) = A   or   op( A ) = A'   or   op( A ) = conjg( A' ).
022: *
023: *  The matrix X is overwritten on B.
024: *
025: *  Arguments
026: *  ==========
027: *
028: *  SIDE   - CHARACTER*1.
029: *           On entry, SIDE specifies whether op( A ) appears on the left
030: *           or right of X as follows:
031: *
032: *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
033: *
034: *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
035: *
036: *           Unchanged on exit.
037: *
038: *  UPLO   - CHARACTER*1.
039: *           On entry, UPLO specifies whether the matrix A is an upper or
040: *           lower triangular matrix as follows:
041: *
042: *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
043: *
044: *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
045: *
046: *           Unchanged on exit.
047: *
048: *  TRANSA - CHARACTER*1.
049: *           On entry, TRANSA specifies the form of op( A ) to be used in
050: *           the matrix multiplication as follows:
051: *
052: *              TRANSA = 'N' or 'n'   op( A ) = A.
053: *
054: *              TRANSA = 'T' or 't'   op( A ) = A'.
055: *
056: *              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ).
057: *
058: *           Unchanged on exit.
059: *
060: *  DIAG   - CHARACTER*1.
061: *           On entry, DIAG specifies whether or not A is unit triangular
062: *           as follows:
063: *
064: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
065: *
066: *              DIAG = 'N' or 'n'   A is not assumed to be unit
067: *                                  triangular.
068: *
069: *           Unchanged on exit.
070: *
071: *  M      - INTEGER.
072: *           On entry, M specifies the number of rows of B. M must be at
073: *           least zero.
074: *           Unchanged on exit.
075: *
076: *  N      - INTEGER.
077: *           On entry, N specifies the number of columns of B.  N must be
078: *           at least zero.
079: *           Unchanged on exit.
080: *
081: *  ALPHA  - COMPLEX         .
082: *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
083: *           zero then  A is not referenced and  B need not be set before
084: *           entry.
085: *           Unchanged on exit.
086: *
087: *  A      - COMPLEX          array of DIMENSION ( LDA, k ), where k is m
088: *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
089: *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
090: *           upper triangular part of the array  A must contain the upper
091: *           triangular matrix  and the strictly lower triangular part of
092: *           A is not referenced.
093: *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
094: *           lower triangular part of the array  A must contain the lower
095: *           triangular matrix  and the strictly upper triangular part of
096: *           A is not referenced.
097: *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
098: *           A  are not referenced either,  but are assumed to be  unity.
099: *           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.0) THEN
178:           INFO = 5
179:       ELSE IF (N.LT.0) THEN
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.0) THEN
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.0) RETURN
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' )*B
249: *           or    B := alpha*inv( conjg( A' ) )*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' )
340: *           or    B := alpha*B*inv( conjg( A' ) ).
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
411: