LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dpftrf.f
Go to the documentation of this file.
1*> \brief \b DPFTRF
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DPFTRF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpftrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpftrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpftrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DPFTRF( TRANSR, UPLO, N, A, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER TRANSR, UPLO
25* INTEGER N, INFO
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION A( 0: * )
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> DPFTRF 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 DOUBLE PRECISION 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 principal minor of order i
95*> is not positive, 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 pftrf
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 dpftrf( 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 DOUBLE PRECISION A( 0: * )
209*
210* =====================================================================
211*
212* .. Parameters ..
213 DOUBLE PRECISION ONE
214 parameter( one = 1.0d+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, dsyrk, dpotrf, dtrsm
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( 'DPFTRF', -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 dpotrf( 'L', n1, a( 0 ), n, info )
291 IF( info.GT.0 )
292 $ RETURN
293 CALL dtrsm( 'R', 'L', 'T', 'N', n2, n1, one, a( 0 ), n,
294 $ a( n1 ), n )
295 CALL dsyrk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
296 $ a( n ), n )
297 CALL dpotrf( '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 dpotrf( 'L', n1, a( n2 ), n, info )
308 IF( info.GT.0 )
309 $ RETURN
310 CALL dtrsm( 'L', 'L', 'N', 'N', n1, n2, one, a( n2 ), n,
311 $ a( 0 ), n )
312 CALL dsyrk( 'U', 'T', n2, n1, -one, a( 0 ), n, one,
313 $ a( n1 ), n )
314 CALL dpotrf( '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 dpotrf( 'U', n1, a( 0 ), n1, info )
331 IF( info.GT.0 )
332 $ RETURN
333 CALL dtrsm( 'L', 'U', 'T', 'N', n1, n2, one, a( 0 ), n1,
334 $ a( n1*n1 ), n1 )
335 CALL dsyrk( 'L', 'T', n2, n1, -one, a( n1*n1 ), n1, one,
336 $ a( 1 ), n1 )
337 CALL dpotrf( '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 dpotrf( 'U', n1, a( n2*n2 ), n2, info )
348 IF( info.GT.0 )
349 $ RETURN
350 CALL dtrsm( 'R', 'U', 'N', 'N', n2, n1, one, a( n2*n2 ),
351 $ n2, a( 0 ), n2 )
352 CALL dsyrk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
353 $ a( n1*n2 ), n2 )
354 CALL dpotrf( '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 dpotrf( 'L', k, a( 1 ), n+1, info )
377 IF( info.GT.0 )
378 $ RETURN
379 CALL dtrsm( 'R', 'L', 'T', 'N', k, k, one, a( 1 ), n+1,
380 $ a( k+1 ), n+1 )
381 CALL dsyrk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
382 $ a( 0 ), n+1 )
383 CALL dpotrf( '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 dpotrf( 'L', k, a( k+1 ), n+1, info )
394 IF( info.GT.0 )
395 $ RETURN
396 CALL dtrsm( 'L', 'L', 'N', 'N', k, k, one, a( k+1 ),
397 $ n+1, a( 0 ), n+1 )
398 CALL dsyrk( 'U', 'T', k, k, -one, a( 0 ), n+1, one,
399 $ a( k ), n+1 )
400 CALL dpotrf( '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 dpotrf( 'U', k, a( 0+k ), k, info )
417 IF( info.GT.0 )
418 $ RETURN
419 CALL dtrsm( 'L', 'U', 'T', 'N', k, k, one, a( k ), n1,
420 $ a( k*( k+1 ) ), k )
421 CALL dsyrk( 'L', 'T', k, k, -one, a( k*( k+1 ) ), k, one,
422 $ a( 0 ), k )
423 CALL dpotrf( '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 dpotrf( 'U', k, a( k*( k+1 ) ), k, info )
434 IF( info.GT.0 )
435 $ RETURN
436 CALL dtrsm( 'R', 'U', 'N', 'N', k, k, one,
437 $ a( k*( k+1 ) ), k, a( 0 ), k )
438 CALL dsyrk( 'L', 'N', k, k, -one, a( 0 ), k, one,
439 $ a( k*k ), k )
440 CALL dpotrf( '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 DPFTRF
453*
454 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
Definition dsyrk.f:169
subroutine dpftrf(transr, uplo, n, a, info)
DPFTRF
Definition dpftrf.f:198
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:107
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181