LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zpftri.f
Go to the documentation of this file.
1 *> \brief \b ZPFTRI
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZPFTRI + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpftri.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpftri.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpftri.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZPFTRI( TRANSR, UPLO, N, A, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER TRANSR, UPLO
25 * INTEGER INFO, N
26 * .. Array Arguments ..
27 * COMPLEX*16 A( 0: * )
28 * ..
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> ZPFTRI 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 ZPFTRF.
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*16 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 complex16OTHERcomputational
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 zpftri( 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*16 A( 0: * )
222 * ..
223 *
224 * =====================================================================
225 *
226 * .. Parameters ..
227  DOUBLE PRECISION ONE
228  COMPLEX*16 CONE
229  parameter( one = 1.d0, cone = ( 1.d0, 0.d0 ) )
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, ztftri, zlauum, ztrmm, zherk
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( 'ZPFTRI', -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 ztftri( 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 zlauum( 'L', n1, a( 0 ), n, info )
313  CALL zherk( 'L', 'C', n1, n2, one, a( n1 ), n, one,
314  $ a( 0 ), n )
315  CALL ztrmm( 'L', 'U', 'N', 'N', n2, n1, cone, a( n ), n,
316  $ a( n1 ), n )
317  CALL zlauum( '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 zlauum( 'L', n1, a( n2 ), n, info )
326  CALL zherk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
327  $ a( n2 ), n )
328  CALL ztrmm( 'R', 'U', 'C', 'N', n1, n2, cone, a( n1 ), n,
329  $ a( 0 ), n )
330  CALL zlauum( '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 zlauum( 'U', n1, a( 0 ), n1, info )
344  CALL zherk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1, one,
345  $ a( 0 ), n1 )
346  CALL ztrmm( 'R', 'L', 'N', 'N', n1, n2, cone, a( 1 ), n1,
347  $ a( n1*n1 ), n1 )
348  CALL zlauum( '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 zlauum( 'U', n1, a( n2*n2 ), n2, info )
356  CALL zherk( 'U', 'C', n1, n2, one, a( 0 ), n2, one,
357  $ a( n2*n2 ), n2 )
358  CALL ztrmm( 'L', 'L', 'C', 'N', n2, n1, cone, a( n1*n2 ),
359  $ n2, a( 0 ), n2 )
360  CALL zlauum( '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 zlauum( 'L', k, a( 1 ), n+1, info )
381  CALL zherk( 'L', 'C', k, k, one, a( k+1 ), n+1, one,
382  $ a( 1 ), n+1 )
383  CALL ztrmm( 'L', 'U', 'N', 'N', k, k, cone, a( 0 ), n+1,
384  $ a( k+1 ), n+1 )
385  CALL zlauum( '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 zlauum( 'L', k, a( k+1 ), n+1, info )
394  CALL zherk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
395  $ a( k+1 ), n+1 )
396  CALL ztrmm( 'R', 'U', 'C', 'N', k, k, cone, a( k ), n+1,
397  $ a( 0 ), n+1 )
398  CALL zlauum( '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 zlauum( 'U', k, a( k ), k, info )
413  CALL zherk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k, one,
414  $ a( k ), k )
415  CALL ztrmm( 'R', 'L', 'N', 'N', k, k, cone, a( 0 ), k,
416  $ a( k*( k+1 ) ), k )
417  CALL zlauum( '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 zlauum( 'U', k, a( k*( k+1 ) ), k, info )
426  CALL zherk( 'U', 'C', k, k, one, a( 0 ), k, one,
427  $ a( k*( k+1 ) ), k )
428  CALL ztrmm( 'L', 'L', 'C', 'N', k, k, cone, a( k*k ), k,
429  $ a( 0 ), k )
430  CALL zlauum( '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 ZPFTRI
441 *
442  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:173
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:177
subroutine zlauum(UPLO, N, A, LDA, INFO)
ZLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition: zlauum.f:102
subroutine ztftri(TRANSR, UPLO, DIAG, N, A, INFO)
ZTFTRI
Definition: ztftri.f:221
subroutine zpftri(TRANSR, UPLO, N, A, INFO)
ZPFTRI
Definition: zpftri.f:212