LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
101 *
102 *> \ingroup complex16OTHERcomputational
103 *
104 *> \par Further Details:
105 * =====================
106 *>
107 *> \verbatim
108 *>
109 *> We first consider Standard Packed Format when N is even.
110 *> We give an example where N = 6.
111 *>
112 *> AP is Upper AP is Lower
113 *>
114 *> 00 01 02 03 04 05 00
115 *> 11 12 13 14 15 10 11
116 *> 22 23 24 25 20 21 22
117 *> 33 34 35 30 31 32 33
118 *> 44 45 40 41 42 43 44
119 *> 55 50 51 52 53 54 55
120 *>
121 *>
122 *> Let TRANSR = 'N'. RFP holds AP as follows:
123 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
124 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
125 *> conjugate-transpose of the first three columns of AP upper.
126 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
127 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
128 *> conjugate-transpose of the last three columns of AP lower.
129 *> To denote conjugate we place -- above the element. This covers the
130 *> case N even and TRANSR = 'N'.
131 *>
132 *> RFP A RFP A
133 *>
134 *> -- -- --
135 *> 03 04 05 33 43 53
136 *> -- --
137 *> 13 14 15 00 44 54
138 *> --
139 *> 23 24 25 10 11 55
140 *>
141 *> 33 34 35 20 21 22
142 *> --
143 *> 00 44 45 30 31 32
144 *> -- --
145 *> 01 11 55 40 41 42
146 *> -- -- --
147 *> 02 12 22 50 51 52
148 *>
149 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
150 *> transpose of RFP A above. One therefore gets:
151 *>
152 *>
153 *> RFP A RFP A
154 *>
155 *> -- -- -- -- -- -- -- -- -- --
156 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
157 *> -- -- -- -- -- -- -- -- -- --
158 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
159 *> -- -- -- -- -- -- -- -- -- --
160 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
161 *>
162 *>
163 *> We next consider Standard Packed Format when N is odd.
164 *> We give an example where N = 5.
165 *>
166 *> AP is Upper AP is Lower
167 *>
168 *> 00 01 02 03 04 00
169 *> 11 12 13 14 10 11
170 *> 22 23 24 20 21 22
171 *> 33 34 30 31 32 33
172 *> 44 40 41 42 43 44
173 *>
174 *>
175 *> Let TRANSR = 'N'. RFP holds AP as follows:
176 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
177 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
178 *> conjugate-transpose of the first two columns of AP upper.
179 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
180 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
181 *> conjugate-transpose of the last two columns of AP lower.
182 *> To denote conjugate we place -- above the element. This covers the
183 *> case N odd and TRANSR = 'N'.
184 *>
185 *> RFP A RFP A
186 *>
187 *> -- --
188 *> 02 03 04 00 33 43
189 *> --
190 *> 12 13 14 10 11 44
191 *>
192 *> 22 23 24 20 21 22
193 *> --
194 *> 00 33 34 30 31 32
195 *> -- --
196 *> 01 11 44 40 41 42
197 *>
198 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
199 *> transpose of RFP A above. One therefore gets:
200 *>
201 *>
202 *> RFP A RFP A
203 *>
204 *> -- -- -- -- -- -- -- -- --
205 *> 02 12 22 00 01 00 10 20 30 40 50
206 *> -- -- -- -- -- -- -- -- --
207 *> 03 13 23 33 11 33 11 21 31 41 51
208 *> -- -- -- -- -- -- -- -- --
209 *> 04 14 24 34 44 43 44 22 32 42 52
210 *> \endverbatim
211 *>
212 * =====================================================================
213  SUBROUTINE zpftri( TRANSR, UPLO, N, A, INFO )
214 *
215 * -- LAPACK computational routine (version 3.4.0) --
216 * -- LAPACK is a software package provided by Univ. of Tennessee, --
217 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218 * November 2011
219 *
220 * .. Scalar Arguments ..
221  CHARACTER transr, uplo
222  INTEGER info, n
223 * .. Array Arguments ..
224  COMPLEX*16 a( 0: * )
225 * ..
226 *
227 * =====================================================================
228 *
229 * .. Parameters ..
230  DOUBLE PRECISION one
231  COMPLEX*16 cone
232  parameter( one = 1.d0, cone = ( 1.d0, 0.d0 ) )
233 * ..
234 * .. Local Scalars ..
235  LOGICAL lower, nisodd, normaltransr
236  INTEGER n1, n2, k
237 * ..
238 * .. External Functions ..
239  LOGICAL lsame
240  EXTERNAL lsame
241 * ..
242 * .. External Subroutines ..
243  EXTERNAL xerbla, ztftri, zlauum, ztrmm, zherk
244 * ..
245 * .. Intrinsic Functions ..
246  INTRINSIC mod
247 * ..
248 * .. Executable Statements ..
249 *
250 * Test the input parameters.
251 *
252  info = 0
253  normaltransr = lsame( transr, 'N' )
254  lower = lsame( uplo, 'L' )
255  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
256  info = -1
257  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
258  info = -2
259  ELSE IF( n.LT.0 ) THEN
260  info = -3
261  END IF
262  IF( info.NE.0 ) THEN
263  CALL xerbla( 'ZPFTRI', -info )
264  return
265  END IF
266 *
267 * Quick return if possible
268 *
269  IF( n.EQ.0 )
270  $ return
271 *
272 * Invert the triangular Cholesky factor U or L.
273 *
274  CALL ztftri( transr, uplo, 'N', n, a, info )
275  IF( info.GT.0 )
276  $ return
277 *
278 * If N is odd, set NISODD = .TRUE.
279 * If N is even, set K = N/2 and NISODD = .FALSE.
280 *
281  IF( mod( n, 2 ).EQ.0 ) THEN
282  k = n / 2
283  nisodd = .false.
284  ELSE
285  nisodd = .true.
286  END IF
287 *
288 * Set N1 and N2 depending on LOWER
289 *
290  IF( lower ) THEN
291  n2 = n / 2
292  n1 = n - n2
293  ELSE
294  n1 = n / 2
295  n2 = n - n1
296  END IF
297 *
298 * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
299 * inv(L)^C*inv(L). There are eight cases.
300 *
301  IF( nisodd ) THEN
302 *
303 * N is odd
304 *
305  IF( normaltransr ) THEN
306 *
307 * N is odd and TRANSR = 'N'
308 *
309  IF( lower ) THEN
310 *
311 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
312 * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
313 * T1 -> a(0), T2 -> a(n), S -> a(N1)
314 *
315  CALL zlauum( 'L', n1, a( 0 ), n, info )
316  CALL zherk( 'L', 'C', n1, n2, one, a( n1 ), n, one,
317  $ a( 0 ), n )
318  CALL ztrmm( 'L', 'U', 'N', 'N', n2, n1, cone, a( n ), n,
319  $ a( n1 ), n )
320  CALL zlauum( 'U', n2, a( n ), n, info )
321 *
322  ELSE
323 *
324 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
325 * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
326 * T1 -> a(N2), T2 -> a(N1), S -> a(0)
327 *
328  CALL zlauum( 'L', n1, a( n2 ), n, info )
329  CALL zherk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
330  $ a( n2 ), n )
331  CALL ztrmm( 'R', 'U', 'C', 'N', n1, n2, cone, a( n1 ), n,
332  $ a( 0 ), n )
333  CALL zlauum( 'U', n2, a( n1 ), n, info )
334 *
335  END IF
336 *
337  ELSE
338 *
339 * N is odd and TRANSR = 'C'
340 *
341  IF( lower ) THEN
342 *
343 * SRPA for LOWER, TRANSPOSE, and N is odd
344 * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
345 *
346  CALL zlauum( 'U', n1, a( 0 ), n1, info )
347  CALL zherk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1, one,
348  $ a( 0 ), n1 )
349  CALL ztrmm( 'R', 'L', 'N', 'N', n1, n2, cone, a( 1 ), n1,
350  $ a( n1*n1 ), n1 )
351  CALL zlauum( 'L', n2, a( 1 ), n1, info )
352 *
353  ELSE
354 *
355 * SRPA for UPPER, TRANSPOSE, and N is odd
356 * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
357 *
358  CALL zlauum( 'U', n1, a( n2*n2 ), n2, info )
359  CALL zherk( 'U', 'C', n1, n2, one, a( 0 ), n2, one,
360  $ a( n2*n2 ), n2 )
361  CALL ztrmm( 'L', 'L', 'C', 'N', n2, n1, cone, a( n1*n2 ),
362  $ n2, a( 0 ), n2 )
363  CALL zlauum( 'L', n2, a( n1*n2 ), n2, info )
364 *
365  END IF
366 *
367  END IF
368 *
369  ELSE
370 *
371 * N is even
372 *
373  IF( normaltransr ) THEN
374 *
375 * N is even and TRANSR = 'N'
376 *
377  IF( lower ) THEN
378 *
379 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
380 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
381 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
382 *
383  CALL zlauum( 'L', k, a( 1 ), n+1, info )
384  CALL zherk( 'L', 'C', k, k, one, a( k+1 ), n+1, one,
385  $ a( 1 ), n+1 )
386  CALL ztrmm( 'L', 'U', 'N', 'N', k, k, cone, a( 0 ), n+1,
387  $ a( k+1 ), n+1 )
388  CALL zlauum( 'U', k, a( 0 ), n+1, info )
389 *
390  ELSE
391 *
392 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
393 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
394 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
395 *
396  CALL zlauum( 'L', k, a( k+1 ), n+1, info )
397  CALL zherk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
398  $ a( k+1 ), n+1 )
399  CALL ztrmm( 'R', 'U', 'C', 'N', k, k, cone, a( k ), n+1,
400  $ a( 0 ), n+1 )
401  CALL zlauum( 'U', k, a( k ), n+1, info )
402 *
403  END IF
404 *
405  ELSE
406 *
407 * N is even and TRANSR = 'C'
408 *
409  IF( lower ) THEN
410 *
411 * SRPA for LOWER, TRANSPOSE, and N is even (see paper)
412 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
413 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
414 *
415  CALL zlauum( 'U', k, a( k ), k, info )
416  CALL zherk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k, one,
417  $ a( k ), k )
418  CALL ztrmm( 'R', 'L', 'N', 'N', k, k, cone, a( 0 ), k,
419  $ a( k*( k+1 ) ), k )
420  CALL zlauum( 'L', k, a( 0 ), k, info )
421 *
422  ELSE
423 *
424 * SRPA for UPPER, TRANSPOSE, and N is even (see paper)
425 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0),
426 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
427 *
428  CALL zlauum( 'U', k, a( k*( k+1 ) ), k, info )
429  CALL zherk( 'U', 'C', k, k, one, a( 0 ), k, one,
430  $ a( k*( k+1 ) ), k )
431  CALL ztrmm( 'L', 'L', 'C', 'N', k, k, cone, a( k*k ), k,
432  $ a( 0 ), k )
433  CALL zlauum( 'L', k, a( k*k ), k, info )
434 *
435  END IF
436 *
437  END IF
438 *
439  END IF
440 *
441  return
442 *
443 * End of ZPFTRI
444 *
445  END