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