1*> \brief \b SPFTRF
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SPFTRF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spftrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spftrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spftrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SPFTRF( TRANSR, UPLO, N, A, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          TRANSR, UPLO
25*       INTEGER            N, INFO
26*       ..
27*       .. Array Arguments ..
28*       REAL               A( 0: * )
29*
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> SPFTRF computes the Cholesky factorization of a real symmetric
37*> positive definite matrix A.
38*>
39*> The factorization has the form
40*>    A = U**T * U,  if UPLO = 'U', or
41*>    A = L  * L**T,  if UPLO = 'L',
42*> where U is an upper triangular matrix and L is lower triangular.
43*>
44*> This is the block version of the algorithm, calling Level 3 BLAS.
45*> \endverbatim
46*
47*  Arguments:
48*  ==========
49*
50*> \param[in] TRANSR
51*> \verbatim
52*>          TRANSR is CHARACTER*1
53*>          = 'N':  The Normal TRANSR of RFP A is stored;
54*>          = 'T':  The Transpose TRANSR of RFP A is stored.
55*> \endverbatim
56*>
57*> \param[in] UPLO
58*> \verbatim
59*>          UPLO is CHARACTER*1
60*>          = 'U':  Upper triangle of RFP A is stored;
61*>          = 'L':  Lower triangle of RFP A is stored.
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*>          N is INTEGER
67*>          The order of the matrix A.  N >= 0.
68*> \endverbatim
69*>
70*> \param[in,out] A
71*> \verbatim
72*>          A is REAL array, dimension ( N*(N+1)/2 );
73*>          On entry, the symmetric matrix A in RFP format. RFP format is
74*>          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
75*>          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
76*>          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
77*>          the transpose of RFP A as defined when
78*>          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
79*>          follows: If UPLO = 'U' the RFP A contains the NT elements of
80*>          upper packed A. If UPLO = 'L' the RFP A contains the elements
81*>          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
82*>          'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
83*>          is odd. See the Note below for more details.
84*>
85*>          On exit, if INFO = 0, the factor U or L from the Cholesky
86*>          factorization RFP A = U**T*U or RFP A = L*L**T.
87*> \endverbatim
88*>
89*> \param[out] INFO
90*> \verbatim
91*>          INFO is INTEGER
92*>          = 0:  successful exit
93*>          < 0:  if INFO = -i, the i-th argument had an illegal value
94*>          > 0:  if INFO = i, the leading minor of order i is not
95*>                positive definite, and the factorization could not be
96*>                completed.
97*> \endverbatim
98*
99*  Authors:
100*  ========
101*
102*> \author Univ. of Tennessee
103*> \author Univ. of California Berkeley
104*> \author Univ. of Colorado Denver
105*> \author NAG Ltd.
106*
107*> \date November 2011
108*
109*> \ingroup realOTHERcomputational
110*
111*> \par Further Details:
112*  =====================
113*>
114*> \verbatim
115*>
116*>  We first consider Rectangular Full Packed (RFP) Format when N is
117*>  even. We give an example where N = 6.
118*>
119*>      AP is Upper             AP is Lower
120*>
121*>   00 01 02 03 04 05       00
122*>      11 12 13 14 15       10 11
123*>         22 23 24 25       20 21 22
124*>            33 34 35       30 31 32 33
125*>               44 45       40 41 42 43 44
126*>                  55       50 51 52 53 54 55
127*>
128*>
129*>  Let TRANSR = 'N'. RFP holds AP as follows:
130*>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
131*>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
132*>  the transpose of the first three columns of AP upper.
133*>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
134*>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
135*>  the transpose of the last three columns of AP lower.
136*>  This covers the case N even and TRANSR = 'N'.
137*>
138*>         RFP A                   RFP A
139*>
140*>        03 04 05                33 43 53
141*>        13 14 15                00 44 54
142*>        23 24 25                10 11 55
143*>        33 34 35                20 21 22
144*>        00 44 45                30 31 32
145*>        01 11 55                40 41 42
146*>        02 12 22                50 51 52
147*>
148*>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
149*>  transpose of RFP A above. One therefore gets:
150*>
151*>
152*>           RFP A                   RFP A
153*>
154*>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
155*>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
156*>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
157*>
158*>
159*>  We then consider Rectangular Full Packed (RFP) Format when N is
160*>  odd. We give an example where N = 5.
161*>
162*>     AP is Upper                 AP is Lower
163*>
164*>   00 01 02 03 04              00
165*>      11 12 13 14              10 11
166*>         22 23 24              20 21 22
167*>            33 34              30 31 32 33
168*>               44              40 41 42 43 44
169*>
170*>
171*>  Let TRANSR = 'N'. RFP holds AP as follows:
172*>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
173*>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
174*>  the transpose of the first two columns of AP upper.
175*>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
176*>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
177*>  the transpose of the last two columns of AP lower.
178*>  This covers the case N odd and TRANSR = 'N'.
179*>
180*>         RFP A                   RFP A
181*>
182*>        02 03 04                00 33 43
183*>        12 13 14                10 11 44
184*>        22 23 24                20 21 22
185*>        00 33 34                30 31 32
186*>        01 11 44                40 41 42
187*>
188*>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
189*>  transpose of RFP A above. One therefore gets:
190*>
191*>           RFP A                   RFP A
192*>
193*>     02 12 22 00 01             00 10 20 30 40 50
194*>     03 13 23 33 11             33 11 21 31 41 51
195*>     04 14 24 34 44             43 44 22 32 42 52
196*> \endverbatim
197*>
198*  =====================================================================
199      SUBROUTINE SPFTRF( TRANSR, UPLO, N, A, INFO )
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*
457      END
458