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

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