Annotation of rpl/lapack/lapack/dpotrf2.f, revision 1.6
1.1 bertrand 1: *> \brief \b DPOTRF2
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.3 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
8: * Definition:
9: * ===========
10: *
11: * RECURSIVE SUBROUTINE DPOTRF2( UPLO, N, A, LDA, INFO )
1.3 bertrand 12: *
1.1 bertrand 13: * .. Scalar Arguments ..
14: * CHARACTER UPLO
15: * INTEGER INFO, LDA, N
16: * ..
17: * .. Array Arguments ..
18: * REAL A( LDA, * )
19: * ..
1.3 bertrand 20: *
1.1 bertrand 21: *
22: *> \par Purpose:
23: * =============
24: *>
25: *> \verbatim
26: *>
27: *> DPOTRF2 computes the Cholesky factorization of a real symmetric
28: *> positive definite matrix A using the recursive algorithm.
29: *>
30: *> The factorization has the form
31: *> A = U**T * U, if UPLO = 'U', or
32: *> A = L * L**T, if UPLO = 'L',
33: *> where U is an upper triangular matrix and L is lower triangular.
34: *>
35: *> This is the recursive version of the algorithm. It divides
36: *> the matrix into four submatrices:
37: *>
38: *> [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2
39: *> A = [ -----|----- ] with n1 = n/2
40: *> [ A21 | A22 ] n2 = n-n1
41: *>
42: *> The subroutine calls itself to factor A11. Update and scale A21
43: *> or A12, update A22 then calls itself to factor A22.
1.3 bertrand 44: *>
1.1 bertrand 45: *> \endverbatim
46: *
47: * Arguments:
48: * ==========
49: *
50: *> \param[in] UPLO
51: *> \verbatim
52: *> UPLO is CHARACTER*1
53: *> = 'U': Upper triangle of A is stored;
54: *> = 'L': Lower triangle of A is stored.
55: *> \endverbatim
56: *>
57: *> \param[in] N
58: *> \verbatim
59: *> N is INTEGER
60: *> The order of the matrix A. N >= 0.
61: *> \endverbatim
62: *>
63: *> \param[in,out] A
64: *> \verbatim
65: *> A is DOUBLE PRECISION array, dimension (LDA,N)
66: *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
67: *> N-by-N upper triangular part of A contains the upper
68: *> triangular part of the matrix A, and the strictly lower
69: *> triangular part of A is not referenced. If UPLO = 'L', the
70: *> leading N-by-N lower triangular part of A contains the lower
71: *> triangular part of the matrix A, and the strictly upper
72: *> triangular part of A is not referenced.
73: *>
74: *> On exit, if INFO = 0, the factor U or L from the Cholesky
75: *> factorization A = U**T*U or A = L*L**T.
76: *> \endverbatim
77: *>
78: *> \param[in] LDA
79: *> \verbatim
80: *> LDA is INTEGER
81: *> The leading dimension of the array A. LDA >= max(1,N).
82: *> \endverbatim
83: *>
84: *> \param[out] INFO
85: *> \verbatim
86: *> INFO is INTEGER
87: *> = 0: successful exit
88: *> < 0: if INFO = -i, the i-th argument had an illegal value
89: *> > 0: if INFO = i, the leading minor of order i is not
90: *> positive definite, and the factorization could not be
91: *> completed.
92: *> \endverbatim
93: *
94: * Authors:
95: * ========
96: *
1.3 bertrand 97: *> \author Univ. of Tennessee
98: *> \author Univ. of California Berkeley
99: *> \author Univ. of Colorado Denver
100: *> \author NAG Ltd.
1.1 bertrand 101: *
102: *> \ingroup doublePOcomputational
103: *
104: * =====================================================================
105: RECURSIVE SUBROUTINE DPOTRF2( UPLO, N, A, LDA, INFO )
106: *
1.6 ! bertrand 107: * -- LAPACK computational routine --
1.1 bertrand 108: * -- LAPACK is a software package provided by Univ. of Tennessee, --
109: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
110: *
111: * .. Scalar Arguments ..
112: CHARACTER UPLO
113: INTEGER INFO, LDA, N
114: * ..
115: * .. Array Arguments ..
116: DOUBLE PRECISION A( LDA, * )
117: * ..
118: *
119: * =====================================================================
120: *
121: * .. Parameters ..
122: DOUBLE PRECISION ONE, ZERO
123: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
124: * ..
125: * .. Local Scalars ..
1.3 bertrand 126: LOGICAL UPPER
1.1 bertrand 127: INTEGER N1, N2, IINFO
128: * ..
129: * .. External Functions ..
130: LOGICAL LSAME, DISNAN
131: EXTERNAL LSAME, DISNAN
132: * ..
133: * .. External Subroutines ..
134: EXTERNAL DSYRK, DTRSM, XERBLA
135: * ..
136: * .. Intrinsic Functions ..
137: INTRINSIC MAX, SQRT
138: * ..
139: * .. Executable Statements ..
140: *
141: * Test the input parameters
142: *
143: INFO = 0
144: UPPER = LSAME( UPLO, 'U' )
145: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
146: INFO = -1
147: ELSE IF( N.LT.0 ) THEN
148: INFO = -2
149: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
150: INFO = -4
151: END IF
152: IF( INFO.NE.0 ) THEN
153: CALL XERBLA( 'DPOTRF2', -INFO )
154: RETURN
155: END IF
156: *
157: * Quick return if possible
158: *
159: IF( N.EQ.0 )
160: $ RETURN
161: *
162: * N=1 case
163: *
164: IF( N.EQ.1 ) THEN
165: *
166: * Test for non-positive-definiteness
167: *
168: IF( A( 1, 1 ).LE.ZERO.OR.DISNAN( A( 1, 1 ) ) ) THEN
169: INFO = 1
170: RETURN
171: END IF
172: *
173: * Factor
174: *
175: A( 1, 1 ) = SQRT( A( 1, 1 ) )
176: *
177: * Use recursive code
178: *
179: ELSE
180: N1 = N/2
181: N2 = N-N1
182: *
183: * Factor A11
184: *
185: CALL DPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO )
186: IF ( IINFO.NE.0 ) THEN
187: INFO = IINFO
188: RETURN
1.3 bertrand 189: END IF
1.1 bertrand 190: *
191: * Compute the Cholesky factorization A = U**T*U
192: *
193: IF( UPPER ) THEN
194: *
195: * Update and scale A12
196: *
197: CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE,
1.3 bertrand 198: $ A( 1, 1 ), LDA, A( 1, N1+1 ), LDA )
1.1 bertrand 199: *
200: * Update and factor A22
1.3 bertrand 201: *
1.1 bertrand 202: CALL DSYRK( UPLO, 'T', N2, N1, -ONE, A( 1, N1+1 ), LDA,
203: $ ONE, A( N1+1, N1+1 ), LDA )
204: CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
205: IF ( IINFO.NE.0 ) THEN
206: INFO = IINFO + N1
207: RETURN
208: END IF
209: *
210: * Compute the Cholesky factorization A = L*L**T
211: *
212: ELSE
213: *
214: * Update and scale A21
215: *
1.3 bertrand 216: CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE,
1.1 bertrand 217: $ A( 1, 1 ), LDA, A( N1+1, 1 ), LDA )
218: *
219: * Update and factor A22
220: *
221: CALL DSYRK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA,
222: $ ONE, A( N1+1, N1+1 ), LDA )
223: CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
224: IF ( IINFO.NE.0 ) THEN
225: INFO = IINFO + N1
226: RETURN
227: END IF
228: END IF
229: END IF
230: RETURN
231: *
232: * End of DPOTRF2
233: *
234: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpotrf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack