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

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