swiftnav: clarfb.c Source File

Go to the documentation of this file.
00001 /* clarfb.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 = {1.f,0.f};
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int clarfb_(char *side, char *trans, char *direct, char *
00022         storev, integer *m, integer *n, integer *k, complex *v, integer *ldv, 
00023         complex *t, integer *ldt, complex *c__, integer *ldc, complex *work, 
00024         integer *ldwork)
00025 {
00026     /* System generated locals */
00027     integer c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1, 
00028             work_offset, i__1, i__2, i__3, i__4, i__5;
00029     complex q__1, q__2;
00030 
00031     /* Builtin functions */
00032     void r_cnjg(complex *, complex *);
00033 
00034     /* Local variables */
00035     integer i__, j;
00036     extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *, 
00037             integer *, complex *, complex *, integer *, complex *, integer *, 
00038             complex *, complex *, integer *);
00039     extern logical lsame_(char *, char *);
00040     integer lastc;
00041     extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
00042             complex *, integer *), ctrmm_(char *, char *, char *, char *, 
00043             integer *, integer *, complex *, complex *, integer *, complex *, 
00044             integer *);
00045     integer lastv;
00046     extern integer ilaclc_(integer *, integer *, complex *, integer *);
00047     extern /* Subroutine */ int clacgv_(integer *, complex *, integer *);
00048     extern integer ilaclr_(integer *, integer *, complex *, integer *);
00049     char transt[1];
00050 
00051 
00052 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00053 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00054 /*     November 2006 */
00055 
00056 /*     .. Scalar Arguments .. */
00057 /*     .. */
00058 /*     .. Array Arguments .. */
00059 /*     .. */
00060 
00061 /*  Purpose */
00062 /*  ======= */
00063 
00064 /*  CLARFB applies a complex block reflector H or its transpose H' to a */
00065 /*  complex M-by-N matrix C, from either the left or the right. */
00066 
00067 /*  Arguments */
00068 /*  ========= */
00069 
00070 /*  SIDE    (input) CHARACTER*1 */
00071 /*          = 'L': apply H or H' from the Left */
00072 /*          = 'R': apply H or H' from the Right */
00073 
00074 /*  TRANS   (input) CHARACTER*1 */
00075 /*          = 'N': apply H (No transpose) */
00076 /*          = 'C': apply H' (Conjugate transpose) */
00077 
00078 /*  DIRECT  (input) CHARACTER*1 */
00079 /*          Indicates how H is formed from a product of elementary */
00080 /*          reflectors */
00081 /*          = 'F': H = H(1) H(2) . . . H(k) (Forward) */
00082 /*          = 'B': H = H(k) . . . H(2) H(1) (Backward) */
00083 
00084 /*  STOREV  (input) CHARACTER*1 */
00085 /*          Indicates how the vectors which define the elementary */
00086 /*          reflectors are stored: */
00087 /*          = 'C': Columnwise */
00088 /*          = 'R': Rowwise */
00089 
00090 /*  M       (input) INTEGER */
00091 /*          The number of rows of the matrix C. */
00092 
00093 /*  N       (input) INTEGER */
00094 /*          The number of columns of the matrix C. */
00095 
00096 /*  K       (input) INTEGER */
00097 /*          The order of the matrix T (= the number of elementary */
00098 /*          reflectors whose product defines the block reflector). */
00099 
00100 /*  V       (input) COMPLEX array, dimension */
00101 /*                                (LDV,K) if STOREV = 'C' */
00102 /*                                (LDV,M) if STOREV = 'R' and SIDE = 'L' */
00103 /*                                (LDV,N) if STOREV = 'R' and SIDE = 'R' */
00104 /*          The matrix V. See further details. */
00105 
00106 /*  LDV     (input) INTEGER */
00107 /*          The leading dimension of the array V. */
00108 /*          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); */
00109 /*          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); */
00110 /*          if STOREV = 'R', LDV >= K. */
00111 
00112 /*  T       (input) COMPLEX array, dimension (LDT,K) */
00113 /*          The triangular K-by-K matrix T in the representation of the */
00114 /*          block reflector. */
00115 
00116 /*  LDT     (input) INTEGER */
00117 /*          The leading dimension of the array T. LDT >= K. */
00118 
00119 /*  C       (input/output) COMPLEX array, dimension (LDC,N) */
00120 /*          On entry, the M-by-N matrix C. */
00121 /*          On exit, C is overwritten by H*C or H'*C or C*H or C*H'. */
00122 
00123 /*  LDC     (input) INTEGER */
00124 /*          The leading dimension of the array C. LDC >= max(1,M). */
00125 
00126 /*  WORK    (workspace) COMPLEX array, dimension (LDWORK,K) */
00127 
00128 /*  LDWORK  (input) INTEGER */
00129 /*          The leading dimension of the array WORK. */
00130 /*          If SIDE = 'L', LDWORK >= max(1,N); */
00131 /*          if SIDE = 'R', LDWORK >= max(1,M). */
00132 
00133 /*  ===================================================================== */
00134 
00135 /*     .. Parameters .. */
00136 /*     .. */
00137 /*     .. Local Scalars .. */
00138 /*     .. */
00139 /*     .. External Functions .. */
00140 /*     .. */
00141 /*     .. External Subroutines .. */
00142 /*     .. */
00143 /*     .. Intrinsic Functions .. */
00144 /*     .. */
00145 /*     .. Executable Statements .. */
00146 
00147 /*     Quick return if possible */
00148 
00149     /* Parameter adjustments */
00150     v_dim1 = *ldv;
00151     v_offset = 1 + v_dim1;
00152     v -= v_offset;
00153     t_dim1 = *ldt;
00154     t_offset = 1 + t_dim1;
00155     t -= t_offset;
00156     c_dim1 = *ldc;
00157     c_offset = 1 + c_dim1;
00158     c__ -= c_offset;
00159     work_dim1 = *ldwork;
00160     work_offset = 1 + work_dim1;
00161     work -= work_offset;
00162 
00163     /* Function Body */
00164     if (*m <= 0 || *n <= 0) {
00165         return 0;
00166     }
00167 
00168     if (lsame_(trans, "N")) {
00169         *(unsigned char *)transt = 'C';
00170     } else {
00171         *(unsigned char *)transt = 'N';
00172     }
00173 
00174     if (lsame_(storev, "C")) {
00175 
00176         if (lsame_(direct, "F")) {
00177 
00178 /*           Let  V =  ( V1 )    (first K rows) */
00179 /*                     ( V2 ) */
00180 /*           where  V1  is unit lower triangular. */
00181 
00182             if (lsame_(side, "L")) {
00183 
00184 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
00185 /*                                                  ( C2 ) */
00186 
00187 /* Computing MAX */
00188                 i__1 = *k, i__2 = ilaclr_(m, k, &v[v_offset], ldv);
00189                 lastv = max(i__1,i__2);
00190                 lastc = ilaclc_(&lastv, n, &c__[c_offset], ldc);
00191 
00192 /*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK) */
00193 
00194 /*              W := C1' */
00195 
00196                 i__1 = *k;
00197                 for (j = 1; j <= i__1; ++j) {
00198                     ccopy_(&lastc, &c__[j + c_dim1], ldc, &work[j * work_dim1 
00199                             + 1], &c__1);
00200                     clacgv_(&lastc, &work[j * work_dim1 + 1], &c__1);
00201 /* L10: */
00202                 }
00203 
00204 /*              W := W * V1 */
00205 
00206                 ctrmm_("Right", "Lower", "No transpose", "Unit", &lastc, k, &
00207                         c_b1, &v[v_offset], ldv, &work[work_offset], ldwork);
00208                 if (lastv > *k) {
00209 
00210 /*                 W := W + C2'*V2 */
00211 
00212                     i__1 = lastv - *k;
00213                     cgemm_("Conjugate transpose", "No transpose", &lastc, k, &
00214                             i__1, &c_b1, &c__[*k + 1 + c_dim1], ldc, &v[*k + 
00215                             1 + v_dim1], ldv, &c_b1, &work[work_offset], 
00216                             ldwork);
00217                 }
00218 
00219 /*              W := W * T'  or  W * T */
00220 
00221                 ctrmm_("Right", "Upper", transt, "Non-unit", &lastc, k, &c_b1, 
00222                          &t[t_offset], ldt, &work[work_offset], ldwork);
00223 
00224 /*              C := C - V * W' */
00225 
00226                 if (*m > *k) {
00227 
00228 /*                 C2 := C2 - V2 * W' */
00229 
00230                     i__1 = lastv - *k;
00231                     q__1.r = -1.f, q__1.i = -0.f;
00232                     cgemm_("No transpose", "Conjugate transpose", &i__1, &
00233                             lastc, k, &q__1, &v[*k + 1 + v_dim1], ldv, &work[
00234                             work_offset], ldwork, &c_b1, &c__[*k + 1 + c_dim1]
00235 , ldc);
00236                 }
00237 
00238 /*              W := W * V1' */
00239 
00240                 ctrmm_("Right", "Lower", "Conjugate transpose", "Unit", &
00241                         lastc, k, &c_b1, &v[v_offset], ldv, &work[work_offset]
00242 , ldwork)
00243                         ;
00244 
00245 /*              C1 := C1 - W' */
00246 
00247                 i__1 = *k;
00248                 for (j = 1; j <= i__1; ++j) {
00249                     i__2 = lastc;
00250                     for (i__ = 1; i__ <= i__2; ++i__) {
00251                         i__3 = j + i__ * c_dim1;
00252                         i__4 = j + i__ * c_dim1;
00253                         r_cnjg(&q__2, &work[i__ + j * work_dim1]);
00254                         q__1.r = c__[i__4].r - q__2.r, q__1.i = c__[i__4].i - 
00255                                 q__2.i;
00256                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00257 /* L20: */
00258                     }
00259 /* L30: */
00260                 }
00261 
00262             } else if (lsame_(side, "R")) {
00263 
00264 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
00265 
00266 /* Computing MAX */
00267                 i__1 = *k, i__2 = ilaclr_(n, k, &v[v_offset], ldv);
00268                 lastv = max(i__1,i__2);
00269                 lastc = ilaclr_(m, &lastv, &c__[c_offset], ldc);
00270 
00271 /*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK) */
00272 
00273 /*              W := C1 */
00274 
00275                 i__1 = *k;
00276                 for (j = 1; j <= i__1; ++j) {
00277                     ccopy_(&lastc, &c__[j * c_dim1 + 1], &c__1, &work[j * 
00278                             work_dim1 + 1], &c__1);
00279 /* L40: */
00280                 }
00281 
00282 /*              W := W * V1 */
00283 
00284                 ctrmm_("Right", "Lower", "No transpose", "Unit", &lastc, k, &
00285                         c_b1, &v[v_offset], ldv, &work[work_offset], ldwork);
00286                 if (lastv > *k) {
00287 
00288 /*                 W := W + C2 * V2 */
00289 
00290                     i__1 = lastv - *k;
00291                     cgemm_("No transpose", "No transpose", &lastc, k, &i__1, &
00292                             c_b1, &c__[(*k + 1) * c_dim1 + 1], ldc, &v[*k + 1 
00293                             + v_dim1], ldv, &c_b1, &work[work_offset], ldwork);
00294                 }
00295 
00296 /*              W := W * T  or  W * T' */
00297 
00298                 ctrmm_("Right", "Upper", trans, "Non-unit", &lastc, k, &c_b1, 
00299                         &t[t_offset], ldt, &work[work_offset], ldwork);
00300 
00301 /*              C := C - W * V' */
00302 
00303                 if (lastv > *k) {
00304 
00305 /*                 C2 := C2 - W * V2' */
00306 
00307                     i__1 = lastv - *k;
00308                     q__1.r = -1.f, q__1.i = -0.f;
00309                     cgemm_("No transpose", "Conjugate transpose", &lastc, &
00310                             i__1, k, &q__1, &work[work_offset], ldwork, &v[*k 
00311                             + 1 + v_dim1], ldv, &c_b1, &c__[(*k + 1) * c_dim1 
00312                             + 1], ldc);
00313                 }
00314 
00315 /*              W := W * V1' */
00316 
00317                 ctrmm_("Right", "Lower", "Conjugate transpose", "Unit", &
00318                         lastc, k, &c_b1, &v[v_offset], ldv, &work[work_offset]
00319 , ldwork)
00320                         ;
00321 
00322 /*              C1 := C1 - W */
00323 
00324                 i__1 = *k;
00325                 for (j = 1; j <= i__1; ++j) {
00326                     i__2 = lastc;
00327                     for (i__ = 1; i__ <= i__2; ++i__) {
00328                         i__3 = i__ + j * c_dim1;
00329                         i__4 = i__ + j * c_dim1;
00330                         i__5 = i__ + j * work_dim1;
00331                         q__1.r = c__[i__4].r - work[i__5].r, q__1.i = c__[
00332                                 i__4].i - work[i__5].i;
00333                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00334 /* L50: */
00335                     }
00336 /* L60: */
00337                 }
00338             }
00339 
00340         } else {
00341 
00342 /*           Let  V =  ( V1 ) */
00343 /*                     ( V2 )    (last K rows) */
00344 /*           where  V2  is unit upper triangular. */
00345 
00346             if (lsame_(side, "L")) {
00347 
00348 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
00349 /*                                                  ( C2 ) */
00350 
00351 /* Computing MAX */
00352                 i__1 = *k, i__2 = ilaclr_(m, k, &v[v_offset], ldv);
00353                 lastv = max(i__1,i__2);
00354                 lastc = ilaclc_(&lastv, n, &c__[c_offset], ldc);
00355 
00356 /*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK) */
00357 
00358 /*              W := C2' */
00359 
00360                 i__1 = *k;
00361                 for (j = 1; j <= i__1; ++j) {
00362                     ccopy_(&lastc, &c__[lastv - *k + j + c_dim1], ldc, &work[
00363                             j * work_dim1 + 1], &c__1);
00364                     clacgv_(&lastc, &work[j * work_dim1 + 1], &c__1);
00365 /* L70: */
00366                 }
00367 
00368 /*              W := W * V2 */
00369 
00370                 ctrmm_("Right", "Upper", "No transpose", "Unit", &lastc, k, &
00371                         c_b1, &v[lastv - *k + 1 + v_dim1], ldv, &work[
00372                         work_offset], ldwork);
00373                 if (lastv > *k) {
00374 
00375 /*                 W := W + C1'*V1 */
00376 
00377                     i__1 = lastv - *k;
00378                     cgemm_("Conjugate transpose", "No transpose", &lastc, k, &
00379                             i__1, &c_b1, &c__[c_offset], ldc, &v[v_offset], 
00380                             ldv, &c_b1, &work[work_offset], ldwork);
00381                 }
00382 
00383 /*              W := W * T'  or  W * T */
00384 
00385                 ctrmm_("Right", "Lower", transt, "Non-unit", &lastc, k, &c_b1, 
00386                          &t[t_offset], ldt, &work[work_offset], ldwork);
00387 
00388 /*              C := C - V * W' */
00389 
00390                 if (lastv > *k) {
00391 
00392 /*                 C1 := C1 - V1 * W' */
00393 
00394                     i__1 = lastv - *k;
00395                     q__1.r = -1.f, q__1.i = -0.f;
00396                     cgemm_("No transpose", "Conjugate transpose", &i__1, &
00397                             lastc, k, &q__1, &v[v_offset], ldv, &work[
00398                             work_offset], ldwork, &c_b1, &c__[c_offset], ldc);
00399                 }
00400 
00401 /*              W := W * V2' */
00402 
00403                 ctrmm_("Right", "Upper", "Conjugate transpose", "Unit", &
00404                         lastc, k, &c_b1, &v[lastv - *k + 1 + v_dim1], ldv, &
00405                         work[work_offset], ldwork);
00406 
00407 /*              C2 := C2 - W' */
00408 
00409                 i__1 = *k;
00410                 for (j = 1; j <= i__1; ++j) {
00411                     i__2 = lastc;
00412                     for (i__ = 1; i__ <= i__2; ++i__) {
00413                         i__3 = lastv - *k + j + i__ * c_dim1;
00414                         i__4 = lastv - *k + j + i__ * c_dim1;
00415                         r_cnjg(&q__2, &work[i__ + j * work_dim1]);
00416                         q__1.r = c__[i__4].r - q__2.r, q__1.i = c__[i__4].i - 
00417                                 q__2.i;
00418                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00419 /* L80: */
00420                     }
00421 /* L90: */
00422                 }
00423 
00424             } else if (lsame_(side, "R")) {
00425 
00426 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
00427 
00428 /* Computing MAX */
00429                 i__1 = *k, i__2 = ilaclr_(n, k, &v[v_offset], ldv);
00430                 lastv = max(i__1,i__2);
00431                 lastc = ilaclr_(m, &lastv, &c__[c_offset], ldc);
00432 
00433 /*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK) */
00434 
00435 /*              W := C2 */
00436 
00437                 i__1 = *k;
00438                 for (j = 1; j <= i__1; ++j) {
00439                     ccopy_(&lastc, &c__[(lastv - *k + j) * c_dim1 + 1], &c__1, 
00440                              &work[j * work_dim1 + 1], &c__1);
00441 /* L100: */
00442                 }
00443 
00444 /*              W := W * V2 */
00445 
00446                 ctrmm_("Right", "Upper", "No transpose", "Unit", &lastc, k, &
00447                         c_b1, &v[lastv - *k + 1 + v_dim1], ldv, &work[
00448                         work_offset], ldwork);
00449                 if (lastv > *k) {
00450 
00451 /*                 W := W + C1 * V1 */
00452 
00453                     i__1 = lastv - *k;
00454                     cgemm_("No transpose", "No transpose", &lastc, k, &i__1, &
00455                             c_b1, &c__[c_offset], ldc, &v[v_offset], ldv, &
00456                             c_b1, &work[work_offset], ldwork);
00457                 }
00458 
00459 /*              W := W * T  or  W * T' */
00460 
00461                 ctrmm_("Right", "Lower", trans, "Non-unit", &lastc, k, &c_b1, 
00462                         &t[t_offset], ldt, &work[work_offset], ldwork);
00463 
00464 /*              C := C - W * V' */
00465 
00466                 if (lastv > *k) {
00467 
00468 /*                 C1 := C1 - W * V1' */
00469 
00470                     i__1 = lastv - *k;
00471                     q__1.r = -1.f, q__1.i = -0.f;
00472                     cgemm_("No transpose", "Conjugate transpose", &lastc, &
00473                             i__1, k, &q__1, &work[work_offset], ldwork, &v[
00474                             v_offset], ldv, &c_b1, &c__[c_offset], ldc);
00475                 }
00476 
00477 /*              W := W * V2' */
00478 
00479                 ctrmm_("Right", "Upper", "Conjugate transpose", "Unit", &
00480                         lastc, k, &c_b1, &v[lastv - *k + 1 + v_dim1], ldv, &
00481                         work[work_offset], ldwork);
00482 
00483 /*              C2 := C2 - W */
00484 
00485                 i__1 = *k;
00486                 for (j = 1; j <= i__1; ++j) {
00487                     i__2 = lastc;
00488                     for (i__ = 1; i__ <= i__2; ++i__) {
00489                         i__3 = i__ + (lastv - *k + j) * c_dim1;
00490                         i__4 = i__ + (lastv - *k + j) * c_dim1;
00491                         i__5 = i__ + j * work_dim1;
00492                         q__1.r = c__[i__4].r - work[i__5].r, q__1.i = c__[
00493                                 i__4].i - work[i__5].i;
00494                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00495 /* L110: */
00496                     }
00497 /* L120: */
00498                 }
00499             }
00500         }
00501 
00502     } else if (lsame_(storev, "R")) {
00503 
00504         if (lsame_(direct, "F")) {
00505 
00506 /*           Let  V =  ( V1  V2 )    (V1: first K columns) */
00507 /*           where  V1  is unit upper triangular. */
00508 
00509             if (lsame_(side, "L")) {
00510 
00511 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
00512 /*                                                  ( C2 ) */
00513 
00514 /* Computing MAX */
00515                 i__1 = *k, i__2 = ilaclc_(k, m, &v[v_offset], ldv);
00516                 lastv = max(i__1,i__2);
00517                 lastc = ilaclc_(&lastv, n, &c__[c_offset], ldc);
00518 
00519 /*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK) */
00520 
00521 /*              W := C1' */
00522 
00523                 i__1 = *k;
00524                 for (j = 1; j <= i__1; ++j) {
00525                     ccopy_(&lastc, &c__[j + c_dim1], ldc, &work[j * work_dim1 
00526                             + 1], &c__1);
00527                     clacgv_(&lastc, &work[j * work_dim1 + 1], &c__1);
00528 /* L130: */
00529                 }
00530 
00531 /*              W := W * V1' */
00532 
00533                 ctrmm_("Right", "Upper", "Conjugate transpose", "Unit", &
00534                         lastc, k, &c_b1, &v[v_offset], ldv, &work[work_offset]
00535 , ldwork)
00536                         ;
00537                 if (lastv > *k) {
00538 
00539 /*                 W := W + C2'*V2' */
00540 
00541                     i__1 = lastv - *k;
00542                     cgemm_("Conjugate transpose", "Conjugate transpose", &
00543                             lastc, k, &i__1, &c_b1, &c__[*k + 1 + c_dim1], 
00544                             ldc, &v[(*k + 1) * v_dim1 + 1], ldv, &c_b1, &work[
00545                             work_offset], ldwork);
00546                 }
00547 
00548 /*              W := W * T'  or  W * T */
00549 
00550                 ctrmm_("Right", "Upper", transt, "Non-unit", &lastc, k, &c_b1, 
00551                          &t[t_offset], ldt, &work[work_offset], ldwork);
00552 
00553 /*              C := C - V' * W' */
00554 
00555                 if (lastv > *k) {
00556 
00557 /*                 C2 := C2 - V2' * W' */
00558 
00559                     i__1 = lastv - *k;
00560                     q__1.r = -1.f, q__1.i = -0.f;
00561                     cgemm_("Conjugate transpose", "Conjugate transpose", &
00562                             i__1, &lastc, k, &q__1, &v[(*k + 1) * v_dim1 + 1], 
00563                              ldv, &work[work_offset], ldwork, &c_b1, &c__[*k 
00564                             + 1 + c_dim1], ldc);
00565                 }
00566 
00567 /*              W := W * V1 */
00568 
00569                 ctrmm_("Right", "Upper", "No transpose", "Unit", &lastc, k, &
00570                         c_b1, &v[v_offset], ldv, &work[work_offset], ldwork);
00571 
00572 /*              C1 := C1 - W' */
00573 
00574                 i__1 = *k;
00575                 for (j = 1; j <= i__1; ++j) {
00576                     i__2 = lastc;
00577                     for (i__ = 1; i__ <= i__2; ++i__) {
00578                         i__3 = j + i__ * c_dim1;
00579                         i__4 = j + i__ * c_dim1;
00580                         r_cnjg(&q__2, &work[i__ + j * work_dim1]);
00581                         q__1.r = c__[i__4].r - q__2.r, q__1.i = c__[i__4].i - 
00582                                 q__2.i;
00583                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00584 /* L140: */
00585                     }
00586 /* L150: */
00587                 }
00588 
00589             } else if (lsame_(side, "R")) {
00590 
00591 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
00592 
00593 /* Computing MAX */
00594                 i__1 = *k, i__2 = ilaclc_(k, n, &v[v_offset], ldv);
00595                 lastv = max(i__1,i__2);
00596                 lastc = ilaclr_(m, &lastv, &c__[c_offset], ldc);
00597 
00598 /*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK) */
00599 
00600 /*              W := C1 */
00601 
00602                 i__1 = *k;
00603                 for (j = 1; j <= i__1; ++j) {
00604                     ccopy_(&lastc, &c__[j * c_dim1 + 1], &c__1, &work[j * 
00605                             work_dim1 + 1], &c__1);
00606 /* L160: */
00607                 }
00608 
00609 /*              W := W * V1' */
00610 
00611                 ctrmm_("Right", "Upper", "Conjugate transpose", "Unit", &
00612                         lastc, k, &c_b1, &v[v_offset], ldv, &work[work_offset]
00613 , ldwork)
00614                         ;
00615                 if (lastv > *k) {
00616 
00617 /*                 W := W + C2 * V2' */
00618 
00619                     i__1 = lastv - *k;
00620                     cgemm_("No transpose", "Conjugate transpose", &lastc, k, &
00621                             i__1, &c_b1, &c__[(*k + 1) * c_dim1 + 1], ldc, &v[
00622                             (*k + 1) * v_dim1 + 1], ldv, &c_b1, &work[
00623                             work_offset], ldwork);
00624                 }
00625 
00626 /*              W := W * T  or  W * T' */
00627 
00628                 ctrmm_("Right", "Upper", trans, "Non-unit", &lastc, k, &c_b1, 
00629                         &t[t_offset], ldt, &work[work_offset], ldwork);
00630 
00631 /*              C := C - W * V */
00632 
00633                 if (lastv > *k) {
00634 
00635 /*                 C2 := C2 - W * V2 */
00636 
00637                     i__1 = lastv - *k;
00638                     q__1.r = -1.f, q__1.i = -0.f;
00639                     cgemm_("No transpose", "No transpose", &lastc, &i__1, k, &
00640                             q__1, &work[work_offset], ldwork, &v[(*k + 1) * 
00641                             v_dim1 + 1], ldv, &c_b1, &c__[(*k + 1) * c_dim1 + 
00642                             1], ldc);
00643                 }
00644 
00645 /*              W := W * V1 */
00646 
00647                 ctrmm_("Right", "Upper", "No transpose", "Unit", &lastc, k, &
00648                         c_b1, &v[v_offset], ldv, &work[work_offset], ldwork);
00649 
00650 /*              C1 := C1 - W */
00651 
00652                 i__1 = *k;
00653                 for (j = 1; j <= i__1; ++j) {
00654                     i__2 = lastc;
00655                     for (i__ = 1; i__ <= i__2; ++i__) {
00656                         i__3 = i__ + j * c_dim1;
00657                         i__4 = i__ + j * c_dim1;
00658                         i__5 = i__ + j * work_dim1;
00659                         q__1.r = c__[i__4].r - work[i__5].r, q__1.i = c__[
00660                                 i__4].i - work[i__5].i;
00661                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00662 /* L170: */
00663                     }
00664 /* L180: */
00665                 }
00666 
00667             }
00668 
00669         } else {
00670 
00671 /*           Let  V =  ( V1  V2 )    (V2: last K columns) */
00672 /*           where  V2  is unit lower triangular. */
00673 
00674             if (lsame_(side, "L")) {
00675 
00676 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
00677 /*                                                  ( C2 ) */
00678 
00679 /* Computing MAX */
00680                 i__1 = *k, i__2 = ilaclc_(k, m, &v[v_offset], ldv);
00681                 lastv = max(i__1,i__2);
00682                 lastc = ilaclc_(&lastv, n, &c__[c_offset], ldc);
00683 
00684 /*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK) */
00685 
00686 /*              W := C2' */
00687 
00688                 i__1 = *k;
00689                 for (j = 1; j <= i__1; ++j) {
00690                     ccopy_(&lastc, &c__[lastv - *k + j + c_dim1], ldc, &work[
00691                             j * work_dim1 + 1], &c__1);
00692                     clacgv_(&lastc, &work[j * work_dim1 + 1], &c__1);
00693 /* L190: */
00694                 }
00695 
00696 /*              W := W * V2' */
00697 
00698                 ctrmm_("Right", "Lower", "Conjugate transpose", "Unit", &
00699                         lastc, k, &c_b1, &v[(lastv - *k + 1) * v_dim1 + 1], 
00700                         ldv, &work[work_offset], ldwork);
00701                 if (lastv > *k) {
00702 
00703 /*                 W := W + C1'*V1' */
00704 
00705                     i__1 = lastv - *k;
00706                     cgemm_("Conjugate transpose", "Conjugate transpose", &
00707                             lastc, k, &i__1, &c_b1, &c__[c_offset], ldc, &v[
00708                             v_offset], ldv, &c_b1, &work[work_offset], ldwork);
00709                 }
00710 
00711 /*              W := W * T'  or  W * T */
00712 
00713                 ctrmm_("Right", "Lower", transt, "Non-unit", &lastc, k, &c_b1, 
00714                          &t[t_offset], ldt, &work[work_offset], ldwork);
00715 
00716 /*              C := C - V' * W' */
00717 
00718                 if (lastv > *k) {
00719 
00720 /*                 C1 := C1 - V1' * W' */
00721 
00722                     i__1 = lastv - *k;
00723                     q__1.r = -1.f, q__1.i = -0.f;
00724                     cgemm_("Conjugate transpose", "Conjugate transpose", &
00725                             i__1, &lastc, k, &q__1, &v[v_offset], ldv, &work[
00726                             work_offset], ldwork, &c_b1, &c__[c_offset], ldc);
00727                 }
00728 
00729 /*              W := W * V2 */
00730 
00731                 ctrmm_("Right", "Lower", "No transpose", "Unit", &lastc, k, &
00732                         c_b1, &v[(lastv - *k + 1) * v_dim1 + 1], ldv, &work[
00733                         work_offset], ldwork);
00734 
00735 /*              C2 := C2 - W' */
00736 
00737                 i__1 = *k;
00738                 for (j = 1; j <= i__1; ++j) {
00739                     i__2 = lastc;
00740                     for (i__ = 1; i__ <= i__2; ++i__) {
00741                         i__3 = lastv - *k + j + i__ * c_dim1;
00742                         i__4 = lastv - *k + j + i__ * c_dim1;
00743                         r_cnjg(&q__2, &work[i__ + j * work_dim1]);
00744                         q__1.r = c__[i__4].r - q__2.r, q__1.i = c__[i__4].i - 
00745                                 q__2.i;
00746                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00747 /* L200: */
00748                     }
00749 /* L210: */
00750                 }
00751 
00752             } else if (lsame_(side, "R")) {
00753 
00754 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
00755 
00756 /* Computing MAX */
00757                 i__1 = *k, i__2 = ilaclc_(k, n, &v[v_offset], ldv);
00758                 lastv = max(i__1,i__2);
00759                 lastc = ilaclr_(m, &lastv, &c__[c_offset], ldc);
00760 
00761 /*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK) */
00762 
00763 /*              W := C2 */
00764 
00765                 i__1 = *k;
00766                 for (j = 1; j <= i__1; ++j) {
00767                     ccopy_(&lastc, &c__[(lastv - *k + j) * c_dim1 + 1], &c__1, 
00768                              &work[j * work_dim1 + 1], &c__1);
00769 /* L220: */
00770                 }
00771 
00772 /*              W := W * V2' */
00773 
00774                 ctrmm_("Right", "Lower", "Conjugate transpose", "Unit", &
00775                         lastc, k, &c_b1, &v[(lastv - *k + 1) * v_dim1 + 1], 
00776                         ldv, &work[work_offset], ldwork);
00777                 if (lastv > *k) {
00778 
00779 /*                 W := W + C1 * V1' */
00780 
00781                     i__1 = lastv - *k;
00782                     cgemm_("No transpose", "Conjugate transpose", &lastc, k, &
00783                             i__1, &c_b1, &c__[c_offset], ldc, &v[v_offset], 
00784                             ldv, &c_b1, &work[work_offset], ldwork);
00785                 }
00786 
00787 /*              W := W * T  or  W * T' */
00788 
00789                 ctrmm_("Right", "Lower", trans, "Non-unit", &lastc, k, &c_b1, 
00790                         &t[t_offset], ldt, &work[work_offset], ldwork);
00791 
00792 /*              C := C - W * V */
00793 
00794                 if (lastv > *k) {
00795 
00796 /*                 C1 := C1 - W * V1 */
00797 
00798                     i__1 = lastv - *k;
00799                     q__1.r = -1.f, q__1.i = -0.f;
00800                     cgemm_("No transpose", "No transpose", &lastc, &i__1, k, &
00801                             q__1, &work[work_offset], ldwork, &v[v_offset], 
00802                             ldv, &c_b1, &c__[c_offset], ldc);
00803                 }
00804 
00805 /*              W := W * V2 */
00806 
00807                 ctrmm_("Right", "Lower", "No transpose", "Unit", &lastc, k, &
00808                         c_b1, &v[(lastv - *k + 1) * v_dim1 + 1], ldv, &work[
00809                         work_offset], ldwork);
00810 
00811 /*              C1 := C1 - W */
00812 
00813                 i__1 = *k;
00814                 for (j = 1; j <= i__1; ++j) {
00815                     i__2 = lastc;
00816                     for (i__ = 1; i__ <= i__2; ++i__) {
00817                         i__3 = i__ + (lastv - *k + j) * c_dim1;
00818                         i__4 = i__ + (lastv - *k + j) * c_dim1;
00819                         i__5 = i__ + j * work_dim1;
00820                         q__1.r = c__[i__4].r - work[i__5].r, q__1.i = c__[
00821                                 i__4].i - work[i__5].i;
00822                         c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00823 /* L230: */
00824                     }
00825 /* L240: */
00826                 }
00827 
00828             }
00829 
00830         }
00831     }
00832 
00833     return 0;
00834 
00835 /*     End of CLARFB */
00836 
00837 } /* clarfb_ */