3 /* Subroutine */ int dsytf2_(char *uplo, integer *n, doublereal *a, integer *
4 lda, integer *ipiv, integer *info)
6 /* -- LAPACK routine (version 3.1) --
7 Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
14 DSYTF2 computes the factorization of a real symmetric matrix A using
15 the Bunch-Kaufman diagonal pivoting method:
17 A = U*D*U' or A = L*D*L'
19 where U (or L) is a product of permutation and unit upper (lower)
20 triangular matrices, U' is the transpose of U, and D is symmetric and
21 block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
23 This is the unblocked version of the algorithm, calling Level 2 BLAS.
28 UPLO (input) CHARACTER*1
29 Specifies whether the upper or lower triangular part of the
30 symmetric matrix A is stored:
31 = 'U': Upper triangular
32 = 'L': Lower triangular
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.
64 < 0: if INFO = -k, the k-th argument had an illegal value
65 > 0: if INFO = k, D(k,k) is exactly zero. The factorization
66 has been completed, but the block diagonal matrix D is
67 exactly singular, and division by zero will occur if it
68 is used to solve a system of equations.
74 Bobby Cheng, MathWorks
76 Replace l.204 and l.372
77 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
79 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
81 01-01-96 - Based on modifications by
82 J. Lewis, Boeing Computer Services Company
83 A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
84 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
87 If UPLO = 'U', then A = U*D*U', where
88 U = P(n)*U(n)* ... *P(k)U(k)* ...,
89 i.e., U is a product of terms P(k)*U(k), where k decreases from n to
90 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
91 and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
92 defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
93 that if the diagonal block D(k) is of order s (s = 1 or 2), then
100 If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
101 If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
102 and A(k,k), and v overwrites A(1:k-2,k-1:k).
104 If UPLO = 'L', then A = L*D*L', where
105 L = P(1)*L(1)* ... *P(k)*L(k)* ...,
106 i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
107 n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
108 and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
109 defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
110 that if the diagonal block D(k) is of order s (s = 1 or 2), then
117 If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
118 If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
119 and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
121 =====================================================================
124 Test the input parameters.
126 Parameter adjustments */
127 /* Table of constant values */
128 static integer c__1 = 1;
130 /* System generated locals */
131 integer a_dim1, a_offset, i__1, i__2;
132 doublereal d__1, d__2, d__3;
133 /* Builtin functions */
134 double sqrt(doublereal);
135 /* Local variables */
136 static integer i__, j, k;
137 static doublereal t, r1, d11, d12, d21, d22;
138 static integer kk, kp;
139 static doublereal wk, wkm1, wkp1;
140 static integer imax, jmax;
141 extern /* Subroutine */ int dsyr_(char *, integer *, doublereal *,
142 doublereal *, integer *, doublereal *, integer *);
143 static doublereal alpha;
144 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
146 extern logical lsame_(char *, char *);
147 extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
148 doublereal *, integer *);
149 static integer kstep;
150 static logical upper;
151 static doublereal absakk;
152 extern integer idamax_(integer *, doublereal *, integer *);
153 extern logical disnan_(doublereal *);
154 extern /* Subroutine */ int xerbla_(char *, integer *);
155 static doublereal colmax, rowmax;
159 a_offset = 1 + a_dim1;
165 upper = lsame_(uplo, "U");
166 if (! upper && ! lsame_(uplo, "L")) {
170 } else if (*lda < max(1,*n)) {
175 xerbla_("DSYTF2", &i__1);
179 /* Initialize ALPHA for use in choosing pivot block size. */
181 alpha = (sqrt(17.) + 1.) / 8.;
185 /* Factorize A as U*D*U' using the upper triangle of A
187 K is the main loop index, decreasing from N to 1 in steps of
193 /* If K < 1, exit from loop */
200 /* Determine rows and columns to be interchanged and whether
201 a 1-by-1 or 2-by-2 pivot block will be used */
203 absakk = (d__1 = a[k + k * a_dim1], abs(d__1));
205 /* IMAX is the row-index of the largest off-diagonal element in
206 column K, and COLMAX is its absolute value */
210 imax = idamax_(&i__1, &a[k * a_dim1 + 1], &c__1);
211 colmax = (d__1 = a[imax + k * a_dim1], abs(d__1));
216 if (max(absakk,colmax) == 0. || disnan_(&absakk)) {
218 /* Column K is zero or contains a NaN: set INFO and continue */
225 if (absakk >= alpha * colmax) {
227 /* no interchange, use 1-by-1 pivot block */
232 /* JMAX is the column-index of the largest off-diagonal
233 element in row IMAX, and ROWMAX is its absolute value */
236 jmax = imax + idamax_(&i__1, &a[imax + (imax + 1) * a_dim1],
238 rowmax = (d__1 = a[imax + jmax * a_dim1], abs(d__1));
241 jmax = idamax_(&i__1, &a[imax * a_dim1 + 1], &c__1);
243 d__2 = rowmax, d__3 = (d__1 = a[jmax + imax * a_dim1],
245 rowmax = max(d__2,d__3);
248 if (absakk >= alpha * colmax * (colmax / rowmax)) {
250 /* no interchange, use 1-by-1 pivot block */
253 } else if ((d__1 = a[imax + imax * a_dim1], abs(d__1)) >=
256 /* interchange rows and columns K and IMAX, use 1-by-1
262 /* interchange rows and columns K-1 and IMAX, use 2-by-2
273 /* Interchange rows and columns KK and KP in the leading
274 submatrix A(1:k,1:k) */
277 dswap_(&i__1, &a[kk * a_dim1 + 1], &c__1, &a[kp * a_dim1 + 1],
280 dswap_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + (kp +
282 t = a[kk + kk * a_dim1];
283 a[kk + kk * a_dim1] = a[kp + kp * a_dim1];
284 a[kp + kp * a_dim1] = t;
286 t = a[k - 1 + k * a_dim1];
287 a[k - 1 + k * a_dim1] = a[kp + k * a_dim1];
288 a[kp + k * a_dim1] = t;
292 /* Update the leading submatrix */
296 /* 1-by-1 pivot block D(k): column k now holds
300 where U(k) is the k-th column of U
302 Perform a rank-1 update of A(1:k-1,1:k-1) as
304 A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)' */
306 r1 = 1. / a[k + k * a_dim1];
309 dsyr_(uplo, &i__1, &d__1, &a[k * a_dim1 + 1], &c__1, &a[
312 /* Store U(k) in column k */
315 dscal_(&i__1, &r1, &a[k * a_dim1 + 1], &c__1);
318 /* 2-by-2 pivot block D(k): columns k and k-1 now hold
320 ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
322 where U(k) and U(k-1) are the k-th and (k-1)-th columns
325 Perform a rank-2 update of A(1:k-2,1:k-2) as
327 A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
328 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )' */
332 d12 = a[k - 1 + k * a_dim1];
333 d22 = a[k - 1 + (k - 1) * a_dim1] / d12;
334 d11 = a[k + k * a_dim1] / d12;
335 t = 1. / (d11 * d22 - 1.);
338 for (j = k - 2; j >= 1; --j) {
339 wkm1 = d12 * (d11 * a[j + (k - 1) * a_dim1] - a[j + k
341 wk = d12 * (d22 * a[j + k * a_dim1] - a[j + (k - 1) *
343 for (i__ = j; i__ >= 1; --i__) {
344 a[i__ + j * a_dim1] = a[i__ + j * a_dim1] - a[i__
345 + k * a_dim1] * wk - a[i__ + (k - 1) *
349 a[j + k * a_dim1] = wk;
350 a[j + (k - 1) * a_dim1] = wkm1;
359 /* Store details of the interchanges in IPIV */
368 /* Decrease K and return to the start of the main loop */
375 /* Factorize A as L*D*L' using the lower triangle of A
377 K is the main loop index, increasing from 1 to N in steps of
383 /* If K > N, exit from loop */
390 /* Determine rows and columns to be interchanged and whether
391 a 1-by-1 or 2-by-2 pivot block will be used */
393 absakk = (d__1 = a[k + k * a_dim1], abs(d__1));
395 /* IMAX is the row-index of the largest off-diagonal element in
396 column K, and COLMAX is its absolute value */
400 imax = k + idamax_(&i__1, &a[k + 1 + k * a_dim1], &c__1);
401 colmax = (d__1 = a[imax + k * a_dim1], abs(d__1));
406 if (max(absakk,colmax) == 0. || disnan_(&absakk)) {
408 /* Column K is zero or contains a NaN: set INFO and continue */
415 if (absakk >= alpha * colmax) {
417 /* no interchange, use 1-by-1 pivot block */
422 /* JMAX is the column-index of the largest off-diagonal
423 element in row IMAX, and ROWMAX is its absolute value */
426 jmax = k - 1 + idamax_(&i__1, &a[imax + k * a_dim1], lda);
427 rowmax = (d__1 = a[imax + jmax * a_dim1], abs(d__1));
430 jmax = imax + idamax_(&i__1, &a[imax + 1 + imax * a_dim1],
433 d__2 = rowmax, d__3 = (d__1 = a[jmax + imax * a_dim1],
435 rowmax = max(d__2,d__3);
438 if (absakk >= alpha * colmax * (colmax / rowmax)) {
440 /* no interchange, use 1-by-1 pivot block */
443 } else if ((d__1 = a[imax + imax * a_dim1], abs(d__1)) >=
446 /* interchange rows and columns K and IMAX, use 1-by-1
452 /* interchange rows and columns K+1 and IMAX, use 2-by-2
463 /* Interchange rows and columns KK and KP in the trailing
464 submatrix A(k:n,k:n) */
468 dswap_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + 1
469 + kp * a_dim1], &c__1);
472 dswap_(&i__1, &a[kk + 1 + kk * a_dim1], &c__1, &a[kp + (kk +
474 t = a[kk + kk * a_dim1];
475 a[kk + kk * a_dim1] = a[kp + kp * a_dim1];
476 a[kp + kp * a_dim1] = t;
478 t = a[k + 1 + k * a_dim1];
479 a[k + 1 + k * a_dim1] = a[kp + k * a_dim1];
480 a[kp + k * a_dim1] = t;
484 /* Update the trailing submatrix */
488 /* 1-by-1 pivot block D(k): column k now holds
492 where L(k) is the k-th column of L */
496 /* Perform a rank-1 update of A(k+1:n,k+1:n) as
498 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)' */
500 d11 = 1. / a[k + k * a_dim1];
503 dsyr_(uplo, &i__1, &d__1, &a[k + 1 + k * a_dim1], &c__1, &
504 a[k + 1 + (k + 1) * a_dim1], lda);
506 /* Store L(k) in column K */
509 dscal_(&i__1, &d11, &a[k + 1 + k * a_dim1], &c__1);
513 /* 2-by-2 pivot block D(k) */
517 /* Perform a rank-2 update of A(k+2:n,k+2:n) as
519 A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))'
521 where L(k) and L(k+1) are the k-th and (k+1)-th
524 d21 = a[k + 1 + k * a_dim1];
525 d11 = a[k + 1 + (k + 1) * a_dim1] / d21;
526 d22 = a[k + k * a_dim1] / d21;
527 t = 1. / (d11 * d22 - 1.);
531 for (j = k + 2; j <= i__1; ++j) {
533 wk = d21 * (d11 * a[j + k * a_dim1] - a[j + (k + 1) *
535 wkp1 = d21 * (d22 * a[j + (k + 1) * a_dim1] - a[j + k
539 for (i__ = j; i__ <= i__2; ++i__) {
540 a[i__ + j * a_dim1] = a[i__ + j * a_dim1] - a[i__
541 + k * a_dim1] * wk - a[i__ + (k + 1) *
546 a[j + k * a_dim1] = wk;
547 a[j + (k + 1) * a_dim1] = wkp1;
555 /* Store details of the interchanges in IPIV */
564 /* Increase K and return to the start of the main loop */