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

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