LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
stfsm.f
Go to the documentation of this file.
1*> \brief \b STFSM solves a matrix equation (one operand is a triangular matrix in RFP format).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download STFSM + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stfsm.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stfsm.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stfsm.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE STFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
20* B, LDB )
21*
22* .. Scalar Arguments ..
23* CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
24* INTEGER LDB, M, N
25* REAL ALPHA
26* ..
27* .. Array Arguments ..
28* REAL A( 0: * ), B( 0: LDB-1, 0: * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> Level 3 BLAS like routine for A in RFP Format.
38*>
39*> STFSM solves the matrix equation
40*>
41*> op( A )*X = alpha*B or X*op( A ) = alpha*B
42*>
43*> where alpha is a scalar, X and B are m by n matrices, A is a unit, or
44*> non-unit, upper or lower triangular matrix and op( A ) is one of
45*>
46*> op( A ) = A or op( A ) = A**T.
47*>
48*> A is in Rectangular Full Packed (RFP) Format.
49*>
50*> The matrix X is overwritten on B.
51*> \endverbatim
52*
53* Arguments:
54* ==========
55*
56*> \param[in] TRANSR
57*> \verbatim
58*> TRANSR is CHARACTER*1
59*> = 'N': The Normal Form of RFP A is stored;
60*> = 'T': The Transpose Form of RFP A is stored.
61*> \endverbatim
62*>
63*> \param[in] SIDE
64*> \verbatim
65*> SIDE is CHARACTER*1
66*> On entry, SIDE specifies whether op( A ) appears on the left
67*> or right of X as follows:
68*>
69*> SIDE = 'L' or 'l' op( A )*X = alpha*B.
70*>
71*> SIDE = 'R' or 'r' X*op( A ) = alpha*B.
72*>
73*> Unchanged on exit.
74*> \endverbatim
75*>
76*> \param[in] UPLO
77*> \verbatim
78*> UPLO is CHARACTER*1
79*> On entry, UPLO specifies whether the RFP matrix A came from
80*> an upper or lower triangular matrix as follows:
81*> UPLO = 'U' or 'u' RFP A came from an upper triangular matrix
82*> UPLO = 'L' or 'l' RFP A came from a lower triangular matrix
83*>
84*> Unchanged on exit.
85*> \endverbatim
86*>
87*> \param[in] TRANS
88*> \verbatim
89*> TRANS is CHARACTER*1
90*> On entry, TRANS specifies the form of op( A ) to be used
91*> in the matrix multiplication as follows:
92*>
93*> TRANS = 'N' or 'n' op( A ) = A.
94*>
95*> TRANS = 'T' or 't' op( A ) = A'.
96*>
97*> Unchanged on exit.
98*> \endverbatim
99*>
100*> \param[in] DIAG
101*> \verbatim
102*> DIAG is CHARACTER*1
103*> On entry, DIAG specifies whether or not RFP A is unit
104*> triangular as follows:
105*>
106*> DIAG = 'U' or 'u' A is assumed to be unit triangular.
107*>
108*> DIAG = 'N' or 'n' A is not assumed to be unit
109*> triangular.
110*>
111*> Unchanged on exit.
112*> \endverbatim
113*>
114*> \param[in] M
115*> \verbatim
116*> M is INTEGER
117*> On entry, M specifies the number of rows of B. M must be at
118*> least zero.
119*> Unchanged on exit.
120*> \endverbatim
121*>
122*> \param[in] N
123*> \verbatim
124*> N is INTEGER
125*> On entry, N specifies the number of columns of B. N must be
126*> at least zero.
127*> Unchanged on exit.
128*> \endverbatim
129*>
130*> \param[in] ALPHA
131*> \verbatim
132*> ALPHA is REAL
133*> On entry, ALPHA specifies the scalar alpha. When alpha is
134*> zero then A is not referenced and B need not be set before
135*> entry.
136*> Unchanged on exit.
137*> \endverbatim
138*>
139*> \param[in] A
140*> \verbatim
141*> A is REAL array, dimension (NT)
142*> NT = N*(N+1)/2 if SIDE='R' and NT = M*(M+1)/2 otherwise.
143*> On entry, the matrix A in RFP Format.
144*> RFP Format is described by TRANSR, UPLO and N as follows:
145*> If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;
146*> K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If
147*> TRANSR = 'T' then RFP is the transpose of RFP A as
148*> defined when TRANSR = 'N'. The contents of RFP A are defined
149*> by UPLO as follows: If UPLO = 'U' the RFP A contains the NT
150*> elements of upper packed A either in normal or
151*> transpose Format. If UPLO = 'L' the RFP A contains
152*> the NT elements of lower packed A either in normal or
153*> transpose Format. The LDA of RFP A is (N+1)/2 when
154*> TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
155*> even and is N when is odd.
156*> See the Note below for more details. Unchanged on exit.
157*> \endverbatim
158*>
159*> \param[in,out] B
160*> \verbatim
161*> B is REAL array, dimension (LDB,N)
162*> Before entry, the leading m by n part of the array B must
163*> contain the right-hand side matrix B, and on exit is
164*> overwritten by the solution matrix X.
165*> \endverbatim
166*>
167*> \param[in] LDB
168*> \verbatim
169*> LDB is INTEGER
170*> On entry, LDB specifies the first dimension of B as declared
171*> in the calling (sub) program. LDB must be at least
172*> max( 1, m ).
173*> Unchanged on exit.
174*> \endverbatim
175*
176* Authors:
177* ========
178*
179*> \author Univ. of Tennessee
180*> \author Univ. of California Berkeley
181*> \author Univ. of Colorado Denver
182*> \author NAG Ltd.
183*
184*> \ingroup tfsm
185*
186*> \par Further Details:
187* =====================
188*>
189*> \verbatim
190*>
191*> We first consider Rectangular Full Packed (RFP) Format when N is
192*> even. We give an example where N = 6.
193*>
194*> AP is Upper AP is Lower
195*>
196*> 00 01 02 03 04 05 00
197*> 11 12 13 14 15 10 11
198*> 22 23 24 25 20 21 22
199*> 33 34 35 30 31 32 33
200*> 44 45 40 41 42 43 44
201*> 55 50 51 52 53 54 55
202*>
203*>
204*> Let TRANSR = 'N'. RFP holds AP as follows:
205*> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
206*> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
207*> the transpose of the first three columns of AP upper.
208*> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
209*> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
210*> the transpose of the last three columns of AP lower.
211*> This covers the case N even and TRANSR = 'N'.
212*>
213*> RFP A RFP A
214*>
215*> 03 04 05 33 43 53
216*> 13 14 15 00 44 54
217*> 23 24 25 10 11 55
218*> 33 34 35 20 21 22
219*> 00 44 45 30 31 32
220*> 01 11 55 40 41 42
221*> 02 12 22 50 51 52
222*>
223*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
224*> transpose of RFP A above. One therefore gets:
225*>
226*>
227*> RFP A RFP A
228*>
229*> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
230*> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
231*> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
232*>
233*>
234*> We then consider Rectangular Full Packed (RFP) Format when N is
235*> odd. We give an example where N = 5.
236*>
237*> AP is Upper AP is Lower
238*>
239*> 00 01 02 03 04 00
240*> 11 12 13 14 10 11
241*> 22 23 24 20 21 22
242*> 33 34 30 31 32 33
243*> 44 40 41 42 43 44
244*>
245*>
246*> Let TRANSR = 'N'. RFP holds AP as follows:
247*> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
248*> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
249*> the transpose of the first two columns of AP upper.
250*> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
251*> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
252*> the transpose of the last two columns of AP lower.
253*> This covers the case N odd and TRANSR = 'N'.
254*>
255*> RFP A RFP A
256*>
257*> 02 03 04 00 33 43
258*> 12 13 14 10 11 44
259*> 22 23 24 20 21 22
260*> 00 33 34 30 31 32
261*> 01 11 44 40 41 42
262*>
263*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
264*> transpose of RFP A above. One therefore gets:
265*>
266*> RFP A RFP A
267*>
268*> 02 12 22 00 01 00 10 20 30 40 50
269*> 03 13 23 33 11 33 11 21 31 41 51
270*> 04 14 24 34 44 43 44 22 32 42 52
271*> \endverbatim
272*
273* =====================================================================
274 SUBROUTINE stfsm( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA,
275 $ A,
276 $ B, LDB )
277*
278* -- LAPACK computational routine --
279* -- LAPACK is a software package provided by Univ. of Tennessee, --
280* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
281*
282* .. Scalar Arguments ..
283 CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
284 INTEGER LDB, M, N
285 REAL ALPHA
286* ..
287* .. Array Arguments ..
288 REAL A( 0: * ), B( 0: LDB-1, 0: * )
289* ..
290*
291* =====================================================================
292*
293* ..
294* .. Parameters ..
295 REAL ONE, ZERO
296 PARAMETER ( ONE = 1.0e+0, zero = 0.0e+0 )
297* ..
298* .. Local Scalars ..
299 LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR,
300 $ NOTRANS
301 INTEGER M1, M2, N1, N2, K, INFO, I, J
302* ..
303* .. External Functions ..
304 LOGICAL LSAME
305 EXTERNAL LSAME
306* ..
307* .. External Subroutines ..
308 EXTERNAL sgemm, strsm, xerbla
309* ..
310* .. Intrinsic Functions ..
311 INTRINSIC max, mod
312* ..
313* .. Executable Statements ..
314*
315* Test the input parameters.
316*
317 info = 0
318 normaltransr = lsame( transr, 'N' )
319 lside = lsame( side, 'L' )
320 lower = lsame( uplo, 'L' )
321 notrans = lsame( trans, 'N' )
322 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
323 info = -1
324 ELSE IF( .NOT.lside .AND. .NOT.lsame( side, 'R' ) ) THEN
325 info = -2
326 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
327 info = -3
328 ELSE IF( .NOT.notrans .AND. .NOT.lsame( trans, 'T' ) ) THEN
329 info = -4
330 ELSE IF( .NOT.lsame( diag, 'N' ) .AND.
331 $ .NOT.lsame( diag, 'U' ) )
332 $ THEN
333 info = -5
334 ELSE IF( m.LT.0 ) THEN
335 info = -6
336 ELSE IF( n.LT.0 ) THEN
337 info = -7
338 ELSE IF( ldb.LT.max( 1, m ) ) THEN
339 info = -11
340 END IF
341 IF( info.NE.0 ) THEN
342 CALL xerbla( 'STFSM ', -info )
343 RETURN
344 END IF
345*
346* Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
347*
348 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )
349 $ RETURN
350*
351* Quick return when ALPHA.EQ.(0D+0)
352*
353 IF( alpha.EQ.zero ) THEN
354 DO 20 j = 0, n - 1
355 DO 10 i = 0, m - 1
356 b( i, j ) = zero
357 10 CONTINUE
358 20 CONTINUE
359 RETURN
360 END IF
361*
362 IF( lside ) THEN
363*
364* SIDE = 'L'
365*
366* A is M-by-M.
367* If M is odd, set NISODD = .TRUE., and M1 and M2.
368* If M is even, NISODD = .FALSE., and M.
369*
370 IF( mod( m, 2 ).EQ.0 ) THEN
371 misodd = .false.
372 k = m / 2
373 ELSE
374 misodd = .true.
375 IF( lower ) THEN
376 m2 = m / 2
377 m1 = m - m2
378 ELSE
379 m1 = m / 2
380 m2 = m - m1
381 END IF
382 END IF
383*
384 IF( misodd ) THEN
385*
386* SIDE = 'L' and N is odd
387*
388 IF( normaltransr ) THEN
389*
390* SIDE = 'L', N is odd, and TRANSR = 'N'
391*
392 IF( lower ) THEN
393*
394* SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
395*
396 IF( notrans ) THEN
397*
398* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
399* TRANS = 'N'
400*
401 IF( m.EQ.1 ) THEN
402 CALL strsm( 'L', 'L', 'N', diag, m1, n,
403 $ alpha,
404 $ a, m, b, ldb )
405 ELSE
406 CALL strsm( 'L', 'L', 'N', diag, m1, n,
407 $ alpha,
408 $ a( 0 ), m, b, ldb )
409 CALL sgemm( 'N', 'N', m2, n, m1, -one,
410 $ a( m1 ),
411 $ m, b, ldb, alpha, b( m1, 0 ), ldb )
412 CALL strsm( 'L', 'U', 'T', diag, m2, n, one,
413 $ a( m ), m, b( m1, 0 ), ldb )
414 END IF
415*
416 ELSE
417*
418* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
419* TRANS = 'T'
420*
421 IF( m.EQ.1 ) THEN
422 CALL strsm( 'L', 'L', 'T', diag, m1, n,
423 $ alpha,
424 $ a( 0 ), m, b, ldb )
425 ELSE
426 CALL strsm( 'L', 'U', 'N', diag, m2, n,
427 $ alpha,
428 $ a( m ), m, b( m1, 0 ), ldb )
429 CALL sgemm( 'T', 'N', m1, n, m2, -one,
430 $ a( m1 ),
431 $ m, b( m1, 0 ), ldb, alpha, b, ldb )
432 CALL strsm( 'L', 'L', 'T', diag, m1, n, one,
433 $ a( 0 ), m, b, ldb )
434 END IF
435*
436 END IF
437*
438 ELSE
439*
440* SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
441*
442 IF( .NOT.notrans ) THEN
443*
444* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
445* TRANS = 'N'
446*
447 CALL strsm( 'L', 'L', 'N', diag, m1, n, alpha,
448 $ a( m2 ), m, b, ldb )
449 CALL sgemm( 'T', 'N', m2, n, m1, -one, a( 0 ),
450 $ m,
451 $ b, ldb, alpha, b( m1, 0 ), ldb )
452 CALL strsm( 'L', 'U', 'T', diag, m2, n, one,
453 $ a( m1 ), m, b( m1, 0 ), ldb )
454*
455 ELSE
456*
457* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
458* TRANS = 'T'
459*
460 CALL strsm( 'L', 'U', 'N', diag, m2, n, alpha,
461 $ a( m1 ), m, b( m1, 0 ), ldb )
462 CALL sgemm( 'N', 'N', m1, n, m2, -one, a( 0 ),
463 $ m,
464 $ b( m1, 0 ), ldb, alpha, b, ldb )
465 CALL strsm( 'L', 'L', 'T', diag, m1, n, one,
466 $ a( m2 ), m, b, ldb )
467*
468 END IF
469*
470 END IF
471*
472 ELSE
473*
474* SIDE = 'L', N is odd, and TRANSR = 'T'
475*
476 IF( lower ) THEN
477*
478* SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
479*
480 IF( notrans ) THEN
481*
482* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
483* TRANS = 'N'
484*
485 IF( m.EQ.1 ) THEN
486 CALL strsm( 'L', 'U', 'T', diag, m1, n,
487 $ alpha,
488 $ a( 0 ), m1, b, ldb )
489 ELSE
490 CALL strsm( 'L', 'U', 'T', diag, m1, n,
491 $ alpha,
492 $ a( 0 ), m1, b, ldb )
493 CALL sgemm( 'T', 'N', m2, n, m1, -one,
494 $ a( m1*m1 ), m1, b, ldb, alpha,
495 $ b( m1, 0 ), ldb )
496 CALL strsm( 'L', 'L', 'N', diag, m2, n, one,
497 $ a( 1 ), m1, b( m1, 0 ), ldb )
498 END IF
499*
500 ELSE
501*
502* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
503* TRANS = 'T'
504*
505 IF( m.EQ.1 ) THEN
506 CALL strsm( 'L', 'U', 'N', diag, m1, n,
507 $ alpha,
508 $ a( 0 ), m1, b, ldb )
509 ELSE
510 CALL strsm( 'L', 'L', 'T', diag, m2, n,
511 $ alpha,
512 $ a( 1 ), m1, b( m1, 0 ), ldb )
513 CALL sgemm( 'N', 'N', m1, n, m2, -one,
514 $ a( m1*m1 ), m1, b( m1, 0 ), ldb,
515 $ alpha, b, ldb )
516 CALL strsm( 'L', 'U', 'N', diag, m1, n, one,
517 $ a( 0 ), m1, b, ldb )
518 END IF
519*
520 END IF
521*
522 ELSE
523*
524* SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
525*
526 IF( .NOT.notrans ) THEN
527*
528* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
529* TRANS = 'N'
530*
531 CALL strsm( 'L', 'U', 'T', diag, m1, n, alpha,
532 $ a( m2*m2 ), m2, b, ldb )
533 CALL sgemm( 'N', 'N', m2, n, m1, -one, a( 0 ),
534 $ m2,
535 $ b, ldb, alpha, b( m1, 0 ), ldb )
536 CALL strsm( 'L', 'L', 'N', diag, m2, n, one,
537 $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
538*
539 ELSE
540*
541* SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
542* TRANS = 'T'
543*
544 CALL strsm( 'L', 'L', 'T', diag, m2, n, alpha,
545 $ a( m1*m2 ), m2, b( m1, 0 ), ldb )
546 CALL sgemm( 'T', 'N', m1, n, m2, -one, a( 0 ),
547 $ m2,
548 $ b( m1, 0 ), ldb, alpha, b, ldb )
549 CALL strsm( 'L', 'U', 'N', diag, m1, n, one,
550 $ a( m2*m2 ), m2, b, ldb )
551*
552 END IF
553*
554 END IF
555*
556 END IF
557*
558 ELSE
559*
560* SIDE = 'L' and N is even
561*
562 IF( normaltransr ) THEN
563*
564* SIDE = 'L', N is even, and TRANSR = 'N'
565*
566 IF( lower ) THEN
567*
568* SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
569*
570 IF( notrans ) THEN
571*
572* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
573* and TRANS = 'N'
574*
575 CALL strsm( 'L', 'L', 'N', diag, k, n, alpha,
576 $ a( 1 ), m+1, b, ldb )
577 CALL sgemm( 'N', 'N', k, n, k, -one, a( k+1 ),
578 $ m+1, b, ldb, alpha, b( k, 0 ), ldb )
579 CALL strsm( 'L', 'U', 'T', diag, k, n, one,
580 $ a( 0 ), m+1, b( k, 0 ), ldb )
581*
582 ELSE
583*
584* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
585* and TRANS = 'T'
586*
587 CALL strsm( 'L', 'U', 'N', diag, k, n, alpha,
588 $ a( 0 ), m+1, b( k, 0 ), ldb )
589 CALL sgemm( 'T', 'N', k, n, k, -one, a( k+1 ),
590 $ m+1, b( k, 0 ), ldb, alpha, b, ldb )
591 CALL strsm( 'L', 'L', 'T', diag, k, n, one,
592 $ a( 1 ), m+1, b, ldb )
593*
594 END IF
595*
596 ELSE
597*
598* SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
599*
600 IF( .NOT.notrans ) THEN
601*
602* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
603* and TRANS = 'N'
604*
605 CALL strsm( 'L', 'L', 'N', diag, k, n, alpha,
606 $ a( k+1 ), m+1, b, ldb )
607 CALL sgemm( 'T', 'N', k, n, k, -one, a( 0 ),
608 $ m+1,
609 $ b, ldb, alpha, b( k, 0 ), ldb )
610 CALL strsm( 'L', 'U', 'T', diag, k, n, one,
611 $ a( k ), m+1, b( k, 0 ), ldb )
612*
613 ELSE
614*
615* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
616* and TRANS = 'T'
617 CALL strsm( 'L', 'U', 'N', diag, k, n, alpha,
618 $ a( k ), m+1, b( k, 0 ), ldb )
619 CALL sgemm( 'N', 'N', k, n, k, -one, a( 0 ),
620 $ m+1,
621 $ b( k, 0 ), ldb, alpha, b, ldb )
622 CALL strsm( 'L', 'L', 'T', diag, k, n, one,
623 $ a( k+1 ), m+1, b, ldb )
624*
625 END IF
626*
627 END IF
628*
629 ELSE
630*
631* SIDE = 'L', N is even, and TRANSR = 'T'
632*
633 IF( lower ) THEN
634*
635* SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L'
636*
637 IF( notrans ) THEN
638*
639* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
640* and TRANS = 'N'
641*
642 CALL strsm( 'L', 'U', 'T', diag, k, n, alpha,
643 $ a( k ), k, b, ldb )
644 CALL sgemm( 'T', 'N', k, n, k, -one,
645 $ a( k*( k+1 ) ), k, b, ldb, alpha,
646 $ b( k, 0 ), ldb )
647 CALL strsm( 'L', 'L', 'N', diag, k, n, one,
648 $ a( 0 ), k, b( k, 0 ), ldb )
649*
650 ELSE
651*
652* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
653* and TRANS = 'T'
654*
655 CALL strsm( 'L', 'L', 'T', diag, k, n, alpha,
656 $ a( 0 ), k, b( k, 0 ), ldb )
657 CALL sgemm( 'N', 'N', k, n, k, -one,
658 $ a( k*( k+1 ) ), k, b( k, 0 ), ldb,
659 $ alpha, b, ldb )
660 CALL strsm( 'L', 'U', 'N', diag, k, n, one,
661 $ a( k ), k, b, ldb )
662*
663 END IF
664*
665 ELSE
666*
667* SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U'
668*
669 IF( .NOT.notrans ) THEN
670*
671* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
672* and TRANS = 'N'
673*
674 CALL strsm( 'L', 'U', 'T', diag, k, n, alpha,
675 $ a( k*( k+1 ) ), k, b, ldb )
676 CALL sgemm( 'N', 'N', k, n, k, -one, a( 0 ), k,
677 $ b,
678 $ ldb, alpha, b( k, 0 ), ldb )
679 CALL strsm( 'L', 'L', 'N', diag, k, n, one,
680 $ a( k*k ), k, b( k, 0 ), ldb )
681*
682 ELSE
683*
684* SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
685* and TRANS = 'T'
686*
687 CALL strsm( 'L', 'L', 'T', diag, k, n, alpha,
688 $ a( k*k ), k, b( k, 0 ), ldb )
689 CALL sgemm( 'T', 'N', k, n, k, -one, a( 0 ), k,
690 $ b( k, 0 ), ldb, alpha, b, ldb )
691 CALL strsm( 'L', 'U', 'N', diag, k, n, one,
692 $ a( k*( k+1 ) ), k, b, ldb )
693*
694 END IF
695*
696 END IF
697*
698 END IF
699*
700 END IF
701*
702 ELSE
703*
704* SIDE = 'R'
705*
706* A is N-by-N.
707* If N is odd, set NISODD = .TRUE., and N1 and N2.
708* If N is even, NISODD = .FALSE., and K.
709*
710 IF( mod( n, 2 ).EQ.0 ) THEN
711 nisodd = .false.
712 k = n / 2
713 ELSE
714 nisodd = .true.
715 IF( lower ) THEN
716 n2 = n / 2
717 n1 = n - n2
718 ELSE
719 n1 = n / 2
720 n2 = n - n1
721 END IF
722 END IF
723*
724 IF( nisodd ) THEN
725*
726* SIDE = 'R' and N is odd
727*
728 IF( normaltransr ) THEN
729*
730* SIDE = 'R', N is odd, and TRANSR = 'N'
731*
732 IF( lower ) THEN
733*
734* SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
735*
736 IF( notrans ) THEN
737*
738* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
739* TRANS = 'N'
740*
741 CALL strsm( 'R', 'U', 'T', diag, m, n2, alpha,
742 $ a( n ), n, b( 0, n1 ), ldb )
743 CALL sgemm( 'N', 'N', m, n1, n2, -one, b( 0,
744 $ n1 ),
745 $ ldb, a( n1 ), n, alpha, b( 0, 0 ),
746 $ ldb )
747 CALL strsm( 'R', 'L', 'N', diag, m, n1, one,
748 $ a( 0 ), n, b( 0, 0 ), ldb )
749*
750 ELSE
751*
752* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
753* TRANS = 'T'
754*
755 CALL strsm( 'R', 'L', 'T', diag, m, n1, alpha,
756 $ a( 0 ), n, b( 0, 0 ), ldb )
757 CALL sgemm( 'N', 'T', m, n2, n1, -one, b( 0,
758 $ 0 ),
759 $ ldb, a( n1 ), n, alpha, b( 0, n1 ),
760 $ ldb )
761 CALL strsm( 'R', 'U', 'N', diag, m, n2, one,
762 $ a( n ), n, b( 0, n1 ), ldb )
763*
764 END IF
765*
766 ELSE
767*
768* SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
769*
770 IF( notrans ) THEN
771*
772* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
773* TRANS = 'N'
774*
775 CALL strsm( 'R', 'L', 'T', diag, m, n1, alpha,
776 $ a( n2 ), n, b( 0, 0 ), ldb )
777 CALL sgemm( 'N', 'N', m, n2, n1, -one, b( 0,
778 $ 0 ),
779 $ ldb, a( 0 ), n, alpha, b( 0, n1 ),
780 $ ldb )
781 CALL strsm( 'R', 'U', 'N', diag, m, n2, one,
782 $ a( n1 ), n, b( 0, n1 ), ldb )
783*
784 ELSE
785*
786* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
787* TRANS = 'T'
788*
789 CALL strsm( 'R', 'U', 'T', diag, m, n2, alpha,
790 $ a( n1 ), n, b( 0, n1 ), ldb )
791 CALL sgemm( 'N', 'T', m, n1, n2, -one, b( 0,
792 $ n1 ),
793 $ ldb, a( 0 ), n, alpha, b( 0, 0 ), ldb )
794 CALL strsm( 'R', 'L', 'N', diag, m, n1, one,
795 $ a( n2 ), n, b( 0, 0 ), ldb )
796*
797 END IF
798*
799 END IF
800*
801 ELSE
802*
803* SIDE = 'R', N is odd, and TRANSR = 'T'
804*
805 IF( lower ) THEN
806*
807* SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
808*
809 IF( notrans ) THEN
810*
811* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
812* TRANS = 'N'
813*
814 CALL strsm( 'R', 'L', 'N', diag, m, n2, alpha,
815 $ a( 1 ), n1, b( 0, n1 ), ldb )
816 CALL sgemm( 'N', 'T', m, n1, n2, -one, b( 0,
817 $ n1 ),
818 $ ldb, a( n1*n1 ), n1, alpha, b( 0, 0 ),
819 $ ldb )
820 CALL strsm( 'R', 'U', 'T', diag, m, n1, one,
821 $ a( 0 ), n1, b( 0, 0 ), ldb )
822*
823 ELSE
824*
825* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
826* TRANS = 'T'
827*
828 CALL strsm( 'R', 'U', 'N', diag, m, n1, alpha,
829 $ a( 0 ), n1, b( 0, 0 ), ldb )
830 CALL sgemm( 'N', 'N', m, n2, n1, -one, b( 0,
831 $ 0 ),
832 $ ldb, a( n1*n1 ), n1, alpha, b( 0, n1 ),
833 $ ldb )
834 CALL strsm( 'R', 'L', 'T', diag, m, n2, one,
835 $ a( 1 ), n1, b( 0, n1 ), ldb )
836*
837 END IF
838*
839 ELSE
840*
841* SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
842*
843 IF( notrans ) THEN
844*
845* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
846* TRANS = 'N'
847*
848 CALL strsm( 'R', 'U', 'N', diag, m, n1, alpha,
849 $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
850 CALL sgemm( 'N', 'T', m, n2, n1, -one, b( 0,
851 $ 0 ),
852 $ ldb, a( 0 ), n2, alpha, b( 0, n1 ),
853 $ ldb )
854 CALL strsm( 'R', 'L', 'T', diag, m, n2, one,
855 $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
856*
857 ELSE
858*
859* SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
860* TRANS = 'T'
861*
862 CALL strsm( 'R', 'L', 'N', diag, m, n2, alpha,
863 $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
864 CALL sgemm( 'N', 'N', m, n1, n2, -one, b( 0,
865 $ n1 ),
866 $ ldb, a( 0 ), n2, alpha, b( 0, 0 ),
867 $ ldb )
868 CALL strsm( 'R', 'U', 'T', diag, m, n1, one,
869 $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
870*
871 END IF
872*
873 END IF
874*
875 END IF
876*
877 ELSE
878*
879* SIDE = 'R' and N is even
880*
881 IF( normaltransr ) THEN
882*
883* SIDE = 'R', N is even, and TRANSR = 'N'
884*
885 IF( lower ) THEN
886*
887* SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
888*
889 IF( notrans ) THEN
890*
891* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
892* and TRANS = 'N'
893*
894 CALL strsm( 'R', 'U', 'T', diag, m, k, alpha,
895 $ a( 0 ), n+1, b( 0, k ), ldb )
896 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
897 $ ldb, a( k+1 ), n+1, alpha, b( 0, 0 ),
898 $ ldb )
899 CALL strsm( 'R', 'L', 'N', diag, m, k, one,
900 $ a( 1 ), n+1, b( 0, 0 ), ldb )
901*
902 ELSE
903*
904* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
905* and TRANS = 'T'
906*
907 CALL strsm( 'R', 'L', 'T', diag, m, k, alpha,
908 $ a( 1 ), n+1, b( 0, 0 ), ldb )
909 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
910 $ ldb, a( k+1 ), n+1, alpha, b( 0, k ),
911 $ ldb )
912 CALL strsm( 'R', 'U', 'N', diag, m, k, one,
913 $ a( 0 ), n+1, b( 0, k ), ldb )
914*
915 END IF
916*
917 ELSE
918*
919* SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
920*
921 IF( notrans ) THEN
922*
923* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
924* and TRANS = 'N'
925*
926 CALL strsm( 'R', 'L', 'T', diag, m, k, alpha,
927 $ a( k+1 ), n+1, b( 0, 0 ), ldb )
928 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
929 $ ldb, a( 0 ), n+1, alpha, b( 0, k ),
930 $ ldb )
931 CALL strsm( 'R', 'U', 'N', diag, m, k, one,
932 $ a( k ), n+1, b( 0, k ), ldb )
933*
934 ELSE
935*
936* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
937* and TRANS = 'T'
938*
939 CALL strsm( 'R', 'U', 'T', diag, m, k, alpha,
940 $ a( k ), n+1, b( 0, k ), ldb )
941 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
942 $ ldb, a( 0 ), n+1, alpha, b( 0, 0 ),
943 $ ldb )
944 CALL strsm( 'R', 'L', 'N', diag, m, k, one,
945 $ a( k+1 ), n+1, b( 0, 0 ), ldb )
946*
947 END IF
948*
949 END IF
950*
951 ELSE
952*
953* SIDE = 'R', N is even, and TRANSR = 'T'
954*
955 IF( lower ) THEN
956*
957* SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L'
958*
959 IF( notrans ) THEN
960*
961* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
962* and TRANS = 'N'
963*
964 CALL strsm( 'R', 'L', 'N', diag, m, k, alpha,
965 $ a( 0 ), k, b( 0, k ), ldb )
966 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, k ),
967 $ ldb, a( ( k+1 )*k ), k, alpha,
968 $ b( 0, 0 ), ldb )
969 CALL strsm( 'R', 'U', 'T', diag, m, k, one,
970 $ a( k ), k, b( 0, 0 ), ldb )
971*
972 ELSE
973*
974* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
975* and TRANS = 'T'
976*
977 CALL strsm( 'R', 'U', 'N', diag, m, k, alpha,
978 $ a( k ), k, b( 0, 0 ), ldb )
979 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, 0 ),
980 $ ldb, a( ( k+1 )*k ), k, alpha,
981 $ b( 0, k ), ldb )
982 CALL strsm( 'R', 'L', 'T', diag, m, k, one,
983 $ a( 0 ), k, b( 0, k ), ldb )
984*
985 END IF
986*
987 ELSE
988*
989* SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U'
990*
991 IF( notrans ) THEN
992*
993* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
994* and TRANS = 'N'
995*
996 CALL strsm( 'R', 'U', 'N', diag, m, k, alpha,
997 $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
998 CALL sgemm( 'N', 'T', m, k, k, -one, b( 0, 0 ),
999 $ ldb, a( 0 ), k, alpha, b( 0, k ), ldb )
1000 CALL strsm( 'R', 'L', 'T', diag, m, k, one,
1001 $ a( k*k ), k, b( 0, k ), ldb )
1002*
1003 ELSE
1004*
1005* SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
1006* and TRANS = 'T'
1007*
1008 CALL strsm( 'R', 'L', 'N', diag, m, k, alpha,
1009 $ a( k*k ), k, b( 0, k ), ldb )
1010 CALL sgemm( 'N', 'N', m, k, k, -one, b( 0, k ),
1011 $ ldb, a( 0 ), k, alpha, b( 0, 0 ), ldb )
1012 CALL strsm( 'R', 'U', 'T', diag, m, k, one,
1013 $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
1014*
1015 END IF
1016*
1017 END IF
1018*
1019 END IF
1020*
1021 END IF
1022 END IF
1023*
1024 RETURN
1025*
1026* End of STFSM
1027*
1028 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
stfsm
subroutine stfsm(transr, side, uplo, trans, diag, m, n, alpha, a, b, ldb)
STFSM solves a matrix equation (one operand is a triangular matrix in RFP format).
Definition stfsm.f:277
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181