00001 /* dsytd2.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 static doublereal c_b8 = 0.; 00020 static doublereal c_b14 = -1.; 00021 00022 /* Subroutine */ int dsytd2_(char *uplo, integer *n, doublereal *a, integer * 00023 lda, doublereal *d__, doublereal *e, doublereal *tau, integer *info) 00024 { 00025 /* System generated locals */ 00026 integer a_dim1, a_offset, i__1, i__2, i__3; 00027 00028 /* Local variables */ 00029 integer i__; 00030 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 00031 integer *); 00032 doublereal taui; 00033 extern /* Subroutine */ int dsyr2_(char *, integer *, doublereal *, 00034 doublereal *, integer *, doublereal *, integer *, doublereal *, 00035 integer *); 00036 doublereal alpha; 00037 extern logical lsame_(char *, char *); 00038 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 00039 integer *, doublereal *, integer *); 00040 logical upper; 00041 extern /* Subroutine */ int dsymv_(char *, integer *, doublereal *, 00042 doublereal *, integer *, doublereal *, integer *, doublereal *, 00043 doublereal *, integer *), dlarfg_(integer *, doublereal *, 00044 doublereal *, integer *, doublereal *), xerbla_(char *, integer * 00045 ); 00046 00047 00048 /* -- LAPACK routine (version 3.2) -- */ 00049 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00050 /* November 2006 */ 00051 00052 /* .. Scalar Arguments .. */ 00053 /* .. */ 00054 /* .. Array Arguments .. */ 00055 /* .. */ 00056 00057 /* Purpose */ 00058 /* ======= */ 00059 00060 /* DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal */ 00061 /* form T by an orthogonal similarity transformation: Q' * A * Q = T. */ 00062 00063 /* Arguments */ 00064 /* ========= */ 00065 00066 /* UPLO (input) CHARACTER*1 */ 00067 /* Specifies whether the upper or lower triangular part of the */ 00068 /* symmetric matrix A is stored: */ 00069 /* = 'U': Upper triangular */ 00070 /* = 'L': Lower triangular */ 00071 00072 /* N (input) INTEGER */ 00073 /* The order of the matrix A. N >= 0. */ 00074 00075 /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ 00076 /* On entry, the symmetric matrix A. If UPLO = 'U', the leading */ 00077 /* n-by-n upper triangular part of A contains the upper */ 00078 /* triangular part of the matrix A, and the strictly lower */ 00079 /* triangular part of A is not referenced. If UPLO = 'L', the */ 00080 /* leading n-by-n lower triangular part of A contains the lower */ 00081 /* triangular part of the matrix A, and the strictly upper */ 00082 /* triangular part of A is not referenced. */ 00083 /* On exit, if UPLO = 'U', the diagonal and first superdiagonal */ 00084 /* of A are overwritten by the corresponding elements of the */ 00085 /* tridiagonal matrix T, and the elements above the first */ 00086 /* superdiagonal, with the array TAU, represent the orthogonal */ 00087 /* matrix Q as a product of elementary reflectors; if UPLO */ 00088 /* = 'L', the diagonal and first subdiagonal of A are over- */ 00089 /* written by the corresponding elements of the tridiagonal */ 00090 /* matrix T, and the elements below the first subdiagonal, with */ 00091 /* the array TAU, represent the orthogonal matrix Q as a product */ 00092 /* of elementary reflectors. See Further Details. */ 00093 00094 /* LDA (input) INTEGER */ 00095 /* The leading dimension of the array A. LDA >= max(1,N). */ 00096 00097 /* D (output) DOUBLE PRECISION array, dimension (N) */ 00098 /* The diagonal elements of the tridiagonal matrix T: */ 00099 /* D(i) = A(i,i). */ 00100 00101 /* E (output) DOUBLE PRECISION array, dimension (N-1) */ 00102 /* The off-diagonal elements of the tridiagonal matrix T: */ 00103 /* E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. */ 00104 00105 /* TAU (output) DOUBLE PRECISION array, dimension (N-1) */ 00106 /* The scalar factors of the elementary reflectors (see Further */ 00107 /* Details). */ 00108 00109 /* INFO (output) INTEGER */ 00110 /* = 0: successful exit */ 00111 /* < 0: if INFO = -i, the i-th argument had an illegal value. */ 00112 00113 /* Further Details */ 00114 /* =============== */ 00115 00116 /* If UPLO = 'U', the matrix Q is represented as a product of elementary */ 00117 /* reflectors */ 00118 00119 /* Q = H(n-1) . . . H(2) H(1). */ 00120 00121 /* Each H(i) has the form */ 00122 00123 /* H(i) = I - tau * v * v' */ 00124 00125 /* where tau is a real scalar, and v is a real vector with */ 00126 /* v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in */ 00127 /* A(1:i-1,i+1), and tau in TAU(i). */ 00128 00129 /* If UPLO = 'L', the matrix Q is represented as a product of elementary */ 00130 /* reflectors */ 00131 00132 /* Q = H(1) H(2) . . . H(n-1). */ 00133 00134 /* Each H(i) has the form */ 00135 00136 /* H(i) = I - tau * v * v' */ 00137 00138 /* where tau is a real scalar, and v is a real vector with */ 00139 /* v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), */ 00140 /* and tau in TAU(i). */ 00141 00142 /* The contents of A on exit are illustrated by the following examples */ 00143 /* with n = 5: */ 00144 00145 /* if UPLO = 'U': if UPLO = 'L': */ 00146 00147 /* ( d e v2 v3 v4 ) ( d ) */ 00148 /* ( d e v3 v4 ) ( e d ) */ 00149 /* ( d e v4 ) ( v1 e d ) */ 00150 /* ( d e ) ( v1 v2 e d ) */ 00151 /* ( d ) ( v1 v2 v3 e d ) */ 00152 00153 /* where d and e denote diagonal and off-diagonal elements of T, and vi */ 00154 /* denotes an element of the vector defining H(i). */ 00155 00156 /* ===================================================================== */ 00157 00158 /* .. Parameters .. */ 00159 /* .. */ 00160 /* .. Local Scalars .. */ 00161 /* .. */ 00162 /* .. External Subroutines .. */ 00163 /* .. */ 00164 /* .. External Functions .. */ 00165 /* .. */ 00166 /* .. Intrinsic Functions .. */ 00167 /* .. */ 00168 /* .. Executable Statements .. */ 00169 00170 /* Test the input parameters */ 00171 00172 /* Parameter adjustments */ 00173 a_dim1 = *lda; 00174 a_offset = 1 + a_dim1; 00175 a -= a_offset; 00176 --d__; 00177 --e; 00178 --tau; 00179 00180 /* Function Body */ 00181 *info = 0; 00182 upper = lsame_(uplo, "U"); 00183 if (! upper && ! lsame_(uplo, "L")) { 00184 *info = -1; 00185 } else if (*n < 0) { 00186 *info = -2; 00187 } else if (*lda < max(1,*n)) { 00188 *info = -4; 00189 } 00190 if (*info != 0) { 00191 i__1 = -(*info); 00192 xerbla_("DSYTD2", &i__1); 00193 return 0; 00194 } 00195 00196 /* Quick return if possible */ 00197 00198 if (*n <= 0) { 00199 return 0; 00200 } 00201 00202 if (upper) { 00203 00204 /* Reduce the upper triangle of A */ 00205 00206 for (i__ = *n - 1; i__ >= 1; --i__) { 00207 00208 /* Generate elementary reflector H(i) = I - tau * v * v' */ 00209 /* to annihilate A(1:i-1,i+1) */ 00210 00211 dlarfg_(&i__, &a[i__ + (i__ + 1) * a_dim1], &a[(i__ + 1) * a_dim1 00212 + 1], &c__1, &taui); 00213 e[i__] = a[i__ + (i__ + 1) * a_dim1]; 00214 00215 if (taui != 0.) { 00216 00217 /* Apply H(i) from both sides to A(1:i,1:i) */ 00218 00219 a[i__ + (i__ + 1) * a_dim1] = 1.; 00220 00221 /* Compute x := tau * A * v storing x in TAU(1:i) */ 00222 00223 dsymv_(uplo, &i__, &taui, &a[a_offset], lda, &a[(i__ + 1) * 00224 a_dim1 + 1], &c__1, &c_b8, &tau[1], &c__1); 00225 00226 /* Compute w := x - 1/2 * tau * (x'*v) * v */ 00227 00228 alpha = taui * -.5 * ddot_(&i__, &tau[1], &c__1, &a[(i__ + 1) 00229 * a_dim1 + 1], &c__1); 00230 daxpy_(&i__, &alpha, &a[(i__ + 1) * a_dim1 + 1], &c__1, &tau[ 00231 1], &c__1); 00232 00233 /* Apply the transformation as a rank-2 update: */ 00234 /* A := A - v * w' - w * v' */ 00235 00236 dsyr2_(uplo, &i__, &c_b14, &a[(i__ + 1) * a_dim1 + 1], &c__1, 00237 &tau[1], &c__1, &a[a_offset], lda); 00238 00239 a[i__ + (i__ + 1) * a_dim1] = e[i__]; 00240 } 00241 d__[i__ + 1] = a[i__ + 1 + (i__ + 1) * a_dim1]; 00242 tau[i__] = taui; 00243 /* L10: */ 00244 } 00245 d__[1] = a[a_dim1 + 1]; 00246 } else { 00247 00248 /* Reduce the lower triangle of A */ 00249 00250 i__1 = *n - 1; 00251 for (i__ = 1; i__ <= i__1; ++i__) { 00252 00253 /* Generate elementary reflector H(i) = I - tau * v * v' */ 00254 /* to annihilate A(i+2:n,i) */ 00255 00256 i__2 = *n - i__; 00257 /* Computing MIN */ 00258 i__3 = i__ + 2; 00259 dlarfg_(&i__2, &a[i__ + 1 + i__ * a_dim1], &a[min(i__3, *n)+ i__ * 00260 a_dim1], &c__1, &taui); 00261 e[i__] = a[i__ + 1 + i__ * a_dim1]; 00262 00263 if (taui != 0.) { 00264 00265 /* Apply H(i) from both sides to A(i+1:n,i+1:n) */ 00266 00267 a[i__ + 1 + i__ * a_dim1] = 1.; 00268 00269 /* Compute x := tau * A * v storing y in TAU(i:n-1) */ 00270 00271 i__2 = *n - i__; 00272 dsymv_(uplo, &i__2, &taui, &a[i__ + 1 + (i__ + 1) * a_dim1], 00273 lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b8, &tau[ 00274 i__], &c__1); 00275 00276 /* Compute w := x - 1/2 * tau * (x'*v) * v */ 00277 00278 i__2 = *n - i__; 00279 alpha = taui * -.5 * ddot_(&i__2, &tau[i__], &c__1, &a[i__ + 00280 1 + i__ * a_dim1], &c__1); 00281 i__2 = *n - i__; 00282 daxpy_(&i__2, &alpha, &a[i__ + 1 + i__ * a_dim1], &c__1, &tau[ 00283 i__], &c__1); 00284 00285 /* Apply the transformation as a rank-2 update: */ 00286 /* A := A - v * w' - w * v' */ 00287 00288 i__2 = *n - i__; 00289 dsyr2_(uplo, &i__2, &c_b14, &a[i__ + 1 + i__ * a_dim1], &c__1, 00290 &tau[i__], &c__1, &a[i__ + 1 + (i__ + 1) * a_dim1], 00291 lda); 00292 00293 a[i__ + 1 + i__ * a_dim1] = e[i__]; 00294 } 00295 d__[i__] = a[i__ + i__ * a_dim1]; 00296 tau[i__] = taui; 00297 /* L20: */ 00298 } 00299 d__[*n] = a[*n + *n * a_dim1]; 00300 } 00301 00302 return 0; 00303 00304 /* End of DSYTD2 */ 00305 00306 } /* dsytd2_ */