LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dpftri()

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.
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 190 of file dpftri.f.

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