LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cpftri ( character  TRANSR,
character  UPLO,
integer  N,
complex, dimension( 0: * )  A,
integer  INFO 
)

CPFTRI

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

Purpose: 
 CPFTRI 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 CPFTRF.
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 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 cpftri.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 a( 0: * )
225 * ..
226 *
227 * =====================================================================
228 *
229 * .. Parameters ..
230  REAL one
231  COMPLEX cone
232  parameter ( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ) )
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, ctftri, clauum, ctrmm, cherk
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( 'CPFTRI', -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 ctftri( 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 clauum( 'L', n1, a( 0 ), n, info )
316  CALL cherk( 'L', 'C', n1, n2, one, a( n1 ), n, one,
317  $ a( 0 ), n )
318  CALL ctrmm( 'L', 'U', 'N', 'N', n2, n1, cone, a( n ), n,
319  $ a( n1 ), n )
320  CALL clauum( '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 clauum( 'L', n1, a( n2 ), n, info )
329  CALL cherk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
330  $ a( n2 ), n )
331  CALL ctrmm( 'R', 'U', 'C', 'N', n1, n2, cone, a( n1 ), n,
332  $ a( 0 ), n )
333  CALL clauum( '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 clauum( 'U', n1, a( 0 ), n1, info )
347  CALL cherk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1, one,
348  $ a( 0 ), n1 )
349  CALL ctrmm( 'R', 'L', 'N', 'N', n1, n2, cone, a( 1 ), n1,
350  $ a( n1*n1 ), n1 )
351  CALL clauum( '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 clauum( 'U', n1, a( n2*n2 ), n2, info )
359  CALL cherk( 'U', 'C', n1, n2, one, a( 0 ), n2, one,
360  $ a( n2*n2 ), n2 )
361  CALL ctrmm( 'L', 'L', 'C', 'N', n2, n1, cone, a( n1*n2 ),
362  $ n2, a( 0 ), n2 )
363  CALL clauum( '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 clauum( 'L', k, a( 1 ), n+1, info )
384  CALL cherk( 'L', 'C', k, k, one, a( k+1 ), n+1, one,
385  $ a( 1 ), n+1 )
386  CALL ctrmm( 'L', 'U', 'N', 'N', k, k, cone, a( 0 ), n+1,
387  $ a( k+1 ), n+1 )
388  CALL clauum( '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 clauum( 'L', k, a( k+1 ), n+1, info )
397  CALL cherk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
398  $ a( k+1 ), n+1 )
399  CALL ctrmm( 'R', 'U', 'C', 'N', k, k, cone, a( k ), n+1,
400  $ a( 0 ), n+1 )
401  CALL clauum( '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 clauum( 'U', k, a( k ), k, info )
416  CALL cherk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k, one,
417  $ a( k ), k )
418  CALL ctrmm( 'R', 'L', 'N', 'N', k, k, cone, a( 0 ), k,
419  $ a( k*( k+1 ) ), k )
420  CALL clauum( '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 clauum( 'U', k, a( k*( k+1 ) ), k, info )
429  CALL cherk( 'U', 'C', k, k, one, a( 0 ), k, one,
430  $ a( k*( k+1 ) ), k )
431  CALL ctrmm( 'L', 'L', 'C', 'N', k, k, cone, a( k*k ), k,
432  $ a( 0 ), k )
433  CALL clauum( '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 CPFTRI
444 *
ctftri
subroutine ctftri(TRANSR, UPLO, DIAG, N, A, INFO)
CTFTRI
Definition: ctftri.f:223
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
ctrmm
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179
clauum
subroutine clauum(UPLO, N, A, LDA, INFO)
CLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition: clauum.f:104
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: