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

DPFTRI

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

Purpose: 
 DPFTRI computes the inverse of a (real) symmetric positive definite
 matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
 computed by DPFTRF.
Parameters
[in]TRANSR
          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The 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 DOUBLE PRECISION array, dimension ( N*(N+1)/2 )
          On entry, the symmetric 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 = 'T' then RFP is
          the 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 =
          'T'. 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 symmetric 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 Rectangular Full Packed (RFP) 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
  the 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
  the transpose of the last three columns of AP lower.
  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 = 'T'. RFP A in both UPLO cases is just the
  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 then consider Rectangular Full Packed (RFP) 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
  the 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
  the transpose of the last two columns of AP lower.
  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 = 'T'. RFP A in both UPLO cases is just the
  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 193 of file dpftri.f.

193 *
194 * -- LAPACK computational routine (version 3.4.0) --
195 * -- LAPACK is a software package provided by Univ. of Tennessee, --
196 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197 * November 2011
198 *
199 * .. Scalar Arguments ..
200  CHARACTER transr, uplo
201  INTEGER info, n
202 * .. Array Arguments ..
203  DOUBLE PRECISION a( 0: * )
204 * ..
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209  DOUBLE PRECISION one
210  parameter ( one = 1.0d+0 )
211 * ..
212 * .. Local Scalars ..
213  LOGICAL lower, nisodd, normaltransr
214  INTEGER n1, n2, k
215 * ..
216 * .. External Functions ..
217  LOGICAL lsame
218  EXTERNAL lsame
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL xerbla, dtftri, dlauum, dtrmm, dsyrk
222 * ..
223 * .. Intrinsic Functions ..
224  INTRINSIC mod
225 * ..
226 * .. Executable Statements ..
227 *
228 * Test the input parameters.
229 *
230  info = 0
231  normaltransr = lsame( transr, 'N' )
232  lower = lsame( uplo, 'L' )
233  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
234  info = -1
235  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
236  info = -2
237  ELSE IF( n.LT.0 ) THEN
238  info = -3
239  END IF
240  IF( info.NE.0 ) THEN
241  CALL xerbla( 'DPFTRI', -info )
242  RETURN
243  END IF
244 *
245 * Quick return if possible
246 *
247  IF( n.EQ.0 )
248  $ RETURN
249 *
250 * Invert the triangular Cholesky factor U or L.
251 *
252  CALL dtftri( transr, uplo, 'N', n, a, info )
253  IF( info.GT.0 )
254  $ RETURN
255 *
256 * If N is odd, set NISODD = .TRUE.
257 * If N is even, set K = N/2 and NISODD = .FALSE.
258 *
259  IF( mod( n, 2 ).EQ.0 ) THEN
260  k = n / 2
261  nisodd = .false.
262  ELSE
263  nisodd = .true.
264  END IF
265 *
266 * Set N1 and N2 depending on LOWER
267 *
268  IF( lower ) THEN
269  n2 = n / 2
270  n1 = n - n2
271  ELSE
272  n1 = n / 2
273  n2 = n - n1
274  END IF
275 *
276 * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
277 * inv(L)^C*inv(L). There are eight cases.
278 *
279  IF( nisodd ) THEN
280 *
281 * N is odd
282 *
283  IF( normaltransr ) THEN
284 *
285 * N is odd and TRANSR = 'N'
286 *
287  IF( lower ) THEN
288 *
289 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
290 * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
291 * T1 -> a(0), T2 -> a(n), S -> a(N1)
292 *
293  CALL dlauum( 'L', n1, a( 0 ), n, info )
294  CALL dsyrk( 'L', 'T', n1, n2, one, a( n1 ), n, one,
295  $ a( 0 ), n )
296  CALL dtrmm( 'L', 'U', 'N', 'N', n2, n1, one, a( n ), n,
297  $ a( n1 ), n )
298  CALL dlauum( 'U', n2, a( n ), n, info )
299 *
300  ELSE
301 *
302 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
303 * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
304 * T1 -> a(N2), T2 -> a(N1), S -> a(0)
305 *
306  CALL dlauum( 'L', n1, a( n2 ), n, info )
307  CALL dsyrk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
308  $ a( n2 ), n )
309  CALL dtrmm( 'R', 'U', 'T', 'N', n1, n2, one, a( n1 ), n,
310  $ a( 0 ), n )
311  CALL dlauum( 'U', n2, a( n1 ), n, info )
312 *
313  END IF
314 *
315  ELSE
316 *
317 * N is odd and TRANSR = 'T'
318 *
319  IF( lower ) THEN
320 *
321 * SRPA for LOWER, TRANSPOSE, and N is odd
322 * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
323 *
324  CALL dlauum( 'U', n1, a( 0 ), n1, info )
325  CALL dsyrk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1, one,
326  $ a( 0 ), n1 )
327  CALL dtrmm( 'R', 'L', 'N', 'N', n1, n2, one, a( 1 ), n1,
328  $ a( n1*n1 ), n1 )
329  CALL dlauum( 'L', n2, a( 1 ), n1, info )
330 *
331  ELSE
332 *
333 * SRPA for UPPER, TRANSPOSE, and N is odd
334 * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
335 *
336  CALL dlauum( 'U', n1, a( n2*n2 ), n2, info )
337  CALL dsyrk( 'U', 'T', n1, n2, one, a( 0 ), n2, one,
338  $ a( n2*n2 ), n2 )
339  CALL dtrmm( 'L', 'L', 'T', 'N', n2, n1, one, a( n1*n2 ),
340  $ n2, a( 0 ), n2 )
341  CALL dlauum( 'L', n2, a( n1*n2 ), n2, info )
342 *
343  END IF
344 *
345  END IF
346 *
347  ELSE
348 *
349 * N is even
350 *
351  IF( normaltransr ) THEN
352 *
353 * N is even and TRANSR = 'N'
354 *
355  IF( lower ) THEN
356 *
357 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
358 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
359 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
360 *
361  CALL dlauum( 'L', k, a( 1 ), n+1, info )
362  CALL dsyrk( 'L', 'T', k, k, one, a( k+1 ), n+1, one,
363  $ a( 1 ), n+1 )
364  CALL dtrmm( 'L', 'U', 'N', 'N', k, k, one, a( 0 ), n+1,
365  $ a( k+1 ), n+1 )
366  CALL dlauum( 'U', k, a( 0 ), n+1, info )
367 *
368  ELSE
369 *
370 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
371 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
372 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
373 *
374  CALL dlauum( 'L', k, a( k+1 ), n+1, info )
375  CALL dsyrk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
376  $ a( k+1 ), n+1 )
377  CALL dtrmm( 'R', 'U', 'T', 'N', k, k, one, a( k ), n+1,
378  $ a( 0 ), n+1 )
379  CALL dlauum( 'U', k, a( k ), n+1, info )
380 *
381  END IF
382 *
383  ELSE
384 *
385 * N is even and TRANSR = 'T'
386 *
387  IF( lower ) THEN
388 *
389 * SRPA for LOWER, TRANSPOSE, and N is even (see paper)
390 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
391 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
392 *
393  CALL dlauum( 'U', k, a( k ), k, info )
394  CALL dsyrk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k, one,
395  $ a( k ), k )
396  CALL dtrmm( 'R', 'L', 'N', 'N', k, k, one, a( 0 ), k,
397  $ a( k*( k+1 ) ), k )
398  CALL dlauum( 'L', k, a( 0 ), k, info )
399 *
400  ELSE
401 *
402 * SRPA for UPPER, TRANSPOSE, and N is even (see paper)
403 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0),
404 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
405 *
406  CALL dlauum( 'U', k, a( k*( k+1 ) ), k, info )
407  CALL dsyrk( 'U', 'T', k, k, one, a( 0 ), k, one,
408  $ a( k*( k+1 ) ), k )
409  CALL dtrmm( 'L', 'L', 'T', 'N', k, k, one, a( k*k ), k,
410  $ a( 0 ), k )
411  CALL dlauum( 'L', k, a( k*k ), k, info )
412 *
413  END IF
414 *
415  END IF
416 *
417  END IF
418 *
419  RETURN
420 *
421 * End of DPFTRI
422 *
dtrmm
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:179
dsyrk
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
Definition: dsyrk.f:171
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dtftri
subroutine dtftri(TRANSR, UPLO, DIAG, N, A, INFO)
DTFTRI
Definition: dtftri.f:203
dlauum
subroutine dlauum(UPLO, N, A, LDA, INFO)
DLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition: dlauum.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: