001:       SUBROUTINE STFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
002:      +                  B, LDB )
003: *
004: *  -- LAPACK routine (version 3.2.1)                                    --
005: *
006: *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
007: *  -- April 2009                                                      --
008: *
009: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
010: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
011: *
012: *     ..
013: *     .. Scalar Arguments ..
014:       CHARACTER          TRANSR, DIAG, SIDE, TRANS, UPLO
015:       INTEGER            LDB, M, N
016:       REAL               ALPHA
017: *     ..
018: *     .. Array Arguments ..
019:       REAL               A( 0: * ), B( 0: LDB-1, 0: * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  Level 3 BLAS like routine for A in RFP Format.
026: *
027: *  STFSM  solves the matrix equation
028: *
029: *     op( A )*X = alpha*B  or  X*op( A ) = alpha*B
030: *
031: *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
032: *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
033: *
034: *     op( A ) = A   or   op( A ) = A'.
035: *
036: *  A is in Rectangular Full Packed (RFP) Format.
037: *
038: *  The matrix X is overwritten on B.
039: *
040: *  Arguments
041: *  ==========
042: *
043: *  TRANSR - (input) CHARACTER
044: *          = 'N':  The Normal Form of RFP A is stored;
045: *          = 'T':  The Transpose Form of RFP A is stored.
046: *
047: *  SIDE   - (input) CHARACTER
048: *           On entry, SIDE specifies whether op( A ) appears on the left
049: *           or right of X as follows:
050: *
051: *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
052: *
053: *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
054: *
055: *           Unchanged on exit.
056: *
057: *  UPLO   - (input) CHARACTER
058: *           On entry, UPLO specifies whether the RFP matrix A came from
059: *           an upper or lower triangular matrix as follows:
060: *           UPLO = 'U' or 'u' RFP A came from an upper triangular matrix
061: *           UPLO = 'L' or 'l' RFP A came from a  lower triangular matrix
062: *
063: *           Unchanged on exit.
064: *
065: *  TRANS  - (input) CHARACTER
066: *           On entry, TRANS  specifies the form of op( A ) to be used
067: *           in the matrix multiplication as follows:
068: *
069: *              TRANS  = 'N' or 'n'   op( A ) = A.
070: *
071: *              TRANS  = 'T' or 't'   op( A ) = A'.
072: *
073: *           Unchanged on exit.
074: *
075: *  DIAG   - (input) CHARACTER
076: *           On entry, DIAG specifies whether or not RFP A is unit
077: *           triangular as follows:
078: *
079: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
080: *
081: *              DIAG = 'N' or 'n'   A is not assumed to be unit
082: *                                  triangular.
083: *
084: *           Unchanged on exit.
085: *
086: *  M      - (input) INTEGER.
087: *           On entry, M specifies the number of rows of B. M must be at
088: *           least zero.
089: *           Unchanged on exit.
090: *
091: *  N      - (input) INTEGER.
092: *           On entry, N specifies the number of columns of B.  N must be
093: *           at least zero.
094: *           Unchanged on exit.
095: *
096: *  ALPHA  - (input) REAL.
097: *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
098: *           zero then  A is not referenced and  B need not be set before
099: *           entry.
100: *           Unchanged on exit.
101: *
102: *  A      - (input) REAL array, dimension (NT);
103: *           NT = N*(N+1)/2. On entry, the matrix A in RFP Format.
104: *           RFP Format is described by TRANSR, UPLO and N as follows:
105: *           If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;
106: *           K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If
107: *           TRANSR = 'T' then RFP is the transpose of RFP A as
108: *           defined when TRANSR = 'N'. The contents of RFP A are defined
109: *           by UPLO as follows: If UPLO = 'U' the RFP A contains the NT
110: *           elements of upper packed A either in normal or
111: *           transpose Format. If UPLO = 'L' the RFP A contains
112: *           the NT elements of lower packed A either in normal or
113: *           transpose Format. The LDA of RFP A is (N+1)/2 when
114: *           TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
115: *           even and is N when is odd.
116: *           See the Note below for more details. Unchanged on exit.
117: *
118: *  B      - (input/ouptut) REAL array,  DIMENSION (LDB,N)
119: *           Before entry,  the leading  m by n part of the array  B must
120: *           contain  the  right-hand  side  matrix  B,  and  on exit  is
121: *           overwritten by the solution matrix  X.
122: *
123: *  LDB    - (input) INTEGER.
124: *           On entry, LDB specifies the first dimension of B as declared
125: *           in  the  calling  (sub)  program.   LDB  must  be  at  least
126: *           max( 1, m ).
127: *           Unchanged on exit.
128: *
129: *  Further Details
130: *  ===============
131: *
132: *  We first consider Rectangular Full Packed (RFP) Format when N is
133: *  even. We give an example where N = 6.
134: *
135: *      AP is Upper             AP is Lower
136: *
137: *   00 01 02 03 04 05       00
138: *      11 12 13 14 15       10 11
139: *         22 23 24 25       20 21 22
140: *            33 34 35       30 31 32 33
141: *               44 45       40 41 42 43 44
142: *                  55       50 51 52 53 54 55
143: *
144: *
145: *  Let TRANSR = 'N'. RFP holds AP as follows:
146: *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
147: *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
148: *  the transpose of the first three columns of AP upper.
149: *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
150: *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
151: *  the transpose of the last three columns of AP lower.
152: *  This covers the case N even and TRANSR = 'N'.
153: *
154: *         RFP A                   RFP A
155: *
156: *        03 04 05                33 43 53
157: *        13 14 15                00 44 54
158: *        23 24 25                10 11 55
159: *        33 34 35                20 21 22
160: *        00 44 45                30 31 32
161: *        01 11 55                40 41 42
162: *        02 12 22                50 51 52
163: *
164: *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
165: *  transpose of RFP A above. One therefore gets:
166: *
167: *
168: *           RFP A                   RFP A
169: *
170: *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
171: *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
172: *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
173: *
174: *
175: *  We first consider Rectangular Full Packed (RFP) Format when N is
176: *  odd. We give an example where N = 5.
177: *
178: *     AP is Upper                 AP is Lower
179: *
180: *   00 01 02 03 04              00
181: *      11 12 13 14              10 11
182: *         22 23 24              20 21 22
183: *            33 34              30 31 32 33
184: *               44              40 41 42 43 44
185: *
186: *
187: *  Let TRANSR = 'N'. RFP holds AP as follows:
188: *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
189: *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
190: *  the transpose of the first two columns of AP upper.
191: *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
192: *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
193: *  the transpose of the last two columns of AP lower.
194: *  This covers the case N odd and TRANSR = 'N'.
195: *
196: *         RFP A                   RFP A
197: *
198: *        02 03 04                00 33 43
199: *        12 13 14                10 11 44
200: *        22 23 24                20 21 22
201: *        00 33 34                30 31 32
202: *        01 11 44                40 41 42
203: *
204: *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
205: *  transpose of RFP A above. One therefore gets:
206: *
207: *           RFP A                   RFP A
208: *
209: *     02 12 22 00 01             00 10 20 30 40 50
210: *     03 13 23 33 11             33 11 21 31 41 51
211: *     04 14 24 34 44             43 44 22 32 42 52
212: *
213: *  Reference
214: *  =========
215: *
216: *  =====================================================================
217: *
218: *     ..
219: *     .. Parameters ..
220:       REAL               ONE, ZERO
221:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
222: *     ..
223: *     .. Local Scalars ..
224:       LOGICAL            LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR,
225:      +                   NOTRANS
226:       INTEGER            M1, M2, N1, N2, K, INFO, I, J
227: *     ..
228: *     .. External Functions ..
229:       LOGICAL            LSAME
230:       EXTERNAL           LSAME
231: *     ..
232: *     .. External Subroutines ..
233:       EXTERNAL           SGEMM, STRSM, XERBLA
234: *     ..
235: *     .. Intrinsic Functions ..
236:       INTRINSIC          MAX, MOD
237: *     ..
238: *     .. Executable Statements ..
239: *
240: *     Test the input parameters.
241: *
242:       INFO = 0
243:       NORMALTRANSR = LSAME( TRANSR, 'N' )
244:       LSIDE = LSAME( SIDE, 'L' )
245:       LOWER = LSAME( UPLO, 'L' )
246:       NOTRANS = LSAME( TRANS, 'N' )
247:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
248:          INFO = -1
249:       ELSE IF( .NOT.LSIDE .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
250:          INFO = -2
251:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
252:          INFO = -3
253:       ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
254:          INFO = -4
255:       ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) )
256:      +         THEN
257:          INFO = -5
258:       ELSE IF( M.LT.0 ) THEN
259:          INFO = -6
260:       ELSE IF( N.LT.0 ) THEN
261:          INFO = -7
262:       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
263:          INFO = -11
264:       END IF
265:       IF( INFO.NE.0 ) THEN
266:          CALL XERBLA( 'STFSM ', -INFO )
267:          RETURN
268:       END IF
269: *
270: *     Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
271: *
272:       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
273:      +   RETURN
274: *
275: *     Quick return when ALPHA.EQ.(0D+0)
276: *
277:       IF( ALPHA.EQ.ZERO ) THEN
278:          DO 20 J = 0, N - 1
279:             DO 10 I = 0, M - 1
280:                B( I, J ) = ZERO
281:    10       CONTINUE
282:    20    CONTINUE
283:          RETURN
284:       END IF
285: *
286:       IF( LSIDE ) THEN
287: *
288: *        SIDE = 'L'
289: *
290: *        A is M-by-M.
291: *        If M is odd, set NISODD = .TRUE., and M1 and M2.
292: *        If M is even, NISODD = .FALSE., and M.
293: *
294:          IF( MOD( M, 2 ).EQ.0 ) THEN
295:             MISODD = .FALSE.
296:             K = M / 2
297:          ELSE
298:             MISODD = .TRUE.
299:             IF( LOWER ) THEN
300:                M2 = M / 2
301:                M1 = M - M2
302:             ELSE
303:                M1 = M / 2
304:                M2 = M - M1
305:             END IF
306:          END IF
307: *
308:          IF( MISODD ) THEN
309: *
310: *           SIDE = 'L' and N is odd
311: *
312:             IF( NORMALTRANSR ) THEN
313: *
314: *              SIDE = 'L', N is odd, and TRANSR = 'N'
315: *
316:                IF( LOWER ) THEN
317: *
318: *                 SIDE  ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
319: *
320:                   IF( NOTRANS ) THEN
321: *
322: *                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
323: *                    TRANS = 'N'
324: *
325:                      IF( M.EQ.1 ) THEN
326:                         CALL STRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
327:      +                              A, M, B, LDB )
328:                      ELSE
329:                         CALL STRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
330:      +                              A( 0 ), M, B, LDB )
331:                         CALL SGEMM( 'N', 'N', M2, N, M1, -ONE, A( M1 ),
332:      +                              M, B, LDB, ALPHA, B( M1, 0 ), LDB )
333:                         CALL STRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
334:      +                              A( M ), M, B( M1, 0 ), LDB )
335:                      END IF
336: *
337:                   ELSE
338: *
339: *                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
340: *                    TRANS = 'T'
341: *
342:                      IF( M.EQ.1 ) THEN
343:                         CALL STRSM( 'L', 'L', 'T', DIAG, M1, N, ALPHA,
344:      +                              A( 0 ), M, B, LDB )
345:                      ELSE
346:                         CALL STRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
347:      +                              A( M ), M, B( M1, 0 ), LDB )
348:                         CALL SGEMM( 'T', 'N', M1, N, M2, -ONE, A( M1 ), M,
349:      +                              M, B( M1, 0 ), LDB, ALPHA, B, LDB )
350:                         CALL STRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
351:      +                              A( 0 ), M, B, LDB )
352:                      END IF
353: *
354:                   END IF
355: *
356:                ELSE
357: *
358: *                 SIDE  ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
359: *
360:                   IF( .NOT.NOTRANS ) THEN
361: *
362: *                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
363: *                    TRANS = 'N'
364: *
365:                      CALL STRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
366:      +                           A( M2 ), M, B, LDB )
367:                      CALL SGEMM( 'T', 'N', M2, N, M1, -ONE, A( 0 ), M,
368:      +                           B, LDB, ALPHA, B( M1, 0 ), LDB )
369:                      CALL STRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
370:      +                           A( M1 ), M, B( M1, 0 ), LDB )
371: *
372:                   ELSE
373: *
374: *                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
375: *                    TRANS = 'T'
376: *
377:                      CALL STRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
378:      +                           A( M1 ), M, B( M1, 0 ), LDB )
379:                      CALL SGEMM( 'N', 'N', M1, N, M2, -ONE, A( 0 ), M,
380:      +                           B( M1, 0 ), LDB, ALPHA, B, LDB )
381:                      CALL STRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
382:      +                           A( M2 ), M, B, LDB )
383: *
384:                   END IF
385: *
386:                END IF
387: *
388:             ELSE
389: *
390: *              SIDE = 'L', N is odd, and TRANSR = 'T'
391: *
392:                IF( LOWER ) THEN
393: *
394: *                 SIDE  ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
395: *
396:                   IF( NOTRANS ) THEN
397: *
398: *                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
399: *                    TRANS = 'N'
400: *
401:                      IF( M.EQ.1 ) THEN
402:                         CALL STRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
403:      +                              A( 0 ), M1, B, LDB )
404:                      ELSE
405:                         CALL STRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
406:      +                              A( 0 ), M1, B, LDB )
407:                         CALL SGEMM( 'T', 'N', M2, N, M1, -ONE,
408:      +                              A( M1*M1 ), M1, B, LDB, ALPHA,
409:      +                              B( M1, 0 ), LDB )
410:                         CALL STRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
411:      +                              A( 1 ), M1, B( M1, 0 ), LDB )
412:                      END IF
413: *
414:                   ELSE
415: *
416: *                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
417: *                    TRANS = 'T'
418: *
419:                      IF( M.EQ.1 ) THEN
420:                         CALL STRSM( 'L', 'U', 'N', DIAG, M1, N, ALPHA,
421:      +                              A( 0 ), M1, B, LDB )
422:                      ELSE
423:                         CALL STRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
424:      +                              A( 1 ), M1, B( M1, 0 ), LDB )
425:                         CALL SGEMM( 'N', 'N', M1, N, M2, -ONE,
426:      +                              A( M1*M1 ), M1, B( M1, 0 ), LDB,
427:      +                              ALPHA, B, LDB )
428:                         CALL STRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
429:      +                              A( 0 ), M1, B, LDB )
430:                      END IF
431: *
432:                   END IF
433: *
434:                ELSE
435: *
436: *                 SIDE  ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
437: *
438:                   IF( .NOT.NOTRANS ) THEN
439: *
440: *                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
441: *                    TRANS = 'N'
442: *
443:                      CALL STRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
444:      +                           A( M2*M2 ), M2, B, LDB )
445:                      CALL SGEMM( 'N', 'N', M2, N, M1, -ONE, A( 0 ), M2,
446:      +                           B, LDB, ALPHA, B( M1, 0 ), LDB )
447:                      CALL STRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
448:      +                           A( M1*M2 ), M2, B( M1, 0 ), LDB )
449: *
450:                   ELSE
451: *
452: *                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
453: *                    TRANS = 'T'
454: *
455:                      CALL STRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
456:      +                           A( M1*M2 ), M2, B( M1, 0 ), LDB )
457:                      CALL SGEMM( 'T', 'N', M1, N, M2, -ONE, A( 0 ), M2,
458:      +                           B( M1, 0 ), LDB, ALPHA, B, LDB )
459:                      CALL STRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
460:      +                           A( M2*M2 ), M2, B, LDB )
461: *
462:                   END IF
463: *
464:                END IF
465: *
466:             END IF
467: *
468:          ELSE
469: *
470: *           SIDE = 'L' and N is even
471: *
472:             IF( NORMALTRANSR ) THEN
473: *
474: *              SIDE = 'L', N is even, and TRANSR = 'N'
475: *
476:                IF( LOWER ) THEN
477: *
478: *                 SIDE  ='L', N is even, TRANSR = 'N', and UPLO = 'L'
479: *
480:                   IF( NOTRANS ) THEN
481: *
482: *                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'L',
483: *                    and TRANS = 'N'
484: *
485:                      CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
486:      +                           A( 1 ), M+1, B, LDB )
487:                      CALL SGEMM( 'N', 'N', K, N, K, -ONE, A( K+1 ),
488:      +                           M+1, B, LDB, ALPHA, B( K, 0 ), LDB )
489:                      CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
490:      +                           A( 0 ), M+1, B( K, 0 ), LDB )
491: *
492:                   ELSE
493: *
494: *                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'L',
495: *                    and TRANS = 'T'
496: *
497:                      CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
498:      +                           A( 0 ), M+1, B( K, 0 ), LDB )
499:                      CALL SGEMM( 'T', 'N', K, N, K, -ONE, A( K+1 ),
500:      +                           M+1, B( K, 0 ), LDB, ALPHA, B, LDB )
501:                      CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
502:      +                           A( 1 ), M+1, B, LDB )
503: *
504:                   END IF
505: *
506:                ELSE
507: *
508: *                 SIDE  ='L', N is even, TRANSR = 'N', and UPLO = 'U'
509: *
510:                   IF( .NOT.NOTRANS ) THEN
511: *
512: *                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'U',
513: *                    and TRANS = 'N'
514: *
515:                      CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
516:      +                           A( K+1 ), M+1, B, LDB )
517:                      CALL SGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), M+1,
518:      +                           B, LDB, ALPHA, B( K, 0 ), LDB )
519:                      CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
520:      +                           A( K ), M+1, B( K, 0 ), LDB )
521: *
522:                   ELSE
523: *
524: *                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'U',
525: *                    and TRANS = 'T'
526:                      CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
527:      +                           A( K ), M+1, B( K, 0 ), LDB )
528:                      CALL SGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), M+1,
529:      +                           B( K, 0 ), LDB, ALPHA, B, LDB )
530:                      CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
531:      +                           A( K+1 ), M+1, B, LDB )
532: *
533:                   END IF
534: *
535:                END IF
536: *
537:             ELSE
538: *
539: *              SIDE = 'L', N is even, and TRANSR = 'T'
540: *
541:                IF( LOWER ) THEN
542: *
543: *                 SIDE  ='L', N is even, TRANSR = 'T', and UPLO = 'L'
544: *
545:                   IF( NOTRANS ) THEN
546: *
547: *                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'L',
548: *                    and TRANS = 'N'
549: *
550:                      CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
551:      +                           A( K ), K, B, LDB )
552:                      CALL SGEMM( 'T', 'N', K, N, K, -ONE,
553:      +                           A( K*( K+1 ) ), K, B, LDB, ALPHA,
554:      +                           B( K, 0 ), LDB )
555:                      CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
556:      +                           A( 0 ), K, B( K, 0 ), LDB )
557: *
558:                   ELSE
559: *
560: *                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'L',
561: *                    and TRANS = 'T'
562: *
563:                      CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
564:      +                           A( 0 ), K, B( K, 0 ), LDB )
565:                      CALL SGEMM( 'N', 'N', K, N, K, -ONE,
566:      +                           A( K*( K+1 ) ), K, B( K, 0 ), LDB,
567:      +                           ALPHA, B, LDB )
568:                      CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
569:      +                           A( K ), K, B, LDB )
570: *
571:                   END IF
572: *
573:                ELSE
574: *
575: *                 SIDE  ='L', N is even, TRANSR = 'T', and UPLO = 'U'
576: *
577:                   IF( .NOT.NOTRANS ) THEN
578: *
579: *                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'U',
580: *                    and TRANS = 'N'
581: *
582:                      CALL STRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
583:      +                           A( K*( K+1 ) ), K, B, LDB )
584:                      CALL SGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), K, B,
585:      +                           LDB, ALPHA, B( K, 0 ), LDB )
586:                      CALL STRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
587:      +                           A( K*K ), K, B( K, 0 ), LDB )
588: *
589:                   ELSE
590: *
591: *                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'U',
592: *                    and TRANS = 'T'
593: *
594:                      CALL STRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
595:      +                           A( K*K ), K, B( K, 0 ), LDB )
596:                      CALL SGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), K,
597:      +                           B( K, 0 ), LDB, ALPHA, B, LDB )
598:                      CALL STRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
599:      +                           A( K*( K+1 ) ), K, B, LDB )
600: *
601:                   END IF
602: *
603:                END IF
604: *
605:             END IF
606: *
607:          END IF
608: *
609:       ELSE
610: *
611: *        SIDE = 'R'
612: *
613: *        A is N-by-N.
614: *        If N is odd, set NISODD = .TRUE., and N1 and N2.
615: *        If N is even, NISODD = .FALSE., and K.
616: *
617:          IF( MOD( N, 2 ).EQ.0 ) THEN
618:             NISODD = .FALSE.
619:             K = N / 2
620:          ELSE
621:             NISODD = .TRUE.
622:             IF( LOWER ) THEN
623:                N2 = N / 2
624:                N1 = N - N2
625:             ELSE
626:                N1 = N / 2
627:                N2 = N - N1
628:             END IF
629:          END IF
630: *
631:          IF( NISODD ) THEN
632: *
633: *           SIDE = 'R' and N is odd
634: *
635:             IF( NORMALTRANSR ) THEN
636: *
637: *              SIDE = 'R', N is odd, and TRANSR = 'N'
638: *
639:                IF( LOWER ) THEN
640: *
641: *                 SIDE  ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
642: *
643:                   IF( NOTRANS ) THEN
644: *
645: *                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
646: *                    TRANS = 'N'
647: *
648:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
649:      +                           A( N ), N, B( 0, N1 ), LDB )
650:                      CALL SGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
651:      +                           LDB, A( N1 ), N, ALPHA, B( 0, 0 ),
652:      +                           LDB )
653:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
654:      +                           A( 0 ), N, B( 0, 0 ), LDB )
655: *
656:                   ELSE
657: *
658: *                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
659: *                    TRANS = 'T'
660: *
661:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
662:      +                           A( 0 ), N, B( 0, 0 ), LDB )
663:                      CALL SGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
664:      +                           LDB, A( N1 ), N, ALPHA, B( 0, N1 ),
665:      +                           LDB )
666:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
667:      +                           A( N ), N, B( 0, N1 ), LDB )
668: *
669:                   END IF
670: *
671:                ELSE
672: *
673: *                 SIDE  ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
674: *
675:                   IF( NOTRANS ) THEN
676: *
677: *                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
678: *                    TRANS = 'N'
679: *
680:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
681:      +                           A( N2 ), N, B( 0, 0 ), LDB )
682:                      CALL SGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
683:      +                           LDB, A( 0 ), N, ALPHA, B( 0, N1 ),
684:      +                           LDB )
685:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
686:      +                           A( N1 ), N, B( 0, N1 ), LDB )
687: *
688:                   ELSE
689: *
690: *                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
691: *                    TRANS = 'T'
692: *
693:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
694:      +                           A( N1 ), N, B( 0, N1 ), LDB )
695:                      CALL SGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
696:      +                           LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB )
697:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
698:      +                           A( N2 ), N, B( 0, 0 ), LDB )
699: *
700:                   END IF
701: *
702:                END IF
703: *
704:             ELSE
705: *
706: *              SIDE = 'R', N is odd, and TRANSR = 'T'
707: *
708:                IF( LOWER ) THEN
709: *
710: *                 SIDE  ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
711: *
712:                   IF( NOTRANS ) THEN
713: *
714: *                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
715: *                    TRANS = 'N'
716: *
717:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
718:      +                           A( 1 ), N1, B( 0, N1 ), LDB )
719:                      CALL SGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
720:      +                           LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ),
721:      +                           LDB )
722:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
723:      +                           A( 0 ), N1, B( 0, 0 ), LDB )
724: *
725:                   ELSE
726: *
727: *                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
728: *                    TRANS = 'T'
729: *
730:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
731:      +                           A( 0 ), N1, B( 0, 0 ), LDB )
732:                      CALL SGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
733:      +                           LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ),
734:      +                           LDB )
735:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
736:      +                           A( 1 ), N1, B( 0, N1 ), LDB )
737: *
738:                   END IF
739: *
740:                ELSE
741: *
742: *                 SIDE  ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
743: *
744:                   IF( NOTRANS ) THEN
745: *
746: *                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
747: *                    TRANS = 'N'
748: *
749:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
750:      +                           A( N2*N2 ), N2, B( 0, 0 ), LDB )
751:                      CALL SGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
752:      +                           LDB, A( 0 ), N2, ALPHA, B( 0, N1 ),
753:      +                           LDB )
754:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
755:      +                           A( N1*N2 ), N2, B( 0, N1 ), LDB )
756: *
757:                   ELSE
758: *
759: *                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
760: *                    TRANS = 'T'
761: *
762:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
763:      +                           A( N1*N2 ), N2, B( 0, N1 ), LDB )
764:                      CALL SGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
765:      +                           LDB, A( 0 ), N2, ALPHA, B( 0, 0 ),
766:      +                           LDB )
767:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
768:      +                           A( N2*N2 ), N2, B( 0, 0 ), LDB )
769: *
770:                   END IF
771: *
772:                END IF
773: *
774:             END IF
775: *
776:          ELSE
777: *
778: *           SIDE = 'R' and N is even
779: *
780:             IF( NORMALTRANSR ) THEN
781: *
782: *              SIDE = 'R', N is even, and TRANSR = 'N'
783: *
784:                IF( LOWER ) THEN
785: *
786: *                 SIDE  ='R', N is even, TRANSR = 'N', and UPLO = 'L'
787: *
788:                   IF( NOTRANS ) THEN
789: *
790: *                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'L',
791: *                    and TRANS = 'N'
792: *
793:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
794:      +                           A( 0 ), N+1, B( 0, K ), LDB )
795:                      CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
796:      +                           LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ),
797:      +                           LDB )
798:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
799:      +                           A( 1 ), N+1, B( 0, 0 ), LDB )
800: *
801:                   ELSE
802: *
803: *                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'L',
804: *                    and TRANS = 'T'
805: *
806:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
807:      +                           A( 1 ), N+1, B( 0, 0 ), LDB )
808:                      CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
809:      +                           LDB, A( K+1 ), N+1, ALPHA, B( 0, K ),
810:      +                           LDB )
811:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
812:      +                           A( 0 ), N+1, B( 0, K ), LDB )
813: *
814:                   END IF
815: *
816:                ELSE
817: *
818: *                 SIDE  ='R', N is even, TRANSR = 'N', and UPLO = 'U'
819: *
820:                   IF( NOTRANS ) THEN
821: *
822: *                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'U',
823: *                    and TRANS = 'N'
824: *
825:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
826:      +                           A( K+1 ), N+1, B( 0, 0 ), LDB )
827:                      CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
828:      +                           LDB, A( 0 ), N+1, ALPHA, B( 0, K ),
829:      +                           LDB )
830:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
831:      +                           A( K ), N+1, B( 0, K ), LDB )
832: *
833:                   ELSE
834: *
835: *                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'U',
836: *                    and TRANS = 'T'
837: *
838:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
839:      +                           A( K ), N+1, B( 0, K ), LDB )
840:                      CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
841:      +                           LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ),
842:      +                           LDB )
843:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
844:      +                           A( K+1 ), N+1, B( 0, 0 ), LDB )
845: *
846:                   END IF
847: *
848:                END IF
849: *
850:             ELSE
851: *
852: *              SIDE = 'R', N is even, and TRANSR = 'T'
853: *
854:                IF( LOWER ) THEN
855: *
856: *                 SIDE  ='R', N is even, TRANSR = 'T', and UPLO = 'L'
857: *
858:                   IF( NOTRANS ) THEN
859: *
860: *                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'L',
861: *                    and TRANS = 'N'
862: *
863:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
864:      +                           A( 0 ), K, B( 0, K ), LDB )
865:                      CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
866:      +                           LDB, A( ( K+1 )*K ), K, ALPHA,
867:      +                           B( 0, 0 ), LDB )
868:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
869:      +                           A( K ), K, B( 0, 0 ), LDB )
870: *
871:                   ELSE
872: *
873: *                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'L',
874: *                    and TRANS = 'T'
875: *
876:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
877:      +                           A( K ), K, B( 0, 0 ), LDB )
878:                      CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
879:      +                           LDB, A( ( K+1 )*K ), K, ALPHA,
880:      +                           B( 0, K ), LDB )
881:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
882:      +                           A( 0 ), K, B( 0, K ), LDB )
883: *
884:                   END IF
885: *
886:                ELSE
887: *
888: *                 SIDE  ='R', N is even, TRANSR = 'T', and UPLO = 'U'
889: *
890:                   IF( NOTRANS ) THEN
891: *
892: *                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'U',
893: *                    and TRANS = 'N'
894: *
895:                      CALL STRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
896:      +                           A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
897:                      CALL SGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
898:      +                           LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB )
899:                      CALL STRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
900:      +                           A( K*K ), K, B( 0, K ), LDB )
901: *
902:                   ELSE
903: *
904: *                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'U',
905: *                    and TRANS = 'T'
906: *
907:                      CALL STRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
908:      +                           A( K*K ), K, B( 0, K ), LDB )
909:                      CALL SGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
910:      +                           LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB )
911:                      CALL STRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
912:      +                           A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
913: *
914:                   END IF
915: *
916:                END IF
917: *
918:             END IF
919: *
920:          END IF
921:       END IF
922: *
923:       RETURN
924: *
925: *     End of STFSM
926: *
927:       END
928: