LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
spftrf.f
Go to the documentation of this file.
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 *> \ingroup realOTHERcomputational
108 *
109 *> \par Further Details:
110 * =====================
111 *>
112 *> \verbatim
113 *>
114 *> We first consider Rectangular Full Packed (RFP) Format when N is
115 *> even. We give an example where N = 6.
116 *>
117 *> AP is Upper AP is Lower
118 *>
119 *> 00 01 02 03 04 05 00
120 *> 11 12 13 14 15 10 11
121 *> 22 23 24 25 20 21 22
122 *> 33 34 35 30 31 32 33
123 *> 44 45 40 41 42 43 44
124 *> 55 50 51 52 53 54 55
125 *>
126 *>
127 *> Let TRANSR = 'N'. RFP holds AP as follows:
128 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
129 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
130 *> the transpose of the first three columns of AP upper.
131 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
132 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
133 *> the transpose of the last three columns of AP lower.
134 *> This covers the case N even and TRANSR = 'N'.
135 *>
136 *> RFP A RFP A
137 *>
138 *> 03 04 05 33 43 53
139 *> 13 14 15 00 44 54
140 *> 23 24 25 10 11 55
141 *> 33 34 35 20 21 22
142 *> 00 44 45 30 31 32
143 *> 01 11 55 40 41 42
144 *> 02 12 22 50 51 52
145 *>
146 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
147 *> transpose of RFP A above. One therefore gets:
148 *>
149 *>
150 *> RFP A RFP A
151 *>
152 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
153 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
154 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
155 *>
156 *>
157 *> We then consider Rectangular Full Packed (RFP) Format when N is
158 *> odd. We give an example where N = 5.
159 *>
160 *> AP is Upper AP is Lower
161 *>
162 *> 00 01 02 03 04 00
163 *> 11 12 13 14 10 11
164 *> 22 23 24 20 21 22
165 *> 33 34 30 31 32 33
166 *> 44 40 41 42 43 44
167 *>
168 *>
169 *> Let TRANSR = 'N'. RFP holds AP as follows:
170 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
171 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
172 *> the transpose of the first two columns of AP upper.
173 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
174 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
175 *> the transpose of the last two columns of AP lower.
176 *> This covers the case N odd and TRANSR = 'N'.
177 *>
178 *> RFP A RFP A
179 *>
180 *> 02 03 04 00 33 43
181 *> 12 13 14 10 11 44
182 *> 22 23 24 20 21 22
183 *> 00 33 34 30 31 32
184 *> 01 11 44 40 41 42
185 *>
186 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
187 *> transpose of RFP A above. One therefore gets:
188 *>
189 *> RFP A RFP A
190 *>
191 *> 02 12 22 00 01 00 10 20 30 40 50
192 *> 03 13 23 33 11 33 11 21 31 41 51
193 *> 04 14 24 34 44 43 44 22 32 42 52
194 *> \endverbatim
195 *>
196 * =====================================================================
197  SUBROUTINE spftrf( TRANSR, UPLO, N, A, INFO )
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 *
454  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine spftrf(TRANSR, UPLO, N, A, INFO)
SPFTRF
Definition: spftrf.f:198
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