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