LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zpftri()

 subroutine zpftri ( character TRANSR, character UPLO, integer N, complex*16, dimension( 0: * ) A, integer INFO )

ZPFTRI

Purpose:
``` ZPFTRI computes the inverse of a complex Hermitian positive definite
matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
computed by ZPFTRF.```
Parameters
 [in] TRANSR ``` TRANSR is CHARACTER*1 = 'N': The Normal TRANSR of RFP A is stored; = 'C': The Conjugate-transpose TRANSR of RFP A is stored.``` [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in,out] A ``` A is COMPLEX*16 array, dimension ( N*(N+1)/2 ); On entry, the Hermitian matrix A in RFP format. RFP format is described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is the Conjugate-transpose of RFP A as defined when TRANSR = 'N'. The contents of RFP A are defined by UPLO as follows: If UPLO = 'U' the RFP A contains the nt elements of upper packed A. If UPLO = 'L' the RFP A contains the elements of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N is odd. See the Note below for more details. On exit, the Hermitian inverse of the original matrix, in the same storage format.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the (i,i) element of the factor U or L is zero, and the inverse could not be computed.```
Further Details:
```  We first consider Standard Packed Format when N is even.
We give an example where N = 6.

AP is Upper             AP is Lower

00 01 02 03 04 05       00
11 12 13 14 15       10 11
22 23 24 25       20 21 22
33 34 35       30 31 32 33
44 45       40 41 42 43 44
55       50 51 52 53 54 55

Let TRANSR = 'N'. RFP holds AP as follows:
For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
three columns of AP upper. The lower triangle A(4:6,0:2) consists of
conjugate-transpose of the first three columns of AP upper.
For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
three columns of AP lower. The upper triangle A(0:2,0:2) consists of
conjugate-transpose of the last three columns of AP lower.
To denote conjugate we place -- above the element. This covers the
case N even and TRANSR = 'N'.

RFP A                   RFP A

-- -- --
03 04 05                33 43 53
-- --
13 14 15                00 44 54
--
23 24 25                10 11 55

33 34 35                20 21 22
--
00 44 45                30 31 32
-- --
01 11 55                40 41 42
-- -- --
02 12 22                50 51 52

Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
transpose of RFP A above. One therefore gets:

RFP A                   RFP A

-- -- -- --                -- -- -- -- -- --
03 13 23 33 00 01 02    33 00 10 20 30 40 50
-- -- -- -- --                -- -- -- -- --
04 14 24 34 44 11 12    43 44 11 21 31 41 51
-- -- -- -- -- --                -- -- -- --
05 15 25 35 45 55 22    53 54 55 22 32 42 52

We next  consider Standard Packed Format when N is odd.
We give an example where N = 5.

AP is Upper                 AP is Lower

00 01 02 03 04              00
11 12 13 14              10 11
22 23 24              20 21 22
33 34              30 31 32 33
44              40 41 42 43 44

Let TRANSR = 'N'. RFP holds AP as follows:
For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
three columns of AP upper. The lower triangle A(3:4,0:1) consists of
conjugate-transpose of the first two   columns of AP upper.
For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
three columns of AP lower. The upper triangle A(0:1,1:2) consists of
conjugate-transpose of the last two   columns of AP lower.
To denote conjugate we place -- above the element. This covers the
case N odd  and TRANSR = 'N'.

RFP A                   RFP A

-- --
02 03 04                00 33 43
--
12 13 14                10 11 44

22 23 24                20 21 22
--
00 33 34                30 31 32
-- --
01 11 44                40 41 42

Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
transpose of RFP A above. One therefore gets:

RFP A                   RFP A

-- -- --                   -- -- -- -- -- --
02 12 22 00 01             00 10 20 30 40 50
-- -- -- --                   -- -- -- -- --
03 13 23 33 11             33 11 21 31 41 51
-- -- -- -- --                   -- -- -- --
04 14 24 34 44             43 44 22 32 42 52```

Definition at line 211 of file zpftri.f.

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*
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
Here is the call graph for this function:
Here is the caller graph for this function: