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

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:32