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

◆ 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 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 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
ssyrk
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
spotrf
subroutine spotrf(uplo, n, a, lda, info)
SPOTRF
Definition spotrf.f:107
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
Here is the call graph for this function:
Here is the caller graph for this function: