1 SUBROUTINE DPPTRF( UPLO, N, AP, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DPPTRF computes the Cholesky factorization of a real symmetric
20 * positive definite matrix A stored in packed format.
21 *
22 * The factorization has the form
23 * A = U**T * U, if UPLO = 'U', or
24 * A = L * L**T, if UPLO = 'L',
25 * where U is an upper triangular matrix and L is lower triangular.
26 *
27 * Arguments
28 * =========
29 *
30 * UPLO (input) CHARACTER*1
31 * = 'U': Upper triangle of A is stored;
32 * = 'L': Lower triangle of A is stored.
33 *
34 * N (input) INTEGER
35 * The order of the matrix A. N >= 0.
36 *
37 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
38 * On entry, the upper or lower triangle of the symmetric matrix
39 * A, packed columnwise in a linear array. The j-th column of A
40 * is stored in the array AP as follows:
41 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
42 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
43 * See below for further details.
44 *
45 * On exit, if INFO = 0, the triangular factor U or L from the
46 * Cholesky factorization A = U**T*U or A = L*L**T, in the same
47 * storage format as A.
48 *
49 * INFO (output) INTEGER
50 * = 0: successful exit
51 * < 0: if INFO = -i, the i-th argument had an illegal value
52 * > 0: if INFO = i, the leading minor of order i is not
53 * positive definite, and the factorization could not be
54 * completed.
55 *
56 * Further Details
57 * ======= =======
58 *
59 * The packed storage scheme is illustrated by the following example
60 * when N = 4, UPLO = 'U':
61 *
62 * Two-dimensional storage of the symmetric matrix A:
63 *
64 * a11 a12 a13 a14
65 * a22 a23 a24
66 * a33 a34 (aij = aji)
67 * a44
68 *
69 * Packed storage of the upper triangle of A:
70 *
71 * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
72 *
73 * =====================================================================
74 *
75 * .. Parameters ..
76 DOUBLE PRECISION ONE, ZERO
77 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
78 * ..
79 * .. Local Scalars ..
80 LOGICAL UPPER
81 INTEGER J, JC, JJ
82 DOUBLE PRECISION AJJ
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 DOUBLE PRECISION DDOT
87 EXTERNAL LSAME, DDOT
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL DSCAL, DSPR, DTPSV, XERBLA
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input parameters.
98 *
99 INFO = 0
100 UPPER = LSAME( UPLO, 'U' )
101 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
102 INFO = -1
103 ELSE IF( N.LT.0 ) THEN
104 INFO = -2
105 END IF
106 IF( INFO.NE.0 ) THEN
107 CALL XERBLA( 'DPPTRF', -INFO )
108 RETURN
109 END IF
110 *
111 * Quick return if possible
112 *
113 IF( N.EQ.0 )
114 $ RETURN
115 *
116 IF( UPPER ) THEN
117 *
118 * Compute the Cholesky factorization A = U**T*U.
119 *
120 JJ = 0
121 DO 10 J = 1, N
122 JC = JJ + 1
123 JJ = JJ + J
124 *
125 * Compute elements 1:J-1 of column J.
126 *
127 IF( J.GT.1 )
128 $ CALL DTPSV( 'Upper', 'Transpose', 'Non-unit', J-1, AP,
129 $ AP( JC ), 1 )
130 *
131 * Compute U(J,J) and test for non-positive-definiteness.
132 *
133 AJJ = AP( JJ ) - DDOT( J-1, AP( JC ), 1, AP( JC ), 1 )
134 IF( AJJ.LE.ZERO ) THEN
135 AP( JJ ) = AJJ
136 GO TO 30
137 END IF
138 AP( JJ ) = SQRT( AJJ )
139 10 CONTINUE
140 ELSE
141 *
142 * Compute the Cholesky factorization A = L*L**T.
143 *
144 JJ = 1
145 DO 20 J = 1, N
146 *
147 * Compute L(J,J) and test for non-positive-definiteness.
148 *
149 AJJ = AP( JJ )
150 IF( AJJ.LE.ZERO ) THEN
151 AP( JJ ) = AJJ
152 GO TO 30
153 END IF
154 AJJ = SQRT( AJJ )
155 AP( JJ ) = AJJ
156 *
157 * Compute elements J+1:N of column J and update the trailing
158 * submatrix.
159 *
160 IF( J.LT.N ) THEN
161 CALL DSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
162 CALL DSPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
163 $ AP( JJ+N-J+1 ) )
164 JJ = JJ + N - J + 1
165 END IF
166 20 CONTINUE
167 END IF
168 GO TO 40
169 *
170 30 CONTINUE
171 INFO = J
172 *
173 40 CONTINUE
174 RETURN
175 *
176 * End of DPPTRF
177 *
178 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DPPTRF computes the Cholesky factorization of a real symmetric
20 * positive definite matrix A stored in packed format.
21 *
22 * The factorization has the form
23 * A = U**T * U, if UPLO = 'U', or
24 * A = L * L**T, if UPLO = 'L',
25 * where U is an upper triangular matrix and L is lower triangular.
26 *
27 * Arguments
28 * =========
29 *
30 * UPLO (input) CHARACTER*1
31 * = 'U': Upper triangle of A is stored;
32 * = 'L': Lower triangle of A is stored.
33 *
34 * N (input) INTEGER
35 * The order of the matrix A. N >= 0.
36 *
37 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
38 * On entry, the upper or lower triangle of the symmetric matrix
39 * A, packed columnwise in a linear array. The j-th column of A
40 * is stored in the array AP as follows:
41 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
42 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
43 * See below for further details.
44 *
45 * On exit, if INFO = 0, the triangular factor U or L from the
46 * Cholesky factorization A = U**T*U or A = L*L**T, in the same
47 * storage format as A.
48 *
49 * INFO (output) INTEGER
50 * = 0: successful exit
51 * < 0: if INFO = -i, the i-th argument had an illegal value
52 * > 0: if INFO = i, the leading minor of order i is not
53 * positive definite, and the factorization could not be
54 * completed.
55 *
56 * Further Details
57 * ======= =======
58 *
59 * The packed storage scheme is illustrated by the following example
60 * when N = 4, UPLO = 'U':
61 *
62 * Two-dimensional storage of the symmetric matrix A:
63 *
64 * a11 a12 a13 a14
65 * a22 a23 a24
66 * a33 a34 (aij = aji)
67 * a44
68 *
69 * Packed storage of the upper triangle of A:
70 *
71 * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
72 *
73 * =====================================================================
74 *
75 * .. Parameters ..
76 DOUBLE PRECISION ONE, ZERO
77 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
78 * ..
79 * .. Local Scalars ..
80 LOGICAL UPPER
81 INTEGER J, JC, JJ
82 DOUBLE PRECISION AJJ
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 DOUBLE PRECISION DDOT
87 EXTERNAL LSAME, DDOT
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL DSCAL, DSPR, DTPSV, XERBLA
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input parameters.
98 *
99 INFO = 0
100 UPPER = LSAME( UPLO, 'U' )
101 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
102 INFO = -1
103 ELSE IF( N.LT.0 ) THEN
104 INFO = -2
105 END IF
106 IF( INFO.NE.0 ) THEN
107 CALL XERBLA( 'DPPTRF', -INFO )
108 RETURN
109 END IF
110 *
111 * Quick return if possible
112 *
113 IF( N.EQ.0 )
114 $ RETURN
115 *
116 IF( UPPER ) THEN
117 *
118 * Compute the Cholesky factorization A = U**T*U.
119 *
120 JJ = 0
121 DO 10 J = 1, N
122 JC = JJ + 1
123 JJ = JJ + J
124 *
125 * Compute elements 1:J-1 of column J.
126 *
127 IF( J.GT.1 )
128 $ CALL DTPSV( 'Upper', 'Transpose', 'Non-unit', J-1, AP,
129 $ AP( JC ), 1 )
130 *
131 * Compute U(J,J) and test for non-positive-definiteness.
132 *
133 AJJ = AP( JJ ) - DDOT( J-1, AP( JC ), 1, AP( JC ), 1 )
134 IF( AJJ.LE.ZERO ) THEN
135 AP( JJ ) = AJJ
136 GO TO 30
137 END IF
138 AP( JJ ) = SQRT( AJJ )
139 10 CONTINUE
140 ELSE
141 *
142 * Compute the Cholesky factorization A = L*L**T.
143 *
144 JJ = 1
145 DO 20 J = 1, N
146 *
147 * Compute L(J,J) and test for non-positive-definiteness.
148 *
149 AJJ = AP( JJ )
150 IF( AJJ.LE.ZERO ) THEN
151 AP( JJ ) = AJJ
152 GO TO 30
153 END IF
154 AJJ = SQRT( AJJ )
155 AP( JJ ) = AJJ
156 *
157 * Compute elements J+1:N of column J and update the trailing
158 * submatrix.
159 *
160 IF( J.LT.N ) THEN
161 CALL DSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
162 CALL DSPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
163 $ AP( JJ+N-J+1 ) )
164 JJ = JJ + N - J + 1
165 END IF
166 20 CONTINUE
167 END IF
168 GO TO 40
169 *
170 30 CONTINUE
171 INFO = J
172 *
173 40 CONTINUE
174 RETURN
175 *
176 * End of DPPTRF
177 *
178 END