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