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