LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 principal minor of order i
                is not positive, 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dsyrk
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
Definition dsyrk.f:169
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dpotrf
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:107
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181
Here is the call graph for this function:
Here is the caller graph for this function: