swiftnav: zlarfb.c Source File

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