swiftnav: clarft.c Source File

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