LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zpftri()

subroutine zpftri ( character  transr,
character  uplo,
integer  n,
complex*16, dimension( 0: * )  a,
integer  info 
)

ZPFTRI

Download ZPFTRI + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zherk
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
Definition zherk.f:173
zlauum
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ztftri
subroutine ztftri(transr, uplo, diag, n, a, info)
ZTFTRI
Definition ztftri.f:221
ztrmm
subroutine ztrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRMM
Definition ztrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: