swiftnav: dlarfb.c Source File

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