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

◆ cpftri()

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.
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 cpftri.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 A( 0: * )
222* ..
223*
224* =====================================================================
225*
226* .. Parameters ..
227 REAL ONE
228 COMPLEX CONE
229 parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ) )
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, ctftri, clauum, ctrmm, cherk
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( 'CPFTRI', -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 ctftri( 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 clauum( 'L', n1, a( 0 ), n, info )
313 CALL cherk( 'L', 'C', n1, n2, one, a( n1 ), n, one,
314 $ a( 0 ), n )
315 CALL ctrmm( 'L', 'U', 'N', 'N', n2, n1, cone, a( n ), n,
316 $ a( n1 ), n )
317 CALL clauum( '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 clauum( 'L', n1, a( n2 ), n, info )
326 CALL cherk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
327 $ a( n2 ), n )
328 CALL ctrmm( 'R', 'U', 'C', 'N', n1, n2, cone, a( n1 ), n,
329 $ a( 0 ), n )
330 CALL clauum( '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 clauum( 'U', n1, a( 0 ), n1, info )
344 CALL cherk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1, one,
345 $ a( 0 ), n1 )
346 CALL ctrmm( 'R', 'L', 'N', 'N', n1, n2, cone, a( 1 ), n1,
347 $ a( n1*n1 ), n1 )
348 CALL clauum( '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 clauum( 'U', n1, a( n2*n2 ), n2, info )
356 CALL cherk( 'U', 'C', n1, n2, one, a( 0 ), n2, one,
357 $ a( n2*n2 ), n2 )
358 CALL ctrmm( 'L', 'L', 'C', 'N', n2, n1, cone, a( n1*n2 ),
359 $ n2, a( 0 ), n2 )
360 CALL clauum( '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 clauum( 'L', k, a( 1 ), n+1, info )
381 CALL cherk( 'L', 'C', k, k, one, a( k+1 ), n+1, one,
382 $ a( 1 ), n+1 )
383 CALL ctrmm( 'L', 'U', 'N', 'N', k, k, cone, a( 0 ), n+1,
384 $ a( k+1 ), n+1 )
385 CALL clauum( '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 clauum( 'L', k, a( k+1 ), n+1, info )
394 CALL cherk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
395 $ a( k+1 ), n+1 )
396 CALL ctrmm( 'R', 'U', 'C', 'N', k, k, cone, a( k ), n+1,
397 $ a( 0 ), n+1 )
398 CALL clauum( '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 clauum( 'U', k, a( k ), k, info )
413 CALL cherk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k, one,
414 $ a( k ), k )
415 CALL ctrmm( 'R', 'L', 'N', 'N', k, k, cone, a( 0 ), k,
416 $ a( k*( k+1 ) ), k )
417 CALL clauum( '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 clauum( 'U', k, a( k*( k+1 ) ), k, info )
426 CALL cherk( 'U', 'C', k, k, one, a( 0 ), k, one,
427 $ a( k*( k+1 ) ), k )
428 CALL ctrmm( 'L', 'L', 'C', 'N', k, k, cone, a( k*k ), k,
429 $ a( 0 ), k )
430 CALL clauum( '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 CPFTRI
441*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cherk
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
Definition cherk.f:173
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:102
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ctftri
subroutine ctftri(transr, uplo, diag, n, a, info)
CTFTRI
Definition ctftri.f:221
ctrmm
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: