00001 /* chetrd.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 integer c_n1 = -1; 00020 static integer c__3 = 3; 00021 static integer c__2 = 2; 00022 static real c_b23 = 1.f; 00023 00024 /* Subroutine */ int chetrd_(char *uplo, integer *n, complex *a, integer *lda, 00025 real *d__, real *e, complex *tau, complex *work, integer *lwork, 00026 integer *info) 00027 { 00028 /* System generated locals */ 00029 integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5; 00030 complex q__1; 00031 00032 /* Local variables */ 00033 integer i__, j, nb, kk, nx, iws; 00034 extern logical lsame_(char *, char *); 00035 integer nbmin, iinfo; 00036 logical upper; 00037 extern /* Subroutine */ int chetd2_(char *, integer *, complex *, integer 00038 *, real *, real *, complex *, integer *), cher2k_(char *, 00039 char *, integer *, integer *, complex *, complex *, integer *, 00040 complex *, integer *, real *, complex *, integer *), clatrd_(char *, integer *, integer *, complex *, integer 00041 *, real *, complex *, complex *, integer *), xerbla_(char 00042 *, integer *); 00043 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 00044 integer *, integer *); 00045 integer ldwork, lwkopt; 00046 logical lquery; 00047 00048 00049 /* -- LAPACK routine (version 3.2) -- */ 00050 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00051 /* November 2006 */ 00052 00053 /* .. Scalar Arguments .. */ 00054 /* .. */ 00055 /* .. Array Arguments .. */ 00056 /* .. */ 00057 00058 /* Purpose */ 00059 /* ======= */ 00060 00061 /* CHETRD reduces a complex Hermitian matrix A to real symmetric */ 00062 /* tridiagonal form T by a unitary similarity transformation: */ 00063 /* Q**H * A * Q = T. */ 00064 00065 /* Arguments */ 00066 /* ========= */ 00067 00068 /* UPLO (input) CHARACTER*1 */ 00069 /* = 'U': Upper triangle of A is stored; */ 00070 /* = 'L': Lower triangle of A is stored. */ 00071 00072 /* N (input) INTEGER */ 00073 /* The order of the matrix A. N >= 0. */ 00074 00075 /* A (input/output) COMPLEX array, dimension (LDA,N) */ 00076 /* On entry, the Hermitian 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 unitary */ 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 unitary 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) REAL array, dimension (N) */ 00098 /* The diagonal elements of the tridiagonal matrix T: */ 00099 /* D(i) = A(i,i). */ 00100 00101 /* E (output) REAL 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) COMPLEX array, dimension (N-1) */ 00106 /* The scalar factors of the elementary reflectors (see Further */ 00107 /* Details). */ 00108 00109 /* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) */ 00110 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ 00111 00112 /* LWORK (input) INTEGER */ 00113 /* The dimension of the array WORK. LWORK >= 1. */ 00114 /* For optimum performance LWORK >= N*NB, where NB is the */ 00115 /* optimal blocksize. */ 00116 00117 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00118 /* only calculates the optimal size of the WORK array, returns */ 00119 /* this value as the first entry of the WORK array, and no error */ 00120 /* message related to LWORK is issued by XERBLA. */ 00121 00122 /* INFO (output) INTEGER */ 00123 /* = 0: successful exit */ 00124 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00125 00126 /* Further Details */ 00127 /* =============== */ 00128 00129 /* If UPLO = 'U', the matrix Q is represented as a product of elementary */ 00130 /* reflectors */ 00131 00132 /* Q = H(n-1) . . . H(2) H(1). */ 00133 00134 /* Each H(i) has the form */ 00135 00136 /* H(i) = I - tau * v * v' */ 00137 00138 /* where tau is a complex scalar, and v is a complex vector with */ 00139 /* v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in */ 00140 /* A(1:i-1,i+1), and tau in TAU(i). */ 00141 00142 /* If UPLO = 'L', the matrix Q is represented as a product of elementary */ 00143 /* reflectors */ 00144 00145 /* Q = H(1) H(2) . . . H(n-1). */ 00146 00147 /* Each H(i) has the form */ 00148 00149 /* H(i) = I - tau * v * v' */ 00150 00151 /* where tau is a complex scalar, and v is a complex vector with */ 00152 /* v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), */ 00153 /* and tau in TAU(i). */ 00154 00155 /* The contents of A on exit are illustrated by the following examples */ 00156 /* with n = 5: */ 00157 00158 /* if UPLO = 'U': if UPLO = 'L': */ 00159 00160 /* ( d e v2 v3 v4 ) ( d ) */ 00161 /* ( d e v3 v4 ) ( e d ) */ 00162 /* ( d e v4 ) ( v1 e d ) */ 00163 /* ( d e ) ( v1 v2 e d ) */ 00164 /* ( d ) ( v1 v2 v3 e d ) */ 00165 00166 /* where d and e denote diagonal and off-diagonal elements of T, and vi */ 00167 /* denotes an element of the vector defining H(i). */ 00168 00169 /* ===================================================================== */ 00170 00171 /* .. Parameters .. */ 00172 /* .. */ 00173 /* .. Local Scalars .. */ 00174 /* .. */ 00175 /* .. External Subroutines .. */ 00176 /* .. */ 00177 /* .. Intrinsic Functions .. */ 00178 /* .. */ 00179 /* .. External Functions .. */ 00180 /* .. */ 00181 /* .. Executable Statements .. */ 00182 00183 /* Test the input parameters */ 00184 00185 /* Parameter adjustments */ 00186 a_dim1 = *lda; 00187 a_offset = 1 + a_dim1; 00188 a -= a_offset; 00189 --d__; 00190 --e; 00191 --tau; 00192 --work; 00193 00194 /* Function Body */ 00195 *info = 0; 00196 upper = lsame_(uplo, "U"); 00197 lquery = *lwork == -1; 00198 if (! upper && ! lsame_(uplo, "L")) { 00199 *info = -1; 00200 } else if (*n < 0) { 00201 *info = -2; 00202 } else if (*lda < max(1,*n)) { 00203 *info = -4; 00204 } else if (*lwork < 1 && ! lquery) { 00205 *info = -9; 00206 } 00207 00208 if (*info == 0) { 00209 00210 /* Determine the block size. */ 00211 00212 nb = ilaenv_(&c__1, "CHETRD", uplo, n, &c_n1, &c_n1, &c_n1); 00213 lwkopt = *n * nb; 00214 work[1].r = (real) lwkopt, work[1].i = 0.f; 00215 } 00216 00217 if (*info != 0) { 00218 i__1 = -(*info); 00219 xerbla_("CHETRD", &i__1); 00220 return 0; 00221 } else if (lquery) { 00222 return 0; 00223 } 00224 00225 /* Quick return if possible */ 00226 00227 if (*n == 0) { 00228 work[1].r = 1.f, work[1].i = 0.f; 00229 return 0; 00230 } 00231 00232 nx = *n; 00233 iws = 1; 00234 if (nb > 1 && nb < *n) { 00235 00236 /* Determine when to cross over from blocked to unblocked code */ 00237 /* (last block is always handled by unblocked code). */ 00238 00239 /* Computing MAX */ 00240 i__1 = nb, i__2 = ilaenv_(&c__3, "CHETRD", uplo, n, &c_n1, &c_n1, & 00241 c_n1); 00242 nx = max(i__1,i__2); 00243 if (nx < *n) { 00244 00245 /* Determine if workspace is large enough for blocked code. */ 00246 00247 ldwork = *n; 00248 iws = ldwork * nb; 00249 if (*lwork < iws) { 00250 00251 /* Not enough workspace to use optimal NB: determine the */ 00252 /* minimum value of NB, and reduce NB or force use of */ 00253 /* unblocked code by setting NX = N. */ 00254 00255 /* Computing MAX */ 00256 i__1 = *lwork / ldwork; 00257 nb = max(i__1,1); 00258 nbmin = ilaenv_(&c__2, "CHETRD", uplo, n, &c_n1, &c_n1, &c_n1); 00259 if (nb < nbmin) { 00260 nx = *n; 00261 } 00262 } 00263 } else { 00264 nx = *n; 00265 } 00266 } else { 00267 nb = 1; 00268 } 00269 00270 if (upper) { 00271 00272 /* Reduce the upper triangle of A. */ 00273 /* Columns 1:kk are handled by the unblocked method. */ 00274 00275 kk = *n - (*n - nx + nb - 1) / nb * nb; 00276 i__1 = kk + 1; 00277 i__2 = -nb; 00278 for (i__ = *n - nb + 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += 00279 i__2) { 00280 00281 /* Reduce columns i:i+nb-1 to tridiagonal form and form the */ 00282 /* matrix W which is needed to update the unreduced part of */ 00283 /* the matrix */ 00284 00285 i__3 = i__ + nb - 1; 00286 clatrd_(uplo, &i__3, &nb, &a[a_offset], lda, &e[1], &tau[1], & 00287 work[1], &ldwork); 00288 00289 /* Update the unreduced submatrix A(1:i-1,1:i-1), using an */ 00290 /* update of the form: A := A - V*W' - W*V' */ 00291 00292 i__3 = i__ - 1; 00293 q__1.r = -1.f, q__1.i = -0.f; 00294 cher2k_(uplo, "No transpose", &i__3, &nb, &q__1, &a[i__ * a_dim1 00295 + 1], lda, &work[1], &ldwork, &c_b23, &a[a_offset], lda); 00296 00297 /* Copy superdiagonal elements back into A, and diagonal */ 00298 /* elements into D */ 00299 00300 i__3 = i__ + nb - 1; 00301 for (j = i__; j <= i__3; ++j) { 00302 i__4 = j - 1 + j * a_dim1; 00303 i__5 = j - 1; 00304 a[i__4].r = e[i__5], a[i__4].i = 0.f; 00305 i__4 = j; 00306 i__5 = j + j * a_dim1; 00307 d__[i__4] = a[i__5].r; 00308 /* L10: */ 00309 } 00310 /* L20: */ 00311 } 00312 00313 /* Use unblocked code to reduce the last or only block */ 00314 00315 chetd2_(uplo, &kk, &a[a_offset], lda, &d__[1], &e[1], &tau[1], &iinfo); 00316 } else { 00317 00318 /* Reduce the lower triangle of A */ 00319 00320 i__2 = *n - nx; 00321 i__1 = nb; 00322 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) { 00323 00324 /* Reduce columns i:i+nb-1 to tridiagonal form and form the */ 00325 /* matrix W which is needed to update the unreduced part of */ 00326 /* the matrix */ 00327 00328 i__3 = *n - i__ + 1; 00329 clatrd_(uplo, &i__3, &nb, &a[i__ + i__ * a_dim1], lda, &e[i__], & 00330 tau[i__], &work[1], &ldwork); 00331 00332 /* Update the unreduced submatrix A(i+nb:n,i+nb:n), using */ 00333 /* an update of the form: A := A - V*W' - W*V' */ 00334 00335 i__3 = *n - i__ - nb + 1; 00336 q__1.r = -1.f, q__1.i = -0.f; 00337 cher2k_(uplo, "No transpose", &i__3, &nb, &q__1, &a[i__ + nb + 00338 i__ * a_dim1], lda, &work[nb + 1], &ldwork, &c_b23, &a[ 00339 i__ + nb + (i__ + nb) * a_dim1], lda); 00340 00341 /* Copy subdiagonal elements back into A, and diagonal */ 00342 /* elements into D */ 00343 00344 i__3 = i__ + nb - 1; 00345 for (j = i__; j <= i__3; ++j) { 00346 i__4 = j + 1 + j * a_dim1; 00347 i__5 = j; 00348 a[i__4].r = e[i__5], a[i__4].i = 0.f; 00349 i__4 = j; 00350 i__5 = j + j * a_dim1; 00351 d__[i__4] = a[i__5].r; 00352 /* L30: */ 00353 } 00354 /* L40: */ 00355 } 00356 00357 /* Use unblocked code to reduce the last or only block */ 00358 00359 i__1 = *n - i__ + 1; 00360 chetd2_(uplo, &i__1, &a[i__ + i__ * a_dim1], lda, &d__[i__], &e[i__], 00361 &tau[i__], &iinfo); 00362 } 00363 00364 work[1].r = (real) lwkopt, work[1].i = 0.f; 00365 return 0; 00366 00367 /* End of CHETRD */ 00368 00369 } /* chetrd_ */