swiftnav: zlahr2.c Source File

Go to the documentation of this file.
00001 /* zlahr2.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 zlahr2_(integer *n, integer *k, integer *nb, 
00023         doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *t, 
00024         integer *ldt, doublecomplex *y, integer *ldy)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, t_dim1, t_offset, y_dim1, y_offset, i__1, i__2, 
00028             i__3;
00029     doublecomplex z__1;
00030 
00031     /* Local variables */
00032     integer i__;
00033     doublecomplex ei;
00034     extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
00035             doublecomplex *, integer *), zgemm_(char *, char *, integer *, 
00036             integer *, integer *, doublecomplex *, doublecomplex *, integer *, 
00037              doublecomplex *, integer *, doublecomplex *, doublecomplex *, 
00038             integer *), zgemv_(char *, integer *, integer *, 
00039             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00040             integer *, doublecomplex *, doublecomplex *, integer *), 
00041             zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 
00042             integer *), ztrmm_(char *, char *, char *, char *, integer *, 
00043             integer *, doublecomplex *, doublecomplex *, integer *, 
00044             doublecomplex *, integer *), 
00045             zaxpy_(integer *, doublecomplex *, doublecomplex *, integer *, 
00046             doublecomplex *, integer *), ztrmv_(char *, char *, char *, 
00047             integer *, doublecomplex *, integer *, doublecomplex *, integer *), zlarfg_(integer *, doublecomplex *, 
00048             doublecomplex *, integer *, doublecomplex *), zlacgv_(integer *, 
00049             doublecomplex *, integer *), zlacpy_(char *, integer *, integer *, 
00050              doublecomplex *, integer *, doublecomplex *, integer *);
00051 
00052 
00053 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00054 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00055 /*     November 2006 */
00056 
00057 /*     .. Scalar Arguments .. */
00058 /*     .. */
00059 /*     .. Array Arguments .. */
00060 /*     .. */
00061 
00062 /*  Purpose */
00063 /*  ======= */
00064 
00065 /*  ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) */
00066 /*  matrix A so that elements below the k-th subdiagonal are zero. The */
00067 /*  reduction is performed by an unitary similarity transformation */
00068 /*  Q' * A * Q. The routine returns the matrices V and T which determine */
00069 /*  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. */
00070 
00071 /*  This is an auxiliary routine called by ZGEHRD. */
00072 
00073 /*  Arguments */
00074 /*  ========= */
00075 
00076 /*  N       (input) INTEGER */
00077 /*          The order of the matrix A. */
00078 
00079 /*  K       (input) INTEGER */
00080 /*          The offset for the reduction. Elements below the k-th */
00081 /*          subdiagonal in the first NB columns are reduced to zero. */
00082 /*          K < N. */
00083 
00084 /*  NB      (input) INTEGER */
00085 /*          The number of columns to be reduced. */
00086 
00087 /*  A       (input/output) COMPLEX*16 array, dimension (LDA,N-K+1) */
00088 /*          On entry, the n-by-(n-k+1) general matrix A. */
00089 /*          On exit, the elements on and above the k-th subdiagonal in */
00090 /*          the first NB columns are overwritten with the corresponding */
00091 /*          elements of the reduced matrix; the elements below the k-th */
00092 /*          subdiagonal, with the array TAU, represent the matrix Q as a */
00093 /*          product of elementary reflectors. The other columns of A are */
00094 /*          unchanged. See Further Details. */
00095 
00096 /*  LDA     (input) INTEGER */
00097 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00098 
00099 /*  TAU     (output) COMPLEX*16 array, dimension (NB) */
00100 /*          The scalar factors of the elementary reflectors. See Further */
00101 /*          Details. */
00102 
00103 /*  T       (output) COMPLEX*16 array, dimension (LDT,NB) */
00104 /*          The upper triangular matrix T. */
00105 
00106 /*  LDT     (input) INTEGER */
00107 /*          The leading dimension of the array T.  LDT >= NB. */
00108 
00109 /*  Y       (output) COMPLEX*16 array, dimension (LDY,NB) */
00110 /*          The n-by-nb matrix Y. */
00111 
00112 /*  LDY     (input) INTEGER */
00113 /*          The leading dimension of the array Y. LDY >= N. */
00114 
00115 /*  Further Details */
00116 /*  =============== */
00117 
00118 /*  The matrix Q is represented as a product of nb elementary reflectors */
00119 
00120 /*     Q = H(1) H(2) . . . H(nb). */
00121 
00122 /*  Each H(i) has the form */
00123 
00124 /*     H(i) = I - tau * v * v' */
00125 
00126 /*  where tau is a complex scalar, and v is a complex vector with */
00127 /*  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in */
00128 /*  A(i+k+1:n,i), and tau in TAU(i). */
00129 
00130 /*  The elements of the vectors v together form the (n-k+1)-by-nb matrix */
00131 /*  V which is needed, with T and Y, to apply the transformation to the */
00132 /*  unreduced part of the matrix, using an update of the form: */
00133 /*  A := (I - V*T*V') * (A - Y*V'). */
00134 
00135 /*  The contents of A on exit are illustrated by the following example */
00136 /*  with n = 7, k = 3 and nb = 2: */
00137 
00138 /*     ( a   a   a   a   a ) */
00139 /*     ( a   a   a   a   a ) */
00140 /*     ( a   a   a   a   a ) */
00141 /*     ( h   h   a   a   a ) */
00142 /*     ( v1  h   a   a   a ) */
00143 /*     ( v1  v2  a   a   a ) */
00144 /*     ( v1  v2  a   a   a ) */
00145 
00146 /*  where a denotes an element of the original matrix A, h denotes a */
00147 /*  modified element of the upper Hessenberg matrix H, and vi denotes an */
00148 /*  element of the vector defining H(i). */
00149 
00150 /*  This file is a slight modification of LAPACK-3.0's ZLAHRD */
00151 /*  incorporating improvements proposed by Quintana-Orti and Van de */
00152 /*  Gejin. Note that the entries of A(1:K,2:NB) differ from those */
00153 /*  returned by the original LAPACK routine. This function is */
00154 /*  not backward compatible with LAPACK3.0. */
00155 
00156 /*  ===================================================================== */
00157 
00158 /*     .. Parameters .. */
00159 /*     .. */
00160 /*     .. Local Scalars .. */
00161 /*     .. */
00162 /*     .. External Subroutines .. */
00163 /*     .. */
00164 /*     .. Intrinsic Functions .. */
00165 /*     .. */
00166 /*     .. Executable Statements .. */
00167 
00168 /*     Quick return if possible */
00169 
00170     /* Parameter adjustments */
00171     --tau;
00172     a_dim1 = *lda;
00173     a_offset = 1 + a_dim1;
00174     a -= a_offset;
00175     t_dim1 = *ldt;
00176     t_offset = 1 + t_dim1;
00177     t -= t_offset;
00178     y_dim1 = *ldy;
00179     y_offset = 1 + y_dim1;
00180     y -= y_offset;
00181 
00182     /* Function Body */
00183     if (*n <= 1) {
00184         return 0;
00185     }
00186 
00187     i__1 = *nb;
00188     for (i__ = 1; i__ <= i__1; ++i__) {
00189         if (i__ > 1) {
00190 
00191 /*           Update A(K+1:N,I) */
00192 
00193 /*           Update I-th column of A - Y * V' */
00194 
00195             i__2 = i__ - 1;
00196             zlacgv_(&i__2, &a[*k + i__ - 1 + a_dim1], lda);
00197             i__2 = *n - *k;
00198             i__3 = i__ - 1;
00199             z__1.r = -1., z__1.i = -0.;
00200             zgemv_("NO TRANSPOSE", &i__2, &i__3, &z__1, &y[*k + 1 + y_dim1], 
00201                     ldy, &a[*k + i__ - 1 + a_dim1], lda, &c_b2, &a[*k + 1 + 
00202                     i__ * a_dim1], &c__1);
00203             i__2 = i__ - 1;
00204             zlacgv_(&i__2, &a[*k + i__ - 1 + a_dim1], lda);
00205 
00206 /*           Apply I - V * T' * V' to this column (call it b) from the */
00207 /*           left, using the last column of T as workspace */
00208 
00209 /*           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows) */
00210 /*                    ( V2 )             ( b2 ) */
00211 
00212 /*           where V1 is unit lower triangular */
00213 
00214 /*           w := V1' * b1 */
00215 
00216             i__2 = i__ - 1;
00217             zcopy_(&i__2, &a[*k + 1 + i__ * a_dim1], &c__1, &t[*nb * t_dim1 + 
00218                     1], &c__1);
00219             i__2 = i__ - 1;
00220             ztrmv_("Lower", "Conjugate transpose", "UNIT", &i__2, &a[*k + 1 + 
00221                     a_dim1], lda, &t[*nb * t_dim1 + 1], &c__1);
00222 
00223 /*           w := w + V2'*b2 */
00224 
00225             i__2 = *n - *k - i__ + 1;
00226             i__3 = i__ - 1;
00227             zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*k + i__ + 
00228                     a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b2, &
00229                     t[*nb * t_dim1 + 1], &c__1);
00230 
00231 /*           w := T'*w */
00232 
00233             i__2 = i__ - 1;
00234             ztrmv_("Upper", "Conjugate transpose", "NON-UNIT", &i__2, &t[
00235                     t_offset], ldt, &t[*nb * t_dim1 + 1], &c__1);
00236 
00237 /*           b2 := b2 - V2*w */
00238 
00239             i__2 = *n - *k - i__ + 1;
00240             i__3 = i__ - 1;
00241             z__1.r = -1., z__1.i = -0.;
00242             zgemv_("NO TRANSPOSE", &i__2, &i__3, &z__1, &a[*k + i__ + a_dim1], 
00243                      lda, &t[*nb * t_dim1 + 1], &c__1, &c_b2, &a[*k + i__ + 
00244                     i__ * a_dim1], &c__1);
00245 
00246 /*           b1 := b1 - V1*w */
00247 
00248             i__2 = i__ - 1;
00249             ztrmv_("Lower", "NO TRANSPOSE", "UNIT", &i__2, &a[*k + 1 + a_dim1]
00250 , lda, &t[*nb * t_dim1 + 1], &c__1);
00251             i__2 = i__ - 1;
00252             z__1.r = -1., z__1.i = -0.;
00253             zaxpy_(&i__2, &z__1, &t[*nb * t_dim1 + 1], &c__1, &a[*k + 1 + i__ 
00254                     * a_dim1], &c__1);
00255 
00256             i__2 = *k + i__ - 1 + (i__ - 1) * a_dim1;
00257             a[i__2].r = ei.r, a[i__2].i = ei.i;
00258         }
00259 
00260 /*        Generate the elementary reflector H(I) to annihilate */
00261 /*        A(K+I+1:N,I) */
00262 
00263         i__2 = *n - *k - i__ + 1;
00264 /* Computing MIN */
00265         i__3 = *k + i__ + 1;
00266         zlarfg_(&i__2, &a[*k + i__ + i__ * a_dim1], &a[min(i__3, *n)+ i__ * 
00267                 a_dim1], &c__1, &tau[i__]);
00268         i__2 = *k + i__ + i__ * a_dim1;
00269         ei.r = a[i__2].r, ei.i = a[i__2].i;
00270         i__2 = *k + i__ + i__ * a_dim1;
00271         a[i__2].r = 1., a[i__2].i = 0.;
00272 
00273 /*        Compute  Y(K+1:N,I) */
00274 
00275         i__2 = *n - *k;
00276         i__3 = *n - *k - i__ + 1;
00277         zgemv_("NO TRANSPOSE", &i__2, &i__3, &c_b2, &a[*k + 1 + (i__ + 1) * 
00278                 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b1, &y[*
00279                 k + 1 + i__ * y_dim1], &c__1);
00280         i__2 = *n - *k - i__ + 1;
00281         i__3 = i__ - 1;
00282         zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*k + i__ + 
00283                 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b1, &t[
00284                 i__ * t_dim1 + 1], &c__1);
00285         i__2 = *n - *k;
00286         i__3 = i__ - 1;
00287         z__1.r = -1., z__1.i = -0.;
00288         zgemv_("NO TRANSPOSE", &i__2, &i__3, &z__1, &y[*k + 1 + y_dim1], ldy, 
00289                 &t[i__ * t_dim1 + 1], &c__1, &c_b2, &y[*k + 1 + i__ * y_dim1], 
00290                  &c__1);
00291         i__2 = *n - *k;
00292         zscal_(&i__2, &tau[i__], &y[*k + 1 + i__ * y_dim1], &c__1);
00293 
00294 /*        Compute T(1:I,I) */
00295 
00296         i__2 = i__ - 1;
00297         i__3 = i__;
00298         z__1.r = -tau[i__3].r, z__1.i = -tau[i__3].i;
00299         zscal_(&i__2, &z__1, &t[i__ * t_dim1 + 1], &c__1);
00300         i__2 = i__ - 1;
00301         ztrmv_("Upper", "No Transpose", "NON-UNIT", &i__2, &t[t_offset], ldt, 
00302                 &t[i__ * t_dim1 + 1], &c__1)
00303                 ;
00304         i__2 = i__ + i__ * t_dim1;
00305         i__3 = i__;
00306         t[i__2].r = tau[i__3].r, t[i__2].i = tau[i__3].i;
00307 
00308 /* L10: */
00309     }
00310     i__1 = *k + *nb + *nb * a_dim1;
00311     a[i__1].r = ei.r, a[i__1].i = ei.i;
00312 
00313 /*     Compute Y(1:K,1:NB) */
00314 
00315     zlacpy_("ALL", k, nb, &a[(a_dim1 << 1) + 1], lda, &y[y_offset], ldy);
00316     ztrmm_("RIGHT", "Lower", "NO TRANSPOSE", "UNIT", k, nb, &c_b2, &a[*k + 1 
00317             + a_dim1], lda, &y[y_offset], ldy);
00318     if (*n > *k + *nb) {
00319         i__1 = *n - *k - *nb;
00320         zgemm_("NO TRANSPOSE", "NO TRANSPOSE", k, nb, &i__1, &c_b2, &a[(*nb + 
00321                 2) * a_dim1 + 1], lda, &a[*k + 1 + *nb + a_dim1], lda, &c_b2, 
00322                 &y[y_offset], ldy);
00323     }
00324     ztrmm_("RIGHT", "Upper", "NO TRANSPOSE", "NON-UNIT", k, nb, &c_b2, &t[
00325             t_offset], ldt, &y[y_offset], ldy);
00326 
00327     return 0;
00328 
00329 /*     End of ZLAHR2 */
00330 
00331 } /* zlahr2_ */