zlarfx.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/* lapack/complex16/zlarfx.f -- translated by f2c (version 20050501).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

/* Disable "global optimizations" to avoid optimizer bugs in this file.
   For GCC the file should be compiled with the -fno-gcse option.  Here
   is a note from the GCC man page:

    Note: When compiling a program using computed gotos, a GCC exten-
    sion, you may get better runtime performance if you disable the
    global common subexpression elimination pass by adding -fno-gcse to
    the command line.
*/
#ifdef _MSC_VER
# pragma optimize ("g",off)
#endif

#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"

/* Table of constant values */

static doublecomplex c_b1 = {0.,0.};
static doublecomplex c_b2 = {1.,0.};
static integer c__1 = 1;

/*<       SUBROUTINE ZLARFX( SIDE, M, N, V, TAU, C, LDC, WORK ) >*/
/* Subroutine */ int zlarfx_(char *side, integer *m, integer *n, 
        doublecomplex *v, doublecomplex *tau, doublecomplex *c__, integer *
        ldc, doublecomplex *work, ftnlen side_len)
{
    /* 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 */
    integer j;
    doublecomplex t1, t2, t3, t4, t5, t6, t7, t8, t9, v1, v2, v3, v4, v5, v6, 
            v7, v8, v9, t10, v10, sum;
    extern logical lsame_(char *, char *, ftnlen, ftnlen);
    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 *, ftnlen);
    (void)side_len;

/*  -- LAPACK auxiliary routine (version 3.0) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/*     Courant Institute, Argonne National Lab, and Rice University */
/*     September 30, 1994 */

/*     .. Scalar Arguments .. */
/*<       CHARACTER          SIDE >*/
/*<       INTEGER            LDC, M, N >*/
/*<       COMPLEX*16         TAU >*/
/*     .. */
/*     .. Array Arguments .. */
/*<       COMPLEX*16         C( LDC, * ), V( * ), WORK( * ) >*/
/*     .. */

/*  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. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*<       COMPLEX*16         ZERO, ONE >*/
/*<    >*/
/*     .. */
/*     .. Local Scalars .. */
/*<       INTEGER            J >*/
/*<    >*/
/*     .. */
/*     .. External Functions .. */
/*<       LOGICAL            LSAME >*/
/*<       EXTERNAL           LSAME >*/
/*     .. */
/*     .. External Subroutines .. */
/*<       EXTERNAL           ZGEMV, ZGERC >*/
/*     .. */
/*     .. Intrinsic Functions .. */
/*<       INTRINSIC          DCONJG >*/
/*     .. */
/*     .. Executable Statements .. */

/*<    >*/
    /* Parameter adjustments */
    --v;
    c_dim1 = *ldc;
    c_offset = 1 + c_dim1;
    c__ -= c_offset;
    --work;

    /* Function Body */
    if (tau->r == 0. && tau->i == 0.) {
        return 0;
    }
/*<       IF( LSAME( SIDE, 'L' ) ) THEN >*/
    if (lsame_(side, "L", (ftnlen)1, (ftnlen)1)) {

/*        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__[c_offset], ldc, &v[1],
                 &c__1, &c_b1, &work[1], &c__1, (ftnlen)19);

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

/*<          CALL ZGERC( M, N, -TAU, V, 1, WORK, 1, C, LDC ) >*/
        z__1.r = -tau->r, z__1.i = -tau->i;
        zgerc_(m, n, &z__1, &v[1], &c__1, &work[1], &c__1, &c__[c_offset], 
                ldc);
/*<          GO TO 410 >*/
        goto L410;
/*<    10    CONTINUE >*/
L10:

/*        Special code for 1 x 1 Householder */

/*<          T1 = ONE - TAU*V( 1 )*DCONJG( V( 1 ) ) >*/
        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;
/*<          DO 20 J = 1, N >*/
        i__1 = *n;
        for (j = 1; j <= i__1; ++j) {
/*<             C( 1, J ) = T1*C( 1, J ) >*/
            i__2 = j * c_dim1 + 1;
            i__3 = j * c_dim1 + 1;
            z__1.r = t1.r * c__[i__3].r - t1.i * c__[i__3].i, z__1.i = t1.r * 
                    c__[i__3].i + t1.i * c__[i__3].r;
            c__[i__2].r = z__1.r, c__[i__2].i = z__1.i;
/*<    20    CONTINUE >*/
/* L20: */
        }
/*<          GO TO 410 >*/
        goto L410;
/*<    30    CONTINUE >*/
L30:

/*        Special code for 2 x 2 Householder */

/*<          V1 = DCONJG( V( 1 ) ) >*/
        d_cnjg(&z__1, &v[1]);
        v1.r = z__1.r, v1.i = z__1.i;
/*<          T1 = TAU*DCONJG( V1 ) >*/
        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;
/*<          V2 = DCONJG( V( 2 ) ) >*/
        d_cnjg(&z__1, &v[2]);
        v2.r = z__1.r, v2.i = z__1.i;
/*<          T2 = TAU*DCONJG( V2 ) >*/
        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;
/*<          DO 40 J = 1, N >*/
        i__1 = *n;
        for (j = 1; j <= i__1; ++j) {
/*<             SUM = V1*C( 1, J ) + V2*C( 2, J ) >*/
            i__2 = j * c_dim1 + 1;
            z__2.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__2.i = v1.r * 
                    c__[i__2].i + v1.i * c__[i__2].r;
            i__3 = j * c_dim1 + 2;
            z__3.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__3.i = v2.r * 
                    c__[i__3].i + v2.i * c__[i__3].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;
/*<             C( 1, J ) = C( 1, J ) - SUM*T1 >*/
            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__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i;
            c__[i__2].r = z__1.r, c__[i__2].i = z__1.i;
/*<             C( 2, J ) = C( 2, J ) - SUM*T2 >*/
            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__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i;
            c__[i__2].r = z__1.r, c__[i__2].i = z__1.i;
/*<    40    CONTINUE >*/
/* L40: */
        }
/*<          GO TO 410 >*/
        goto L410;
/*<    50    CONTINUE >*/
L50:

/*        Special code for 3 x 3 Householder */

/*<          V1 = DCONJG( V( 1 ) ) >*/
        d_cnjg(&z__1, &v[1]);
        v1.r = z__1.r, v1.i = z__1.i;
/*<          T1 = TAU*DCONJG( V1 ) >*/
        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;
/*<          V2 = DCONJG( V( 2 ) ) >*/
        d_cnjg(&z__1, &v[2]);
        v2.r = z__1.r, v2.i = z__1.i;
/*<          T2 = TAU*DCONJG( V2 ) >*/
        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;
/*<          V3 = DCONJG( V( 3 ) ) >*/
        d_cnjg(&z__1, &v[3]);
        v3.r = z__1.r, v3.i = z__1.i;
/*<          T3 = TAU*DCONJG( V3 ) >*/
        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;
/*<          DO 60 J = 1, N >*/