zlarfx.c

来自「算断裂的」· C语言 代码 · 共 1,690 行 · 第 1/5 页

#include "f2c.h"

/* Subroutine */ int zlarfx_(char *side, integer *m, integer *n, 
	doublecomplex *v, doublecomplex *tau, doublecomplex *c, integer *ldc, 
	doublecomplex *work)
{
/*  -- 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.   

    ===================================================================== 
  


    
   Parameter adjustments   
       Function Body */
    /* Table of constant values */
    static doublecomplex c_b1 = {0.,0.};
    static doublecomplex c_b2 = {1.,0.};
    static integer c__1 = 1;
    
    /* System generated locals */
    integer c_dim1, c_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, 
	    i__9, i__10, i__11;
    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6, z__7, z__8, z__9, z__10,
	     z__11, z__12, z__13, z__14, z__15, z__16, z__17, z__18, z__19;
    /* Builtin functions */
    void d_cnjg(doublecomplex *, doublecomplex *);
    /* Local variables */
    static integer j;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *), zgemv_(char *, integer *, integer *, 
	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
	    integer *, doublecomplex *, doublecomplex *, integer *);
    static doublecomplex t1, t2, t3, t4, t5, t6, t7, t8, t9, v1, v2, v3, v4, 
	    v5, v6, v7, v8, v9, t10, v10, sum;



#define V(I) v[(I)-1]
#define WORK(I) work[(I)-1]

#define C(I,J) c[(I)-1 + ((J)-1)* ( *ldc)]

    if (tau->r == 0. && tau->i == 0.) {
	return 0;
    }
    if (lsame_(side, "L")) {

/*        Form  H * C, where H has order m. */

	switch (*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_b2, &C(1,1), ldc, &V(1), &
		c__1, &c_b1, &WORK(1), &c__1);

/*        C := C - tau * v * w' */

	z__1.r = -tau->r, z__1.i = -tau->i;
	zgerc_(m, n, &z__1, &V(1), &c__1, &WORK(1), &c__1, &C(1,1), ldc);
	goto L410;
L10:

/*        Special code for 1 x 1 Householder */

	z__3.r = tau->r * V(1).r - tau->i * V(1).i, z__3.i = tau->r * V(1).i 
		+ tau->i * V(1).r;
	d_cnjg(&z__4, &V(1));
	z__2.r = z__3.r * z__4.r - z__3.i * z__4.i, z__2.i = z__3.r * z__4.i 
		+ z__3.i * z__4.r;
	z__1.r = 1. - z__2.r, z__1.i = 0. - z__2.i;
	t1.r = z__1.r, t1.i = z__1.i;
	i__1 = *n;
	for (j = 1; j <= *n; ++j) {
	    i__2 = j * c_dim1 + 1;
	    i__3 = j * c_dim1 + 1;
	    z__1.r = t1.r * C(1,j).r - t1.i * C(1,j).i, z__1.i = t1.r * C(1,j).i + t1.i * C(1,j).r;
	    C(1,j).r = z__1.r, C(1,j).i = z__1.i;
/* L20: */
	}
	goto L410;
L30:

/*        Special code for 2 x 2 Householder */

	d_cnjg(&z__1, &V(1));
	v1.r = z__1.r, v1.i = z__1.i;
	d_cnjg(&z__2, &v1);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t1.r = z__1.r, t1.i = z__1.i;
	d_cnjg(&z__1, &V(2));
	v2.r = z__1.r, v2.i = z__1.i;
	d_cnjg(&z__2, &v2);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t2.r = z__1.r, t2.i = z__1.i;
	i__1 = *n;
	for (j = 1; j <= *n; ++j) {
	    i__2 = j * c_dim1 + 1;
	    z__2.r = v1.r * C(1,j).r - v1.i * C(1,j).i, z__2.i = v1.r * C(1,j).i + v1.i * C(1,j).r;
	    i__3 = j * c_dim1 + 2;
	    z__3.r = v2.r * C(2,j).r - v2.i * C(2,j).i, z__3.i = v2.r * C(2,j).i + v2.i * C(2,j).r;
	    z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
	    sum.r = z__1.r, sum.i = z__1.i;
	    i__2 = j * c_dim1 + 1;
	    i__3 = j * c_dim1 + 1;
	    z__2.r = sum.r * t1.r - sum.i * t1.i, z__2.i = sum.r * t1.i + 
		    sum.i * t1.r;
	    z__1.r = C(1,j).r - z__2.r, z__1.i = C(1,j).i - z__2.i;
	    C(1,j).r = z__1.r, C(1,j).i = z__1.i;
	    i__2 = j * c_dim1 + 2;
	    i__3 = j * c_dim1 + 2;
	    z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + 
		    sum.i * t2.r;
	    z__1.r = C(2,j).r - z__2.r, z__1.i = C(2,j).i - z__2.i;
	    C(2,j).r = z__1.r, C(2,j).i = z__1.i;
/* L40: */
	}
	goto L410;
L50:

/*        Special code for 3 x 3 Householder */

	d_cnjg(&z__1, &V(1));
	v1.r = z__1.r, v1.i = z__1.i;
	d_cnjg(&z__2, &v1);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t1.r = z__1.r, t1.i = z__1.i;
	d_cnjg(&z__1, &V(2));
	v2.r = z__1.r, v2.i = z__1.i;
	d_cnjg(&z__2, &v2);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t2.r = z__1.r, t2.i = z__1.i;
	d_cnjg(&z__1, &V(3));
	v3.r = z__1.r, v3.i = z__1.i;
	d_cnjg(&z__2, &v3);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t3.r = z__1.r, t3.i = z__1.i;
	i__1 = *n;
	for (j = 1; j <= *n; ++j) {
	    i__2 = j * c_dim1 + 1;
	    z__3.r = v1.r * C(1,j).r - v1.i * C(1,j).i, z__3.i = v1.r * C(1,j).i + v1.i * C(1,j).r;
	    i__3 = j * c_dim1 + 2;
	    z__4.r = v2.r * C(2,j).r - v2.i * C(2,j).i, z__4.i = v2.r * C(2,j).i + v2.i * C(2,j).r;
	    z__2.r = z__3.r + z__4.r, z__2.i = z__3.i + z__4.i;
	    i__4 = j * c_dim1 + 3;
	    z__5.r = v3.r * C(3,j).r - v3.i * C(3,j).i, z__5.i = v3.r * C(3,j).i + v3.i * C(3,j).r;
	    z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
	    sum.r = z__1.r, sum.i = z__1.i;
	    i__2 = j * c_dim1 + 1;
	    i__3 = j * c_dim1 + 1;
	    z__2.r = sum.r * t1.r - sum.i * t1.i, z__2.i = sum.r * t1.i + 
		    sum.i * t1.r;
	    z__1.r = C(1,j).r - z__2.r, z__1.i = C(1,j).i - z__2.i;
	    C(1,j).r = z__1.r, C(1,j).i = z__1.i;
	    i__2 = j * c_dim1 + 2;
	    i__3 = j * c_dim1 + 2;
	    z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + 
		    sum.i * t2.r;
	    z__1.r = C(2,j).r - z__2.r, z__1.i = C(2,j).i - z__2.i;
	    C(2,j).r = z__1.r, C(2,j).i = z__1.i;
	    i__2 = j * c_dim1 + 3;
	    i__3 = j * c_dim1 + 3;
	    z__2.r = sum.r * t3.r - sum.i * t3.i, z__2.i = sum.r * t3.i + 
		    sum.i * t3.r;
	    z__1.r = C(3,j).r - z__2.r, z__1.i = C(3,j).i - z__2.i;
	    C(3,j).r = z__1.r, C(3,j).i = z__1.i;
/* L60: */
	}
	goto L410;
L70:

/*        Special code for 4 x 4 Householder */

	d_cnjg(&z__1, &V(1));
	v1.r = z__1.r, v1.i = z__1.i;
	d_cnjg(&z__2, &v1);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t1.r = z__1.r, t1.i = z__1.i;
	d_cnjg(&z__1, &V(2));
	v2.r = z__1.r, v2.i = z__1.i;
	d_cnjg(&z__2, &v2);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t2.r = z__1.r, t2.i = z__1.i;
	d_cnjg(&z__1, &V(3));
	v3.r = z__1.r, v3.i = z__1.i;
	d_cnjg(&z__2, &v3);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t3.r = z__1.r, t3.i = z__1.i;
	d_cnjg(&z__1, &V(4));
	v4.r = z__1.r, v4.i = z__1.i;
	d_cnjg(&z__2, &v4);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t4.r = z__1.r, t4.i = z__1.i;
	i__1 = *n;
	for (j = 1; j <= *n; ++j) {
	    i__2 = j * c_dim1 + 1;
	    z__4.r = v1.r * C(1,j).r - v1.i * C(1,j).i, z__4.i = v1.r * C(1,j).i + v1.i * C(1,j).r;
	    i__3 = j * c_dim1 + 2;
	    z__5.r = v2.r * C(2,j).r - v2.i * C(2,j).i, z__5.i = v2.r * C(2,j).i + v2.i * C(2,j).r;
	    z__3.r = z__4.r + z__5.r, z__3.i = z__4.i + z__5.i;
	    i__4 = j * c_dim1 + 3;
	    z__6.r = v3.r * C(3,j).r - v3.i * C(3,j).i, z__6.i = v3.r * C(3,j).i + v3.i * C(3,j).r;
	    z__2.r = z__3.r + z__6.r, z__2.i = z__3.i + z__6.i;
	    i__5 = j * c_dim1 + 4;
	    z__7.r = v4.r * C(4,j).r - v4.i * C(4,j).i, z__7.i = v4.r * C(4,j).i + v4.i * C(4,j).r;
	    z__1.r = z__2.r + z__7.r, z__1.i = z__2.i + z__7.i;
	    sum.r = z__1.r, sum.i = z__1.i;
	    i__2 = j * c_dim1 + 1;
	    i__3 = j * c_dim1 + 1;
	    z__2.r = sum.r * t1.r - sum.i * t1.i, z__2.i = sum.r * t1.i + 
		    sum.i * t1.r;
	    z__1.r = C(1,j).r - z__2.r, z__1.i = C(1,j).i - z__2.i;
	    C(1,j).r = z__1.r, C(1,j).i = z__1.i;
	    i__2 = j * c_dim1 + 2;
	    i__3 = j * c_dim1 + 2;
	    z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + 
		    sum.i * t2.r;
	    z__1.r = C(2,j).r - z__2.r, z__1.i = C(2,j).i - z__2.i;
	    C(2,j).r = z__1.r, C(2,j).i = z__1.i;
	    i__2 = j * c_dim1 + 3;
	    i__3 = j * c_dim1 + 3;
	    z__2.r = sum.r * t3.r - sum.i * t3.i, z__2.i = sum.r * t3.i + 
		    sum.i * t3.r;
	    z__1.r = C(3,j).r - z__2.r, z__1.i = C(3,j).i - z__2.i;
	    C(3,j).r = z__1.r, C(3,j).i = z__1.i;
	    i__2 = j * c_dim1 + 4;
	    i__3 = j * c_dim1 + 4;
	    z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + 
		    sum.i * t4.r;
	    z__1.r = C(4,j).r - z__2.r, z__1.i = C(4,j).i - z__2.i;
	    C(4,j).r = z__1.r, C(4,j).i = z__1.i;
/* L80: */
	}
	goto L410;
L90:

/*        Special code for 5 x 5 Householder */

	d_cnjg(&z__1, &V(1));
	v1.r = z__1.r, v1.i = z__1.i;
	d_cnjg(&z__2, &v1);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t1.r = z__1.r, t1.i = z__1.i;
	d_cnjg(&z__1, &V(2));
	v2.r = z__1.r, v2.i = z__1.i;
	d_cnjg(&z__2, &v2);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t2.r = z__1.r, t2.i = z__1.i;
	d_cnjg(&z__1, &V(3));
	v3.r = z__1.r, v3.i = z__1.i;
	d_cnjg(&z__2, &v3);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i 
		+ tau->i * z__2.r;
	t3.r = z__1.r, t3.i = z__1.i;
	d_cnjg(&z__1, &V(4));
	v4.r = z__1.r, v4.i = z__1.i;
	d_cnjg(&z__2, &v4);
	z__1.r = tau->r * z__2.r - tau->i * z__2.i, z__1.i = tau->r * z__2.i