3 /* Subroutine */ int dsytrf_(char *uplo, integer *n, doublereal *a, integer *
4 lda, integer *ipiv, doublereal *work, integer *lwork, integer *info)
6 /* -- LAPACK routine (version 3.0) --
7 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
8 Courant Institute, Argonne National Lab, and Rice University
15 DSYTRF computes the factorization of a real symmetric matrix A using
16 the Bunch-Kaufman diagonal pivoting method. The form of the
19 A = U*D*U**T or A = L*D*L**T
21 where U (or L) is a product of permutation and unit upper (lower)
22 triangular matrices, and D is symmetric and block diagonal with
23 1-by-1 and 2-by-2 diagonal blocks.
25 This is the blocked version of the algorithm, calling Level 3 BLAS.
30 UPLO (input) CHARACTER*1
31 = 'U': Upper triangle of A is stored;
32 = 'L': Lower triangle of A is stored.
35 The order of the matrix A. N >= 0.
37 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
38 On entry, the symmetric matrix A. If UPLO = 'U', the leading
39 N-by-N upper triangular part of A contains the upper
40 triangular part of the matrix A, and the strictly lower
41 triangular part of A is not referenced. If UPLO = 'L', the
42 leading N-by-N lower triangular part of A contains the lower
43 triangular part of the matrix A, and the strictly upper
44 triangular part of A is not referenced.
46 On exit, the block diagonal matrix D and the multipliers used
47 to obtain the factor U or L (see below for further details).
50 The leading dimension of the array A. LDA >= max(1,N).
52 IPIV (output) INTEGER array, dimension (N)
53 Details of the interchanges and the block structure of D.
54 If IPIV(k) > 0, then rows and columns k and IPIV(k) were
55 interchanged and D(k,k) is a 1-by-1 diagonal block.
56 If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
57 columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
58 is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
59 IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
60 interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
62 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
63 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
66 The length of WORK. LWORK >=1. For best performance
67 LWORK >= N*NB, where NB is the block size returned by ILAENV.
69 If LWORK = -1, then a workspace query is assumed; the routine
70 only calculates the optimal size of the WORK array, returns
71 this value as the first entry of the WORK array, and no error
72 message related to LWORK is issued by XERBLA.
76 < 0: if INFO = -i, the i-th argument had an illegal value
77 > 0: if INFO = i, D(i,i) is exactly zero. The factorization
78 has been completed, but the block diagonal matrix D is
79 exactly singular, and division by zero will occur if it
80 is used to solve a system of equations.
85 If UPLO = 'U', then A = U*D*U', where
86 U = P(n)*U(n)* ... *P(k)U(k)* ...,
87 i.e., U is a product of terms P(k)*U(k), where k decreases from n to
88 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
89 and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
90 defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
91 that if the diagonal block D(k) is of order s (s = 1 or 2), then
98 If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
99 If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
100 and A(k,k), and v overwrites A(1:k-2,k-1:k).
102 If UPLO = 'L', then A = L*D*L', where
103 L = P(1)*L(1)* ... *P(k)*L(k)* ...,
104 i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
105 n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
106 and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
107 defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
108 that if the diagonal block D(k) is of order s (s = 1 or 2), then
115 If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
116 If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
117 and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
119 =====================================================================
122 Test the input parameters.
124 Parameter adjustments */
125 /* Table of constant values */
126 static integer c__1 = 1;
127 static integer c_n1 = -1;
128 static integer c__2 = 2;
130 /* System generated locals */
131 integer a_dim1, a_offset, i__1, i__2;
132 /* Local variables */
134 extern logical lsame_(char *, char *);
135 static integer nbmin, iinfo;
136 static logical upper;
137 extern /* Subroutine */ int dsytf2_(char *, integer *, doublereal *,
138 integer *, integer *, integer *);
139 static integer kb, nb;
140 extern /* Subroutine */ int xerbla_(char *, integer *);
141 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
142 integer *, integer *, ftnlen, ftnlen);
143 extern /* Subroutine */ int dlasyf_(char *, integer *, integer *, integer
144 *, doublereal *, integer *, integer *, doublereal *, integer *,
146 static integer ldwork, lwkopt;
147 static logical lquery;
149 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
153 a_offset = 1 + a_dim1 * 1;
160 upper = lsame_(uplo, "U");
161 lquery = *lwork == -1;
162 if (! upper && ! lsame_(uplo, "L")) {
166 } else if (*lda < max(1,*n)) {
168 } else if (*lwork < 1 && ! lquery) {
174 /* Determine the block size */
176 nb = ilaenv_(&c__1, "DSYTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
179 work[1] = (doublereal) lwkopt;
184 xerbla_("DSYTRF", &i__1);
192 if (nb > 1 && nb < *n) {
196 i__1 = *lwork / ldwork;
199 i__1 = 2, i__2 = ilaenv_(&c__2, "DSYTRF", uplo, n, &c_n1, &c_n1, &
200 c_n1, (ftnlen)6, (ftnlen)1);
201 nbmin = max(i__1,i__2);
212 /* Factorize A as U*D*U' using the upper triangle of A
214 K is the main loop index, decreasing from N to 1 in steps of
215 KB, where KB is the number of columns factorized by DLASYF;
216 KB is either NB or NB-1, or K for the last block */
221 /* If K < 1, exit from loop */
229 /* Factorize columns k-kb+1:k of A and use blocked code to
230 update columns 1:k-kb */
232 dlasyf_(uplo, &k, &nb, &kb, &a[a_offset], lda, &ipiv[1], &work[1],
236 /* Use unblocked code to factorize columns 1:k of A */
238 dsytf2_(uplo, &k, &a[a_offset], lda, &ipiv[1], &iinfo);
242 /* Set INFO on the first occurrence of a zero pivot */
244 if (*info == 0 && iinfo > 0) {
248 /* Decrease K and return to the start of the main loop */
255 /* Factorize A as L*D*L' using the lower triangle of A
257 K is the main loop index, increasing from 1 to N in steps of
258 KB, where KB is the number of columns factorized by DLASYF;
259 KB is either NB or NB-1, or N-K+1 for the last block */
264 /* If K > N, exit from loop */
272 /* Factorize columns k:k+kb-1 of A and use blocked code to
273 update columns k+kb:n */
276 dlasyf_(uplo, &i__1, &nb, &kb, &a_ref(k, k), lda, &ipiv[k], &work[
277 1], &ldwork, &iinfo);
280 /* Use unblocked code to factorize columns k:n of A */
283 dsytf2_(uplo, &i__1, &a_ref(k, k), lda, &ipiv[k], &iinfo);
287 /* Set INFO on the first occurrence of a zero pivot */
289 if (*info == 0 && iinfo > 0) {
290 *info = iinfo + k - 1;
296 for (j = k; j <= i__1; ++j) {
298 ipiv[j] = ipiv[j] + k - 1;
300 ipiv[j] = ipiv[j] - k + 1;
305 /* Increase K and return to the start of the main loop */
313 work[1] = (doublereal) lwkopt;