LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011
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 214 of file zpftri.f.

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 *
ztftri
subroutine ztftri(TRANSR, UPLO, DIAG, N, A, INFO)
ZTFTRI
Definition: ztftri.f:223
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:104
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
ztrmm
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
zherk
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:175
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: