swiftnav: zlarft.c Source File

Go to the documentation of this file.
00001 /* zlarft.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_b2 = {0.,0.};
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int zlarft_(char *direct, char *storev, integer *n, integer *
00022         k, doublecomplex *v, integer *ldv, doublecomplex *tau, doublecomplex *
00023         t, integer *ldt)
00024 {
00025     /* System generated locals */
00026     integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;
00027     doublecomplex z__1;
00028 
00029     /* Local variables */
00030     integer i__, j, prevlastv;
00031     doublecomplex vii;
00032     extern logical lsame_(char *, char *);
00033     extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 
00034             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00035             integer *, doublecomplex *, doublecomplex *, integer *);
00036     integer lastv;
00037     extern /* Subroutine */ int ztrmv_(char *, char *, char *, integer *, 
00038             doublecomplex *, integer *, doublecomplex *, integer *), zlacgv_(integer *, doublecomplex *, integer *);
00039 
00040 
00041 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00042 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00043 /*     November 2006 */
00044 
00045 /*     .. Scalar Arguments .. */
00046 /*     .. */
00047 /*     .. Array Arguments .. */
00048 /*     .. */
00049 
00050 /*  Purpose */
00051 /*  ======= */
00052 
00053 /*  ZLARFT forms the triangular factor T of a complex block reflector H */
00054 /*  of order n, which is defined as a product of k elementary reflectors. */
00055 
00056 /*  If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; */
00057 
00058 /*  If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. */
00059 
00060 /*  If STOREV = 'C', the vector which defines the elementary reflector */
00061 /*  H(i) is stored in the i-th column of the array V, and */
00062 
00063 /*     H  =  I - V * T * V' */
00064 
00065 /*  If STOREV = 'R', the vector which defines the elementary reflector */
00066 /*  H(i) is stored in the i-th row of the array V, and */
00067 
00068 /*     H  =  I - V' * T * V */
00069 
00070 /*  Arguments */
00071 /*  ========= */
00072 
00073 /*  DIRECT  (input) CHARACTER*1 */
00074 /*          Specifies the order in which the elementary reflectors are */
00075 /*          multiplied to form the block reflector: */
00076 /*          = 'F': H = H(1) H(2) . . . H(k) (Forward) */
00077 /*          = 'B': H = H(k) . . . H(2) H(1) (Backward) */
00078 
00079 /*  STOREV  (input) CHARACTER*1 */
00080 /*          Specifies how the vectors which define the elementary */
00081 /*          reflectors are stored (see also Further Details): */
00082 /*          = 'C': columnwise */
00083 /*          = 'R': rowwise */
00084 
00085 /*  N       (input) INTEGER */
00086 /*          The order of the block reflector H. N >= 0. */
00087 
00088 /*  K       (input) INTEGER */
00089 /*          The order of the triangular factor T (= the number of */
00090 /*          elementary reflectors). K >= 1. */
00091 
00092 /*  V       (input/output) COMPLEX*16 array, dimension */
00093 /*                               (LDV,K) if STOREV = 'C' */
00094 /*                               (LDV,N) if STOREV = 'R' */
00095 /*          The matrix V. See further details. */
00096 
00097 /*  LDV     (input) INTEGER */
00098 /*          The leading dimension of the array V. */
00099 /*          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. */
00100 
00101 /*  TAU     (input) COMPLEX*16 array, dimension (K) */
00102 /*          TAU(i) must contain the scalar factor of the elementary */
00103 /*          reflector H(i). */
00104 
00105 /*  T       (output) COMPLEX*16 array, dimension (LDT,K) */
00106 /*          The k by k triangular factor T of the block reflector. */
00107 /*          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is */
00108 /*          lower triangular. The rest of the array is not used. */
00109 
00110 /*  LDT     (input) INTEGER */
00111 /*          The leading dimension of the array T. LDT >= K. */
00112 
00113 /*  Further Details */
00114 /*  =============== */
00115 
00116 /*  The shape of the matrix V and the storage of the vectors which define */
00117 /*  the H(i) is best illustrated by the following example with n = 5 and */
00118 /*  k = 3. The elements equal to 1 are not stored; the corresponding */
00119 /*  array elements are modified but restored on exit. The rest of the */
00120 /*  array is not used. */
00121 
00122 /*  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R': */
00123 
00124 /*               V = (  1       )                 V = (  1 v1 v1 v1 v1 ) */
00125 /*                   ( v1  1    )                     (     1 v2 v2 v2 ) */
00126 /*                   ( v1 v2  1 )                     (        1 v3 v3 ) */
00127 /*                   ( v1 v2 v3 ) */
00128 /*                   ( v1 v2 v3 ) */
00129 
00130 /*  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R': */
00131 
00132 /*               V = ( v1 v2 v3 )                 V = ( v1 v1  1       ) */
00133 /*                   ( v1 v2 v3 )                     ( v2 v2 v2  1    ) */
00134 /*                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 ) */
00135 /*                   (     1 v3 ) */
00136 /*                   (        1 ) */
00137 
00138 /*  ===================================================================== */
00139 
00140 /*     .. Parameters .. */
00141 /*     .. */
00142 /*     .. Local Scalars .. */
00143 /*     .. */
00144 /*     .. External Subroutines .. */
00145 /*     .. */
00146 /*     .. External Functions .. */
00147 /*     .. */
00148 /*     .. Executable Statements .. */
00149 
00150 /*     Quick return if possible */
00151 
00152     /* Parameter adjustments */
00153     v_dim1 = *ldv;
00154     v_offset = 1 + v_dim1;
00155     v -= v_offset;
00156     --tau;
00157     t_dim1 = *ldt;
00158     t_offset = 1 + t_dim1;
00159     t -= t_offset;
00160 
00161     /* Function Body */
00162     if (*n == 0) {
00163         return 0;
00164     }
00165 
00166     if (lsame_(direct, "F")) {
00167         prevlastv = *n;
00168         i__1 = *k;
00169         for (i__ = 1; i__ <= i__1; ++i__) {
00170             prevlastv = max(prevlastv,i__);
00171             i__2 = i__;
00172             if (tau[i__2].r == 0. && tau[i__2].i == 0.) {
00173 
00174 /*              H(i)  =  I */
00175 
00176                 i__2 = i__;
00177                 for (j = 1; j <= i__2; ++j) {
00178                     i__3 = j + i__ * t_dim1;
00179                     t[i__3].r = 0., t[i__3].i = 0.;
00180 /* L10: */
00181                 }
00182             } else {
00183 
00184 /*              general case */
00185 
00186                 i__2 = i__ + i__ * v_dim1;
00187                 vii.r = v[i__2].r, vii.i = v[i__2].i;
00188                 i__2 = i__ + i__ * v_dim1;
00189                 v[i__2].r = 1., v[i__2].i = 0.;
00190                 if (lsame_(storev, "C")) {
00191 /*                 Skip any trailing zeros. */
00192                     i__2 = i__ + 1;
00193                     for (lastv = *n; lastv >= i__2; --lastv) {
00194                         i__3 = lastv + i__ * v_dim1;
00195                         if (v[i__3].r != 0. || v[i__3].i != 0.) {
00196                             break;
00197                         }
00198                     }
00199                     j = min(lastv,prevlastv);
00200 
00201 /*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)' * V(i:j,i) */
00202 
00203                     i__2 = j - i__ + 1;
00204                     i__3 = i__ - 1;
00205                     i__4 = i__;
00206                     z__1.r = -tau[i__4].r, z__1.i = -tau[i__4].i;
00207                     zgemv_("Conjugate transpose", &i__2, &i__3, &z__1, &v[i__ 
00208                             + v_dim1], ldv, &v[i__ + i__ * v_dim1], &c__1, &
00209                             c_b2, &t[i__ * t_dim1 + 1], &c__1);
00210                 } else {
00211 /*                 Skip any trailing zeros. */
00212                     i__2 = i__ + 1;
00213                     for (lastv = *n; lastv >= i__2; --lastv) {
00214                         i__3 = i__ + lastv * v_dim1;
00215                         if (v[i__3].r != 0. || v[i__3].i != 0.) {
00216                             break;
00217                         }
00218                     }
00219                     j = min(lastv,prevlastv);
00220 
00221 /*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)' */
00222 
00223                     if (i__ < j) {
00224                         i__2 = j - i__;
00225                         zlacgv_(&i__2, &v[i__ + (i__ + 1) * v_dim1], ldv);
00226                     }
00227                     i__2 = i__ - 1;
00228                     i__3 = j - i__ + 1;
00229                     i__4 = i__;
00230                     z__1.r = -tau[i__4].r, z__1.i = -tau[i__4].i;
00231                     zgemv_("No transpose", &i__2, &i__3, &z__1, &v[i__ * 
00232                             v_dim1 + 1], ldv, &v[i__ + i__ * v_dim1], ldv, &
00233                             c_b2, &t[i__ * t_dim1 + 1], &c__1);
00234                     if (i__ < j) {
00235                         i__2 = j - i__;
00236                         zlacgv_(&i__2, &v[i__ + (i__ + 1) * v_dim1], ldv);
00237                     }
00238                 }
00239                 i__2 = i__ + i__ * v_dim1;
00240                 v[i__2].r = vii.r, v[i__2].i = vii.i;
00241 
00242 /*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) */
00243 
00244                 i__2 = i__ - 1;
00245                 ztrmv_("Upper", "No transpose", "Non-unit", &i__2, &t[
00246                         t_offset], ldt, &t[i__ * t_dim1 + 1], &c__1);
00247                 i__2 = i__ + i__ * t_dim1;
00248                 i__3 = i__;
00249                 t[i__2].r = tau[i__3].r, t[i__2].i = tau[i__3].i;
00250                 if (i__ > 1) {
00251                     prevlastv = max(prevlastv,lastv);
00252                 } else {
00253                     prevlastv = lastv;
00254                 }
00255             }
00256 /* L20: */
00257         }
00258     } else {
00259         prevlastv = 1;
00260         for (i__ = *k; i__ >= 1; --i__) {
00261             i__1 = i__;
00262             if (tau[i__1].r == 0. && tau[i__1].i == 0.) {
00263 
00264 /*              H(i)  =  I */
00265 
00266                 i__1 = *k;
00267                 for (j = i__; j <= i__1; ++j) {
00268                     i__2 = j + i__ * t_dim1;
00269                     t[i__2].r = 0., t[i__2].i = 0.;
00270 /* L30: */
00271                 }
00272             } else {
00273 
00274 /*              general case */
00275 
00276                 if (i__ < *k) {
00277                     if (lsame_(storev, "C")) {
00278                         i__1 = *n - *k + i__ + i__ * v_dim1;
00279                         vii.r = v[i__1].r, vii.i = v[i__1].i;
00280                         i__1 = *n - *k + i__ + i__ * v_dim1;
00281                         v[i__1].r = 1., v[i__1].i = 0.;
00282 /*                    Skip any leading zeros. */
00283                         i__1 = i__ - 1;
00284                         for (lastv = 1; lastv <= i__1; ++lastv) {
00285                             i__2 = lastv + i__ * v_dim1;
00286                             if (v[i__2].r != 0. || v[i__2].i != 0.) {
00287                                 break;
00288                             }
00289                         }
00290                         j = max(lastv,prevlastv);
00291 
00292 /*                    T(i+1:k,i) := */
00293 /*                            - tau(i) * V(j:n-k+i,i+1:k)' * V(j:n-k+i,i) */
00294 
00295                         i__1 = *n - *k + i__ - j + 1;
00296                         i__2 = *k - i__;
00297                         i__3 = i__;
00298                         z__1.r = -tau[i__3].r, z__1.i = -tau[i__3].i;
00299                         zgemv_("Conjugate transpose", &i__1, &i__2, &z__1, &v[
00300                                 j + (i__ + 1) * v_dim1], ldv, &v[j + i__ * 
00301                                 v_dim1], &c__1, &c_b2, &t[i__ + 1 + i__ * 
00302                                 t_dim1], &c__1);
00303                         i__1 = *n - *k + i__ + i__ * v_dim1;
00304                         v[i__1].r = vii.r, v[i__1].i = vii.i;
00305                     } else {
00306                         i__1 = i__ + (*n - *k + i__) * v_dim1;
00307                         vii.r = v[i__1].r, vii.i = v[i__1].i;
00308                         i__1 = i__ + (*n - *k + i__) * v_dim1;
00309                         v[i__1].r = 1., v[i__1].i = 0.;
00310 /*                    Skip any leading zeros. */
00311                         i__1 = i__ - 1;
00312                         for (lastv = 1; lastv <= i__1; ++lastv) {
00313                             i__2 = i__ + lastv * v_dim1;
00314                             if (v[i__2].r != 0. || v[i__2].i != 0.) {
00315                                 break;
00316                             }
00317                         }
00318                         j = max(lastv,prevlastv);
00319 
00320 /*                    T(i+1:k,i) := */
00321 /*                            - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)' */
00322 
00323                         i__1 = *n - *k + i__ - 1 - j + 1;
00324                         zlacgv_(&i__1, &v[i__ + j * v_dim1], ldv);
00325                         i__1 = *k - i__;
00326                         i__2 = *n - *k + i__ - j + 1;
00327                         i__3 = i__;
00328                         z__1.r = -tau[i__3].r, z__1.i = -tau[i__3].i;
00329                         zgemv_("No transpose", &i__1, &i__2, &z__1, &v[i__ + 
00330                                 1 + j * v_dim1], ldv, &v[i__ + j * v_dim1], 
00331                                 ldv, &c_b2, &t[i__ + 1 + i__ * t_dim1], &c__1);
00332                         i__1 = *n - *k + i__ - 1 - j + 1;
00333                         zlacgv_(&i__1, &v[i__ + j * v_dim1], ldv);
00334                         i__1 = i__ + (*n - *k + i__) * v_dim1;
00335                         v[i__1].r = vii.r, v[i__1].i = vii.i;
00336                     }
00337 
00338 /*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) */
00339 
00340                     i__1 = *k - i__;
00341                     ztrmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__ 
00342                             + 1 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ *
00343                              t_dim1], &c__1)
00344                             ;
00345                     if (i__ > 1) {
00346                         prevlastv = min(prevlastv,lastv);
00347                     } else {
00348                         prevlastv = lastv;
00349                     }
00350                 }
00351                 i__1 = i__ + i__ * t_dim1;
00352                 i__2 = i__;
00353                 t[i__1].r = tau[i__2].r, t[i__1].i = tau[i__2].i;
00354             }
00355 /* L40: */
00356         }
00357     }
00358     return 0;
00359 
00360 /*     End of ZLARFT */
00361 
00362 } /* zlarft_ */