clarzt.c
Go to the documentation of this file.
00001 /* clarzt.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_b1 = {0.f,0.f};
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int clarzt_(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;
00026     complex q__1;
00027 
00028     /* Local variables */
00029     integer i__, j, info;
00030     extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
00031 , complex *, integer *, complex *, integer *, complex *, complex *
00032 , integer *);
00033     extern logical lsame_(char *, char *);
00034     extern /* Subroutine */ int ctrmv_(char *, char *, char *, integer *, 
00035             complex *, integer *, complex *, integer *), clacgv_(integer *, complex *, integer *), xerbla_(char *, 
00036              integer *);
00037 
00038 
00039 /*  -- LAPACK routine (version 3.2) -- */
00040 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00041 /*     November 2006 */
00042 
00043 /*     .. Scalar Arguments .. */
00044 /*     .. */
00045 /*     .. Array Arguments .. */
00046 /*     .. */
00047 
00048 /*  Purpose */
00049 /*  ======= */
00050 
00051 /*  CLARZT forms the triangular factor T of a complex block reflector */
00052 /*  H of order > n, which is defined as a product of k elementary */
00053 /*  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 /*  Currently, only STOREV = 'R' and DIRECT = 'B' are supported. */
00070 
00071 /*  Arguments */
00072 /*  ========= */
00073 
00074 /*  DIRECT  (input) CHARACTER*1 */
00075 /*          Specifies the order in which the elementary reflectors are */
00076 /*          multiplied to form the block reflector: */
00077 /*          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet) */
00078 /*          = 'B': H = H(k) . . . H(2) H(1) (Backward) */
00079 
00080 /*  STOREV  (input) CHARACTER*1 */
00081 /*          Specifies how the vectors which define the elementary */
00082 /*          reflectors are stored (see also Further Details): */
00083 /*          = 'C': columnwise                        (not supported yet) */
00084 /*          = 'R': rowwise */
00085 
00086 /*  N       (input) INTEGER */
00087 /*          The order of the block reflector H. N >= 0. */
00088 
00089 /*  K       (input) INTEGER */
00090 /*          The order of the triangular factor T (= the number of */
00091 /*          elementary reflectors). K >= 1. */
00092 
00093 /*  V       (input/output) COMPLEX array, dimension */
00094 /*                               (LDV,K) if STOREV = 'C' */
00095 /*                               (LDV,N) if STOREV = 'R' */
00096 /*          The matrix V. See further details. */
00097 
00098 /*  LDV     (input) INTEGER */
00099 /*          The leading dimension of the array V. */
00100 /*          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. */
00101 
00102 /*  TAU     (input) COMPLEX array, dimension (K) */
00103 /*          TAU(i) must contain the scalar factor of the elementary */
00104 /*          reflector H(i). */
00105 
00106 /*  T       (output) COMPLEX array, dimension (LDT,K) */
00107 /*          The k by k triangular factor T of the block reflector. */
00108 /*          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is */
00109 /*          lower triangular. The rest of the array is not used. */
00110 
00111 /*  LDT     (input) INTEGER */
00112 /*          The leading dimension of the array T. LDT >= K. */
00113 
00114 /*  Further Details */
00115 /*  =============== */
00116 
00117 /*  Based on contributions by */
00118 /*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA */
00119 
00120 /*  The shape of the matrix V and the storage of the vectors which define */
00121 /*  the H(i) is best illustrated by the following example with n = 5 and */
00122 /*  k = 3. The elements equal to 1 are not stored; the corresponding */
00123 /*  array elements are modified but restored on exit. The rest of the */
00124 /*  array is not used. */
00125 
00126 /*  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R': */
00127 
00128 /*                                              ______V_____ */
00129 /*         ( v1 v2 v3 )                        /            \ */
00130 /*         ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 ) */
00131 /*     V = ( v1 v2 v3 )                      ( v2 v2 v2 v2 v2 . . . 1   ) */
00132 /*         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     ) */
00133 /*         ( v1 v2 v3 ) */
00134 /*            .  .  . */
00135 /*            .  .  . */
00136 /*            1  .  . */
00137 /*               1  . */
00138 /*                  1 */
00139 
00140 /*  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R': */
00141 
00142 /*                                                        ______V_____ */
00143 /*            1                                          /            \ */
00144 /*            .  1                           ( 1 . . . . v1 v1 v1 v1 v1 ) */
00145 /*            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 ) */
00146 /*            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 ) */
00147 /*            .  .  . */
00148 /*         ( v1 v2 v3 ) */
00149 /*         ( v1 v2 v3 ) */
00150 /*     V = ( v1 v2 v3 ) */
00151 /*         ( v1 v2 v3 ) */
00152 /*         ( v1 v2 v3 ) */
00153 
00154 /*  ===================================================================== */
00155 
00156 /*     .. Parameters .. */
00157 /*     .. */
00158 /*     .. Local Scalars .. */
00159 /*     .. */
00160 /*     .. External Subroutines .. */
00161 /*     .. */
00162 /*     .. External Functions .. */
00163 /*     .. */
00164 /*     .. Executable Statements .. */
00165 
00166 /*     Check for currently supported options */
00167 
00168     /* Parameter adjustments */
00169     v_dim1 = *ldv;
00170     v_offset = 1 + v_dim1;
00171     v -= v_offset;
00172     --tau;
00173     t_dim1 = *ldt;
00174     t_offset = 1 + t_dim1;
00175     t -= t_offset;
00176 
00177     /* Function Body */
00178     info = 0;
00179     if (! lsame_(direct, "B")) {
00180         info = -1;
00181     } else if (! lsame_(storev, "R")) {
00182         info = -2;
00183     }
00184     if (info != 0) {
00185         i__1 = -info;
00186         xerbla_("CLARZT", &i__1);
00187         return 0;
00188     }
00189 
00190     for (i__ = *k; i__ >= 1; --i__) {
00191         i__1 = i__;
00192         if (tau[i__1].r == 0.f && tau[i__1].i == 0.f) {
00193 
00194 /*           H(i)  =  I */
00195 
00196             i__1 = *k;
00197             for (j = i__; j <= i__1; ++j) {
00198                 i__2 = j + i__ * t_dim1;
00199                 t[i__2].r = 0.f, t[i__2].i = 0.f;
00200 /* L10: */
00201             }
00202         } else {
00203 
00204 /*           general case */
00205 
00206             if (i__ < *k) {
00207 
00208 /*              T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)' */
00209 
00210                 clacgv_(n, &v[i__ + v_dim1], ldv);
00211                 i__1 = *k - i__;
00212                 i__2 = i__;
00213                 q__1.r = -tau[i__2].r, q__1.i = -tau[i__2].i;
00214                 cgemv_("No transpose", &i__1, n, &q__1, &v[i__ + 1 + v_dim1], 
00215                         ldv, &v[i__ + v_dim1], ldv, &c_b1, &t[i__ + 1 + i__ * 
00216                         t_dim1], &c__1);
00217                 clacgv_(n, &v[i__ + v_dim1], ldv);
00218 
00219 /*              T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i) */
00220 
00221                 i__1 = *k - i__;
00222                 ctrmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__ + 1 
00223                         + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ * t_dim1]
00224 , &c__1);
00225             }
00226             i__1 = i__ + i__ * t_dim1;
00227             i__2 = i__;
00228             t[i__1].r = tau[i__2].r, t[i__1].i = tau[i__2].i;
00229         }
00230 /* L20: */
00231     }
00232     return 0;
00233 
00234 /*     End of CLARZT */
00235 
00236 } /* clarzt_ */


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