swiftnav: slarfb.c Source File

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