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