LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ spftrf()

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.
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 spftrf.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  REAL A( 0: * )
209 *
210 * =====================================================================
211 *
212 * .. Parameters ..
213  REAL ONE
214  parameter( one = 1.0e+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, ssyrk, spotrf, strsm
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( 'SPFTRF', -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 spotrf( 'L', n1, a( 0 ), n, info )
291  IF( info.GT.0 )
292  $ RETURN
293  CALL strsm( 'R', 'L', 'T', 'N', n2, n1, one, a( 0 ), n,
294  $ a( n1 ), n )
295  CALL ssyrk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
296  $ a( n ), n )
297  CALL spotrf( '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 spotrf( 'L', n1, a( n2 ), n, info )
308  IF( info.GT.0 )
309  $ RETURN
310  CALL strsm( 'L', 'L', 'N', 'N', n1, n2, one, a( n2 ), n,
311  $ a( 0 ), n )
312  CALL ssyrk( 'U', 'T', n2, n1, -one, a( 0 ), n, one,
313  $ a( n1 ), n )
314  CALL spotrf( '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 spotrf( 'U', n1, a( 0 ), n1, info )
331  IF( info.GT.0 )
332  $ RETURN
333  CALL strsm( 'L', 'U', 'T', 'N', n1, n2, one, a( 0 ), n1,
334  $ a( n1*n1 ), n1 )
335  CALL ssyrk( 'L', 'T', n2, n1, -one, a( n1*n1 ), n1, one,
336  $ a( 1 ), n1 )
337  CALL spotrf( '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 spotrf( 'U', n1, a( n2*n2 ), n2, info )
348  IF( info.GT.0 )
349  $ RETURN
350  CALL strsm( 'R', 'U', 'N', 'N', n2, n1, one, a( n2*n2 ),
351  $ n2, a( 0 ), n2 )
352  CALL ssyrk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
353  $ a( n1*n2 ), n2 )
354  CALL spotrf( '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 spotrf( 'L', k, a( 1 ), n+1, info )
377  IF( info.GT.0 )
378  $ RETURN
379  CALL strsm( 'R', 'L', 'T', 'N', k, k, one, a( 1 ), n+1,
380  $ a( k+1 ), n+1 )
381  CALL ssyrk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
382  $ a( 0 ), n+1 )
383  CALL spotrf( '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 spotrf( 'L', k, a( k+1 ), n+1, info )
394  IF( info.GT.0 )
395  $ RETURN
396  CALL strsm( 'L', 'L', 'N', 'N', k, k, one, a( k+1 ),
397  $ n+1, a( 0 ), n+1 )
398  CALL ssyrk( 'U', 'T', k, k, -one, a( 0 ), n+1, one,
399  $ a( k ), n+1 )
400  CALL spotrf( '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 spotrf( 'U', k, a( 0+k ), k, info )
417  IF( info.GT.0 )
418  $ RETURN
419  CALL strsm( 'L', 'U', 'T', 'N', k, k, one, a( k ), n1,
420  $ a( k*( k+1 ) ), k )
421  CALL ssyrk( 'L', 'T', k, k, -one, a( k*( k+1 ) ), k, one,
422  $ a( 0 ), k )
423  CALL spotrf( '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 spotrf( 'U', k, a( k*( k+1 ) ), k, info )
434  IF( info.GT.0 )
435  $ RETURN
436  CALL strsm( 'R', 'U', 'N', 'N', k, k, one,
437  $ a( k*( k+1 ) ), k, a( 0 ), k )
438  CALL ssyrk( 'L', 'N', k, k, -one, a( 0 ), k, one,
439  $ a( k*k ), k )
440  CALL spotrf( '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 SPFTRF
453 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine spotrf(UPLO, N, A, LDA, INFO)
SPOTRF
Definition: spotrf.f:107
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:181
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169
Here is the call graph for this function:
Here is the caller graph for this function: