LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dpftri()

 subroutine dpftri ( character TRANSR, character UPLO, integer N, double precision, dimension( 0: * ) A, integer INFO )

DPFTRI

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.```
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: