001:       SUBROUTINE STRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
002: *     .. Scalar Arguments ..
003:       REAL ALPHA
004:       INTEGER LDA,LDB,M,N
005:       CHARACTER DIAG,SIDE,TRANSA,UPLO
006: *     ..
007: *     .. Array Arguments ..
008:       REAL A(LDA,*),B(LDB,*)
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  STRSM  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'.
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 ) = 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  - REAL            .
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      - REAL             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      - REAL             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: *     ..
134: *     .. External Subroutines ..
135:       EXTERNAL XERBLA
136: *     ..
137: *     .. Intrinsic Functions ..
138:       INTRINSIC MAX
139: *     ..
140: *     .. Local Scalars ..
141:       REAL TEMP
142:       INTEGER I,INFO,J,K,NROWA
143:       LOGICAL LSIDE,NOUNIT,UPPER
144: *     ..
145: *     .. Parameters ..
146:       REAL ONE,ZERO
147:       PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
148: *     ..
149: *
150: *     Test the input parameters.
151: *
152:       LSIDE = LSAME(SIDE,'L')
153:       IF (LSIDE) THEN
154:           NROWA = M
155:       ELSE
156:           NROWA = N
157:       END IF
158:       NOUNIT = LSAME(DIAG,'N')
159:       UPPER = LSAME(UPLO,'U')
160: *
161:       INFO = 0
162:       IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
163:           INFO = 1
164:       ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
165:           INFO = 2
166:       ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
167:      +         (.NOT.LSAME(TRANSA,'T')) .AND.
168:      +         (.NOT.LSAME(TRANSA,'C'))) THEN
169:           INFO = 3
170:       ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.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('STRSM ',INFO)
183:           RETURN
184:       END IF
185: *
186: *     Quick return if possible.
187: *
188:       IF (M.EQ.0 .OR. N.EQ.0) RETURN
189: *
190: *     And when  alpha.eq.zero.
191: *
192:       IF (ALPHA.EQ.ZERO) THEN
193:           DO 20 J = 1,N
194:               DO 10 I = 1,M
195:                   B(I,J) = ZERO
196:    10         CONTINUE
197:    20     CONTINUE
198:           RETURN
199:       END IF
200: *
201: *     Start the operations.
202: *
203:       IF (LSIDE) THEN
204:           IF (LSAME(TRANSA,'N')) THEN
205: *
206: *           Form  B := alpha*inv( A )*B.
207: *
208:               IF (UPPER) THEN
209:                   DO 60 J = 1,N
210:                       IF (ALPHA.NE.ONE) THEN
211:                           DO 30 I = 1,M
212:                               B(I,J) = ALPHA*B(I,J)
213:    30                     CONTINUE
214:                       END IF
215:                       DO 50 K = M,1,-1
216:                           IF (B(K,J).NE.ZERO) THEN
217:                               IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
218:                               DO 40 I = 1,K - 1
219:                                   B(I,J) = B(I,J) - B(K,J)*A(I,K)
220:    40                         CONTINUE
221:                           END IF
222:    50                 CONTINUE
223:    60             CONTINUE
224:               ELSE
225:                   DO 100 J = 1,N
226:                       IF (ALPHA.NE.ONE) THEN
227:                           DO 70 I = 1,M
228:                               B(I,J) = ALPHA*B(I,J)
229:    70                     CONTINUE
230:                       END IF
231:                       DO 90 K = 1,M
232:                           IF (B(K,J).NE.ZERO) THEN
233:                               IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
234:                               DO 80 I = K + 1,M
235:                                   B(I,J) = B(I,J) - B(K,J)*A(I,K)
236:    80                         CONTINUE
237:                           END IF
238:    90                 CONTINUE
239:   100             CONTINUE
240:               END IF
241:           ELSE
242: *
243: *           Form  B := alpha*inv( A' )*B.
244: *
245:               IF (UPPER) THEN
246:                   DO 130 J = 1,N
247:                       DO 120 I = 1,M
248:                           TEMP = ALPHA*B(I,J)
249:                           DO 110 K = 1,I - 1
250:                               TEMP = TEMP - A(K,I)*B(K,J)
251:   110                     CONTINUE
252:                           IF (NOUNIT) TEMP = TEMP/A(I,I)
253:                           B(I,J) = TEMP
254:   120                 CONTINUE
255:   130             CONTINUE
256:               ELSE
257:                   DO 160 J = 1,N
258:                       DO 150 I = M,1,-1
259:                           TEMP = ALPHA*B(I,J)
260:                           DO 140 K = I + 1,M
261:                               TEMP = TEMP - A(K,I)*B(K,J)
262:   140                     CONTINUE
263:                           IF (NOUNIT) TEMP = TEMP/A(I,I)
264:                           B(I,J) = TEMP
265:   150                 CONTINUE
266:   160             CONTINUE
267:               END IF
268:           END IF
269:       ELSE
270:           IF (LSAME(TRANSA,'N')) THEN
271: *
272: *           Form  B := alpha*B*inv( A ).
273: *
274:               IF (UPPER) THEN
275:                   DO 210 J = 1,N
276:                       IF (ALPHA.NE.ONE) THEN
277:                           DO 170 I = 1,M
278:                               B(I,J) = ALPHA*B(I,J)
279:   170                     CONTINUE
280:                       END IF
281:                       DO 190 K = 1,J - 1
282:                           IF (A(K,J).NE.ZERO) THEN
283:                               DO 180 I = 1,M
284:                                   B(I,J) = B(I,J) - A(K,J)*B(I,K)
285:   180                         CONTINUE
286:                           END IF
287:   190                 CONTINUE
288:                       IF (NOUNIT) THEN
289:                           TEMP = ONE/A(J,J)
290:                           DO 200 I = 1,M
291:                               B(I,J) = TEMP*B(I,J)
292:   200                     CONTINUE
293:                       END IF
294:   210             CONTINUE
295:               ELSE
296:                   DO 260 J = N,1,-1
297:                       IF (ALPHA.NE.ONE) THEN
298:                           DO 220 I = 1,M
299:                               B(I,J) = ALPHA*B(I,J)
300:   220                     CONTINUE
301:                       END IF
302:                       DO 240 K = J + 1,N
303:                           IF (A(K,J).NE.ZERO) THEN
304:                               DO 230 I = 1,M
305:                                   B(I,J) = B(I,J) - A(K,J)*B(I,K)
306:   230                         CONTINUE
307:                           END IF
308:   240                 CONTINUE
309:                       IF (NOUNIT) THEN
310:                           TEMP = ONE/A(J,J)
311:                           DO 250 I = 1,M
312:                               B(I,J) = TEMP*B(I,J)
313:   250                     CONTINUE
314:                       END IF
315:   260             CONTINUE
316:               END IF
317:           ELSE
318: *
319: *           Form  B := alpha*B*inv( A' ).
320: *
321:               IF (UPPER) THEN
322:                   DO 310 K = N,1,-1
323:                       IF (NOUNIT) THEN
324:                           TEMP = ONE/A(K,K)
325:                           DO 270 I = 1,M
326:                               B(I,K) = TEMP*B(I,K)
327:   270                     CONTINUE
328:                       END IF
329:                       DO 290 J = 1,K - 1
330:                           IF (A(J,K).NE.ZERO) THEN
331:                               TEMP = A(J,K)
332:                               DO 280 I = 1,M
333:                                   B(I,J) = B(I,J) - TEMP*B(I,K)
334:   280                         CONTINUE
335:                           END IF
336:   290                 CONTINUE
337:                       IF (ALPHA.NE.ONE) THEN
338:                           DO 300 I = 1,M
339:                               B(I,K) = ALPHA*B(I,K)
340:   300                     CONTINUE
341:                       END IF
342:   310             CONTINUE
343:               ELSE
344:                   DO 360 K = 1,N
345:                       IF (NOUNIT) THEN
346:                           TEMP = ONE/A(K,K)
347:                           DO 320 I = 1,M
348:                               B(I,K) = TEMP*B(I,K)
349:   320                     CONTINUE
350:                       END IF
351:                       DO 340 J = K + 1,N
352:                           IF (A(J,K).NE.ZERO) THEN
353:                               TEMP = A(J,K)
354:                               DO 330 I = 1,M
355:                                   B(I,J) = B(I,J) - TEMP*B(I,K)
356:   330                         CONTINUE
357:                           END IF
358:   340                 CONTINUE
359:                       IF (ALPHA.NE.ONE) THEN
360:                           DO 350 I = 1,M
361:                               B(I,K) = ALPHA*B(I,K)
362:   350                     CONTINUE
363:                       END IF
364:   360             CONTINUE
365:               END IF
366:           END IF
367:       END IF
368: *
369:       RETURN
370: *
371: *     End of STRSM .
372: *
373:       END
374: