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

DPFTRF

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

Purpose: 
 DPFTRF computes the Cholesky factorization of a real symmetric
 positive definite matrix A.

 The factorization has the form
    A = U**T * U,  if UPLO = 'U', or
    A = L  * L**T,  if UPLO = 'L',
 where U is an upper triangular matrix and L is lower triangular.

 This is the block version of the algorithm, calling Level 3 BLAS.
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 RFP A is stored;
          = 'L':  Lower triangle of RFP 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, if INFO = 0, the factor U or L from the Cholesky
          factorization RFP A = U**T*U or RFP A = L*L**T.
[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 leading minor of order i is not
                positive definite, and the factorization could not be
                completed.
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 200 of file dpftrf.f.

200 *
201 * -- LAPACK computational routine (version 3.4.0) --
202 * -- LAPACK is a software package provided by Univ. of Tennessee, --
203 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
204 * November 2011
205 *
206 * .. Scalar Arguments ..
207  CHARACTER transr, uplo
208  INTEGER n, info
209 * ..
210 * .. Array Arguments ..
211  DOUBLE PRECISION a( 0: * )
212 *
213 * =====================================================================
214 *
215 * .. Parameters ..
216  DOUBLE PRECISION one
217  parameter ( one = 1.0d+0 )
218 * ..
219 * .. Local Scalars ..
220  LOGICAL lower, nisodd, normaltransr
221  INTEGER n1, n2, k
222 * ..
223 * .. External Functions ..
224  LOGICAL lsame
225  EXTERNAL lsame
226 * ..
227 * .. External Subroutines ..
228  EXTERNAL xerbla, dsyrk, dpotrf, dtrsm
229 * ..
230 * .. Intrinsic Functions ..
231  INTRINSIC mod
232 * ..
233 * .. Executable Statements ..
234 *
235 * Test the input parameters.
236 *
237  info = 0
238  normaltransr = lsame( transr, 'N' )
239  lower = lsame( uplo, 'L' )
240  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
241  info = -1
242  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
243  info = -2
244  ELSE IF( n.LT.0 ) THEN
245  info = -3
246  END IF
247  IF( info.NE.0 ) THEN
248  CALL xerbla( 'DPFTRF', -info )
249  RETURN
250  END IF
251 *
252 * Quick return if possible
253 *
254  IF( n.EQ.0 )
255  $ RETURN
256 *
257 * If N is odd, set NISODD = .TRUE.
258 * If N is even, set K = N/2 and NISODD = .FALSE.
259 *
260  IF( mod( n, 2 ).EQ.0 ) THEN
261  k = n / 2
262  nisodd = .false.
263  ELSE
264  nisodd = .true.
265  END IF
266 *
267 * Set N1 and N2 depending on LOWER
268 *
269  IF( lower ) THEN
270  n2 = n / 2
271  n1 = n - n2
272  ELSE
273  n1 = n / 2
274  n2 = n - n1
275  END IF
276 *
277 * start execution: 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 dpotrf( 'L', n1, a( 0 ), n, info )
294  IF( info.GT.0 )
295  $ RETURN
296  CALL dtrsm( 'R', 'L', 'T', 'N', n2, n1, one, a( 0 ), n,
297  $ a( n1 ), n )
298  CALL dsyrk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
299  $ a( n ), n )
300  CALL dpotrf( 'U', n2, a( n ), n, info )
301  IF( info.GT.0 )
302  $ info = info + n1
303 *
304  ELSE
305 *
306 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
307 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
308 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
309 *
310  CALL dpotrf( 'L', n1, a( n2 ), n, info )
311  IF( info.GT.0 )
312  $ RETURN
313  CALL dtrsm( 'L', 'L', 'N', 'N', n1, n2, one, a( n2 ), n,
314  $ a( 0 ), n )
315  CALL dsyrk( 'U', 'T', n2, n1, -one, a( 0 ), n, one,
316  $ a( n1 ), n )
317  CALL dpotrf( 'U', n2, a( n1 ), n, info )
318  IF( info.GT.0 )
319  $ info = info + n1
320 *
321  END IF
322 *
323  ELSE
324 *
325 * N is odd and TRANSR = 'T'
326 *
327  IF( lower ) THEN
328 *
329 * SRPA for LOWER, TRANSPOSE and N is odd
330 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
331 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
332 *
333  CALL dpotrf( 'U', n1, a( 0 ), n1, info )
334  IF( info.GT.0 )
335  $ RETURN
336  CALL dtrsm( 'L', 'U', 'T', 'N', n1, n2, one, a( 0 ), n1,
337  $ a( n1*n1 ), n1 )
338  CALL dsyrk( 'L', 'T', n2, n1, -one, a( n1*n1 ), n1, one,
339  $ a( 1 ), n1 )
340  CALL dpotrf( 'L', n2, a( 1 ), n1, info )
341  IF( info.GT.0 )
342  $ info = info + n1
343 *
344  ELSE
345 *
346 * SRPA for UPPER, TRANSPOSE and N is odd
347 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
348 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
349 *
350  CALL dpotrf( 'U', n1, a( n2*n2 ), n2, info )
351  IF( info.GT.0 )
352  $ RETURN
353  CALL dtrsm( 'R', 'U', 'N', 'N', n2, n1, one, a( n2*n2 ),
354  $ n2, a( 0 ), n2 )
355  CALL dsyrk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
356  $ a( n1*n2 ), n2 )
357  CALL dpotrf( 'L', n2, a( n1*n2 ), n2, info )
358  IF( info.GT.0 )
359  $ info = info + n1
360 *
361  END IF
362 *
363  END IF
364 *
365  ELSE
366 *
367 * N is even
368 *
369  IF( normaltransr ) THEN
370 *
371 * N is even and TRANSR = 'N'
372 *
373  IF( lower ) THEN
374 *
375 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
376 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
377 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
378 *
379  CALL dpotrf( 'L', k, a( 1 ), n+1, info )
380  IF( info.GT.0 )
381  $ RETURN
382  CALL dtrsm( 'R', 'L', 'T', 'N', k, k, one, a( 1 ), n+1,
383  $ a( k+1 ), n+1 )
384  CALL dsyrk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
385  $ a( 0 ), n+1 )
386  CALL dpotrf( 'U', k, a( 0 ), n+1, info )
387  IF( info.GT.0 )
388  $ info = info + k
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 dpotrf( 'L', k, a( k+1 ), n+1, info )
397  IF( info.GT.0 )
398  $ RETURN
399  CALL dtrsm( 'L', 'L', 'N', 'N', k, k, one, a( k+1 ),
400  $ n+1, a( 0 ), n+1 )
401  CALL dsyrk( 'U', 'T', k, k, -one, a( 0 ), n+1, one,
402  $ a( k ), n+1 )
403  CALL dpotrf( 'U', k, a( k ), n+1, info )
404  IF( info.GT.0 )
405  $ info = info + k
406 *
407  END IF
408 *
409  ELSE
410 *
411 * N is even and TRANSR = 'T'
412 *
413  IF( lower ) THEN
414 *
415 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
416 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
417 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
418 *
419  CALL dpotrf( 'U', k, a( 0+k ), k, info )
420  IF( info.GT.0 )
421  $ RETURN
422  CALL dtrsm( 'L', 'U', 'T', 'N', k, k, one, a( k ), n1,
423  $ a( k*( k+1 ) ), k )
424  CALL dsyrk( 'L', 'T', k, k, -one, a( k*( k+1 ) ), k, one,
425  $ a( 0 ), k )
426  CALL dpotrf( 'L', k, a( 0 ), k, info )
427  IF( info.GT.0 )
428  $ info = info + k
429 *
430  ELSE
431 *
432 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
433 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
434 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
435 *
436  CALL dpotrf( 'U', k, a( k*( k+1 ) ), k, info )
437  IF( info.GT.0 )
438  $ RETURN
439  CALL dtrsm( 'R', 'U', 'N', 'N', k, k, one,
440  $ a( k*( k+1 ) ), k, a( 0 ), k )
441  CALL dsyrk( 'L', 'N', k, k, -one, a( 0 ), k, one,
442  $ a( k*k ), k )
443  CALL dpotrf( 'L', k, a( k*k ), k, info )
444  IF( info.GT.0 )
445  $ info = info + k
446 *
447  END IF
448 *
449  END IF
450 *
451  END IF
452 *
453  RETURN
454 *
455 * End of DPFTRF
456 *
dtrsm
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:183
dpotrf
subroutine dpotrf(UPLO, N, A, LDA, INFO)
DPOTRF
Definition: dpotrf.f:109
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
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: