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

SPFTRF

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

Purpose: 
 SPFTRF 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 REAL 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 spftrf.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  REAL a( 0: * )
212 *
213 * =====================================================================
214 *
215 * .. Parameters ..
216  REAL one
217  parameter ( one = 1.0e+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, ssyrk, spotrf, strsm
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( 'SPFTRF', -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 spotrf( 'L', n1, a( 0 ), n, info )
294  IF( info.GT.0 )
295  $ RETURN
296  CALL strsm( 'R', 'L', 'T', 'N', n2, n1, one, a( 0 ), n,
297  $ a( n1 ), n )
298  CALL ssyrk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
299  $ a( n ), n )
300  CALL spotrf( '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 spotrf( 'L', n1, a( n2 ), n, info )
311  IF( info.GT.0 )
312  $ RETURN
313  CALL strsm( 'L', 'L', 'N', 'N', n1, n2, one, a( n2 ), n,
314  $ a( 0 ), n )
315  CALL ssyrk( 'U', 'T', n2, n1, -one, a( 0 ), n, one,
316  $ a( n1 ), n )
317  CALL spotrf( '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 spotrf( 'U', n1, a( 0 ), n1, info )
334  IF( info.GT.0 )
335  $ RETURN
336  CALL strsm( 'L', 'U', 'T', 'N', n1, n2, one, a( 0 ), n1,
337  $ a( n1*n1 ), n1 )
338  CALL ssyrk( 'L', 'T', n2, n1, -one, a( n1*n1 ), n1, one,
339  $ a( 1 ), n1 )
340  CALL spotrf( '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 spotrf( 'U', n1, a( n2*n2 ), n2, info )
351  IF( info.GT.0 )
352  $ RETURN
353  CALL strsm( 'R', 'U', 'N', 'N', n2, n1, one, a( n2*n2 ),
354  $ n2, a( 0 ), n2 )
355  CALL ssyrk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
356  $ a( n1*n2 ), n2 )
357  CALL spotrf( '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 spotrf( 'L', k, a( 1 ), n+1, info )
380  IF( info.GT.0 )
381  $ RETURN
382  CALL strsm( 'R', 'L', 'T', 'N', k, k, one, a( 1 ), n+1,
383  $ a( k+1 ), n+1 )
384  CALL ssyrk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
385  $ a( 0 ), n+1 )
386  CALL spotrf( '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 spotrf( 'L', k, a( k+1 ), n+1, info )
397  IF( info.GT.0 )
398  $ RETURN
399  CALL strsm( 'L', 'L', 'N', 'N', k, k, one, a( k+1 ),
400  $ n+1, a( 0 ), n+1 )
401  CALL ssyrk( 'U', 'T', k, k, -one, a( 0 ), n+1, one,
402  $ a( k ), n+1 )
403  CALL spotrf( '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 spotrf( 'U', k, a( 0+k ), k, info )
420  IF( info.GT.0 )
421  $ RETURN
422  CALL strsm( 'L', 'U', 'T', 'N', k, k, one, a( k ), n1,
423  $ a( k*( k+1 ) ), k )
424  CALL ssyrk( 'L', 'T', k, k, -one, a( k*( k+1 ) ), k, one,
425  $ a( 0 ), k )
426  CALL spotrf( '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 spotrf( 'U', k, a( k*( k+1 ) ), k, info )
437  IF( info.GT.0 )
438  $ RETURN
439  CALL strsm( 'R', 'U', 'N', 'N', k, k, one,
440  $ a( k*( k+1 ) ), k, a( 0 ), k )
441  CALL ssyrk( 'L', 'N', k, k, -one, a( 0 ), k, one,
442  $ a( k*k ), k )
443  CALL spotrf( '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 SPFTRF
456 *
ssyrk
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:171
strsm
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
spotrf
subroutine spotrf(UPLO, N, A, LDA, INFO)
SPOTRF
Definition: spotrf.f:109
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: