📄 zlarfx.c
字号:
#include "f2c.h"
#include "netlib.h"
/* Modified by Peter Vanroose, June 2001: manual optimisation and clean-up */
/* Table of constant values */
static doublecomplex c_b2 = {0.,0.};
static doublecomplex c_b15 = {1.,0.};
static integer c__1 = 1;
/* Subroutine */ void zlarfx_(side, m, n, v, tau, c, ldc, work)
const char *side;
const integer *m, *n;
doublecomplex *v, *tau, *c;
const integer *ldc;
doublecomplex *work;
{
/* System generated locals */
integer i__1;
doublecomplex z__1;
/* Local variables */
static integer j;
static doublecomplex t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, sum;
/* -- LAPACK auxiliary routine (version 2.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* September 30, 1994 */
/* ===================================================================== */
/* */
/* Purpose */
/* ======= */
/* */
/* ZLARFX applies a complex elementary reflector H to a complex m by n */
/* matrix C, from either the left or the right. H is represented in the */
/* form */
/* */
/* H = I - tau * v * v' */
/* */
/* where tau is a complex scalar and v is a complex vector. */
/* */
/* If tau = 0, then H is taken to be the unit matrix */
/* */
/* This version uses inline code if H has order < 11. */
/* */
/* Arguments */
/* ========= */
/* */
/* SIDE (input) CHARACTER*1 */
/* = 'L': form H * C */
/* = 'R': form C * H */
/* */
/* M (input) INTEGER */
/* The number of rows of the matrix C. */
/* */
/* N (input) INTEGER */
/* The number of columns of the matrix C. */
/* */
/* V (input) COMPLEX*16 array, dimension (M) if SIDE = 'L' */
/* or (N) if SIDE = 'R' */
/* The vector v in the representation of H. */
/* */
/* TAU (input) COMPLEX*16 */
/* The value tau in the representation of H. */
/* */
/* C (input/output) COMPLEX*16 array, dimension (LDC,N) */
/* On entry, the m by n matrix C. */
/* On exit, C is overwritten by the matrix H * C if SIDE = 'L', */
/* or C * H if SIDE = 'R'. */
/* */
/* LDC (input) INTEGER */
/* The leading dimension of the array C. LDA >= max(1,M). */
/* */
/* WORK (workspace) COMPLEX*16 array, dimension (N) if SIDE = 'L' */
/* or (M) if SIDE = 'R' */
/* WORK is not referenced if H has order < 11. */
/* */
/* ===================================================================== */
/* Quick return if possible */
if (tau->r == 0. && tau->i == 0.) {
return;
}
if (lsame_(side, "L")) {
/* Form H * C, where H has order m. */
switch ((int)*m) {
case 1: goto L10;
case 2: goto L30;
case 3: goto L50;
case 4: goto L70;
case 5: goto L90;
case 6: goto L110;
case 7: goto L130;
case 8: goto L150;
case 9: goto L170;
case 10: goto L190;
}
/* Code for general M */
/* w := C'*v */
zgemv_("Conjugate transpose", m, n, &c_b15, c, ldc, v, &c__1, &c_b2, work, &c__1);
/* C := C - tau * v * w' */
z__1.r = -tau->r, z__1.i = -tau->i;
zgerc_(m, n, &z__1, v, &c__1, work, &c__1, c, ldc);
return; /* exit zlarfx */
L10:
/* Special code for 1 x 1 Householder */
z__1.r = tau->r * v[0].r - tau->i * v[0].i,
z__1.i = tau->r * v[0].i + tau->i * v[0].r;
t1.r = 1. - z__1.r * v[0].r - z__1.i * v[0].i,
t1.i = 0. + z__1.r * v[0].i - z__1.i * v[0].r;
for (j = 0; j < *n; ++j) {
i__1 = j * *ldc;
z__1.r = t1.r * c[i__1].r - t1.i * c[i__1].i,
z__1.i = t1.r * c[i__1].i + t1.i * c[i__1].r;
c[i__1].r = z__1.r, c[i__1].i = z__1.i;
}
return; /* exit zlarfx */
L30:
/* Special code for 2 x 2 Householder */
t1.r = tau->r * v[0].r - tau->i * v[0].i,
t1.i = tau->r * v[0].i + tau->i * v[0].r;
t2.r = tau->r * v[1].r - tau->i * v[1].i,
t2.i = tau->r * v[1].i + tau->i * v[1].r;
for (j = 0; j < *n; ++j) {
i__1 = j * *ldc;
sum.r = v[0].r * c[i__1].r + v[0].i * c[i__1].i,
sum.i = v[0].r * c[i__1].i - v[0].i * c[i__1].r;
++i__1;
sum.r += v[1].r * c[i__1].r + v[1].i * c[i__1].i,
sum.i += v[1].r * c[i__1].i - v[1].i * c[i__1].r;
i__1 = j * *ldc;
c[i__1].r -= sum.r * t1.r - sum.i * t1.i,
c[i__1].i -= sum.r * t1.i + sum.i * t1.r;
++i__1;
c[i__1].r -= sum.r * t2.r - sum.i * t2.i,
c[i__1].i -= sum.r * t2.i + sum.i * t2.r;
}
return; /* exit zlarfx */
L50:
/* Special code for 3 x 3 Householder */
t1.r = tau->r * v[0].r - tau->i * v[0].i,
t1.i = tau->r * v[0].i + tau->i * v[0].r;
t2.r = tau->r * v[1].r - tau->i * v[1].i,
t2.i = tau->r * v[1].i + tau->i * v[1].r;
t3.r = tau->r * v[2].r - tau->i * v[2].i,
t3.i = tau->r * v[2].i + tau->i * v[2].r;
for (j = 0; j < *n; ++j) {
i__1 = j * *ldc;
sum.r = v[0].r * c[i__1].r + v[0].i * c[i__1].i,
sum.i = v[0].r * c[i__1].i - v[0].i * c[i__1].r;
++i__1;
sum.r += v[1].r * c[i__1].r + v[1].i * c[i__1].i,
sum.i += v[1].r * c[i__1].i - v[1].i * c[i__1].r;
++i__1;
sum.r += v[2].r * c[i__1].r + v[2].i * c[i__1].i,
sum.i += v[2].r * c[i__1].i - v[2].i * c[i__1].r;
i__1 = j * *ldc;
c[i__1].r -= sum.r * t1.r - sum.i * t1.i,
c[i__1].i -= sum.r * t1.i + sum.i * t1.r;
++i__1;
c[i__1].r -= sum.r * t2.r - sum.i * t2.i,
c[i__1].i -= sum.r * t2.i + sum.i * t2.r;
++i__1;
c[i__1].r -= sum.r * t3.r - sum.i * t3.i,
c[i__1].i -= sum.r * t3.i + sum.i * t3.r;
}
return; /* exit zlarfx */
L70:
/* Special code for 4 x 4 Householder */
t1.r = tau->r * v[0].r - tau->i * v[0].i,
t1.i = tau->r * v[0].i + tau->i * v[0].r;
t2.r = tau->r * v[1].r - tau->i * v[1].i,
t2.i = tau->r * v[1].i + tau->i * v[1].r;
t3.r = tau->r * v[2].r - tau->i * v[2].i,
t3.i = tau->r * v[2].i + tau->i * v[2].r;
t4.r = tau->r * v[3].r - tau->i * v[3].i,
t4.i = tau->r * v[3].i + tau->i * v[3].r;
for (j = 0; j < *n; ++j) {
i__1 = j * *ldc;
sum.r = v[0].r * c[i__1].r + v[0].i * c[i__1].i,
sum.i = v[0].r * c[i__1].i - v[0].i * c[i__1].r;
++i__1;
sum.r += v[1].r * c[i__1].r + v[1].i * c[i__1].i,
sum.i += v[1].r * c[i__1].i - v[1].i * c[i__1].r;
++i__1;
sum.r += v[2].r * c[i__1].r + v[2].i * c[i__1].i,
sum.i += v[2].r * c[i__1].i - v[2].i * c[i__1].r;
++i__1;
sum.r += v[3].r * c[i__1].r + v[3].i * c[i__1].i,
sum.i += v[3].r * c[i__1].i - v[3].i * c[i__1].r;
i__1 = j * *ldc;
c[i__1].r -= sum.r * t1.r - sum.i * t1.i,
c[i__1].i -= sum.r * t1.i + sum.i * t1.r;
++i__1;
c[i__1].r -= sum.r * t2.r - sum.i * t2.i,
c[i__1].i -= sum.r * t2.i + sum.i * t2.r;
++i__1;
c[i__1].r -= sum.r * t3.r - sum.i * t3.i,
c[i__1].i -= sum.r * t3.i + sum.i * t3.r;
++i__1;
c[i__1].r -= sum.r * t4.r - sum.i * t4.i,
c[i__1].i -= sum.r * t4.i + sum.i * t4.r;
}
return; /* exit zlarfx */
L90:
/* Special code for 5 x 5 Householder */
t1.r = tau->r * v[0].r - tau->i * v[0].i,
t1.i = tau->r * v[0].i + tau->i * v[0].r;
t2.r = tau->r * v[1].r - tau->i * v[1].i,
t2.i = tau->r * v[1].i + tau->i * v[1].r;
t3.r = tau->r * v[2].r - tau->i * v[2].i,
t3.i = tau->r * v[2].i + tau->i * v[2].r;
t4.r = tau->r * v[3].r - tau->i * v[3].i,
t4.i = tau->r * v[3].i + tau->i * v[3].r;
t5.r = tau->r * v[4].r - tau->i * v[4].i,
t5.i = tau->r * v[4].i + tau->i * v[4].r;
for (j = 0; j < *n; ++j) {
i__1 = j * *ldc;
sum.r = v[0].r * c[i__1].r + v[0].i * c[i__1].i,
sum.i = v[0].r * c[i__1].i - v[0].i * c[i__1].r;
++i__1;
sum.r += v[1].r * c[i__1].r + v[1].i * c[i__1].i,
sum.i += v[1].r * c[i__1].i - v[1].i * c[i__1].r;
++i__1;
sum.r += v[2].r * c[i__1].r + v[2].i * c[i__1].i,
sum.i += v[2].r * c[i__1].i - v[2].i * c[i__1].r;
++i__1;
sum.r += v[3].r * c[i__1].r + v[3].i * c[i__1].i,
sum.i += v[3].r * c[i__1].i - v[3].i * c[i__1].r;
++i__1;
sum.r += v[4].r * c[i__1].r + v[4].i * c[i__1].i,
sum.i += v[4].r * c[i__1].i - v[4].i * c[i__1].r;
i__1 = j * *ldc;
c[i__1].r -= sum.r * t1.r - sum.i * t1.i,
c[i__1].i -= sum.r * t1.i + sum.i * t1.r;
++i__1;
c[i__1].r -= sum.r * t2.r - sum.i * t2.i,
c[i__1].i -= sum.r * t2.i + sum.i * t2.r;
++i__1;
c[i__1].r -= sum.r * t3.r - sum.i * t3.i,
c[i__1].i -= sum.r * t3.i + sum.i * t3.r;
++i__1;
c[i__1].r -= sum.r * t4.r - sum.i * t4.i,
c[i__1].i -= sum.r * t4.i + sum.i * t4.r;
++i__1;
c[i__1].r -= sum.r * t5.r - sum.i * t5.i,
c[i__1].i -= sum.r * t5.i + sum.i * t5.r;
}
return; /* exit zlarfx */
L110:
/* Special code for 6 x 6 Householder */
t1.r = tau->r * v[0].r - tau->i * v[0].i,
t1.i = tau->r * v[0].i + tau->i * v[0].r;
t2.r = tau->r * v[1].r - tau->i * v[1].i,
t2.i = tau->r * v[1].i + tau->i * v[1].r;
t3.r = tau->r * v[2].r - tau->i * v[2].i,
t3.i = tau->r * v[2].i + tau->i * v[2].r;
t4.r = tau->r * v[3].r - tau->i * v[3].i,
t4.i = tau->r * v[3].i + tau->i * v[3].r;
t5.r = tau->r * v[4].r - tau->i * v[4].i,
t5.i = tau->r * v[4].i + tau->i * v[4].r;
t6.r = tau->r * v[5].r - tau->i * v[5].i,
t6.i = tau->r * v[5].i + tau->i * v[5].r;
for (j = 0; j < *n; ++j) {
i__1 = j * *ldc;
sum.r = v[0].r * c[i__1].r + v[0].i * c[i__1].i,
sum.i = v[0].r * c[i__1].i - v[0].i * c[i__1].r;
++i__1;
sum.r += v[1].r * c[i__1].r + v[1].i * c[i__1].i,
sum.i += v[1].r * c[i__1].i - v[1].i * c[i__1].r;
++i__1;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -