OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ctfsm.f
Go to the documentation of this file.
1*> \brief \b CTFSM 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 CTFSM + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctfsm.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctfsm.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctfsm.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CTFSM( 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* COMPLEX ALPHA
28* ..
29* .. Array Arguments ..
30* COMPLEX 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*> CTFSM 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**H.
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*> = 'C': The Conjugate-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 = 'C' or 'c' op( A ) = conjg( 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 COMPLEX
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 COMPLEX array, dimension (N*(N+1)/2)
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 = 'C' then RFP is the Conjugate-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*> conjugate-transpose Format. If UPLO = 'L' the RFP A contains
153*> the NT elements of lower packed A either in normal or
154*> conjugate-transpose Format. The LDA of RFP A is (N+1)/2 when
155*> TRANSR = 'C'. 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 COMPLEX 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 complexOTHERcomputational
186*
187*> \par Further Details:
188* =====================
189*>
190*> \verbatim
191*>
192*> We first consider Standard Packed Format when N is even.
193*> 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*> conjugate-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*> conjugate-transpose of the last three columns of AP lower.
212*> To denote conjugate we place -- above the element. This covers the
213*> case N even and TRANSR = 'N'.
214*>
215*> RFP A RFP A
216*>
217*> -- -- --
218*> 03 04 05 33 43 53
219*> -- --
220*> 13 14 15 00 44 54
221*> --
222*> 23 24 25 10 11 55
223*>
224*> 33 34 35 20 21 22
225*> --
226*> 00 44 45 30 31 32
227*> -- --
228*> 01 11 55 40 41 42
229*> -- -- --
230*> 02 12 22 50 51 52
231*>
232*> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
233*> transpose of RFP A above. One therefore gets:
234*>
235*>
236*> RFP A RFP A
237*>
238*> -- -- -- -- -- -- -- -- -- --
239*> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
240*> -- -- -- -- -- -- -- -- -- --
241*> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
242*> -- -- -- -- -- -- -- -- -- --
243*> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
244*>
245*>
246*> We next consider Standard Packed Format when N is odd.
247*> We give an example where N = 5.
248*>
249*> AP is Upper AP is Lower
250*>
251*> 00 01 02 03 04 00
252*> 11 12 13 14 10 11
253*> 22 23 24 20 21 22
254*> 33 34 30 31 32 33
255*> 44 40 41 42 43 44
256*>
257*>
258*> Let TRANSR = 'N'. RFP holds AP as follows:
259*> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
260*> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
261*> conjugate-transpose of the first two columns of AP upper.
262*> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
263*> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
264*> conjugate-transpose of the last two columns of AP lower.
265*> To denote conjugate we place -- above the element. This covers the
266*> case N odd and TRANSR = 'N'.
267*>
268*> RFP A RFP A
269*>
270*> -- --
271*> 02 03 04 00 33 43
272*> --
273*> 12 13 14 10 11 44
274*>
275*> 22 23 24 20 21 22
276*> --
277*> 00 33 34 30 31 32
278*> -- --
279*> 01 11 44 40 41 42
280*>
281*> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
282*> transpose of RFP A above. One therefore gets:
283*>
284*>
285*> RFP A RFP A
286*>
287*> -- -- -- -- -- -- -- -- --
288*> 02 12 22 00 01 00 10 20 30 40 50
289*> -- -- -- -- -- -- -- -- --
290*> 03 13 23 33 11 33 11 21 31 41 51
291*> -- -- -- -- -- -- -- -- --
292*> 04 14 24 34 44 43 44 22 32 42 52
293*> \endverbatim
294*>
295* =====================================================================
296 SUBROUTINE ctfsm( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
297 $ B, LDB )
298*
299* -- LAPACK computational routine --
300* -- LAPACK is a software package provided by Univ. of Tennessee, --
301* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
302*
303* .. Scalar Arguments ..
304 CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
305 INTEGER LDB, M, N
306 COMPLEX ALPHA
307* ..
308* .. Array Arguments ..
309 COMPLEX A( 0: * ), B( 0: LDB-1, 0: * )
310* ..
311*
312* =====================================================================
313* ..
314* .. Parameters ..
315 COMPLEX CONE, CZERO
316 parameter( cone = ( 1.0e+0, 0.0e+0 ),
317 $ czero = ( 0.0e+0, 0.0e+0 ) )
318* ..
319* .. Local Scalars ..
320 LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR,
321 $ notrans
322 INTEGER M1, M2, N1, N2, K, INFO, I, J
323* ..
324* .. External Functions ..
325 LOGICAL LSAME
326 EXTERNAL lsame
327* ..
328* .. External Subroutines ..
329 EXTERNAL xerbla, cgemm, ctrsm
330* ..
331* .. Intrinsic Functions ..
332 INTRINSIC max, mod
333* ..
334* .. Executable Statements ..
335*
336* Test the input parameters.
337*
338 info = 0
339 normaltransr = lsame( transr, 'N' )
340 lside = lsame( side, 'L' )
341 lower = lsame( uplo, 'L' )
342 notrans = lsame( trans, 'N' )
343 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
344 info = -1
345 ELSE IF( .NOT.lside .AND. .NOT.lsame( side, 'R' ) ) THEN
346 info = -2
347 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
348 info = -3
349 ELSE IF( .NOT.notrans .AND. .NOT.lsame( trans, 'C' ) ) THEN
350 info = -4
351 ELSE IF( .NOT.lsame( diag, 'N' ) .AND. .NOT.lsame( diag, 'U' ) )
352 $ THEN
353 info = -5
354 ELSE IF( m.LT.0 ) THEN
355 info = -6
356 ELSE IF( n.LT.0 ) THEN
357 info = -7
358 ELSE IF( ldb.LT.max( 1, m ) ) THEN
359 info = -11
360 END IF
361 IF( info.NE.0 ) THEN
362 CALL xerbla( 'CTFSM ', -info )
363 RETURN
364 END IF
365*
366* Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
367*
368 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )
369 $ RETURN
370*
371* Quick return when ALPHA.EQ.(0E+0,0E+0)
372*
373 IF( alpha.EQ.czero ) THEN
374 DO 20 j = 0, n - 1
375 DO 10 i = 0, m - 1
376 b( i, j ) = czero
377 10 CONTINUE
378 20 CONTINUE
379 RETURN
380 END IF
381*
382 IF( lside ) THEN
383*
384* SIDE = 'L'
385*
386* A is M-by-M.
387* If M is odd, set NISODD = .TRUE., and M1 and M2.
388* If M is even, NISODD = .FALSE., and M.
389*
390 IF( mod( m, 2 ).EQ.0 ) THEN
391 misodd = .false.
392 k = m / 2
393 ELSE
394 misodd = .true.
395 IF( lower ) THEN
396 m2 = m / 2
397 m1 = m - m2
398 ELSE
399 m1 = m / 2
400 m2 = m - m1
401 END IF
402 END IF
403*
404 IF( misodd ) THEN
405*
406* SIDE = 'L' and N is odd
407*
408 IF( normaltransr ) THEN
409*
410* SIDE = 'L', N is odd, and TRANSR = 'N'
411*
412 IF( lower ) THEN
413*
414* SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
415*
416 IF( notrans ) THEN
417*
418* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
419* TRANS = 'N'
420*
421 IF( m.EQ.1 ) THEN
422 CALL ctrsm( 'L', 'L', 'N', diag, m1, n, alpha,
423 $ a, m, b, ldb )
424 ELSE
425 CALL ctrsm( 'L', 'L', 'N', diag, m1, n, alpha,
426 $ a( 0 ), m, b, ldb )
427 CALL cgemm( 'N', 'N', m2, n, m1, -cone, a( m1 ),
428 $ m, b, ldb, alpha, b( m1, 0 ), ldb )
429 CALL ctrsm( 'L', 'U', 'C', diag, m2, n, cone,
430 $ a( m ), m, b( m1, 0 ), ldb )
431 END IF
432*
433 ELSE
434*
435* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
436* TRANS = 'C'
437*
438 IF( m.EQ.1 ) THEN
439 CALL ctrsm( 'L', 'L', 'C', diag, m1, n, alpha,
440 $ a( 0 ), m, b, ldb )
441 ELSE
442 CALL ctrsm( 'L', 'U', 'N', diag, m2, n, alpha,
443 $ a( m ), m, b( m1, 0 ), ldb )
444 CALL cgemm( 'C', 'N', m1, n, m2, -cone, a( m1 ),
445 $ m, b( m1, 0 ), ldb, alpha, b, ldb )
446 CALL ctrsm( 'L', 'L', 'C', diag, m1, n, cone,
447 $ a( 0 ), m, b, ldb )
448 END IF
449*
450 END IF
451*
452 ELSE
453*
454* SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
455*
456 IF( .NOT.notrans ) THEN
457*
458* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
459* TRANS = 'N'
460*
461 CALL ctrsm( 'L', 'L', 'N', diag, m1, n, alpha,
462 $ a( m2 ), m, b, ldb )
463 CALL cgemm( 'c', 'n', M2, N, M1, -CONE, A( 0 ), M,
464 $ B, LDB, ALPHA, B( M1, 0 ), LDB )
465 CALL CTRSM( 'l', 'u', 'c', DIAG, M2, N, CONE,
466 $ A( M1 ), M, B( M1, 0 ), LDB )
467*
468 ELSE
469*
470* SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
471* TRANS = 'C'
472*
473 CALL CTRSM( 'l', 'u', 'n', DIAG, M2, N, ALPHA,
474 $ A( M1 ), M, B( M1, 0 ), LDB )
475 CALL CGEMM( 'n', 'n', M1, N, M2, -CONE, A( 0 ), M,
476 $ B( M1, 0 ), LDB, ALPHA, B, LDB )
477 CALL CTRSM( 'l', 'l', 'c', DIAG, M1, N, CONE,
478 $ A( M2 ), M, B, LDB )
479*
480 END IF
481*
482 END IF
483*
484 ELSE
485*
486* SIDE = 'L', N is odd, and TRANSR = 'C'
487*
488 IF( LOWER ) THEN
489*
490* SIDE ='L', N is odd, TRANSR = 'C', and UPLO = 'L'
491*
492 IF( NOTRANS ) THEN
493*
494* SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'L', and
495* TRANS = 'N'
496*
497.EQ. IF( M1 ) THEN
498 CALL CTRSM( 'l', 'u', 'c', DIAG, M1, N, ALPHA,
499 $ A( 0 ), M1, B, LDB )
500 ELSE
501 CALL CTRSM( 'l', 'u', 'c', DIAG, M1, N, ALPHA,
502 $ A( 0 ), M1, B, LDB )
503 CALL CGEMM( 'c', 'n', M2, N, M1, -CONE,
504 $ A( M1*M1 ), M1, B, LDB, ALPHA,
505 $ B( M1, 0 ), LDB )
506 CALL CTRSM( 'l', 'l', 'n', DIAG, M2, N, CONE,
507 $ A( 1 ), M1, B( M1, 0 ), LDB )
508 END IF
509*
510 ELSE
511*
512* SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'L', and
513* TRANS = 'C'
514*
515.EQ. IF( M1 ) THEN
516 CALL CTRSM( 'l', 'u', 'n', DIAG, M1, N, ALPHA,
517 $ A( 0 ), M1, B, LDB )
518 ELSE
519 CALL CTRSM( 'l', 'l', 'c', DIAG, M2, N, ALPHA,
520 $ A( 1 ), M1, B( M1, 0 ), LDB )
521 CALL CGEMM( 'n', 'n', M1, N, M2, -CONE,
522 $ A( M1*M1 ), M1, B( M1, 0 ), LDB,
523 $ ALPHA, B, LDB )
524 CALL CTRSM( 'l', 'u', 'n', DIAG, M1, N, CONE,
525 $ A( 0 ), M1, B, LDB )
526 END IF
527*
528 END IF
529*
530 ELSE
531*
532* SIDE ='L', N is odd, TRANSR = 'C', and UPLO = 'U'
533*
534.NOT. IF( NOTRANS ) THEN
535*
536* SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'U', and
537* TRANS = 'N'
538*
539 CALL CTRSM( 'l', 'u', 'c', DIAG, M1, N, ALPHA,
540 $ A( M2*M2 ), M2, B, LDB )
541 CALL CGEMM( 'n', 'n', M2, N, M1, -CONE, A( 0 ), M2,
542 $ B, LDB, ALPHA, B( M1, 0 ), LDB )
543 CALL CTRSM( 'l', 'l', 'n', DIAG, M2, N, CONE,
544 $ A( M1*M2 ), M2, B( M1, 0 ), LDB )
545*
546 ELSE
547*
548* SIDE ='L', N is odd, TRANSR = 'C', UPLO = 'U', and
549* TRANS = 'C'
550*
551 CALL CTRSM( 'l', 'l', 'c', DIAG, M2, N, ALPHA,
552 $ A( M1*M2 ), M2, B( M1, 0 ), LDB )
553 CALL CGEMM( 'c', 'n', M1, N, M2, -CONE, A( 0 ), M2,
554 $ B( M1, 0 ), LDB, ALPHA, B, LDB )
555 CALL CTRSM( 'l', 'u', 'n', DIAG, M1, N, CONE,
556 $ A( M2*M2 ), M2, B, LDB )
557*
558 END IF
559*
560 END IF
561*
562 END IF
563*
564 ELSE
565*
566* SIDE = 'L' and N is even
567*
568 IF( NORMALTRANSR ) THEN
569*
570* SIDE = 'L', N is even, and TRANSR = 'N'
571*
572 IF( LOWER ) THEN
573*
574* SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
575*
576 IF( NOTRANS ) THEN
577*
578* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
579* and TRANS = 'N'
580*
581 CALL CTRSM( 'l', 'l', 'n', DIAG, K, N, ALPHA,
582 $ A( 1 ), M+1, B, LDB )
583 CALL CGEMM( 'n', 'n', K, N, K, -CONE, A( K+1 ),
584 $ M+1, B, LDB, ALPHA, B( K, 0 ), LDB )
585 CALL CTRSM( 'l', 'U', 'C', diag, k, n, cone,
586 $ a( 0 ), m+1, b( k, 0 ), ldb )
587*
588 ELSE
589*
590* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
591* and TRANS = 'C'
592*
593 CALL ctrsm( 'L', 'U', 'N', diag, k, n, alpha,
594 $ a( 0 ), m+1, b( k, 0 ), ldb )
595 CALL cgemm( 'C', 'N', k, n, k, -cone, a( k+1 ),
596 $ m+1, b( k, 0 ), ldb, alpha, b, ldb )
597 CALL ctrsm( 'L', 'L', 'C', diag, k, n, cone,
598 $ a( 1 ), m+1, b, ldb )
599*
600 END IF
601*
602 ELSE
603*
604* SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
605*
606 IF( .NOT.notrans ) THEN
607*
608* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
609* and TRANS = 'N'
610*
611 CALL ctrsm( 'L', 'L', 'N', diag, k, n, alpha,
612 $ a( k+1 ), m+1, b, ldb )
613 CALL cgemm( 'C', 'N', k, n, k, -cone, a( 0 ), m+1,
614 $ b, ldb, alpha, b( k, 0 ), ldb )
615 CALL ctrsm( 'L', 'U', 'C', diag, k, n, cone,
616 $ a( k ), m+1, b( k, 0 ), ldb )
617*
618 ELSE
619*
620* SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
621* and TRANS = 'C'
622 CALL ctrsm( 'L', 'U', 'N', diag, k, n, alpha,
623 $ a( k ), m+1, b( k, 0 ), ldb )
624 CALL cgemm( 'N', 'N', k, n, k, -cone, a( 0 ), m+1,
625 $ b( k, 0 ), ldb, alpha, b, ldb )
626 CALL ctrsm( 'L', 'L', 'C', diag, k, n, cone,
627 $ a( k+1 ), m+1, b, ldb )
628*
629 END IF
630*
631 END IF
632*
633 ELSE
634*
635* SIDE = 'L', N is even, and TRANSR = 'C'
636*
637 IF( lower ) THEN
638*
639* SIDE ='L', N is even, TRANSR = 'C', and UPLO = 'L'
640*
641 IF( notrans ) THEN
642*
643* SIDE ='L', N is even, TRANSR = 'C', UPLO = 'L',
644* and TRANS = 'N'
645*
646 CALL ctrsm( 'L', 'U', 'C', diag, k, n, alpha,
647 $ a( k ), k, b, ldb )
648 CALL cgemm( 'C', 'N', k, n, k, -cone,
649 $ a( k*( k+1 ) ), k, b, ldb, alpha,
650 $ b( k, 0 ), ldb )
651 CALL ctrsm( 'L', 'L', 'N', diag, k, n, cone,
652 $ a( 0 ), k, b( k, 0 ), ldb )
653*
654 ELSE
655*
656* SIDE ='L', N is even, TRANSR = 'C', UPLO = 'L',
657* and TRANS = 'C'
658*
659 CALL ctrsm( 'L', 'L', 'C', diag, k, n, alpha,
660 $ a( 0 ), k, b( k, 0 ), ldb )
661 CALL cgemm( 'N', 'N', k, n, k, -cone,
662 $ a( k*( k+1 ) ), k, b( k, 0 ), ldb,
663 $ alpha, b, ldb )
664 CALL ctrsm( 'L', 'U', 'N', diag, k, n, cone,
665 $ a( k ), k, b, ldb )
666*
667 END IF
668*
669 ELSE
670*
671* SIDE ='L', N is even, TRANSR = 'C', and UPLO = 'U'
672*
673 IF( .NOT.notrans ) THEN
674*
675* SIDE ='L', N is even, TRANSR = 'C', UPLO = 'U',
676* and TRANS = 'N'
677*
678 CALL ctrsm( 'L', 'U', 'C', diag, k, n, alpha,
679 $ a( k*( k+1 ) ), k, b, ldb )
680 CALL cgemm( 'N', 'N', k, n, k, -cone, a( 0 ), k, b,
681 $ ldb, alpha, b( k, 0 ), ldb )
682 CALL ctrsm( 'L', 'L', 'N', diag, k, n, cone,
683 $ a( k*k ), k, b( k, 0 ), ldb )
684*
685 ELSE
686*
687* SIDE ='L', N is even, TRANSR = 'C', UPLO = 'U',
688* and TRANS = 'C'
689*
690 CALL ctrsm( 'L', 'L', 'C', diag, k, n, alpha,
691 $ a( k*k ), k, b( k, 0 ), ldb )
692 CALL cgemm( 'C', 'N', k, n, k, -cone, a( 0 ), k,
693 $ b( k, 0 ), ldb, alpha, b, ldb )
694 CALL ctrsm( 'L', 'U', 'N', diag, k, n, cone,
695 $ a( k*( k+1 ) ), k, b, ldb )
696*
697 END IF
698*
699 END IF
700*
701 END IF
702*
703 END IF
704*
705 ELSE
706*
707* SIDE = 'R'
708*
709* A is N-by-N.
710* If N is odd, set NISODD = .TRUE., and N1 and N2.
711* If N is even, NISODD = .FALSE., and K.
712*
713 IF( mod( n, 2 ).EQ.0 ) THEN
714 nisodd = .false.
715 k = n / 2
716 ELSE
717 nisodd = .true.
718 IF( lower ) THEN
719 n2 = n / 2
720 n1 = n - n2
721 ELSE
722 n1 = n / 2
723 n2 = n - n1
724 END IF
725 END IF
726*
727 IF( nisodd ) THEN
728*
729* SIDE = 'R' and N is odd
730*
731 IF( normaltransr ) THEN
732*
733* SIDE = 'R', N is odd, and TRANSR = 'N'
734*
735 IF( lower ) THEN
736*
737* SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
738*
739 IF( notrans ) THEN
740*
741* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
742* TRANS = 'N'
743*
744 CALL ctrsm( 'R', 'U', 'C', diag, m, n2, alpha,
745 $ a( n ), n, b( 0, n1 ), ldb )
746 CALL cgemm( 'N', 'N', m, n1, n2, -cone, b( 0, n1 ),
747 $ ldb, a( n1 ), n, alpha, b( 0, 0 ),
748 $ ldb )
749 CALL ctrsm( 'R', 'L', 'N', diag, m, n1, cone,
750 $ a( 0 ), n, b( 0, 0 ), ldb )
751*
752 ELSE
753*
754* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
755* TRANS = 'C'
756*
757 CALL ctrsm( 'R', 'L', 'C', diag, m, n1, alpha,
758 $ a( 0 ), n, b( 0, 0 ), ldb )
759 CALL cgemm( 'N', 'C', m, n2, n1, -cone, b( 0, 0 ),
760 $ ldb, a( n1 ), n, alpha, b( 0, n1 ),
761 $ ldb )
762 CALL ctrsm( 'R', 'U', 'N', diag, m, n2, cone,
763 $ a( n ), n, b( 0, n1 ), ldb )
764*
765 END IF
766*
767 ELSE
768*
769* SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
770*
771 IF( notrans ) THEN
772*
773* SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
774* TRANS = 'N'
775*
776 CALL ctrsm( 'R', 'L', 'C', diag, m, n1, alpha,
777 $ a( n2 ), n, b( 0, 0 ), ldb )
778 CALL cgemm( 'N', 'N', m, n2, n1, -cone, b( 0, 0 ),
779 $ ldb, a( 0 ), n, alpha, b( 0, n1 ),
780 $ ldb )
781 CALL ctrsm( 'R', 'U', 'N', diag, m, n2, cone,
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 = 'C'
788*
789 CALL ctrsm( 'R', 'U', 'C', diag, m, n2, alpha,
790 $ a( n1 ), n, b( 0, n1 ), ldb )
791 CALL cgemm( 'N', 'C', m, n1, n2, -cone, b( 0, n1 ),
792 $ ldb, a( 0 ), n, alpha, b( 0, 0 ), ldb )
793 CALL ctrsm( 'R', 'L', 'N', diag, m, n1, cone,
794 $ a( n2 ), n, b( 0, 0 ), ldb )
795*
796 END IF
797*
798 END IF
799*
800 ELSE
801*
802* SIDE = 'R', N is odd, and TRANSR = 'C'
803*
804 IF( lower ) THEN
805*
806* SIDE ='R', N is odd, TRANSR = 'C', and UPLO = 'L'
807*
808 IF( notrans ) THEN
809*
810* SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'L', and
811* TRANS = 'N'
812*
813 CALL ctrsm( 'R', 'L', 'N', diag, m, n2, alpha,
814 $ a( 1 ), n1, b( 0, n1 ), ldb )
815 CALL cgemm( 'N', 'C', m, n1, n2, -cone, b( 0, n1 ),
816 $ ldb, a( n1*n1 ), n1, alpha, b( 0, 0 ),
817 $ ldb )
818 CALL ctrsm( 'R', 'U', 'C', diag, m, n1, cone,
819 $ a( 0 ), n1, b( 0, 0 ), ldb )
820*
821 ELSE
822*
823* SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'L', and
824* TRANS = 'C'
825*
826 CALL ctrsm( 'R', 'U', 'N', diag, m, n1, alpha,
827 $ a( 0 ), n1, b( 0, 0 ), ldb )
828 CALL cgemm( 'N', 'N', m, n2, n1, -cone, b( 0, 0 ),
829 $ ldb, a( n1*n1 ), n1, alpha, b( 0, n1 ),
830 $ ldb )
831 CALL ctrsm( 'R', 'L', 'C', diag, m, n2, cone,
832 $ a( 1 ), n1, b( 0, n1 ), ldb )
833*
834 END IF
835*
836 ELSE
837*
838* SIDE ='R', N is odd, TRANSR = 'C', and UPLO = 'U'
839*
840 IF( notrans ) THEN
841*
842* SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'U', and
843* TRANS = 'N'
844*
845 CALL ctrsm( 'R', 'U', 'N', diag, m, n1, alpha,
846 $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
847 CALL cgemm( 'N', 'C', m, n2, n1, -cone, b( 0, 0 ),
848 $ ldb, a( 0 ), n2, alpha, b( 0, n1 ),
849 $ ldb )
850 CALL ctrsm( 'R', 'L', 'C', diag, m, n2, cone,
851 $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
852*
853 ELSE
854*
855* SIDE ='R', N is odd, TRANSR = 'C', UPLO = 'U', and
856* TRANS = 'C'
857*
858 CALL ctrsm( 'R', 'L', 'N', diag, m, n2, alpha,
859 $ a( n1*n2 ), n2, b( 0, n1 ), ldb )
860 CALL cgemm( 'N', 'N', m, n1, n2, -cone, b( 0, n1 ),
861 $ ldb, a( 0 ), n2, alpha, b( 0, 0 ),
862 $ ldb )
863 CALL ctrsm( 'R', 'U', 'C', diag, m, n1, cone,
864 $ a( n2*n2 ), n2, b( 0, 0 ), ldb )
865*
866 END IF
867*
868 END IF
869*
870 END IF
871*
872 ELSE
873*
874* SIDE = 'R' and N is even
875*
876 IF( normaltransr ) THEN
877*
878* SIDE = 'R', N is even, and TRANSR = 'N'
879*
880 IF( lower ) THEN
881*
882* SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
883*
884 IF( notrans ) THEN
885*
886* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
887* and TRANS = 'N'
888*
889 CALL ctrsm( 'R', 'U', 'C', diag, m, k, alpha,
890 $ a( 0 ), n+1, b( 0, k ), ldb )
891 CALL cgemm( 'N', 'N', m, k, k, -cone, b( 0, k ),
892 $ ldb, a( k+1 ), n+1, alpha, b( 0, 0 ),
893 $ ldb )
894 CALL ctrsm( 'R', 'L', 'N', diag, m, k, cone,
895 $ a( 1 ), n+1, b( 0, 0 ), ldb )
896*
897 ELSE
898*
899* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
900* and TRANS = 'C'
901*
902 CALL ctrsm( 'R', 'L', 'C', diag, m, k, alpha,
903 $ a( 1 ), n+1, b( 0, 0 ), ldb )
904 CALL cgemm( 'N', 'C', m, k, k, -cone, b( 0, 0 ),
905 $ ldb, a( k+1 ), n+1, alpha, b( 0, k ),
906 $ ldb )
907 CALL ctrsm( 'R', 'U', 'N', diag, m, k, cone,
908 $ a( 0 ), n+1, b( 0, k ), ldb )
909*
910 END IF
911*
912 ELSE
913*
914* SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
915*
916 IF( notrans ) THEN
917*
918* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
919* and TRANS = 'N'
920*
921 CALL ctrsm( 'R', 'L', 'C', diag, m, k, alpha,
922 $ a( k+1 ), n+1, b( 0, 0 ), ldb )
923 CALL cgemm( 'N', 'N', m, k, k, -cone, b( 0, 0 ),
924 $ ldb, a( 0 ), n+1, alpha, b( 0, k ),
925 $ ldb )
926 CALL ctrsm( 'R', 'U', 'N', diag, m, k, cone,
927 $ a( k ), n+1, b( 0, k ), ldb )
928*
929 ELSE
930*
931* SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
932* and TRANS = 'C'
933*
934 CALL ctrsm( 'R', 'U', 'C', diag, m, k, alpha,
935 $ a( k ), n+1, b( 0, k ), ldb )
936 CALL cgemm( 'N', 'C', m, k, k, -cone, b( 0, k ),
937 $ ldb, a( 0 ), n+1, alpha, b( 0, 0 ),
938 $ ldb )
939 CALL ctrsm( 'R', 'L', 'N', diag, m, k, cone,
940 $ a( k+1 ), n+1, b( 0, 0 ), ldb )
941*
942 END IF
943*
944 END IF
945*
946 ELSE
947*
948* SIDE = 'R', N is even, and TRANSR = 'C'
949*
950 IF( lower ) THEN
951*
952* SIDE ='R', N is even, TRANSR = 'C', and UPLO = 'L'
953*
954 IF( notrans ) THEN
955*
956* SIDE ='R', N is even, TRANSR = 'C', UPLO = 'L',
957* and TRANS = 'N'
958*
959 CALL ctrsm( 'R', 'L', 'N', diag, m, k, alpha,
960 $ a( 0 ), k, b( 0, k ), ldb )
961 CALL cgemm( 'N', 'C', m, k, k, -cone, b( 0, k ),
962 $ ldb, a( ( k+1 )*k ), k, alpha,
963 $ b( 0, 0 ), ldb )
964 CALL ctrsm( 'R', 'U', 'C', diag, m, k, cone,
965 $ a( k ), k, b( 0, 0 ), ldb )
966*
967 ELSE
968*
969* SIDE ='R', N is even, TRANSR = 'C', UPLO = 'L',
970* and TRANS = 'C'
971*
972 CALL ctrsm( 'R', 'U', 'N', diag, m, k, alpha,
973 $ a( k ), k, b( 0, 0 ), ldb )
974 CALL cgemm( 'N', 'N', m, k, k, -cone, b( 0, 0 ),
975 $ ldb, a( ( k+1 )*k ), k, alpha,
976 $ b( 0, k ), ldb )
977 CALL ctrsm( 'R', 'L', 'C', diag, m, k, cone,
978 $ a( 0 ), k, b( 0, k ), ldb )
979*
980 END IF
981*
982 ELSE
983*
984* SIDE ='R', N is even, TRANSR = 'C', and UPLO = 'U'
985*
986 IF( notrans ) THEN
987*
988* SIDE ='R', N is even, TRANSR = 'C', UPLO = 'U',
989* and TRANS = 'N'
990*
991 CALL ctrsm( 'R', 'U', 'N', diag, m, k, alpha,
992 $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
993 CALL cgemm( 'N', 'C', m, k, k, -cone, b( 0, 0 ),
994 $ ldb, a( 0 ), k, alpha, b( 0, k ), ldb )
995 CALL ctrsm( 'R', 'L', 'C', diag, m, k, cone,
996 $ a( k*k ), k, b( 0, k ), ldb )
997*
998 ELSE
999*
1000* SIDE ='R', N is even, TRANSR = 'C', UPLO = 'U',
1001* and TRANS = 'C'
1002*
1003 CALL ctrsm( 'R', 'L', 'N', diag, m, k, alpha,
1004 $ a( k*k ), k, b( 0, k ), ldb )
1005 CALL cgemm( 'N', 'N', m, k, k, -cone, b( 0, k ),
1006 $ ldb, a( 0 ), k, alpha, b( 0, 0 ), ldb )
1007 CALL ctrsm( 'R', 'U', 'C', diag, m, k, cone,
1008 $ a( ( k+1 )*k ), k, b( 0, 0 ), ldb )
1009*
1010 END IF
1011*
1012 END IF
1013*
1014 END IF
1015*
1016 END IF
1017 END IF
1018*
1019 RETURN
1020*
1021* End of CTFSM
1022*
1023 END
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
ctfsm
subroutine ctfsm(transr, side, uplo, trans, diag, m, n, alpha, a, b, ldb)
CTFSM solves a matrix equation (one operand is a triangular matrix in RFP format).
Definition ctfsm.f:298
ctrsm
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187
max
#define max(a, b)
Definition macros.h:21