LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dpftrf()

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.
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 197 of file dpftrf.f.

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