LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ctftri.f
Go to the documentation of this file.
1*> \brief \b CTFTRI
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CTFTRI + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctftri.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctftri.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctftri.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CTFTRI( TRANSR, UPLO, DIAG, N, A, INFO )
20*
21* .. Scalar Arguments ..
22* CHARACTER TRANSR, UPLO, DIAG
23* INTEGER INFO, N
24* ..
25* .. Array Arguments ..
26* COMPLEX A( 0: * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> CTFTRI computes the inverse of a triangular matrix A stored in RFP
36*> format.
37*>
38*> This is a Level 3 BLAS version of the algorithm.
39*> \endverbatim
40*
41* Arguments:
42* ==========
43*
44*> \param[in] TRANSR
45*> \verbatim
46*> TRANSR is CHARACTER*1
47*> = 'N': The Normal TRANSR of RFP A is stored;
48*> = 'C': The Conjugate-transpose TRANSR of RFP A is stored.
49*> \endverbatim
50*>
51*> \param[in] UPLO
52*> \verbatim
53*> UPLO is CHARACTER*1
54*> = 'U': A is upper triangular;
55*> = 'L': A is lower triangular.
56*> \endverbatim
57*>
58*> \param[in] DIAG
59*> \verbatim
60*> DIAG is CHARACTER*1
61*> = 'N': A is non-unit triangular;
62*> = 'U': A is unit triangular.
63*> \endverbatim
64*>
65*> \param[in] N
66*> \verbatim
67*> N is INTEGER
68*> The order of the matrix A. N >= 0.
69*> \endverbatim
70*>
71*> \param[in,out] A
72*> \verbatim
73*> A is COMPLEX array, dimension ( N*(N+1)/2 );
74*> On entry, the triangular matrix A in RFP format. RFP format
75*> is described by TRANSR, UPLO, and N as follows: If TRANSR =
76*> 'N' then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
77*> (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
78*> the Conjugate-transpose of RFP A as defined when
79*> TRANSR = 'N'. The contents of RFP A are defined by UPLO as
80*> follows: If UPLO = 'U' the RFP A contains the nt elements of
81*> upper packed A; If UPLO = 'L' the RFP A contains the nt
82*> elements of lower packed A. The LDA of RFP A is (N+1)/2 when
83*> TRANSR = 'C'. When TRANSR is 'N' the LDA is N+1 when N is
84*> even and N is odd. See the Note below for more details.
85*>
86*> On exit, the (triangular) inverse of the original matrix, in
87*> the same storage format.
88*> \endverbatim
89*>
90*> \param[out] INFO
91*> \verbatim
92*> INFO is INTEGER
93*> = 0: successful exit
94*> < 0: if INFO = -i, the i-th argument had an illegal value
95*> > 0: if INFO = i, A(i,i) is exactly zero. The triangular
96*> matrix is singular and its inverse can not be computed.
97*> \endverbatim
98*
99* Authors:
100* ========
101*
102*> \author Univ. of Tennessee
103*> \author Univ. of California Berkeley
104*> \author Univ. of Colorado Denver
105*> \author NAG Ltd.
106*
107*> \ingroup tftri
108*
109*> \par Further Details:
110* =====================
111*>
112*> \verbatim
113*>
114*> We first consider Standard Packed Format when N is even.
115*> We give an example where N = 6.
116*>
117*> AP is Upper AP is Lower
118*>
119*> 00 01 02 03 04 05 00
120*> 11 12 13 14 15 10 11
121*> 22 23 24 25 20 21 22
122*> 33 34 35 30 31 32 33
123*> 44 45 40 41 42 43 44
124*> 55 50 51 52 53 54 55
125*>
126*>
127*> Let TRANSR = 'N'. RFP holds AP as follows:
128*> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
129*> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
130*> conjugate-transpose of the first three columns of AP upper.
131*> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
132*> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
133*> conjugate-transpose of the last three columns of AP lower.
134*> To denote conjugate we place -- above the element. This covers the
135*> case N even and TRANSR = 'N'.
136*>
137*> RFP A RFP A
138*>
139*> -- -- --
140*> 03 04 05 33 43 53
141*> -- --
142*> 13 14 15 00 44 54
143*> --
144*> 23 24 25 10 11 55
145*>
146*> 33 34 35 20 21 22
147*> --
148*> 00 44 45 30 31 32
149*> -- --
150*> 01 11 55 40 41 42
151*> -- -- --
152*> 02 12 22 50 51 52
153*>
154*> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
155*> transpose of RFP A above. One therefore gets:
156*>
157*>
158*> RFP A RFP A
159*>
160*> -- -- -- -- -- -- -- -- -- --
161*> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
162*> -- -- -- -- -- -- -- -- -- --
163*> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
164*> -- -- -- -- -- -- -- -- -- --
165*> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
166*>
167*>
168*> We next consider Standard Packed Format when N is odd.
169*> We give an example where N = 5.
170*>
171*> AP is Upper AP is Lower
172*>
173*> 00 01 02 03 04 00
174*> 11 12 13 14 10 11
175*> 22 23 24 20 21 22
176*> 33 34 30 31 32 33
177*> 44 40 41 42 43 44
178*>
179*>
180*> Let TRANSR = 'N'. RFP holds AP as follows:
181*> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
182*> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
183*> conjugate-transpose of the first two columns of AP upper.
184*> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
185*> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
186*> conjugate-transpose of the last two columns of AP lower.
187*> To denote conjugate we place -- above the element. This covers the
188*> case N odd and TRANSR = 'N'.
189*>
190*> RFP A RFP A
191*>
192*> -- --
193*> 02 03 04 00 33 43
194*> --
195*> 12 13 14 10 11 44
196*>
197*> 22 23 24 20 21 22
198*> --
199*> 00 33 34 30 31 32
200*> -- --
201*> 01 11 44 40 41 42
202*>
203*> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
204*> transpose of RFP A above. One therefore gets:
205*>
206*>
207*> RFP A RFP A
208*>
209*> -- -- -- -- -- -- -- -- --
210*> 02 12 22 00 01 00 10 20 30 40 50
211*> -- -- -- -- -- -- -- -- --
212*> 03 13 23 33 11 33 11 21 31 41 51
213*> -- -- -- -- -- -- -- -- --
214*> 04 14 24 34 44 43 44 22 32 42 52
215*> \endverbatim
216*>
217* =====================================================================
218 SUBROUTINE ctftri( TRANSR, UPLO, DIAG, N, A, INFO )
219*
220* -- LAPACK computational routine --
221* -- LAPACK is a software package provided by Univ. of Tennessee, --
222* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223*
224* .. Scalar Arguments ..
225 CHARACTER TRANSR, UPLO, DIAG
226 INTEGER INFO, N
227* ..
228* .. Array Arguments ..
229 COMPLEX A( 0: * )
230* ..
231*
232* =====================================================================
233*
234* .. Parameters ..
235 COMPLEX CONE
236 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
237* ..
238* .. Local Scalars ..
239 LOGICAL LOWER, NISODD, NORMALTRANSR
240 INTEGER N1, N2, K
241* ..
242* .. External Functions ..
243 LOGICAL LSAME
244 EXTERNAL lsame
245* ..
246* .. External Subroutines ..
247 EXTERNAL xerbla, ctrmm, ctrtri
248* ..
249* .. Intrinsic Functions ..
250 INTRINSIC mod
251* ..
252* .. Executable Statements ..
253*
254* Test the input parameters.
255*
256 info = 0
257 normaltransr = lsame( transr, 'N' )
258 lower = lsame( uplo, 'L' )
259 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
260 info = -1
261 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
262 info = -2
263 ELSE IF( .NOT.lsame( diag, 'N' ) .AND.
264 $ .NOT.lsame( diag, 'U' ) )
265 $ THEN
266 info = -3
267 ELSE IF( n.LT.0 ) THEN
268 info = -4
269 END IF
270 IF( info.NE.0 ) THEN
271 CALL xerbla( 'CTFTRI', -info )
272 RETURN
273 END IF
274*
275* Quick return if possible
276*
277 IF( n.EQ.0 )
278 $ RETURN
279*
280* If N is odd, set NISODD = .TRUE.
281* If N is even, set K = N/2 and NISODD = .FALSE.
282*
283 IF( mod( n, 2 ).EQ.0 ) THEN
284 k = n / 2
285 nisodd = .false.
286 ELSE
287 nisodd = .true.
288 END IF
289*
290* Set N1 and N2 depending on LOWER
291*
292 IF( lower ) THEN
293 n2 = n / 2
294 n1 = n - n2
295 ELSE
296 n1 = n / 2
297 n2 = n - n1
298 END IF
299*
300*
301* start execution: there are eight cases
302*
303 IF( nisodd ) THEN
304*
305* N is odd
306*
307 IF( normaltransr ) THEN
308*
309* N is odd and TRANSR = 'N'
310*
311 IF( lower ) THEN
312*
313* SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
314* T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
315* T1 -> a(0), T2 -> a(n), S -> a(n1)
316*
317 CALL ctrtri( 'L', diag, n1, a( 0 ), n, info )
318 IF( info.GT.0 )
319 $ RETURN
320 CALL ctrmm( 'R', 'L', 'N', diag, n2, n1, -cone,
321 $ a( 0 ),
322 $ n, a( n1 ), n )
323 CALL ctrtri( 'U', diag, n2, a( n ), n, info )
324 IF( info.GT.0 )
325 $ info = info + n1
326 IF( info.GT.0 )
327 $ RETURN
328 CALL ctrmm( 'L', 'U', 'C', diag, n2, n1, cone, a( n ),
329 $ n,
330 $ a( n1 ), n )
331*
332 ELSE
333*
334* SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
335* T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
336* T1 -> a(n2), T2 -> a(n1), S -> a(0)
337*
338 CALL ctrtri( 'L', diag, n1, a( n2 ), n, info )
339 IF( info.GT.0 )
340 $ RETURN
341 CALL ctrmm( 'L', 'L', 'C', diag, n1, n2, -cone,
342 $ a( n2 ),
343 $ n, a( 0 ), n )
344 CALL ctrtri( 'U', diag, n2, a( n1 ), n, info )
345 IF( info.GT.0 )
346 $ info = info + n1
347 IF( info.GT.0 )
348 $ RETURN
349 CALL ctrmm( 'R', 'U', 'N', diag, n1, n2, cone,
350 $ a( n1 ),
351 $ n, a( 0 ), n )
352*
353 END IF
354*
355 ELSE
356*
357* N is odd and TRANSR = 'C'
358*
359 IF( lower ) THEN
360*
361* SRPA for LOWER, TRANSPOSE and N is odd
362* T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
363*
364 CALL ctrtri( 'U', diag, n1, a( 0 ), n1, info )
365 IF( info.GT.0 )
366 $ RETURN
367 CALL ctrmm( 'L', 'U', 'N', diag, n1, n2, -cone,
368 $ a( 0 ),
369 $ n1, a( n1*n1 ), n1 )
370 CALL ctrtri( 'L', diag, n2, a( 1 ), n1, info )
371 IF( info.GT.0 )
372 $ info = info + n1
373 IF( info.GT.0 )
374 $ RETURN
375 CALL ctrmm( 'R', 'L', 'C', diag, n1, n2, cone, a( 1 ),
376 $ n1, a( n1*n1 ), n1 )
377*
378 ELSE
379*
380* SRPA for UPPER, TRANSPOSE and N is odd
381* T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
382*
383 CALL ctrtri( 'U', diag, n1, a( n2*n2 ), n2, info )
384 IF( info.GT.0 )
385 $ RETURN
386 CALL ctrmm( 'R', 'U', 'C', diag, n2, n1, -cone,
387 $ a( n2*n2 ), n2, a( 0 ), n2 )
388 CALL ctrtri( 'L', diag, n2, a( n1*n2 ), n2, info )
389 IF( info.GT.0 )
390 $ info = info + n1
391 IF( info.GT.0 )
392 $ RETURN
393 CALL ctrmm( 'L', 'L', 'N', diag, n2, n1, cone,
394 $ a( n1*n2 ), n2, a( 0 ), n2 )
395 END IF
396*
397 END IF
398*
399 ELSE
400*
401* N is even
402*
403 IF( normaltransr ) THEN
404*
405* N is even and TRANSR = 'N'
406*
407 IF( lower ) THEN
408*
409* SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
410* T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
411* T1 -> a(1), T2 -> a(0), S -> a(k+1)
412*
413 CALL ctrtri( 'L', diag, k, a( 1 ), n+1, info )
414 IF( info.GT.0 )
415 $ RETURN
416 CALL ctrmm( 'R', 'L', 'N', diag, k, k, -cone, a( 1 ),
417 $ n+1, a( k+1 ), n+1 )
418 CALL ctrtri( 'U', diag, k, a( 0 ), n+1, info )
419 IF( info.GT.0 )
420 $ info = info + k
421 IF( info.GT.0 )
422 $ RETURN
423 CALL ctrmm( 'L', 'U', 'C', diag, k, k, cone, a( 0 ),
424 $ n+1,
425 $ a( k+1 ), n+1 )
426*
427 ELSE
428*
429* SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
430* T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
431* T1 -> a(k+1), T2 -> a(k), S -> a(0)
432*
433 CALL ctrtri( 'L', diag, k, a( k+1 ), n+1, info )
434 IF( info.GT.0 )
435 $ RETURN
436 CALL ctrmm( 'L', 'L', 'C', diag, k, k, -cone,
437 $ a( k+1 ),
438 $ n+1, a( 0 ), n+1 )
439 CALL ctrtri( 'U', diag, k, a( k ), n+1, info )
440 IF( info.GT.0 )
441 $ info = info + k
442 IF( info.GT.0 )
443 $ RETURN
444 CALL ctrmm( 'R', 'U', 'N', diag, k, k, cone, a( k ),
445 $ n+1,
446 $ a( 0 ), n+1 )
447 END IF
448 ELSE
449*
450* N is even and TRANSR = 'C'
451*
452 IF( lower ) THEN
453*
454* SRPA for LOWER, TRANSPOSE and N is even (see paper)
455* T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
456* T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
457*
458 CALL ctrtri( 'U', diag, k, a( k ), k, info )
459 IF( info.GT.0 )
460 $ RETURN
461 CALL ctrmm( 'L', 'U', 'N', diag, k, k, -cone, a( k ),
462 $ k,
463 $ a( k*( k+1 ) ), k )
464 CALL ctrtri( 'L', diag, k, a( 0 ), k, info )
465 IF( info.GT.0 )
466 $ info = info + k
467 IF( info.GT.0 )
468 $ RETURN
469 CALL ctrmm( 'R', 'L', 'C', diag, k, k, cone, a( 0 ),
470 $ k,
471 $ a( k*( k+1 ) ), k )
472 ELSE
473*
474* SRPA for UPPER, TRANSPOSE and N is even (see paper)
475* T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
476* T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
477*
478 CALL ctrtri( 'U', diag, k, a( k*( k+1 ) ), k, info )
479 IF( info.GT.0 )
480 $ RETURN
481 CALL ctrmm( 'R', 'U', 'C', diag, k, k, -cone,
482 $ a( k*( k+1 ) ), k, a( 0 ), k )
483 CALL ctrtri( 'L', diag, k, a( k*k ), k, info )
484 IF( info.GT.0 )
485 $ info = info + k
486 IF( info.GT.0 )
487 $ RETURN
488 CALL ctrmm( 'L', 'L', 'N', diag, k, k, cone, a( k*k ),
489 $ k,
490 $ a( 0 ), k )
491 END IF
492 END IF
493 END IF
494*
495 RETURN
496*
497* End of CTFTRI
498*
499 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
ctftri
subroutine ctftri(transr, uplo, diag, n, a, info)
CTFTRI
Definition ctftri.f:219
ctrmm
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177
ctrtri
subroutine ctrtri(uplo, diag, n, a, lda, info)
CTRTRI
Definition ctrtri.f:107