zlatrd.c
Go to the documentation of this file.
00001 /* zlatrd.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 doublecomplex c_b1 = {0.,0.};
00019 static doublecomplex c_b2 = {1.,0.};
00020 static integer c__1 = 1;
00021 
00022 /* Subroutine */ int zlatrd_(char *uplo, integer *n, integer *nb, 
00023         doublecomplex *a, integer *lda, doublereal *e, doublecomplex *tau, 
00024         doublecomplex *w, integer *ldw)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3;
00028     doublereal d__1;
00029     doublecomplex z__1, z__2, z__3, z__4;
00030 
00031     /* Local variables */
00032     integer i__, iw;
00033     doublecomplex alpha;
00034     extern logical lsame_(char *, char *);
00035     extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
00036             doublecomplex *, integer *);
00037     extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 
00038             doublecomplex *, integer *, doublecomplex *, integer *);
00039     extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 
00040             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00041             integer *, doublecomplex *, doublecomplex *, integer *), 
00042             zhemv_(char *, integer *, doublecomplex *, doublecomplex *, 
00043             integer *, doublecomplex *, integer *, doublecomplex *, 
00044             doublecomplex *, integer *), zaxpy_(integer *, 
00045             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00046             integer *), zlarfg_(integer *, doublecomplex *, doublecomplex *, 
00047             integer *, doublecomplex *), zlacgv_(integer *, doublecomplex *, 
00048             integer *);
00049 
00050 
00051 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00052 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00053 /*     November 2006 */
00054 
00055 /*     .. Scalar Arguments .. */
00056 /*     .. */
00057 /*     .. Array Arguments .. */
00058 /*     .. */
00059 
00060 /*  Purpose */
00061 /*  ======= */
00062 
00063 /*  ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to */
00064 /*  Hermitian tridiagonal form by a unitary similarity */
00065 /*  transformation Q' * A * Q, and returns the matrices V and W which are */
00066 /*  needed to apply the transformation to the unreduced part of A. */
00067 
00068 /*  If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a */
00069 /*  matrix, of which the upper triangle is supplied; */
00070 /*  if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a */
00071 /*  matrix, of which the lower triangle is supplied. */
00072 
00073 /*  This is an auxiliary routine called by ZHETRD. */
00074 
00075 /*  Arguments */
00076 /*  ========= */
00077 
00078 /*  UPLO    (input) CHARACTER*1 */
00079 /*          Specifies whether the upper or lower triangular part of the */
00080 /*          Hermitian matrix A is stored: */
00081 /*          = 'U': Upper triangular */
00082 /*          = 'L': Lower triangular */
00083 
00084 /*  N       (input) INTEGER */
00085 /*          The order of the matrix A. */
00086 
00087 /*  NB      (input) INTEGER */
00088 /*          The number of rows and columns to be reduced. */
00089 
00090 /*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) */
00091 /*          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading */
00092 /*          n-by-n upper triangular part of A contains the upper */
00093 /*          triangular part of the matrix A, and the strictly lower */
00094 /*          triangular part of A is not referenced.  If UPLO = 'L', the */
00095 /*          leading n-by-n lower triangular part of A contains the lower */
00096 /*          triangular part of the matrix A, and the strictly upper */
00097 /*          triangular part of A is not referenced. */
00098 /*          On exit: */
00099 /*          if UPLO = 'U', the last NB columns have been reduced to */
00100 /*            tridiagonal form, with the diagonal elements overwriting */
00101 /*            the diagonal elements of A; the elements above the diagonal */
00102 /*            with the array TAU, represent the unitary matrix Q as a */
00103 /*            product of elementary reflectors; */
00104 /*          if UPLO = 'L', the first NB columns have been reduced to */
00105 /*            tridiagonal form, with the diagonal elements overwriting */
00106 /*            the diagonal elements of A; the elements below the diagonal */
00107 /*            with the array TAU, represent the  unitary matrix Q as a */
00108 /*            product of elementary reflectors. */
00109 /*          See Further Details. */
00110 
00111 /*  LDA     (input) INTEGER */
00112 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00113 
00114 /*  E       (output) DOUBLE PRECISION array, dimension (N-1) */
00115 /*          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal */
00116 /*          elements of the last NB columns of the reduced matrix; */
00117 /*          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of */
00118 /*          the first NB columns of the reduced matrix. */
00119 
00120 /*  TAU     (output) COMPLEX*16 array, dimension (N-1) */
00121 /*          The scalar factors of the elementary reflectors, stored in */
00122 /*          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. */
00123 /*          See Further Details. */
00124 
00125 /*  W       (output) COMPLEX*16 array, dimension (LDW,NB) */
00126 /*          The n-by-nb matrix W required to update the unreduced part */
00127 /*          of A. */
00128 
00129 /*  LDW     (input) INTEGER */
00130 /*          The leading dimension of the array W. LDW >= max(1,N). */
00131 
00132 /*  Further Details */
00133 /*  =============== */
00134 
00135 /*  If UPLO = 'U', the matrix Q is represented as a product of elementary */
00136 /*  reflectors */
00137 
00138 /*     Q = H(n) H(n-1) . . . H(n-nb+1). */
00139 
00140 /*  Each H(i) has the form */
00141 
00142 /*     H(i) = I - tau * v * v' */
00143 
00144 /*  where tau is a complex scalar, and v is a complex vector with */
00145 /*  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), */
00146 /*  and tau in TAU(i-1). */
00147 
00148 /*  If UPLO = 'L', the matrix Q is represented as a product of elementary */
00149 /*  reflectors */
00150 
00151 /*     Q = H(1) H(2) . . . H(nb). */
00152 
00153 /*  Each H(i) has the form */
00154 
00155 /*     H(i) = I - tau * v * v' */
00156 
00157 /*  where tau is a complex scalar, and v is a complex vector with */
00158 /*  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), */
00159 /*  and tau in TAU(i). */
00160 
00161 /*  The elements of the vectors v together form the n-by-nb matrix V */
00162 /*  which is needed, with W, to apply the transformation to the unreduced */
00163 /*  part of the matrix, using a Hermitian rank-2k update of the form: */
00164 /*  A := A - V*W' - W*V'. */
00165 
00166 /*  The contents of A on exit are illustrated by the following examples */
00167 /*  with n = 5 and nb = 2: */
00168 
00169 /*  if UPLO = 'U':                       if UPLO = 'L': */
00170 
00171 /*    (  a   a   a   v4  v5 )              (  d                  ) */
00172 /*    (      a   a   v4  v5 )              (  1   d              ) */
00173 /*    (          a   1   v5 )              (  v1  1   a          ) */
00174 /*    (              d   1  )              (  v1  v2  a   a      ) */
00175 /*    (                  d  )              (  v1  v2  a   a   a  ) */
00176 
00177 /*  where d denotes a diagonal element of the reduced matrix, a denotes */
00178 /*  an element of the original matrix that is unchanged, and vi denotes */
00179 /*  an element of the vector defining H(i). */
00180 
00181 /*  ===================================================================== */
00182 
00183 /*     .. Parameters .. */
00184 /*     .. */
00185 /*     .. Local Scalars .. */
00186 /*     .. */
00187 /*     .. External Subroutines .. */
00188 /*     .. */
00189 /*     .. External Functions .. */
00190 /*     .. */
00191 /*     .. Intrinsic Functions .. */
00192 /*     .. */
00193 /*     .. Executable Statements .. */
00194 
00195 /*     Quick return if possible */
00196 
00197     /* Parameter adjustments */
00198     a_dim1 = *lda;
00199     a_offset = 1 + a_dim1;
00200     a -= a_offset;
00201     --e;
00202     --tau;
00203     w_dim1 = *ldw;
00204     w_offset = 1 + w_dim1;
00205     w -= w_offset;
00206 
00207     /* Function Body */
00208     if (*n <= 0) {
00209         return 0;
00210     }
00211 
00212     if (lsame_(uplo, "U")) {
00213 
00214 /*        Reduce last NB columns of upper triangle */
00215 
00216         i__1 = *n - *nb + 1;
00217         for (i__ = *n; i__ >= i__1; --i__) {
00218             iw = i__ - *n + *nb;
00219             if (i__ < *n) {
00220 
00221 /*              Update A(1:i,i) */
00222 
00223                 i__2 = i__ + i__ * a_dim1;
00224                 i__3 = i__ + i__ * a_dim1;
00225                 d__1 = a[i__3].r;
00226                 a[i__2].r = d__1, a[i__2].i = 0.;
00227                 i__2 = *n - i__;
00228                 zlacgv_(&i__2, &w[i__ + (iw + 1) * w_dim1], ldw);
00229                 i__2 = *n - i__;
00230                 z__1.r = -1., z__1.i = -0.;
00231                 zgemv_("No transpose", &i__, &i__2, &z__1, &a[(i__ + 1) * 
00232                         a_dim1 + 1], lda, &w[i__ + (iw + 1) * w_dim1], ldw, &
00233                         c_b2, &a[i__ * a_dim1 + 1], &c__1);
00234                 i__2 = *n - i__;
00235                 zlacgv_(&i__2, &w[i__ + (iw + 1) * w_dim1], ldw);
00236                 i__2 = *n - i__;
00237                 zlacgv_(&i__2, &a[i__ + (i__ + 1) * a_dim1], lda);
00238                 i__2 = *n - i__;
00239                 z__1.r = -1., z__1.i = -0.;
00240                 zgemv_("No transpose", &i__, &i__2, &z__1, &w[(iw + 1) * 
00241                         w_dim1 + 1], ldw, &a[i__ + (i__ + 1) * a_dim1], lda, &
00242                         c_b2, &a[i__ * a_dim1 + 1], &c__1);
00243                 i__2 = *n - i__;
00244                 zlacgv_(&i__2, &a[i__ + (i__ + 1) * a_dim1], lda);
00245                 i__2 = i__ + i__ * a_dim1;
00246                 i__3 = i__ + i__ * a_dim1;
00247                 d__1 = a[i__3].r;
00248                 a[i__2].r = d__1, a[i__2].i = 0.;
00249             }
00250             if (i__ > 1) {
00251 
00252 /*              Generate elementary reflector H(i) to annihilate */
00253 /*              A(1:i-2,i) */
00254 
00255                 i__2 = i__ - 1 + i__ * a_dim1;
00256                 alpha.r = a[i__2].r, alpha.i = a[i__2].i;
00257                 i__2 = i__ - 1;
00258                 zlarfg_(&i__2, &alpha, &a[i__ * a_dim1 + 1], &c__1, &tau[i__ 
00259                         - 1]);
00260                 i__2 = i__ - 1;
00261                 e[i__2] = alpha.r;
00262                 i__2 = i__ - 1 + i__ * a_dim1;
00263                 a[i__2].r = 1., a[i__2].i = 0.;
00264 
00265 /*              Compute W(1:i-1,i) */
00266 
00267                 i__2 = i__ - 1;
00268                 zhemv_("Upper", &i__2, &c_b2, &a[a_offset], lda, &a[i__ * 
00269                         a_dim1 + 1], &c__1, &c_b1, &w[iw * w_dim1 + 1], &c__1);
00270                 if (i__ < *n) {
00271                     i__2 = i__ - 1;
00272                     i__3 = *n - i__;
00273                     zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &w[(iw 
00274                             + 1) * w_dim1 + 1], ldw, &a[i__ * a_dim1 + 1], &
00275                             c__1, &c_b1, &w[i__ + 1 + iw * w_dim1], &c__1);
00276                     i__2 = i__ - 1;
00277                     i__3 = *n - i__;
00278                     z__1.r = -1., z__1.i = -0.;
00279                     zgemv_("No transpose", &i__2, &i__3, &z__1, &a[(i__ + 1) *
00280                              a_dim1 + 1], lda, &w[i__ + 1 + iw * w_dim1], &
00281                             c__1, &c_b2, &w[iw * w_dim1 + 1], &c__1);
00282                     i__2 = i__ - 1;
00283                     i__3 = *n - i__;
00284                     zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[(
00285                             i__ + 1) * a_dim1 + 1], lda, &a[i__ * a_dim1 + 1], 
00286                              &c__1, &c_b1, &w[i__ + 1 + iw * w_dim1], &c__1);
00287                     i__2 = i__ - 1;
00288                     i__3 = *n - i__;
00289                     z__1.r = -1., z__1.i = -0.;
00290                     zgemv_("No transpose", &i__2, &i__3, &z__1, &w[(iw + 1) * 
00291                             w_dim1 + 1], ldw, &w[i__ + 1 + iw * w_dim1], &
00292                             c__1, &c_b2, &w[iw * w_dim1 + 1], &c__1);
00293                 }
00294                 i__2 = i__ - 1;
00295                 zscal_(&i__2, &tau[i__ - 1], &w[iw * w_dim1 + 1], &c__1);
00296                 z__3.r = -.5, z__3.i = -0.;
00297                 i__2 = i__ - 1;
00298                 z__2.r = z__3.r * tau[i__2].r - z__3.i * tau[i__2].i, z__2.i =
00299                          z__3.r * tau[i__2].i + z__3.i * tau[i__2].r;
00300                 i__3 = i__ - 1;
00301                 zdotc_(&z__4, &i__3, &w[iw * w_dim1 + 1], &c__1, &a[i__ * 
00302                         a_dim1 + 1], &c__1);
00303                 z__1.r = z__2.r * z__4.r - z__2.i * z__4.i, z__1.i = z__2.r * 
00304                         z__4.i + z__2.i * z__4.r;
00305                 alpha.r = z__1.r, alpha.i = z__1.i;
00306                 i__2 = i__ - 1;
00307                 zaxpy_(&i__2, &alpha, &a[i__ * a_dim1 + 1], &c__1, &w[iw * 
00308                         w_dim1 + 1], &c__1);
00309             }
00310 
00311 /* L10: */
00312         }
00313     } else {
00314 
00315 /*        Reduce first NB columns of lower triangle */
00316 
00317         i__1 = *nb;
00318         for (i__ = 1; i__ <= i__1; ++i__) {
00319 
00320 /*           Update A(i:n,i) */
00321 
00322             i__2 = i__ + i__ * a_dim1;
00323             i__3 = i__ + i__ * a_dim1;
00324             d__1 = a[i__3].r;
00325             a[i__2].r = d__1, a[i__2].i = 0.;
00326             i__2 = i__ - 1;
00327             zlacgv_(&i__2, &w[i__ + w_dim1], ldw);
00328             i__2 = *n - i__ + 1;
00329             i__3 = i__ - 1;
00330             z__1.r = -1., z__1.i = -0.;
00331             zgemv_("No transpose", &i__2, &i__3, &z__1, &a[i__ + a_dim1], lda, 
00332                      &w[i__ + w_dim1], ldw, &c_b2, &a[i__ + i__ * a_dim1], &
00333                     c__1);
00334             i__2 = i__ - 1;
00335             zlacgv_(&i__2, &w[i__ + w_dim1], ldw);
00336             i__2 = i__ - 1;
00337             zlacgv_(&i__2, &a[i__ + a_dim1], lda);
00338             i__2 = *n - i__ + 1;
00339             i__3 = i__ - 1;
00340             z__1.r = -1., z__1.i = -0.;
00341             zgemv_("No transpose", &i__2, &i__3, &z__1, &w[i__ + w_dim1], ldw, 
00342                      &a[i__ + a_dim1], lda, &c_b2, &a[i__ + i__ * a_dim1], &
00343                     c__1);
00344             i__2 = i__ - 1;
00345             zlacgv_(&i__2, &a[i__ + a_dim1], lda);
00346             i__2 = i__ + i__ * a_dim1;
00347             i__3 = i__ + i__ * a_dim1;
00348             d__1 = a[i__3].r;
00349             a[i__2].r = d__1, a[i__2].i = 0.;
00350             if (i__ < *n) {
00351 
00352 /*              Generate elementary reflector H(i) to annihilate */
00353 /*              A(i+2:n,i) */
00354 
00355                 i__2 = i__ + 1 + i__ * a_dim1;
00356                 alpha.r = a[i__2].r, alpha.i = a[i__2].i;
00357                 i__2 = *n - i__;
00358 /* Computing MIN */
00359                 i__3 = i__ + 2;
00360                 zlarfg_(&i__2, &alpha, &a[min(i__3, *n)+ i__ * a_dim1], &c__1, 
00361                          &tau[i__]);
00362                 i__2 = i__;
00363                 e[i__2] = alpha.r;
00364                 i__2 = i__ + 1 + i__ * a_dim1;
00365                 a[i__2].r = 1., a[i__2].i = 0.;
00366 
00367 /*              Compute W(i+1:n,i) */
00368 
00369                 i__2 = *n - i__;
00370                 zhemv_("Lower", &i__2, &c_b2, &a[i__ + 1 + (i__ + 1) * a_dim1]
00371 , lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b1, &w[
00372                         i__ + 1 + i__ * w_dim1], &c__1);
00373                 i__2 = *n - i__;
00374                 i__3 = i__ - 1;
00375                 zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &w[i__ + 1 
00376                         + w_dim1], ldw, &a[i__ + 1 + i__ * a_dim1], &c__1, &
00377                         c_b1, &w[i__ * w_dim1 + 1], &c__1);
00378                 i__2 = *n - i__;
00379                 i__3 = i__ - 1;
00380                 z__1.r = -1., z__1.i = -0.;
00381                 zgemv_("No transpose", &i__2, &i__3, &z__1, &a[i__ + 1 + 
00382                         a_dim1], lda, &w[i__ * w_dim1 + 1], &c__1, &c_b2, &w[
00383                         i__ + 1 + i__ * w_dim1], &c__1);
00384                 i__2 = *n - i__;
00385                 i__3 = i__ - 1;
00386                 zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[i__ + 1 
00387                         + a_dim1], lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &
00388                         c_b1, &w[i__ * w_dim1 + 1], &c__1);
00389                 i__2 = *n - i__;
00390                 i__3 = i__ - 1;
00391                 z__1.r = -1., z__1.i = -0.;
00392                 zgemv_("No transpose", &i__2, &i__3, &z__1, &w[i__ + 1 + 
00393                         w_dim1], ldw, &w[i__ * w_dim1 + 1], &c__1, &c_b2, &w[
00394                         i__ + 1 + i__ * w_dim1], &c__1);
00395                 i__2 = *n - i__;
00396                 zscal_(&i__2, &tau[i__], &w[i__ + 1 + i__ * w_dim1], &c__1);
00397                 z__3.r = -.5, z__3.i = -0.;
00398                 i__2 = i__;
00399                 z__2.r = z__3.r * tau[i__2].r - z__3.i * tau[i__2].i, z__2.i =
00400                          z__3.r * tau[i__2].i + z__3.i * tau[i__2].r;
00401                 i__3 = *n - i__;
00402                 zdotc_(&z__4, &i__3, &w[i__ + 1 + i__ * w_dim1], &c__1, &a[
00403                         i__ + 1 + i__ * a_dim1], &c__1);
00404                 z__1.r = z__2.r * z__4.r - z__2.i * z__4.i, z__1.i = z__2.r * 
00405                         z__4.i + z__2.i * z__4.r;
00406                 alpha.r = z__1.r, alpha.i = z__1.i;
00407                 i__2 = *n - i__;
00408                 zaxpy_(&i__2, &alpha, &a[i__ + 1 + i__ * a_dim1], &c__1, &w[
00409                         i__ + 1 + i__ * w_dim1], &c__1);
00410             }
00411 
00412 /* L20: */
00413         }
00414     }
00415 
00416     return 0;
00417 
00418 /*     End of ZLATRD */
00419 
00420 } /* zlatrd_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:42