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