自学教程：C++ zgemv_函数代码示例

/*! _zrovector*zgematrix operator */inline _zrovector operator*(const _zrovector& vec, const zgematrix& mat){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] operator*(const _zrovector&, const zgematrix&)"            << std::endl;#endif//CPPL_VERBOSE  #ifdef  CPPL_DEBUG  if(vec.L!=mat.M){    std::cerr << "[ERROR] operator*(const _zrovector&, const zgematrix&)"              << std::endl              << "These vector and matrix can not make a product."              << std::endl              << "Your input was (" << vec.L << ") * ("              << mat.M << "x" << mat.N << ")." << std::endl;    exit(1);  }#endif//CPPL_DEBUG    zrovector newvec(mat.N);  zgemv_( 'T', mat.M, mat.N, std::complex<double>(1.0,0.0), mat.Array, mat.M,          vec.Array, 1, std::complex<double>(0.0,0.0), newvec.array, 1 );    vec.destroy();  return _(newvec);}

PyObject* gemv(PyObject *self, PyObject *args){  Py_complex alpha;  PyArrayObject* a;  PyArrayObject* x;  Py_complex beta;  PyArrayObject* y;  char trans = 't';  if (!PyArg_ParseTuple(args, "DOODO|c", &alpha, &a, &x, &beta, &y, &trans))    return NULL;  int m, n, lda, itemsize, incx, incy;  if (trans == 'n')    {      m = PyArray_DIMS(a)[1];      for (int i = 2; i < PyArray_NDIM(a); i++)        m *= PyArray_DIMS(a)[i];      n = PyArray_DIMS(a)[0];      lda = MAX(1, m);    }  else    {      n = PyArray_DIMS(a)[0];      for (int i = 1; i < PyArray_NDIM(a)-1; i++)        n *= PyArray_DIMS(a)[i];      m = PyArray_DIMS(a)[PyArray_NDIM(a)-1];      lda = MAX(1, m);    }  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)    itemsize = sizeof(double);  else    itemsize = sizeof(double_complex);  incx = PyArray_STRIDES(x)[0]/itemsize;  incy = 1;  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)    dgemv_(&trans, &m, &n,           &(alpha.real),           DOUBLEP(a), &lda,           DOUBLEP(x), &incx,           &(beta.real),           DOUBLEP(y), &incy);  else    zgemv_(&trans, &m, &n,           &alpha,           (void*)COMPLEXP(a), &lda,           (void*)COMPLEXP(x), &incx,           &beta,           (void*)COMPLEXP(y), &incy);  Py_RETURN_NONE;}

intf2c_zgemv(char* trans, integer* M, integer* N,          doublecomplex* alpha,          doublecomplex* A, integer* lda,          doublecomplex* X, integer* incX,          doublecomplex* beta,          doublecomplex* Y, integer* incY){    zgemv_(trans, M, N,           alpha, A, lda, X, incX, beta, Y, incY);    return 0;}

/*! zrovector*zgematrix operator */inline _zrovector operator*(const zrovector& vec, const zgematrix& mat){VERBOSE_REPORT;#ifdef  CPPL_DEBUG  if(vec.l!=mat.m){    ERROR_REPORT;    std::cerr << "These vector and matrix can not make a product." << std::endl              << "Your input was (" << vec.l << ") * (" << mat.m << "x" << mat.n << ")." << std::endl;    exit(1);  }#endif//CPPL_DEBUG    zrovector newvec(mat.n);  zgemv_( 'T', mat.m, mat.n, comple(1.0,0.0), mat.array, mat.m,          vec.array, 1, comple(0.0,0.0), newvec.array, 1 );    return _(newvec);}

/*! _zgematrix*zcovector operator */inline _zcovector operator*(const _zgematrix& mat, const zcovector& vec){VERBOSE_REPORT;#ifdef  CPPL_DEBUG  if(mat.n!=vec.l){    ERROR_REPORT;    std::cerr << "These matrix and vector can not make a product." << std::endl              << "Your input was (" << mat.m << "x" << mat.n << ") * (" << vec.l << ")." << std::endl;    exit(1);  }#endif//CPPL_DEBUG    zcovector newvec(mat.m);  zgemv_( 'n', mat.m, mat.n, comple(1.0,0.0), mat.array, mat.m,          vec.array, 1, comple(0.0,0.0), newvec.array, 1 );    mat.destroy();  return _(newvec);}

/*! zgematrix*zcovector operator */inline _zcovector operator*(const zgematrix& mat, const zcovector& vec){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] operator*(const zgematrix&, const zcovector&)"            << std::endl;#endif//CPPL_VERBOSE  #ifdef  CPPL_DEBUG  if(mat.N!=vec.L){    std::cerr << "[ERROR] operator*(const zgematrix&, const zcovector&)"              << std::endl              << "These matrix and vector can not make a product." << std::endl              << "Your input was (" << mat.M << "x" << mat.N << ") * ("              << vec.L << ")." << std::endl;    exit(1);  }#endif//CPPL_DEBUG    zcovector newvec(mat.M);  zgemv_( 'N', mat.M, mat.N, std::complex<double>(1.0,0.0), mat.Array, mat.M,          vec.Array, 1, std::complex<double>(0.0,0.0), newvec.array, 1 );    return _(newvec);}

//.........这里部分代码省略........./*          ZTZRZF. V is not used if TAU = 0. *//*  INCV    (input) INTEGER *//*          The increment between elements of v. INCV <> 0. *//*  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. LDC >= max(1,M). *//*  WORK    (workspace) COMPLEX*16 array, dimension *//*                         (N) if SIDE = 'L' *//*                      or (M) if SIDE = 'R' *//*  Further Details *//*  =============== *//*  Based on contributions by *//*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA *//*  ===================================================================== */    /* Parameter adjustments */    --v;    c_dim1 = *ldc;    c_offset = 1 + c_dim1;    c__ -= c_offset;    --work;    /* Function Body */    if (lsame_(side, "L")) {/*        Form  H * C */	if (tau->r != 0. || tau->i != 0.) {/*           w( 1:n ) = conjg( C( 1, 1:n ) ) */	    zcopy_(n, &c__[c_offset], ldc, &work[1], &c__1);	    zlacgv_(n, &work[1], &c__1);/*           w( 1:n ) = conjg( w( 1:n ) + C( m-l+1:m, 1:n )' * v( 1:l ) ) */	    zgemv_("Conjugate transpose", l, n, &c_b1, &c__[*m - *l + 1 + 		    c_dim1], ldc, &v[1], incv, &c_b1, &work[1], &c__1);	    zlacgv_(n, &work[1], &c__1);/*           C( 1, 1:n ) = C( 1, 1:n ) - tau * w( 1:n ) */	    z__1.r = -tau->r, z__1.i = -tau->i;	    zaxpy_(n, &z__1, &work[1], &c__1, &c__[c_offset], ldc);/*                               tau * v( 1:l ) * conjg( w( 1:n )' ) */	    z__1.r = -tau->r, z__1.i = -tau->i;	    zgeru_(l, n, &z__1, &v[1], incv, &work[1], &c__1, &c__[*m - *l + 		    1 + c_dim1], ldc);	}    } else {/*        Form  C * H */	if (tau->r != 0. || tau->i != 0.) {/*           w( 1:m ) = C( 1:m, 1 ) */	    zcopy_(m, &c__[c_offset], &c__1, &work[1], &c__1);/*           w( 1:m ) = w( 1:m ) + C( 1:m, n-l+1:n, 1:n ) * v( 1:l ) */	    zgemv_("No transpose", m, l, &c_b1, &c__[(*n - *l + 1) * c_dim1 + 		    1], ldc, &v[1], incv, &c_b1, &work[1], &c__1);/*           C( 1:m, 1 ) = C( 1:m, 1 ) - tau * w( 1:m ) */	    z__1.r = -tau->r, z__1.i = -tau->i;	    zaxpy_(m, &z__1, &work[1], &c__1, &c__[c_offset], &c__1);/*                               tau * w( 1:m ) * v( 1:l )' */	    z__1.r = -tau->r, z__1.i = -tau->i;	    zgerc_(m, l, &z__1, &work[1], &c__1, &v[1], incv, &c__[(*n - *l + 		    1) * c_dim1 + 1], ldc);	}    }    return 0;/*     End of ZLARZ */} /* zlarz_ */

//.........这里部分代码省略.........            On exit, X and Y are the solutions of the GLM problem.       WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)               On exit, if INFO = 0, WORK(1) returns the optimal LWORK.       LWORK   (input) INTEGER               The dimension of the array WORK. LWORK >= max(1,N+M+P).               For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,               where NB is an upper bound for the optimal blocksizes for               ZGEQRF, CGERQF, ZUNMQR and CUNMRQ.       INFO    (output) INTEGER               = 0:  successful exit.               < 0:  if INFO = -i, the i-th argument had an illegal value.       ===================================================================          Test the input parameters          Parameter adjustments          Function Body */    /* Table of constant values */    static doublecomplex c_b2 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;    doublereal d__1;    doublecomplex z__1;    /* Local variables */    static integer lopt, i;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *), ztrsv_(char *, char *, char *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static integer np;    extern /* Subroutine */ int xerbla_(char *, integer *), zggqrf_(	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *, integer *), zunmqr_(char *, char *, 	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, integer *), zunmrq_(char *, char *, 	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, integer *);#define D(I) d[(I)-1]#define X(I) x[(I)-1]#define Y(I) y[(I)-1]#define WORK(I) work[(I)-1]#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]#define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)]    *info = 0;    np = min(*n,*p);    if (*n < 0) {	*info = -1;    } else if (*m < 0 || *m > *n) {

//.........这里部分代码省略.........    possible overflow.       Each eigenvector is normalized so that the element of largest       magnitude has magnitude 1; here the magnitude of a complex number       (x,y) is taken to be |x| + |y|.       =====================================================================          Decode and test the input parameters          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b2 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, 	    i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static logical allv;    static doublereal unfl, ovfl, smin;    static logical over;    static integer i__, j, k;    static doublereal scale;    extern logical lsame_(char *, char *);    static doublereal remax;    static logical leftv, bothv;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static logical somev;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), dlabad_(doublereal *, doublereal *);    static integer ii, ki;    extern doublereal dlamch_(char *);    static integer is;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);    static logical rightv;    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    static doublereal smlnum;    extern /* Subroutine */ int zlatrs_(char *, char *, char *, char *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublereal *, doublereal *, integer *);    static doublereal ulp;#define t_subscr(a_1,a_2) (a_2)*t_dim1 + a_1#define t_ref(a_1,a_2) t[t_subscr(a_1,a_2)]#define vl_subscr(a_1,a_2) (a_2)*vl_dim1 + a_1#define vl_ref(a_1,a_2) vl[vl_subscr(a_1,a_2)]#define vr_subscr(a_1,a_2) (a_2)*vr_dim1 + a_1#define vr_ref(a_1,a_2) vr[vr_subscr(a_1,a_2)]    --select;    t_dim1 = *ldt;    t_offset = 1 + t_dim1 * 1;    t -= t_offset;    vl_dim1 = *ldvl;    vl_offset = 1 + vl_dim1 * 1;

/* Subroutine */ int znaitr_(integer *ido, char *bmat, integer *n, integer *k,	 integer *np, integer *nb, doublecomplex *resid, doublereal *rnorm, 	doublecomplex *v, integer *ldv, doublecomplex *h__, integer *ldh, 	integer *ipntr, doublecomplex *workd, integer *info, ftnlen bmat_len){    /* Initialized data */    static logical first = TRUE_;    /* System generated locals */    integer h_dim1, h_offset, v_dim1, v_offset, i__1, i__2, i__3;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1;    /* Builtin functions */    double d_imag(doublecomplex *), sqrt(doublereal);    /* Local variables */    static integer i__, j;    static real t0, t1, t2, t3, t4, t5;    static integer jj, ipj, irj, ivj;    static doublereal ulp, tst1;    static integer ierr, iter;    static doublereal unfl, ovfl;    static integer itry;    static doublereal temp1;    static logical orth1, orth2, step3, step4;    static doublereal betaj;    static integer infol;    static doublecomplex cnorm;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static doublereal rtemp[2];    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, ftnlen);    static doublereal wnorm;    extern /* Subroutine */ int dvout_(integer *, integer *, doublereal *, 	    integer *, char *, ftnlen), zcopy_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), ivout_(integer *, integer 	    *, integer *, integer *, char *, ftnlen), zaxpy_(integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zmout_(integer *, integer *, integer *, doublecomplex 	    *, integer *, integer *, char *, ftnlen), zvout_(integer *, 	    integer *, doublecomplex *, integer *, char *, ftnlen);    extern doublereal dlapy2_(doublereal *, doublereal *);    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    static doublereal rnorm1;    extern /* Subroutine */ int zgetv0_(integer *, char *, integer *, logical 	    *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublereal *, integer *, doublecomplex *, 	    integer *, ftnlen);    extern doublereal dlamch_(char *, ftnlen);    extern /* Subroutine */ int second_(real *), zdscal_(integer *, 	    doublereal *, doublecomplex *, integer *);    static logical rstart;    static integer msglvl;    static doublereal smlnum;    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 	    doublecomplex *, ftnlen);    extern /* Subroutine */ int zlascl_(char *, integer *, integer *, 	    doublereal *, doublereal *, integer *, integer *, doublecomplex *,	     integer *, integer *, ftnlen);/*     %----------------------------------------------------% *//*     | Include files for debugging and timing information | *//*     %----------------------------------------------------% *//* /SCCS Information: @(#) *//* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %---------------------------------% *//*     | See debug.doc for documentation | *//*     %---------------------------------% *//*     %------------------% *//*     | Scalar Arguments | *//*     %------------------% *//*     %--------------------------------% *//*     | See stat.doc for documentation | *//*     %--------------------------------% *//* /SCCS Information: @(#) *//* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %-----------------% *//*     | Array Arguments | *//*     %-----------------% *//*     %------------% *//*     | Parameters | *//*     %------------% *///.........这里部分代码省略.........

/* Subroutine */ int zlaqps_(integer *m, integer *n, integer *offset, integer 	*nb, integer *kb, doublecomplex *a, integer *lda, integer *jpvt, 	doublecomplex *tau, doublereal *vn1, doublereal *vn2, doublecomplex *	auxv, doublecomplex *f, integer *ldf){    /* System generated locals */    integer a_dim1, a_offset, f_dim1, f_offset, i__1, i__2, i__3;    doublereal d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double sqrt(doublereal);    void d_cnjg(doublecomplex *, doublecomplex *);    double z_abs(doublecomplex *);    integer i_dnnt(doublereal *);    /* Local variables */    integer j, k, rk;    doublecomplex akk;    integer pvt;    doublereal temp, temp2, tol3z;    integer itemp;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zswap_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(	    char *);    extern integer idamax_(integer *, doublereal *, integer *);    integer lsticc;    extern /* Subroutine */ int zlarfp_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *);    integer lastrk;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLAQPS computes a step of QR factorization with column pivoting *//*  of a complex M-by-N matrix A by using Blas-3.  It tries to factorize *//*  NB columns from A starting from the row OFFSET+1, and updates all *//*  of the matrix with Blas-3 xGEMM. *//*  In some cases, due to catastrophic cancellations, it cannot *//*  factorize NB columns.  Hence, the actual number of factorized *//*  columns is returned in KB. *//*  Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A. M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A. N >= 0 *//*  OFFSET  (input) INTEGER *//*          The number of rows of A that have been factorized in *//*          previous steps. *//*  NB      (input) INTEGER *//*          The number of columns to factorize. *//*  KB      (output) INTEGER *//*          The number of columns actually factorized. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the M-by-N matrix A. *//*          On exit, block A(OFFSET+1:M,1:KB) is the triangular *//*          factor obtained and block A(1:OFFSET,1:N) has been *//*          accordingly pivoted, but no factorized. *//*          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has *//*          been updated. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= max(1,M). *//*  JPVT    (input/output) INTEGER array, dimension (N) *//*          JPVT(I) = K <==> Column K of the full matrix A has been *//*          permuted into position I in AP. *//*  TAU     (output) COMPLEX*16 array, dimension (KB) *//*          The scalar factors of the elementary reflectors. *//*  VN1     (input/output) DOUBLE PRECISION array, dimension (N) *///.........这里部分代码省略.........

/* Subroutine */ int zlauu2_(char *uplo, integer *n, doublecomplex *a, 	integer *lda, integer *info){/*  -- 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       Purpose       =======       ZLAUU2 computes the product U * U' or L' * L, where the triangular       factor U or L is stored in the upper or lower triangular part of       the array A.       If UPLO = 'U' or 'u' then the upper triangle of the result is stored,       overwriting the factor U in A.       If UPLO = 'L' or 'l' then the lower triangle of the result is stored,       overwriting the factor L in A.       This is the unblocked form of the algorithm, calling Level 2 BLAS.       Arguments       =========       UPLO    (input) CHARACTER*1               Specifies whether the triangular factor stored in the array A               is upper or lower triangular:               = 'U':  Upper triangular               = 'L':  Lower triangular       N       (input) INTEGER               The order of the triangular factor U or L.  N >= 0.       A       (input/output) COMPLEX*16 array, dimension (LDA,N)               On entry, the triangular factor U or L.               On exit, if UPLO = 'U', the upper triangle of A is               overwritten with the upper triangle of the product U * U';               if UPLO = 'L', the lower triangle of A is overwritten with               the lower triangle of the product L' * L.       LDA     (input) INTEGER               The leading dimension of the array A.  LDA >= max(1,N).       INFO    (output) INTEGER               = 0: successful exit               < 0: if INFO = -k, the k-th argument had an illegal value       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1;    /* Local variables */    static integer i__;    extern logical lsame_(char *, char *);    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static logical upper;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *), zlacgv_(	    integer *, doublecomplex *, integer *);    static doublereal aii;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    /* Function Body */    *info = 0;    upper = lsame_(uplo, "U");    if (! upper && ! lsame_(uplo, "L")) {	*info = -1;    } else if (*n < 0) {	*info = -2;    } else if (*lda < max(1,*n)) {	*info = -4;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("ZLAUU2", &i__1);	return 0;    }//.........这里部分代码省略.........

/* Subroutine */ int zlaghe_(integer *n, integer *k, doublereal *d, 	doublecomplex *a, integer *lda, integer *iseed, doublecomplex *work, 	integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double z_abs(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(	    doublecomplex *, doublecomplex *);    /* Local variables */    extern /* Subroutine */ int zher2_(char *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    static integer i, j;    static doublecomplex alpha;    extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zhemv_(char *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *), zaxpy_(integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    static doublecomplex wa, wb;    static doublereal wn;    extern /* Subroutine */ int xerbla_(char *, integer *), zlarnv_(	    integer *, integer *, integer *, doublecomplex *);    static doublecomplex tau;/*  -- LAPACK auxiliary test 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       =======       ZLAGHE generates a complex hermitian matrix A, by pre- and post-       multiplying a real diagonal matrix D with a random unitary matrix:       A = U*D*U'. The semi-bandwidth may then be reduced to k by additional       unitary transformations.       Arguments       =========       N       (input) INTEGER               The order of the matrix A.  N >= 0.       K       (input) INTEGER               The number of nonzero subdiagonals within the band of A.               0 <= K <= N-1.       D       (input) DOUBLE PRECISION array, dimension (N)               The diagonal elements of the diagonal matrix D.       A       (output) COMPLEX*16 array, dimension (LDA,N)               The generated n by n hermitian matrix A (the full matrix is               stored).       LDA     (input) INTEGER               The leading dimension of the array A.  LDA >= N.       ISEED   (input/output) INTEGER array, dimension (4)               On entry, the seed of the random number generator; the array               elements must be between 0 and 4095, and ISEED(4) must be               odd.               On exit, the seed is updated.       WORK    (workspace) COMPLEX*16 array, dimension (2*N)       INFO    (output) INTEGER               = 0: successful exit               < 0: if INFO = -i, the i-th argument had an illegal value       =====================================================================          Test the input arguments          Parameter adjustments */    --d;    a_dim1 = *lda;    a_offset = a_dim1 + 1;//.........这里部分代码省略.........

//.........这里部分代码省略.........            The leading dimension of the array W.  LDW >= max(1,N).       INFO    (output) INTEGER               = 0: successful exit               > 0: if INFO = k, D(k,k) is exactly zero.  The factorization                    has been completed, but the block diagonal matrix D is                    exactly singular.       =====================================================================          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double sqrt(doublereal), d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *), z_div(doublecomplex *, 	    doublecomplex *, doublecomplex *);    /* Local variables */    static integer imax, jmax, j, k;    static doublereal t, alpha;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    static integer kstep;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static doublereal r1;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zswap_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *);    static doublecomplex d11, d21, d22;    static integer jb, jj, kk, jp, kp;    static doublereal absakk;    static integer kw;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *);    static doublereal colmax;    extern /* Subroutine */ int zlacgv_(integer *, doublecomplex *, integer *)	    ;    extern integer izamax_(integer *, doublecomplex *, integer *);    static doublereal rowmax;    static integer kkw;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]#define w_subscr(a_1,a_2) (a_2)*w_dim1 + a_1#define w_ref(a_1,a_2) w[w_subscr(a_1,a_2)]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    --ipiv;    w_dim1 = *ldw;    w_offset = 1 + w_dim1 * 1;    w -= w_offset;

/* ----------------------------------------------------------------------| *//* Subroutine */ int zgexpv(integer *n, integer *m, doublereal *t, 	doublecomplex *v, doublecomplex *w, doublereal *tol, doublereal *	anorm, doublecomplex *wsp, integer *lwsp, integer *iwsp, integer *	liwsp, S_fp matvec, void *matvecdata, integer *itrace, integer *iflag){    /* System generated locals */    integer i__1, i__2, i__3;    doublereal d__1;    complex q__1;    doublecomplex z__1;    /* Builtin functions */    /* Subroutine */ int s_stop(char *, ftnlen);    double sqrt(doublereal), d_sign(doublereal *, doublereal *), pow_di(	    doublereal *, integer *), pow_dd(doublereal *, doublereal *), 	    d_lg10(doublereal *);    integer i_dnnt(doublereal *);    double d_int(doublereal *);    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 	    e_wsle();    double z_abs(doublecomplex *);    /* Local variables */    static integer ibrkflag;    static doublereal step_min__, step_max__;    static integer i__, j;    static doublereal break_tol__;    static integer k1;    static doublereal p1, p2, p3;    static integer ih, mh, iv, ns, mx;    static doublereal xm;    static integer j1v;    static doublecomplex hij;    static doublereal sgn, eps, hj1j, sqr1, beta, hump;    static integer ifree, lfree;    static doublereal t_old__;    static integer iexph;    static doublereal t_new__;    static integer nexph;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static doublereal t_now__;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, ftnlen);    static integer nstep;    static doublereal t_out__;    static integer nmult;    static doublereal vnorm;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zaxpy_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    static integer nscale;    static doublereal rndoff;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), zgpadm_(integer *, integer *, 	    doublereal *, doublecomplex *, integer *, doublecomplex *, 	    integer *, integer *, integer *, integer *, integer *), znchbv_(	    integer *, doublereal *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *);    static doublereal t_step__, avnorm;    static integer ireject;    static doublereal err_loc__;    static integer nreject, mbrkdwn;    static doublereal tbrkdwn, s_error__, x_error__;    /* Fortran I/O blocks */    static cilist io___40 = { 0, 6, 0, 0, 0 };    static cilist io___48 = { 0, 6, 0, 0, 0 };    static cilist io___49 = { 0, 6, 0, 0, 0 };    static cilist io___50 = { 0, 6, 0, 0, 0 };    static cilist io___51 = { 0, 6, 0, 0, 0 };    static cilist io___52 = { 0, 6, 0, 0, 0 };    static cilist io___53 = { 0, 6, 0, 0, 0 };    static cilist io___54 = { 0, 6, 0, 0, 0 };    static cilist io___55 = { 0, 6, 0, 0, 0 };    static cilist io___56 = { 0, 6, 0, 0, 0 };    static cilist io___57 = { 0, 6, 0, 0, 0 };    static cilist io___58 = { 0, 6, 0, 0, 0 };    static cilist io___59 = { 0, 6, 0, 0, 0 };/* -----Purpose----------------------------------------------------------| *//* ---  ZGEXPV computes w = exp(t*A)*v *//*     for a Zomplex (i.e., complex double precision) matrix A *//*     It does not compute the matrix exponential in isolation but *//*     instead, it computes directly the action of the exponential *//*     operator on the operand vector. This way of doing so allows *//*     for addressing large sparse problems. *//*     The method used is based on Krylov subspace projection *//*     techniques and the matrix under consideration interacts only *//*     via the external routine `matvec' performing the matrix-vector *//*     product (matrix-free method). *//* -----Arguments--------------------------------------------------------| *///.........这里部分代码省略.........

 int zlabrd_(int *m, int *n, int *nb, 	doublecomplex *a, int *lda, double *d__, double *e, 	doublecomplex *tauq, doublecomplex *taup, doublecomplex *x, int *	ldx, doublecomplex *y, int *ldy){    /* System generated locals */    int a_dim1, a_offset, x_dim1, x_offset, y_dim1, y_offset, i__1, i__2, 	    i__3;    doublecomplex z__1;    /* Local variables */    int i__;    doublecomplex alpha;    extern  int zscal_(int *, doublecomplex *, 	    doublecomplex *, int *), zgemv_(char *, int *, int *, 	    doublecomplex *, doublecomplex *, int *, doublecomplex *, 	    int *, doublecomplex *, doublecomplex *, int *), 	    zlarfg_(int *, doublecomplex *, doublecomplex *, int *, 	    doublecomplex *), zlacgv_(int *, doublecomplex *, int *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLABRD reduces the first NB rows and columns of a complex general *//*  m by n matrix A to upper or lower float bidiagonal form by a unitary *//*  transformation Q' * A * P, and returns the matrices X and Y which *//*  are needed to apply the transformation to the unreduced part of A. *//*  If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower *//*  bidiagonal form. *//*  This is an auxiliary routine called by ZGEBRD *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows in the matrix A. *//*  N       (input) INTEGER *//*          The number of columns in the matrix A. *//*  NB      (input) INTEGER *//*          The number of leading rows and columns of A to be reduced. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the m by n general matrix to be reduced. *//*          On exit, the first NB rows and columns of the matrix are *//*          overwritten; the rest of the array is unchanged. *//*          If m >= n, elements on and below the diagonal in the first NB *//*            columns, with the array TAUQ, represent the unitary *//*            matrix Q as a product of elementary reflectors; and *//*            elements above the diagonal in the first NB rows, with the *//*            array TAUP, represent the unitary matrix P as a product *//*            of elementary reflectors. *//*          If m < n, elements below the diagonal in the first NB *//*            columns, with the array TAUQ, represent the unitary *//*            matrix Q as a product of elementary reflectors, and *//*            elements on and above the diagonal in the first NB rows, *//*            with the array TAUP, represent the unitary matrix P as *//*            a product of elementary reflectors. *//*          See Further Details. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= MAX(1,M). *//*  D       (output) DOUBLE PRECISION array, dimension (NB) *//*          The diagonal elements of the first NB rows and columns of *//*          the reduced matrix.  D(i) = A(i,i). *//*  E       (output) DOUBLE PRECISION array, dimension (NB) *//*          The off-diagonal elements of the first NB rows and columns of *//*          the reduced matrix. *//*  TAUQ    (output) COMPLEX*16 array dimension (NB) *//*          The scalar factors of the elementary reflectors which *//*          represent the unitary matrix Q. See Further Details. *//*  TAUP    (output) COMPLEX*16 array, dimension (NB) *//*          The scalar factors of the elementary reflectors which *//*          represent the unitary matrix P. See Further Details. *//*  X       (output) COMPLEX*16 array, dimension (LDX,NB) *//*          The m-by-nb matrix X required to update the unreduced part *//*          of A. *//*  LDX     (input) INTEGER *//*          The leading dimension of the array X. LDX >= MAX(1,M). *//*  Y       (output) COMPLEX*16 array, dimension (LDY,NB) *///.........这里部分代码省略.........

/* Subroutine */ int zgetv0_(integer *ido, char *bmat, integer *itry, logical 	*initv, integer *n, integer *j, doublecomplex *v, integer *ldv, 	doublecomplex *resid, doublereal *rnorm, integer *ipntr, 	doublecomplex *workd, integer *ierr, ftnlen bmat_len){    /* Initialized data */    static logical inits = TRUE_;    /* System generated locals */    integer v_dim1, v_offset, i__1, i__2;    doublereal d__1, d__2;    doublecomplex z__1;    /* Local variables */    static real t0, t1, t2, t3;    static integer jj, iter;    static logical orth;    static integer iseed[4], idist;    static doublecomplex cnorm;    extern /* Double Complex */ void zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static logical first;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, ftnlen), 	    dvout_(integer *, integer *, doublereal *, integer *, char *, 	    ftnlen), zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zvout_(integer *, integer *, 	    doublecomplex *, integer *, char *, ftnlen);    extern doublereal dlapy2_(doublereal *, doublereal *), dznrm2_(integer *, 	    doublecomplex *, integer *);    static doublereal rnorm0;    extern /* Subroutine */ int arscnd_(real *);    static integer msglvl;    extern /* Subroutine */ int zlarnv_(integer *, integer *, integer *, 	    doublecomplex *);/*     %----------------------------------------------------% *//*     | Include files for debugging and timing information | *//*     %----------------------------------------------------% *//* /SCCS Information: @(#) *//* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %---------------------------------% *//*     | See debug.doc for documentation | *//*     %---------------------------------% *//*     %------------------% *//*     | Scalar Arguments | *//*     %------------------% *//*     %--------------------------------% *//*     | See stat.doc for documentation | *//*     %--------------------------------% *//* /SCCS Information: @(#) *//* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %-----------------% *//*     | Array Arguments | *//*     %-----------------% *//*     %------------% *//*     | Parameters | *//*     %------------% *//*     %------------------------% *//*     | Local Scalars & Arrays | *//*     %------------------------% *//*     %----------------------% *//*     | External Subroutines | *//*     %----------------------% *//*     %--------------------% *//*     | External Functions | *//*     %--------------------% *//*     %-----------------% *//*     | Data Statements | *//*     %-----------------% */    /* Parameter adjustments */    --workd;    --resid;    v_dim1 = *ldv;    v_offset = 1 + v_dim1;    v -= v_offset;    --ipntr;//.........这里部分代码省略.........

//.........这里部分代码省略.........    =====================================================================          Decode and test the input parameters          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;    static doublecomplex c_b19 = {1.,0.};    static doublecomplex c_b20 = {0.,0.};    static logical c_false = FALSE_;    static integer c__3 = 3;        /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1, 	    vr_offset, i__1, i__2;    doublereal d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double z_abs(doublecomplex *);    /* Local variables */    static doublereal cond;    static integer ierr, ifst;    static doublereal lnrm;    static doublecomplex yhax, yhbx;    static integer ilst;    static doublereal rnrm;    static integer i__, k;    static doublereal scale;    extern logical lsame_(char *, char *);    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static integer lwmin;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static logical wants;    static integer llwrk, n1, n2;    static doublecomplex dummy[1];    extern doublereal dlapy2_(doublereal *, doublereal *);    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);    static doublecomplex dummy1[1];    extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(	    char *);    static integer ks;    extern /* Subroutine */ int xerbla_(char *, integer *);    static doublereal bignum;    static logical wantbh, wantdf, somcon;    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    ztgexc_(logical *, logical *, integer *, doublecomplex *, integer 	    *, doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, integer *, integer *, integer *);    static doublereal smlnum;    static logical lquery;    extern /* Subroutine */ int ztgsyl_(char *, integer *, integer *, integer 	    *, doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *, doublecomplex *, integer *, integer *,	     integer *);    static doublereal eps;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]#define b_subscr(a_1,a_2) (a_2)*b_dim1 + a_1#define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]

//.........这里部分代码省略.........            On entry, the (m - 1) x n matrix C2 if SIDE = 'L', or the               m x (n - 1) matrix C2 if SIDE = 'R'.               On exit, rows 2:m of P*C if SIDE = 'L', or columns 2:m of C*P               if SIDE = 'R'.       LDC     (input) INTEGER               The leading dimension of the arrays C1 and C2.               LDC >= max(1,M).       WORK    (workspace) COMPLEX*16 array, dimension                           (N) if SIDE = 'L'                           (M) if SIDE = 'R'       =====================================================================          Parameter adjustments          Function Body */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer c1_dim1, c1_offset, c2_dim1, c2_offset, i__1;    doublecomplex z__1;    /* Local variables */    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 *), 	    zgeru_(integer *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *)	    , zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zaxpy_(integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *), zlacgv_(integer *, 	    doublecomplex *, integer *);#define V(I) v[(I)-1]#define WORK(I) work[(I)-1]#define C2(I,J) c2[(I)-1 + ((J)-1)* ( *ldc)]#define C1(I,J) c1[(I)-1 + ((J)-1)* ( *ldc)]    if (min(*m,*n) == 0 || tau->r == 0. && tau->i == 0.) {	return 0;    }    if (lsame_(side, "L")) {/*        w :=  conjg( C1 + v' * C2 ) */	zcopy_(n, &C1(1,1), ldc, &WORK(1), &c__1);	zlacgv_(n, &WORK(1), &c__1);	i__1 = *m - 1;	zgemv_("Conjugate transpose", &i__1, n, &c_b1, &C2(1,1), ldc, &		V(1), incv, &c_b1, &WORK(1), &c__1);/*        [ C1 ] := [ C1 ] - tau* [ 1 ] * w'   

/* Subroutine */ int zlarf_(char *side, integer *m, integer *n, doublecomplex 	*v, integer *incv, doublecomplex *tau, doublecomplex *c__, integer *	ldc, doublecomplex *work){    /* System generated locals */    integer c_dim1, c_offset, i__1;    doublecomplex z__1;    /* Local variables */    integer i__;    logical applyleft;    extern logical lsame_(char *, char *);    integer lastc;    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 *);    integer lastv;    extern integer ilazlc_(integer *, integer *, doublecomplex *, integer *), 	    ilazlr_(integer *, integer *, doublecomplex *, integer *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLARF 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. *//*  To apply H' (the conjugate transpose of H), supply conjg(tau) instead *//*  tau. *//*  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 *//*                     (1 + (M-1)*abs(INCV)) if SIDE = 'L' *//*                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R' *//*          The vector v in the representation of H. V is not used if *//*          TAU = 0. *//*  INCV    (input) INTEGER *//*          The increment between elements of v. INCV <> 0. *//*  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. LDC >= lmax(1,M). *//*  WORK    (workspace) COMPLEX*16 array, dimension *//*                         (N) if SIDE = 'L' *//*                      or (M) if SIDE = 'R' *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. Executable Statements .. */    /* Parameter adjustments */    --v;    c_dim1 = *ldc;//.........这里部分代码省略.........

 int zlatrd_(char *uplo, int *n, int *nb, 	doublecomplex *a, int *lda, double *e, doublecomplex *tau, 	doublecomplex *w, int *ldw){    /* System generated locals */    int a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3;    double d__1;    doublecomplex z__1, z__2, z__3, z__4;    /* Local variables */    int i__, iw;    doublecomplex alpha;    extern int lsame_(char *, char *);    extern  int zscal_(int *, doublecomplex *, 	    doublecomplex *, int *);    extern /* Double Complex */ VOID zdotc_(doublecomplex *, int *, 	    doublecomplex *, int *, doublecomplex *, int *);    extern  int zgemv_(char *, int *, int *, 	    doublecomplex *, doublecomplex *, int *, doublecomplex *, 	    int *, doublecomplex *, doublecomplex *, int *), 	    zhemv_(char *, int *, doublecomplex *, doublecomplex *, 	    int *, doublecomplex *, int *, doublecomplex *, 	    doublecomplex *, int *), zaxpy_(int *, 	    doublecomplex *, doublecomplex *, int *, doublecomplex *, 	    int *), zlarfg_(int *, doublecomplex *, doublecomplex *, 	    int *, doublecomplex *), zlacgv_(int *, doublecomplex *, 	    int *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to *//*  Hermitian tridiagonal form by a unitary similarity *//*  transformation Q' * A * Q, and returns the matrices V and W which are *//*  needed to apply the transformation to the unreduced part of A. *//*  If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a *//*  matrix, of which the upper triangle is supplied; *//*  if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a *//*  matrix, of which the lower triangle is supplied. *//*  This is an auxiliary routine called by ZHETRD. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          Specifies whether the upper or lower triangular part of the *//*          Hermitian matrix A is stored: *//*          = 'U': Upper triangular *//*          = 'L': Lower triangular *//*  N       (input) INTEGER *//*          The order of the matrix A. *//*  NB      (input) INTEGER *//*          The number of rows and columns to be reduced. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading *//*          n-by-n upper triangular part of A contains the upper *//*          triangular part of the matrix A, and the strictly lower *//*          triangular part of A is not referenced.  If UPLO = 'L', the *//*          leading n-by-n lower triangular part of A contains the lower *//*          triangular part of the matrix A, and the strictly upper *//*          triangular part of A is not referenced. *//*          On exit: *//*          if UPLO = 'U', the last NB columns have been reduced to *//*            tridiagonal form, with the diagonal elements overwriting *//*            the diagonal elements of A; the elements above the diagonal *//*            with the array TAU, represent the unitary matrix Q as a *//*            product of elementary reflectors; *//*          if UPLO = 'L', the first NB columns have been reduced to *//*            tridiagonal form, with the diagonal elements overwriting *//*            the diagonal elements of A; the elements below the diagonal *//*            with the array TAU, represent the  unitary matrix Q as a *//*            product of elementary reflectors. *//*          See Further Details. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= MAX(1,N). *//*  E       (output) DOUBLE PRECISION array, dimension (N-1) *//*          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal *//*          elements of the last NB columns of the reduced matrix; *//*          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of *//*          the first NB columns of the reduced matrix. *//*  TAU     (output) COMPLEX*16 array, dimension (N-1) *//*          The scalar factors of the elementary reflectors, stored in *///.........这里部分代码省略.........

//.........这里部分代码省略.........                        if ((d__1 = sum.r, abs(d__1)) + (d__2 = d_imag(&sum),                                                         abs(d__2)) >= bignum * ((d__3 = d__.r, abs(                                                                 d__3)) + (d__4 = d_imag(&d__), abs(d__4)))) {                            temp = 1. / ((d__1 = sum.r, abs(d__1)) + (d__2 =                                             d_imag(&sum), abs(d__2)));                            i__3 = j - 1;                            for (jr = je; jr <= i__3; ++jr) {                                i__4 = jr;                                i__5 = jr;                                z__1.r = temp * work[i__5].r, z__1.i = temp *                                                                       work[i__5].i;                                work[i__4].r = z__1.r, work[i__4].i = z__1.i;                            }                            xmax = temp * xmax;                            z__1.r = temp * sum.r, z__1.i = temp * sum.i;                            sum.r = z__1.r, sum.i = z__1.i;                        }                    }                    i__3 = j;                    z__2.r = -sum.r, z__2.i = -sum.i;                    zladiv_(&z__1, &z__2, &d__);                    work[i__3].r = z__1.r, work[i__3].i = z__1.i;                    /* Computing MAX */                    i__3 = j;                    d__3 = xmax, d__4 = (d__1 = work[i__3].r, abs(d__1)) + (                                            d__2 = d_imag(&work[j]), abs(d__2));                    xmax = max(d__3,d__4);                }                /*              Back transform eigenvector if HOWMNY='B'. */                if (ilback) {                    i__2 = *n + 1 - je;                    zgemv_("N", n, &i__2, &c_b2, &vl[je * vl_dim1 + 1], ldvl,                           &work[je], &c__1, &c_b1, &work[*n + 1], &c__1);                    isrc = 2;                    ibeg = 1;                } else {                    isrc = 1;                    ibeg = je;                }                /*              Copy and scale eigenvector into column of VL */                xmax = 0.;                i__2 = *n;                for (jr = ibeg; jr <= i__2; ++jr) {                    /* Computing MAX */                    i__3 = (isrc - 1) * *n + jr;                    d__3 = xmax, d__4 = (d__1 = work[i__3].r, abs(d__1)) + (                                            d__2 = d_imag(&work[(isrc - 1) * *n + jr]), abs(                                                d__2));                    xmax = max(d__3,d__4);                }                if (xmax > safmin) {                    temp = 1. / xmax;                    i__2 = *n;                    for (jr = ibeg; jr <= i__2; ++jr) {                        i__3 = jr + ieig * vl_dim1;                        i__4 = (isrc - 1) * *n + jr;                        z__1.r = temp * work[i__4].r, z__1.i = temp * work[                                i__4].i;                        vl[i__3].r = z__1.r, vl[i__3].i = z__1.i;                    }                } else {

//.........这里部分代码省略.........		 * Case: sup-col update		 * Perform a triangular solve and block update,		 * then scatter the result of sup-col update to dense		 */		no_zeros = kfnz - fst_col;	        /* Copy U[*,j] segment from dense[*] to tempv[*] */	        isub = lptr + no_zeros;	        for (i = 0; i < segsze; i++) {	  	    irow = lsub[isub];		    tempv[i] = dense[irow];		    ++isub; 	        }	        /* Dense triangular solve -- start effective triangle */		luptr += nsupr * no_zeros + no_zeros; #ifdef USE_VENDOR_BLAS#if ( MACH==CRAY_PVP )		CTRSV( ftcs1, ftcs2, ftcs3, &segsze, &lusup[luptr], 		       &nsupr, tempv, &incx );#else		ztrsv_( "L", "N", "U", &segsze, &lusup[luptr], 		       &nsupr, tempv, &incx );#endif		 		luptr += segsze;  /* Dense matrix-vector */		tempv1 = &tempv[segsze];		alpha = one;		beta = zero;#if ( MACH==CRAY_PVP )		CGEMV( ftcs2, &nrow, &segsze, &alpha, &lusup[luptr], 		       &nsupr, tempv, &incx, &beta, tempv1, &incy );#else		zgemv_( "N", &nrow, &segsze, &alpha, &lusup[luptr], 		       &nsupr, tempv, &incx, &beta, tempv1, &incy );#endif#else		zlsolve ( nsupr, segsze, &lusup[luptr], tempv ); 		luptr += segsze;  /* Dense matrix-vector */		tempv1 = &tempv[segsze];		zmatvec (nsupr, nrow , segsze, &lusup[luptr], tempv, tempv1);#endif                /* Scatter tempv[] into SPA dense[*] */                isub = lptr + no_zeros;                for (i = 0; i < segsze; i++) {                    irow = lsub[isub];                    dense[irow] = tempv[i]; /* Scatter */                    tempv[i] = zero;                    isub++;                }		/* Scatter tempv1[] into SPA dense[*] */		for (i = 0; i < nrow; i++) {		    irow = lsub[isub];                    z_sub(&dense[irow], &dense[irow], &tempv1[i]);		    tempv1[i] = zero;		    ++isub;		}	    } /* else segsze >= 4 */	    	} /* if jsupno ... */    } /* for each segment... */    

voidzgemv(char transa, int m, int n, doublecomplex *alpha, doublecomplex *a, int lda, doublecomplex *x, int incx, doublecomplex *beta, doublecomplex *y, int incy){    zgemv_(&transa, &m, &n, alpha, a, &lda, x, &incx, beta, y, &incy);}

 int zggglm_(int *n, int *m, int *p, 	doublecomplex *a, int *lda, doublecomplex *b, int *ldb, 	doublecomplex *d__, doublecomplex *x, doublecomplex *y, doublecomplex 	*work, int *lwork, int *info){    /* System generated locals */    int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;    doublecomplex z__1;    /* Local variables */    int i__, nb, np, nb1, nb2, nb3, nb4, lopt;    extern  int zgemv_(char *, int *, int *, 	    doublecomplex *, doublecomplex *, int *, doublecomplex *, 	    int *, doublecomplex *, doublecomplex *, int *), 	    zcopy_(int *, doublecomplex *, int *, doublecomplex *, 	    int *), xerbla_(char *, int *);    extern int ilaenv_(int *, char *, char *, int *, int *, 	    int *, int *);    extern  int zggqrf_(int *, int *, int *, 	    doublecomplex *, int *, doublecomplex *, doublecomplex *, 	    int *, doublecomplex *, doublecomplex *, int *, int *)	    ;    int lwkmin, lwkopt;    int lquery;    extern  int zunmqr_(char *, char *, int *, int *, 	    int *, doublecomplex *, int *, doublecomplex *, 	    doublecomplex *, int *, doublecomplex *, int *, int *), zunmrq_(char *, char *, int *, int *, 	    int *, doublecomplex *, int *, doublecomplex *, 	    doublecomplex *, int *, doublecomplex *, int *, int *), ztrtrs_(char *, char *, char *, int *, 	    int *, doublecomplex *, int *, doublecomplex *, int *, 	     int *);/*  -- LAPACK driver routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGGGLM solves a general Gauss-Markov linear model (GLM) problem: *//*          minimize || y ||_2   subject to   d = A*x + B*y *//*              x *//*  where A is an N-by-M matrix, B is an N-by-P matrix, and d is a *//*  given N-vector. It is assumed that M <= N <= M+P, and *//*             rank(A) = M    and    rank( A B ) = N. *//*  Under these assumptions, the constrained equation is always *//*  consistent, and there is a unique solution x and a minimal 2-norm *//*  solution y, which is obtained using a generalized QR factorization *//*  of the matrices (A, B) given by *//*     A = Q*(R),   B = Q*T*Z. *//*           (0) *//*  In particular, if matrix B is square nonsingular, then the problem *//*  GLM is equivalent to the following weighted linear least squares *//*  problem *//*               minimize || inv(B)*(d-A*x) ||_2 *//*                   x *//*  where inv(B) denotes the inverse of B. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The number of rows of the matrices A and B.  N >= 0. *//*  M       (input) INTEGER *//*          The number of columns of the matrix A.  0 <= M <= N. *//*  P       (input) INTEGER *//*          The number of columns of the matrix B.  P >= N-M. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,M) *//*          On entry, the N-by-M matrix A. *//*          On exit, the upper triangular part of the array A contains *//*          the M-by-M upper triangular matrix R. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= MAX(1,N). *//*  B       (input/output) COMPLEX*16 array, dimension (LDB,P) *//*          On entry, the N-by-P matrix B. *//*          On exit, if N <= P, the upper triangle of the subarray *//*          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; *//*          if N > P, the elements on and above the (N-P)th subdiagonal *//*          contain the N-by-P upper trapezoidal matrix T. *//*  LDB     (input) INTEGER *///.........这里部分代码省略.........

/* Subroutine */int ztgevc_(char *side, char *howmny, logical *select, integer *n, doublecomplex *s, integer *lds, doublecomplex *p, integer *ldp, doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer * ldvr, integer *mm, integer *m, doublecomplex *work, doublereal *rwork, integer *info){    /* System generated locals */    integer p_dim1, p_offset, s_dim1, s_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3, d__4, d__5, d__6;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    doublecomplex d__;    integer i__, j;    doublecomplex ca, cb;    integer je, im, jr;    doublereal big;    logical lsa, lsb;    doublereal ulp;    doublecomplex sum;    integer ibeg, ieig, iend;    doublereal dmin__;    integer isrc;    doublereal temp;    doublecomplex suma, sumb;    doublereal xmax, scale;    logical ilall;    integer iside;    doublereal sbeta;    extern logical lsame_(char *, char *);    doublereal small;    logical compl;    doublereal anorm, bnorm;    logical compr;    extern /* Subroutine */    int zgemv_(char *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *), dlabad_(doublereal *, doublereal *);    logical ilbbad;    doublereal acoefa, bcoefa, acoeff;    doublecomplex bcoeff;    logical ilback;    doublereal ascale, bscale;    extern doublereal dlamch_(char *);    doublecomplex salpha;    doublereal safmin;    extern /* Subroutine */    int xerbla_(char *, integer *);    doublereal bignum;    logical ilcomp;    extern /* Double Complex */    VOID zladiv_(doublecomplex *, doublecomplex *, doublecomplex *);    integer ihwmny;    /* -- LAPACK computational routine (version 3.4.0) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* November 2011 */    /* .. Scalar Arguments .. */    /* .. */    /* .. Array Arguments .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. Statement Functions .. */    /* .. */    /* .. Statement Function definitions .. */    /* .. */    /* .. Executable Statements .. */    /* Decode and Test the input parameters */    /* Parameter adjustments */    --select;    s_dim1 = *lds;    s_offset = 1 + s_dim1;    s -= s_offset;    p_dim1 = *ldp;    p_offset = 1 + p_dim1;    p -= p_offset;    vl_dim1 = *ldvl;    vl_offset = 1 + vl_dim1;    vl -= vl_offset;    vr_dim1 = *ldvr;    vr_offset = 1 + vr_dim1;    vr -= vr_offset;    --work;    --rwork;    /* Function Body */    if (lsame_(howmny, "A"))    {        ihwmny = 1;        ilall = TRUE_;        ilback = FALSE_;    }    else if (lsame_(howmny, "S"))    {//.........这里部分代码省略.........

/* Subroutine */ int zhptrs_(char *uplo, integer *n, integer *nrhs, 	doublecomplex *ap, integer *ipiv, doublecomplex *b, integer *ldb, 	integer *info){    /* System generated locals */    integer b_dim1, b_offset, i__1, i__2;    doublecomplex z__1, z__2, z__3;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(	    doublecomplex *, doublecomplex *);    /* Local variables */    integer j, k;    doublereal s;    doublecomplex ak, bk;    integer kc, kp;    doublecomplex akm1, bkm1, akm1k;    extern logical lsame_(char *, char *);    doublecomplex denom;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    logical upper;    extern /* Subroutine */ int zgeru_(integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zswap_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), xerbla_(char *, integer *), zdscal_(integer *, doublereal *, doublecomplex *, 	    integer *), zlacgv_(integer *, doublecomplex *, integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZHPTRS solves a system of linear equations A*X = B with a complex *//*  Hermitian matrix A stored in packed format using the factorization *//*  A = U*D*U**H or A = L*D*L**H computed by ZHPTRF. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          Specifies whether the details of the factorization are stored *//*          as an upper or lower triangular matrix. *//*          = 'U':  Upper triangular, form is A = U*D*U**H; *//*          = 'L':  Lower triangular, form is A = L*D*L**H. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  NRHS    (input) INTEGER *//*          The number of right hand sides, i.e., the number of columns *//*          of the matrix B.  NRHS >= 0. *//*  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2) *//*          The block diagonal matrix D and the multipliers used to *//*          obtain the factor U or L as computed by ZHPTRF, stored as a *//*          packed triangular matrix. *//*  IPIV    (input) INTEGER array, dimension (N) *//*          Details of the interchanges and the block structure of D *//*          as determined by ZHPTRF. *//*  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS) *//*          On entry, the right hand side matrix B. *//*          On exit, the solution matrix X. *//*  LDB     (input) INTEGER *//*          The leading dimension of the array B.  LDB >= max(1,N). *//*  INFO    (output) INTEGER *//*          = 0:  successful exit *//*          < 0: if INFO = -i, the i-th argument had an illegal value *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. */    /* Parameter adjustments */    --ap;    --ipiv;//.........这里部分代码省略.........

//.........这里部分代码省略.........	    /* Form x := inv(L)*x */    	    if ( L->nrow == 0 ) return 0; /* Quick return */	    	    for (k = 0; k <= Lstore->nsuper; k++) {		fsupc = L_FST_SUPC(k);		istart = L_SUB_START(fsupc);		nsupr = L_SUB_START(fsupc+1) - istart;		nsupc = L_FST_SUPC(k+1) - fsupc;		luptr = L_NZ_START(fsupc);		nrow = nsupr - nsupc;	        solve_ops += 4 * nsupc * (nsupc - 1);	        solve_ops += 8 * nrow * nsupc;		if ( nsupc == 1 ) {		    for (iptr=istart+1; iptr < L_SUB_START(fsupc+1); ++iptr) {			irow = L_SUB(iptr);			++luptr;			zz_mult(&comp_zero, &x[fsupc], &Lval[luptr]);			z_sub(&x[irow], &x[irow], &comp_zero);		    }		} else {#ifdef USE_VENDOR_BLAS#ifdef _CRAY		    CTRSV(ftcs1, ftcs2, ftcs3, &nsupc, &Lval[luptr], &nsupr,		       	&x[fsupc], &incx);				    CGEMV(ftcs2, &nrow, &nsupc, &alpha, &Lval[luptr+nsupc], 		       	&nsupr, &x[fsupc], &incx, &beta, &work[0], &incy);#else		    ztrsv_("L", "N", "U", &nsupc, &Lval[luptr], &nsupr,		       	&x[fsupc], &incx);				    zgemv_("N", &nrow, &nsupc, &alpha, &Lval[luptr+nsupc], 		       	&nsupr, &x[fsupc], &incx, &beta, &work[0], &incy);#endif#else		    zlsolve ( nsupr, nsupc, &Lval[luptr], &x[fsupc]);				    zmatvec ( nsupr, nsupr-nsupc, nsupc, &Lval[luptr+nsupc],			&x[fsupc], &work[0] );#endif						    iptr = istart + nsupc;		    for (i = 0; i < nrow; ++i, ++iptr) {			irow = L_SUB(iptr);			z_sub(&x[irow], &x[irow], &work[i]); /* Scatter */			work[i] = comp_zero;		    }	 	}	    } /* for k ... */	    	} else {	    /* Form x := inv(U)*x */	    	    if ( U->nrow == 0 ) return 0; /* Quick return */	    	    for (k = Lstore->nsuper; k >= 0; k--) {	    	fsupc = L_FST_SUPC(k);	    	nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);	    	nsupc = L_FST_SUPC(k+1) - fsupc;	    	luptr = L_NZ_START(fsupc);		    	        solve_ops += 4 * nsupc * (nsupc + 1);

/* Subroutine */ int zlaror_slu(char *side, char *init, integer *m, integer *n, 	doublecomplex *a, integer *lda, integer *iseed, doublecomplex *x, 	integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublecomplex z__1, z__2;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static integer kbeg, jcol;    static doublereal xabs;    static integer irow, j;    static doublecomplex csign;    extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static integer ixfrm;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static integer itype, nxfrm;    static doublereal xnorm;    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    extern int input_error(char *, int *);    static doublereal factor;    extern /* Subroutine */ int zlacgv_slu(integer *, doublecomplex *, integer *)	    ;    extern /* Double Complex */ VOID zlarnd_slu(doublecomplex *, integer *, 	    integer *);    extern /* Subroutine */ int zlaset_slu(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    static doublecomplex xnorms;/*  -- LAPACK auxiliary test 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       =======          ZLAROR pre- or post-multiplies an M by N matrix A by a random          unitary matrix U, overwriting A. A may optionally be          initialized to the identity matrix before multiplying by U.          U is generated using the method of G.W. Stewart          ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ).          (BLAS-2 version)       Arguments       =========       SIDE   - CHARACTER*1                SIDE specifies whether A is multiplied on the left or right                by U.            SIDE = 'L'   Multiply A on the left (premultiply) by U            SIDE = 'R'   Multiply A on the right (postmultiply) by U*            SIDE = 'C'   Multiply A on the left by U and the right by U*            SIDE = 'T'   Multiply A on the left by U and the right by U'                Not modified.       INIT   - CHARACTER*1                INIT specifies whether or not A should be initialized to                the identity matrix.                   INIT = 'I'   Initialize A to (a section of) the                                identity matrix before applying U.                   INIT = 'N'   No initialization.  Apply U to the                                input matrix A.                INIT = 'I' may be used to generate square (i.e., unitary)                or rectangular orthogonal matrices (orthogonality being                in the sense of ZDOTC):                For square matrices, M=N, and SIDE many be either 'L' or                'R'; the rows will be orthogonal to each other, as will the                columns.                For rectangular matrices where M < N, SIDE = 'R' will                produce a dense matrix whose rows will be orthogonal and                whose columns will not, while SIDE = 'L' will produce a                matrix whose rows will be orthogonal, and whose first M                columns will be orthogonal, the remaining columns being                zero.                For matrices where M > N, just use the previous                explaination, interchanging 'L' and 'R' and "rows" and                "columns".                Not modified.       M      - INTEGER                Number of rows of A. Not modified.       N      - INTEGER   //.........这里部分代码省略.........

//.........这里部分代码省略.........                     ( v1 v2 v3 )       DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':                    V = ( v1 v2 v3 )                 V = ( v1 v1  1       )                        ( v1 v2 v3 )                     ( v2 v2 v2  1    )                        (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )                        (     1 v3 )                        (        1 )       =====================================================================          Quick return if possible          Parameter adjustments          Function Body */    /* Table of constant values */    static doublecomplex c_b2 = {0.,0.};    static integer c__1 = 1;        /* System generated locals */    integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;    doublecomplex z__1;    /* Local variables */    static integer i, j;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    ztrmv_(char *, char *, char *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), 	    zlacgv_(integer *, doublecomplex *, integer *);    static doublecomplex vii;#define TAU(I) tau[(I)-1]#define V(I,J) v[(I)-1 + ((J)-1)* ( *ldv)]#define T(I,J) t[(I)-1 + ((J)-1)* ( *ldt)]    if (*n == 0) {	return 0;    }    if (lsame_(direct, "F")) {	i__1 = *k;	for (i = 1; i <= *k; ++i) {	    i__2 = i;	    if (TAU(i).r == 0. && TAU(i).i == 0.) {/*              H(i)  =  I */		i__2 = i;		for (j = 1; j <= i; ++j) {		    i__3 = j + i * t_dim1;		    T(j,i).r = 0., T(j,i).i = 0.;/* L10: */		}