git.maemo.org Git - opencv/blob - 3rdparty/lapack/slarfb.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static real c_b14 = 1.f;
   7 static real c_b25 = -1.f;
   8
   9 /* Subroutine */ int slarfb_(char *side, char *trans, char *direct, char *
  10         storev, integer *m, integer *n, integer *k, real *v, integer *ldv,
  11         real *t, integer *ldt, real *c__, integer *ldc, real *work, integer *
  12         ldwork)
  13 {
  14     /* System generated locals */
  15     integer c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1,
  16             work_offset, i__1, i__2;
  17
  18     /* Local variables */
  19     integer i__, j;
  20     extern logical lsame_(char *, char *);
  21     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
  22             integer *, real *, real *, integer *, real *, integer *, real *,
  23             real *, integer *), scopy_(integer *, real *,
  24             integer *, real *, integer *), strmm_(char *, char *, char *,
  25             char *, integer *, integer *, real *, real *, integer *, real *,
  26             integer *);
  27     char transt[1];
  28
  29
  30 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  31 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  32 /*     November 2006 */
  33
  34 /*     .. Scalar Arguments .. */
  35 /*     .. */
  36 /*     .. Array Arguments .. */
  37 /*     .. */
  38
  39 /*  Purpose */
  40 /*  ======= */
  41
  42 /*  SLARFB applies a real block reflector H or its transpose H' to a */
  43 /*  real m by n matrix C, from either the left or the right. */
  44
  45 /*  Arguments */
  46 /*  ========= */
  47
  48 /*  SIDE    (input) CHARACTER*1 */
  49 /*          = 'L': apply H or H' from the Left */
  50 /*          = 'R': apply H or H' from the Right */
  51
  52 /*  TRANS   (input) CHARACTER*1 */
  53 /*          = 'N': apply H (No transpose) */
  54 /*          = 'T': apply H' (Transpose) */
  55
  56 /*  DIRECT  (input) CHARACTER*1 */
  57 /*          Indicates how H is formed from a product of elementary */
  58 /*          reflectors */
  59 /*          = 'F': H = H(1) H(2) . . . H(k) (Forward) */
  60 /*          = 'B': H = H(k) . . . H(2) H(1) (Backward) */
  61
  62 /*  STOREV  (input) CHARACTER*1 */
  63 /*          Indicates how the vectors which define the elementary */
  64 /*          reflectors are stored: */
  65 /*          = 'C': Columnwise */
  66 /*          = 'R': Rowwise */
  67
  68 /*  M       (input) INTEGER */
  69 /*          The number of rows of the matrix C. */
  70
  71 /*  N       (input) INTEGER */
  72 /*          The number of columns of the matrix C. */
  73
  74 /*  K       (input) INTEGER */
  75 /*          The order of the matrix T (= the number of elementary */
  76 /*          reflectors whose product defines the block reflector). */
  77
  78 /*  V       (input) REAL array, dimension */
  79 /*                                (LDV,K) if STOREV = 'C' */
  80 /*                                (LDV,M) if STOREV = 'R' and SIDE = 'L' */
  81 /*                                (LDV,N) if STOREV = 'R' and SIDE = 'R' */
  82 /*          The matrix V. See further details. */
  83
  84 /*  LDV     (input) INTEGER */
  85 /*          The leading dimension of the array V. */
  86 /*          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); */
  87 /*          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); */
  88 /*          if STOREV = 'R', LDV >= K. */
  89
  90 /*  T       (input) REAL array, dimension (LDT,K) */
  91 /*          The triangular k by k matrix T in the representation of the */
  92 /*          block reflector. */
  93
  94 /*  LDT     (input) INTEGER */
  95 /*          The leading dimension of the array T. LDT >= K. */
  96
  97 /*  C       (input/output) REAL array, dimension (LDC,N) */
  98 /*          On entry, the m by n matrix C. */
  99 /*          On exit, C is overwritten by H*C or H'*C or C*H or C*H'. */
 100
 101 /*  LDC     (input) INTEGER */
 102 /*          The leading dimension of the array C. LDA >= max(1,M). */
 103
 104 /*  WORK    (workspace) REAL array, dimension (LDWORK,K) */
 105
 106 /*  LDWORK  (input) INTEGER */
 107 /*          The leading dimension of the array WORK. */
 108 /*          If SIDE = 'L', LDWORK >= max(1,N); */
 109 /*          if SIDE = 'R', LDWORK >= max(1,M). */
 110
 111 /*  ===================================================================== */
 112
 113 /*     .. Parameters .. */
 114 /*     .. */
 115 /*     .. Local Scalars .. */
 116 /*     .. */
 117 /*     .. External Functions .. */
 118 /*     .. */
 119 /*     .. External Subroutines .. */
 120 /*     .. */
 121 /*     .. Executable Statements .. */
 122
 123 /*     Quick return if possible */
 124
 125     /* Parameter adjustments */
 126     v_dim1 = *ldv;
 127     v_offset = 1 + v_dim1;
 128     v -= v_offset;
 129     t_dim1 = *ldt;
 130     t_offset = 1 + t_dim1;
 131     t -= t_offset;
 132     c_dim1 = *ldc;
 133     c_offset = 1 + c_dim1;
 134     c__ -= c_offset;
 135     work_dim1 = *ldwork;
 136     work_offset = 1 + work_dim1;
 137     work -= work_offset;
 138
 139     /* Function Body */
 140     if (*m <= 0 || *n <= 0) {
 141         return 0;
 142     }
 143
 144     if (lsame_(trans, "N")) {
 145         *(unsigned char *)transt = 'T';
 146     } else {
 147         *(unsigned char *)transt = 'N';
 148     }
 149
 150     if (lsame_(storev, "C")) {
 151
 152         if (lsame_(direct, "F")) {
 153
 154 /*           Let  V =  ( V1 )    (first K rows) */
 155 /*                     ( V2 ) */
 156 /*           where  V1  is unit lower triangular. */
 157
 158             if (lsame_(side, "L")) {
 159
 160 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
 161 /*                                                  ( C2 ) */
 162
 163 /*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK) */
 164
 165 /*              W := C1' */
 166
 167                 i__1 = *k;
 168                 for (j = 1; j <= i__1; ++j) {
 169                     scopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
 170                              &c__1);
 171 /* L10: */
 172                 }
 173
 174 /*              W := W * V1 */
 175
 176                 strmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
 177                          &v[v_offset], ldv, &work[work_offset], ldwork);
 178                 if (*m > *k) {
 179
 180 /*                 W := W + C2'*V2 */
 181
 182                     i__1 = *m - *k;
 183                     sgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
 184                             c__[*k + 1 + c_dim1], ldc, &v[*k + 1 + v_dim1],
 185                             ldv, &c_b14, &work[work_offset], ldwork);
 186                 }
 187
 188 /*              W := W * T'  or  W * T */
 189
 190                 strmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
 191                         t_offset], ldt, &work[work_offset], ldwork);
 192
 193 /*              C := C - V * W' */
 194
 195                 if (*m > *k) {
 196
 197 /*                 C2 := C2 - V2 * W' */
 198
 199                     i__1 = *m - *k;
 200                     sgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
 201                             v[*k + 1 + v_dim1], ldv, &work[work_offset],
 202                             ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
 203                 }
 204
 205 /*              W := W * V1' */
 206
 207                 strmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
 208                         v[v_offset], ldv, &work[work_offset], ldwork);
 209
 210 /*              C1 := C1 - W' */
 211
 212                 i__1 = *k;
 213                 for (j = 1; j <= i__1; ++j) {
 214                     i__2 = *n;
 215                     for (i__ = 1; i__ <= i__2; ++i__) {
 216                         c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
 217 /* L20: */
 218                     }
 219 /* L30: */
 220                 }
 221
 222             } else if (lsame_(side, "R")) {
 223
 224 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
 225
 226 /*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK) */
 227
 228 /*              W := C1 */
 229
 230                 i__1 = *k;
 231                 for (j = 1; j <= i__1; ++j) {
 232                     scopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
 233                             work_dim1 + 1], &c__1);
 234 /* L40: */
 235                 }
 236
 237 /*              W := W * V1 */
 238
 239                 strmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
 240                          &v[v_offset], ldv, &work[work_offset], ldwork);
 241                 if (*n > *k) {
 242
 243 /*                 W := W + C2 * V2 */
 244
 245                     i__1 = *n - *k;
 246                     sgemm_("No transpose", "No transpose", m, k, &i__1, &
 247                             c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc, &v[*k +
 248                             1 + v_dim1], ldv, &c_b14, &work[work_offset],
 249                             ldwork);
 250                 }
 251
 252 /*              W := W * T  or  W * T' */
 253
 254                 strmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
 255                         t_offset], ldt, &work[work_offset], ldwork);
 256
 257 /*              C := C - W * V' */
 258
 259                 if (*n > *k) {
 260
 261 /*                 C2 := C2 - W * V2' */
 262
 263                     i__1 = *n - *k;
 264                     sgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
 265                             work[work_offset], ldwork, &v[*k + 1 + v_dim1],
 266                             ldv, &c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc);
 267                 }
 268
 269 /*              W := W * V1' */
 270
 271                 strmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
 272                         v[v_offset], ldv, &work[work_offset], ldwork);
 273
 274 /*              C1 := C1 - W */
 275
 276                 i__1 = *k;
 277                 for (j = 1; j <= i__1; ++j) {
 278                     i__2 = *m;
 279                     for (i__ = 1; i__ <= i__2; ++i__) {
 280                         c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
 281 /* L50: */
 282                     }
 283 /* L60: */
 284                 }
 285             }
 286
 287         } else {
 288
 289 /*           Let  V =  ( V1 ) */
 290 /*                     ( V2 )    (last K rows) */
 291 /*           where  V2  is unit upper triangular. */
 292
 293             if (lsame_(side, "L")) {
 294
 295 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
 296 /*                                                  ( C2 ) */
 297
 298 /*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK) */
 299
 300 /*              W := C2' */
 301
 302                 i__1 = *k;
 303                 for (j = 1; j <= i__1; ++j) {
 304                     scopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
 305                             work_dim1 + 1], &c__1);
 306 /* L70: */
 307                 }
 308
 309 /*              W := W * V2 */
 310
 311                 strmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
 312                          &v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
 313                         ldwork);
 314                 if (*m > *k) {
 315
 316 /*                 W := W + C1'*V1 */
 317
 318                     i__1 = *m - *k;
 319                     sgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
 320                             c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
 321                             work[work_offset], ldwork);
 322                 }
 323
 324 /*              W := W * T'  or  W * T */
 325
 326                 strmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
 327                         t_offset], ldt, &work[work_offset], ldwork);
 328
 329 /*              C := C - V * W' */
 330
 331                 if (*m > *k) {
 332
 333 /*                 C1 := C1 - V1 * W' */
 334
 335                     i__1 = *m - *k;
 336                     sgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
 337                             v[v_offset], ldv, &work[work_offset], ldwork, &
 338                             c_b14, &c__[c_offset], ldc)
 339                             ;
 340                 }
 341
 342 /*              W := W * V2' */
 343
 344                 strmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
 345                         v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
 346                         ldwork);
 347
 348 /*              C2 := C2 - W' */
 349
 350                 i__1 = *k;
 351                 for (j = 1; j <= i__1; ++j) {
 352                     i__2 = *n;
 353                     for (i__ = 1; i__ <= i__2; ++i__) {
 354                         c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
 355                                 work_dim1];
 356 /* L80: */
 357                     }
 358 /* L90: */
 359                 }
 360
 361             } else if (lsame_(side, "R")) {
 362
 363 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
 364
 365 /*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK) */
 366
 367 /*              W := C2 */
 368
 369                 i__1 = *k;
 370                 for (j = 1; j <= i__1; ++j) {
 371                     scopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
 372                             j * work_dim1 + 1], &c__1);
 373 /* L100: */
 374                 }
 375
 376 /*              W := W * V2 */
 377
 378                 strmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
 379                          &v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
 380                         ldwork);
 381                 if (*n > *k) {
 382
 383 /*                 W := W + C1 * V1 */
 384
 385                     i__1 = *n - *k;
 386                     sgemm_("No transpose", "No transpose", m, k, &i__1, &
 387                             c_b14, &c__[c_offset], ldc, &v[v_offset], ldv, &
 388                             c_b14, &work[work_offset], ldwork);
 389                 }
 390
 391 /*              W := W * T  or  W * T' */
 392
 393                 strmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
 394                         t_offset], ldt, &work[work_offset], ldwork);
 395
 396 /*              C := C - W * V' */
 397
 398                 if (*n > *k) {
 399
 400 /*                 C1 := C1 - W * V1' */
 401
 402                     i__1 = *n - *k;
 403                     sgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
 404                             work[work_offset], ldwork, &v[v_offset], ldv, &
 405                             c_b14, &c__[c_offset], ldc)
 406                             ;
 407                 }
 408
 409 /*              W := W * V2' */
 410
 411                 strmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
 412                         v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
 413                         ldwork);
 414
 415 /*              C2 := C2 - W */
 416
 417                 i__1 = *k;
 418                 for (j = 1; j <= i__1; ++j) {
 419                     i__2 = *m;
 420                     for (i__ = 1; i__ <= i__2; ++i__) {
 421                         c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
 422                                 work_dim1];
 423 /* L110: */
 424                     }
 425 /* L120: */
 426                 }
 427             }
 428         }
 429
 430     } else if (lsame_(storev, "R")) {
 431
 432         if (lsame_(direct, "F")) {
 433
 434 /*           Let  V =  ( V1  V2 )    (V1: first K columns) */
 435 /*           where  V1  is unit upper triangular. */
 436
 437             if (lsame_(side, "L")) {
 438
 439 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
 440 /*                                                  ( C2 ) */
 441
 442 /*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK) */
 443
 444 /*              W := C1' */
 445
 446                 i__1 = *k;
 447                 for (j = 1; j <= i__1; ++j) {
 448                     scopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
 449                              &c__1);
 450 /* L130: */
 451                 }
 452
 453 /*              W := W * V1' */
 454
 455                 strmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
 456                         v[v_offset], ldv, &work[work_offset], ldwork);
 457                 if (*m > *k) {
 458
 459 /*                 W := W + C2'*V2' */
 460
 461                     i__1 = *m - *k;
 462                     sgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
 463                             c__[*k + 1 + c_dim1], ldc, &v[(*k + 1) * v_dim1 +
 464                             1], ldv, &c_b14, &work[work_offset], ldwork);
 465                 }
 466
 467 /*              W := W * T'  or  W * T */
 468
 469                 strmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
 470                         t_offset], ldt, &work[work_offset], ldwork);
 471
 472 /*              C := C - V' * W' */
 473
 474                 if (*m > *k) {
 475
 476 /*                 C2 := C2 - V2' * W' */
 477
 478                     i__1 = *m - *k;
 479                     sgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[(
 480                             *k + 1) * v_dim1 + 1], ldv, &work[work_offset],
 481                             ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
 482                 }
 483
 484 /*              W := W * V1 */
 485
 486                 strmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
 487                          &v[v_offset], ldv, &work[work_offset], ldwork);
 488
 489 /*              C1 := C1 - W' */
 490
 491                 i__1 = *k;
 492                 for (j = 1; j <= i__1; ++j) {
 493                     i__2 = *n;
 494                     for (i__ = 1; i__ <= i__2; ++i__) {
 495                         c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
 496 /* L140: */
 497                     }
 498 /* L150: */
 499                 }
 500
 501             } else if (lsame_(side, "R")) {
 502
 503 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
 504
 505 /*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK) */
 506
 507 /*              W := C1 */
 508
 509                 i__1 = *k;
 510                 for (j = 1; j <= i__1; ++j) {
 511                     scopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
 512                             work_dim1 + 1], &c__1);
 513 /* L160: */
 514                 }
 515
 516 /*              W := W * V1' */
 517
 518                 strmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
 519                         v[v_offset], ldv, &work[work_offset], ldwork);
 520                 if (*n > *k) {
 521
 522 /*                 W := W + C2 * V2' */
 523
 524                     i__1 = *n - *k;
 525                     sgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
 526                             c__[(*k + 1) * c_dim1 + 1], ldc, &v[(*k + 1) *
 527                             v_dim1 + 1], ldv, &c_b14, &work[work_offset],
 528                             ldwork);
 529                 }
 530
 531 /*              W := W * T  or  W * T' */
 532
 533                 strmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
 534                         t_offset], ldt, &work[work_offset], ldwork);
 535
 536 /*              C := C - W * V */
 537
 538                 if (*n > *k) {
 539
 540 /*                 C2 := C2 - W * V2 */
 541
 542                     i__1 = *n - *k;
 543                     sgemm_("No transpose", "No transpose", m, &i__1, k, &
 544                             c_b25, &work[work_offset], ldwork, &v[(*k + 1) *
 545                             v_dim1 + 1], ldv, &c_b14, &c__[(*k + 1) * c_dim1
 546                             + 1], ldc);
 547                 }
 548
 549 /*              W := W * V1 */
 550
 551                 strmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
 552                          &v[v_offset], ldv, &work[work_offset], ldwork);
 553
 554 /*              C1 := C1 - W */
 555
 556                 i__1 = *k;
 557                 for (j = 1; j <= i__1; ++j) {
 558                     i__2 = *m;
 559                     for (i__ = 1; i__ <= i__2; ++i__) {
 560                         c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
 561 /* L170: */
 562                     }
 563 /* L180: */
 564                 }
 565
 566             }
 567
 568         } else {
 569
 570 /*           Let  V =  ( V1  V2 )    (V2: last K columns) */
 571 /*           where  V2  is unit lower triangular. */
 572
 573             if (lsame_(side, "L")) {
 574
 575 /*              Form  H * C  or  H' * C  where  C = ( C1 ) */
 576 /*                                                  ( C2 ) */
 577
 578 /*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK) */
 579
 580 /*              W := C2' */
 581
 582                 i__1 = *k;
 583                 for (j = 1; j <= i__1; ++j) {
 584                     scopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
 585                             work_dim1 + 1], &c__1);
 586 /* L190: */
 587                 }
 588
 589 /*              W := W * V2' */
 590
 591                 strmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
 592                         v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
 593 , ldwork);
 594                 if (*m > *k) {
 595
 596 /*                 W := W + C1'*V1' */
 597
 598                     i__1 = *m - *k;
 599                     sgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
 600                             c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
 601                             work[work_offset], ldwork);
 602                 }
 603
 604 /*              W := W * T'  or  W * T */
 605
 606                 strmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
 607                         t_offset], ldt, &work[work_offset], ldwork);
 608
 609 /*              C := C - V' * W' */
 610
 611                 if (*m > *k) {
 612
 613 /*                 C1 := C1 - V1' * W' */
 614
 615                     i__1 = *m - *k;
 616                     sgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[
 617                             v_offset], ldv, &work[work_offset], ldwork, &
 618                             c_b14, &c__[c_offset], ldc);
 619                 }
 620
 621 /*              W := W * V2 */
 622
 623                 strmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
 624                          &v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[
 625                         work_offset], ldwork);
 626
 627 /*              C2 := C2 - W' */
 628
 629                 i__1 = *k;
 630                 for (j = 1; j <= i__1; ++j) {
 631                     i__2 = *n;
 632                     for (i__ = 1; i__ <= i__2; ++i__) {
 633                         c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
 634                                 work_dim1];
 635 /* L200: */
 636                     }
 637 /* L210: */
 638                 }
 639
 640             } else if (lsame_(side, "R")) {
 641
 642 /*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */
 643
 644 /*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK) */
 645
 646 /*              W := C2 */
 647
 648                 i__1 = *k;
 649                 for (j = 1; j <= i__1; ++j) {
 650                     scopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
 651                             j * work_dim1 + 1], &c__1);
 652 /* L220: */
 653                 }
 654
 655 /*              W := W * V2' */
 656
 657                 strmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
 658                         v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
 659 , ldwork);
 660                 if (*n > *k) {
 661
 662 /*                 W := W + C1 * V1' */
 663
 664                     i__1 = *n - *k;
 665                     sgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
 666                             c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
 667                             work[work_offset], ldwork);
 668                 }
 669
 670 /*              W := W * T  or  W * T' */
 671
 672                 strmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
 673                         t_offset], ldt, &work[work_offset], ldwork);
 674
 675 /*              C := C - W * V */
 676
 677                 if (*n > *k) {
 678
 679 /*                 C1 := C1 - W * V1 */
 680
 681                     i__1 = *n - *k;
 682                     sgemm_("No transpose", "No transpose", m, &i__1, k, &
 683                             c_b25, &work[work_offset], ldwork, &v[v_offset],
 684                             ldv, &c_b14, &c__[c_offset], ldc);
 685                 }
 686
 687 /*              W := W * V2 */
 688
 689                 strmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
 690                          &v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[
 691                         work_offset], ldwork);
 692
 693 /*              C1 := C1 - W */
 694
 695                 i__1 = *k;
 696                 for (j = 1; j <= i__1; ++j) {
 697                     i__2 = *m;
 698                     for (i__ = 1; i__ <= i__2; ++i__) {
 699                         c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
 700                                 work_dim1];
 701 /* L230: */
 702                     }
 703 /* L240: */
 704                 }
 705
 706             }
 707
 708         }
 709     }
 710
 711     return 0;
 712
 713 /*     End of SLARFB */
 714
 715 } /* slarfb_ */