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

intf2c_zscal(integer* N,          doublecomplex* alpha,          doublecomplex* X, integer* incX){    zscal_(N, alpha, X, incX);    return 0;}

/*! _zgbmatrix*std::complex<double> operator */inline _zgbmatrix operator*(const _zgbmatrix& mat, const std::complex<double>& d){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] operator*(const _zgbmatrix&, const std::complex<double>&)"            << std::endl;#endif//CPPL_VERBOSE    zscal_((mat.KL+mat.KU+1)*mat.N, d, mat.Array, 1);  return mat;}

/*! zssmatrix*=std::complex<double> operator */inline zssmatrix& zssmatrix::operator*=(const std::complex<double>& d){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] zssmatrix::operator*=(const std::complex<double>&)"            << std::endl;#endif//CPPL_VERBOSE    zscal_(VOL, d, Array, 1);  return *this;}

/*! zcovector*=std::complex<double> operator */inline zcovector& zcovector::operator*=(const std::complex<double>& d){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] zcovector::operator*=(const std::complex<double>&)"            << std::endl;#endif//CPPL_VERBOSE    zscal_(L, d, Array, 1);  return *this;}

/*! std::complex<double>*_zrovector operator */inline _zrovector operator*(const std::complex<double>& d, const _zrovector& vec){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] operator*(const std::complex<double>&, const _zrovector&)"            << std::endl;#endif//CPPL_VERBOSE    zscal_(vec.L, d, vec.Array, 1);    return vec;}

PyObject* scal(PyObject *self, PyObject *args){  Py_complex alpha;  PyArrayObject* x;  if (!PyArg_ParseTuple(args, "DO", &alpha, &x))    return NULL;  int n = PyArray_DIMS(x)[0];  for (int d = 1; d < PyArray_NDIM(x); d++)    n *= PyArray_DIMS(x)[d];  int incx = 1;  if (PyArray_DESCR(x)->type_num == NPY_DOUBLE)    dscal_(&n, &(alpha.real), DOUBLEP(x), &incx);  else    zscal_(&n, &alpha, (void*)COMPLEXP(x), &incx);  Py_RETURN_NONE;}

/* Subroutine */int zlarfgp_(integer *n, doublecomplex *alpha, doublecomplex *x, integer *incx, doublecomplex *tau){    /* System generated locals */    integer i__1, i__2;    doublereal d__1, d__2;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *), d_sign(doublereal *, doublereal *), z_abs( doublecomplex *);    /* Local variables */    integer j;    doublecomplex savealpha;    integer knt;    doublereal beta, alphi, alphr;    extern /* Subroutine */    int zscal_(integer *, doublecomplex *, doublecomplex *, integer *);    doublereal xnorm;    extern doublereal dlapy2_(doublereal *, doublereal *), dlapy3_(doublereal *, doublereal *, doublereal *), dznrm2_(integer *, doublecomplex * , integer *), dlamch_(char *);    extern /* Subroutine */    int zdscal_(integer *, doublereal *, doublecomplex *, integer *);    doublereal bignum;    extern /* Double Complex */    VOID zladiv_(doublecomplex *, doublecomplex *, doublecomplex *);    doublereal smlnum;    /* -- LAPACK auxiliary routine (version 3.4.2) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* September 2012 */    /* .. Scalar Arguments .. */    /* .. */    /* .. Array Arguments .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Functions .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Executable Statements .. */    /* Parameter adjustments */    --x;    /* Function Body */    if (*n <= 0)    {        tau->r = 0., tau->i = 0.;        return 0;    }    i__1 = *n - 1;    xnorm = dznrm2_(&i__1, &x[1], incx);    alphr = alpha->r;    alphi = d_imag(alpha);    if (xnorm == 0.)    {        /* H = [1-alpha/abs(alpha) 0;        0 I], sign chosen so ALPHA >= 0. */        if (alphi == 0.)        {            if (alphr >= 0.)            {                /* When TAU.eq.ZERO, the vector is special-cased to be */                /* all zeros in the application routines. We do not need */                /* to clear it. */                tau->r = 0., tau->i = 0.;            }            else            {                /* However, the application routines rely on explicit */                /* zero checks when TAU.ne.ZERO, and we must clear X. */                tau->r = 2., tau->i = 0.;                i__1 = *n - 1;                for (j = 1;                        j <= i__1;                        ++j)                {                    i__2 = (j - 1) * *incx + 1;                    x[i__2].r = 0.;                    x[i__2].i = 0.; // , expr subst                }                z__1.r = -alpha->r;                z__1.i = -alpha->i; // , expr subst                alpha->r = z__1.r, alpha->i = z__1.i;            }        }        else        {            /* Only "reflecting" the diagonal entry to be real and non-negative. */            xnorm = dlapy2_(&alphr, &alphi);            d__1 = 1. - alphr / xnorm;            d__2 = -alphi / xnorm;            z__1.r = d__1;            z__1.i = d__2; // , expr subst            tau->r = z__1.r, tau->i = z__1.i;            i__1 = *n - 1;            for (j = 1;                    j <= i__1;//.........这里部分代码省略.........

/* Subroutine */ int zlarfg_(integer *n, doublecomplex *alpha, doublecomplex *	x, integer *incx, doublecomplex *tau){/*  -- 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       =======       ZLARFG generates a complex elementary reflector H of order n, such       that             H' * ( alpha ) = ( beta ),   H' * H = I.                  (   x   )   (   0  )       where alpha and beta are scalars, with beta real, and x is an       (n-1)-element complex vector. H is represented in the form             H = I - tau * ( 1 ) * ( 1 v' ) ,                           ( v )       where tau is a complex scalar and v is a complex (n-1)-element       vector. Note that H is not hermitian.       If the elements of x are all zero and alpha is real, then tau = 0       and H is taken to be the unit matrix.       Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .       Arguments       =========       N       (input) INTEGER               The order of the elementary reflector.       ALPHA   (input/output) COMPLEX*16               On entry, the value alpha.               On exit, it is overwritten with the value beta.       X       (input/output) COMPLEX*16 array, dimension                              (1+(N-2)*abs(INCX))               On entry, the vector x.               On exit, it is overwritten with the vector v.       INCX    (input) INTEGER               The increment between elements of X. INCX > 0.       TAU     (output) COMPLEX*16               The value tau.       =====================================================================          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b5 = {1.,0.};        /* System generated locals */    integer i__1;    doublereal d__1, d__2;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *), d_sign(doublereal *, doublereal *);    /* Local variables */    static doublereal beta;    static integer j;    static doublereal alphi, alphr;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static doublereal xnorm;    extern doublereal dlapy3_(doublereal *, doublereal *, doublereal *), 	    dznrm2_(integer *, doublecomplex *, integer *), dlamch_(char *);    static doublereal safmin;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *);    static doublereal rsafmn;    extern /* Double Complex */ VOID zladiv_(doublecomplex *, doublecomplex *,	     doublecomplex *);    static integer knt;    --x;    /* Function Body */    if (*n <= 0) {	tau->r = 0., tau->i = 0.;	return 0;    }    i__1 = *n - 1;    xnorm = dznrm2_(&i__1, &x[1], incx);    alphr = alpha->r;    alphi = d_imag(alpha);    if (xnorm == 0. && alphi == 0.) {/*        H  =  I *///.........这里部分代码省略.........

 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 *///.........这里部分代码省略.........

/* Subroutine */ int zhgeqz_(char *job, char *compq, char *compz, integer *n, 	integer *ilo, integer *ihi, doublecomplex *a, integer *lda, 	doublecomplex *b, integer *ldb, doublecomplex *alpha, doublecomplex *	beta, doublecomplex *q, integer *ldq, doublecomplex *z__, integer *	ldz, doublecomplex *work, integer *lwork, doublereal *rwork, integer *	info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 	    z_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1, d__2, d__3, d__4, d__5, d__6;    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    double d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), pow_zi(	    doublecomplex *, doublecomplex *, integer *), z_sqrt(	    doublecomplex *, doublecomplex *);    /* Local variables */    static doublereal absb, atol, btol, temp, opst;    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *);    static doublereal temp2, c__;    static integer j;    static doublecomplex s, t;    extern logical lsame_(char *, char *);    static doublecomplex ctemp;    static integer iiter, ilast, jiter;    static doublereal anorm;    static integer maxit;    static doublereal bnorm;    static doublecomplex shift;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static doublereal tempr;    static doublecomplex ctemp2, ctemp3;    static logical ilazr2;    static integer jc, in;    static doublereal ascale, bscale;    static doublecomplex u12;    extern doublereal dlamch_(char *);    static integer jr, nq;    static doublecomplex signbc;    static integer nz;    static doublereal safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    static doublecomplex eshift;    static logical ilschr;    static integer icompq, ilastm;    static doublecomplex rtdisc;    static integer ischur;    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 	    doublereal *);    static logical ilazro;    static integer icompz, ifirst;    extern /* Subroutine */ int zlartg_(doublecomplex *, doublecomplex *, 	    doublereal *, doublecomplex *, doublecomplex *);    static integer ifrstm;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    static integer istart;    static logical lquery;    static doublecomplex ad11, ad12, ad21, ad22;    static integer jch;    static logical ilq, ilz;    static doublereal ulp;    static doublecomplex abi22;#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)]#define q_subscr(a_1,a_2) (a_2)*q_dim1 + a_1#define q_ref(a_1,a_2) q[q_subscr(a_1,a_2)]#define z___subscr(a_1,a_2) (a_2)*z_dim1 + a_1#define z___ref(a_1,a_2) z__[z___subscr(a_1,a_2)]/*  -- LAPACK routine (instrumented to count operations, version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          June 30, 1999          ----------------------- Begin Timing Code ------------------------          Common block to return operation count and iteration count          ITCNT is initialized to 0, OPS is only incremented          OPST is used to accumulate small contributions to OPS          to avoid roundoff error          ------------------------ End Timing Code -------------------------       Purpose       =======       ZHGEQZ implements a single-shift version of the QZ   //.........这里部分代码省略.........

/* Subroutine */ int zlatdf_(integer *ijob, integer *n, doublecomplex *z__, 	integer *ldz, doublecomplex *rhs, doublereal *rdsum, doublereal *	rdscal, integer *ipiv, integer *jpiv){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;    doublecomplex z__1, z__2, z__3;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    double z_abs(doublecomplex *);    void z_sqrt(doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, k;    doublecomplex bm, bp, xm[2], xp[2];    integer info;    doublecomplex temp, work[8];    doublereal scale;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    doublecomplex pmone;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    doublereal rtemp, sminu, rwork[2];    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    doublereal splus;    extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *), zgesc2_(	    integer *, doublecomplex *, integer *, doublecomplex *, integer *, 	     integer *, doublereal *), zgecon_(char *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublecomplex *, doublereal *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *, 	     doublereal *, doublereal *), zlaswp_(integer *, doublecomplex *, 	    integer *, integer *, integer *, integer *, integer *);/*  -- LAPACK auxiliary routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLATDF computes the contribution to the reciprocal Dif-estimate *//*  by solving for x in Z * x = b, where b is chosen such that the norm *//*  of x is as large as possible. It is assumed that LU decomposition *//*  of Z has been computed by ZGETC2. On entry RHS = f holds the *//*  contribution from earlier solved sub-systems, and on return RHS = x. *//*  The factorization of Z returned by ZGETC2 has the form *//*  Z = P * L * U * Q, where P and Q are permutation matrices. L is lower *//*  triangular with unit diagonal elements and U is upper triangular. *//*  Arguments *//*  ========= *//*  IJOB    (input) INTEGER *//*          IJOB = 2: First compute an approximative null-vector e *//*              of Z using ZGECON, e is normalized and solve for *//*              Zx = +-e - f with the sign giving the greater value of *//*              2-norm(x).  About 5 times as expensive as Default. *//*          IJOB .ne. 2: Local look ahead strategy where *//*              all entries of the r.h.s. b is choosen as either +1 or *//*              -1.  Default. *//*  N       (input) INTEGER *//*          The number of columns of the matrix Z. *//*  Z       (input) DOUBLE PRECISION array, dimension (LDZ, N) *//*          On entry, the LU part of the factorization of the n-by-n *//*          matrix Z computed by ZGETC2:  Z = P * L * U * Q *//*  LDZ     (input) INTEGER *//*          The leading dimension of the array Z.  LDA >= max(1, N). *//*  RHS     (input/output) DOUBLE PRECISION array, dimension (N). *//*          On entry, RHS contains contributions from other subsystems. *//*          On exit, RHS contains the solution of the subsystem with *//*          entries according to the value of IJOB (see above). *//*  RDSUM   (input/output) DOUBLE PRECISION *//*          On entry, the sum of squares of computed contributions to *//*          the Dif-estimate under computation by ZTGSYL, where the *//*          scaling factor RDSCAL (see below) has been factored out. *//*          On exit, the corresponding sum of squares updated with the *//*          contributions from the current sub-system. *//*          If TRANS = 'T' RDSUM is not touched. *//*          NOTE: RDSUM only makes sense when ZTGSY2 is called by CTGSYL. *//*  RDSCAL  (input/output) DOUBLE PRECISION *//*          On entry, scaling factor used to prevent overflow in RDSUM. *///.........这里部分代码省略.........

 int ztgsy2_(char *trans, int *ijob, int *m, int *	n, doublecomplex *a, int *lda, doublecomplex *b, int *ldb, 	doublecomplex *c__, int *ldc, doublecomplex *d__, int *ldd, 	doublecomplex *e, int *lde, doublecomplex *f, int *ldf, 	double *scale, double *rdsum, double *rdscal, int *	info){    /* System generated locals */    int a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, d_dim1, 	    d_offset, e_dim1, e_offset, f_dim1, f_offset, i__1, i__2, i__3, 	    i__4;    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    int i__, j, k;    doublecomplex z__[4]	/* was [2][2] */, rhs[2];    int ierr, ipiv[2], jpiv[2];    doublecomplex alpha;    extern int lsame_(char *, char *);    extern  int zscal_(int *, doublecomplex *, 	    doublecomplex *, int *), zaxpy_(int *, doublecomplex *, 	    doublecomplex *, int *, doublecomplex *, int *), zgesc2_(	    int *, doublecomplex *, int *, doublecomplex *, int *, 	     int *, double *), zgetc2_(int *, doublecomplex *, 	    int *, int *, int *, int *);    double scaloc;    extern  int xerbla_(char *, int *), zlatdf_(	    int *, int *, doublecomplex *, int *, doublecomplex *, 	     double *, double *, int *, int *);    int notran;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZTGSY2 solves the generalized Sylvester equation *//*              A * R - L * B = scale *   C               (1) *//*              D * R - L * E = scale * F *//*  using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices, *//*  (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, *//*  N-by-N and M-by-N, respectively. A, B, D and E are upper triangular *//*  (i.e., (A,D) and (B,E) in generalized Schur form). *//*  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output *//*  scaling factor chosen to avoid overflow. *//*  In matrix notation solving equation (1) corresponds to solve *//*  Zx = scale * b, where Z is defined as *//*         Z = [ kron(In, A)  -kron(B', Im) ]             (2) *//*             [ kron(In, D)  -kron(E', Im) ], *//*  Ik is the identity matrix of size k and X' is the transpose of X. *//*  kron(X, Y) is the Kronecker product between the matrices X and Y. *//*  If TRANS = 'C', y in the conjugate transposed system Z'y = scale*b *//*  is solved for, which is equivalent to solve for R and L in *//*              A' * R  + D' * L   = scale *  C           (3) *//*              R  * B' + L  * E'  = scale * -F *//*  This case is used to compute an estimate of Dif[(A, D), (B, E)] = *//*  = sigma_MIN(Z) using reverse communicaton with ZLACON. *//*  ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL *//*  of an upper bound on the separation between to matrix pairs. Then *//*  the input (A, D), (B, E) are sub-pencils of two matrix pairs in *//*  ZTGSYL. *//*  Arguments *//*  ========= *//*  TRANS   (input) CHARACTER*1 *//*          = 'N', solve the generalized Sylvester equation (1). *//*          = 'T': solve the 'transposed' system (3). *//*  IJOB    (input) INTEGER *//*          Specifies what kind of functionality to be performed. *//*          =0: solve (1) only. *//*          =1: A contribution from this subsystem to a Frobenius *//*              norm-based estimate of the separation between two matrix *//*              pairs is computed. (look ahead strategy is used). *//*          =2: A contribution from this subsystem to a Frobenius *//*              norm-based estimate of the separation between two matrix *//*              pairs is computed. (DGECON on sub-systems is used.) *//*          Not referenced if TRANS = 'T'. *///.........这里部分代码省略.........

//.........这里部分代码省略.........        },{7.,8.},{7.,8.},{7.,8.},{7.,8.},{7.,8.},{.3,.1},{.1,.4},{            .4,            .1        },{.1,.2},{2.,3.},{2.,3.},{2.,3.},{2.,3.},{.1,.1},{4.,5.},{            4.,            5.        },{4.,5.},{4.,5.},{4.,5.},{4.,5.},{4.,5.},{.3,-.4},{6.,7.},{            6.,            7.        },{6.,7.},{6.,7.},{6.,7.},{6.,7.},{6.,7.},{.1,-.3},{8.,9.},{            .5,            -.1        },{2.,5.},{2.,5.},{2.,5.},{2.,5.},{2.,5.},{.1,.1},{3.,6.},{            -.6,            .1        },{4.,7.},{.1,-.3},{7.,2.},{7.,2.},{7.,2.},{.3,.1},{5.,8.},{            .1,            .4        },{6.,9.},{.4,.1},{8.,3.},{.1,.2},{9.,4.}    };    /* System generated locals */    integer i__1, i__2, i__3;    doublereal d__1;    doublecomplex z__1;    /* Builtin functions */    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),            e_wsle(void);    /* Subroutine */ int s_stop(char *, ftnlen);    /* Local variables */    static integer i__;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *,                                       doublecomplex *, integer *), ctest_(integer *, doublecomplex *,                                               doublecomplex *, doublecomplex *, doublereal *);    static doublecomplex mwpcs[5], mwpct[5];    extern /* Subroutine */ int itest1_(integer *, integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int stest1_(doublereal *, doublereal *,                                        doublereal *, doublereal *);    static doublecomplex cx[8];    extern /* Subroutine */ int zdscal_(integer *, doublereal *,                                        doublecomplex *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    static integer np1, len;    /* Fortran I/O blocks */    static cilist io___19 = { 0, 6, 0, 0, 0 };#define ctrue5_subscr(a_1,a_2,a_3) ((a_3)*5 + (a_2))*8 + a_1 - 49#define ctrue5_ref(a_1,a_2,a_3) ctrue5[ctrue5_subscr(a_1,a_2,a_3)]#define ctrue6_subscr(a_1,a_2,a_3) ((a_3)*5 + (a_2))*8 + a_1 - 49#define ctrue6_ref(a_1,a_2,a_3) ctrue6[ctrue6_subscr(a_1,a_2,a_3)]#define cv_subscr(a_1,a_2,a_3) ((a_3)*5 + (a_2))*8 + a_1 - 49#define cv_ref(a_1,a_2,a_3) cv[cv_subscr(a_1,a_2,a_3)]    for (combla_1.incx = 1; combla_1.incx <= 2; ++combla_1.incx) {        for (np1 = 1; np1 <= 5; ++np1) {            combla_1.n = np1 - 1;            len = max(combla_1.n,1) << 1;            i__1 = len;            for (i__ = 1; i__ <= i__1; ++i__) {

/* Subroutine */ int zlahr2_(integer *n, integer *k, integer *nb, 	doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *t, 	integer *ldt, doublecomplex *y, integer *ldy){    /* System generated locals */    integer a_dim1, a_offset, t_dim1, t_offset, y_dim1, y_offset, i__1, i__2, 	    i__3;    doublecomplex z__1;    /* Local variables */    integer i__;    doublecomplex ei;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), 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 *), 	    zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *), ztrmm_(char *, char *, char *, char *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *), 	    zaxpy_(integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *), ztrmv_(char *, char *, char *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *), zlarfg_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *), zlacgv_(integer *, 	    doublecomplex *, integer *), zlacpy_(char *, integer *, integer *, 	     doublecomplex *, 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 *//*  ======= *//*  ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) *//*  matrix A so that elements below the k-th subdiagonal are zero. The *//*  reduction is performed by an unitary similarity transformation *//*  Q' * A * Q. The routine returns the matrices V and T which determine *//*  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. *//*  This is an auxiliary routine called by ZGEHRD. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The order of the matrix A. *//*  K       (input) INTEGER *//*          The offset for the reduction. Elements below the k-th *//*          subdiagonal in the first NB columns are reduced to zero. *//*          K < N. *//*  NB      (input) INTEGER *//*          The number of columns to be reduced. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N-K+1) *//*          On entry, the n-by-(n-k+1) general matrix A. *//*          On exit, the elements on and above the k-th subdiagonal in *//*          the first NB columns are overwritten with the corresponding *//*          elements of the reduced matrix; the elements below the k-th *//*          subdiagonal, with the array TAU, represent the matrix Q as a *//*          product of elementary reflectors. The other columns of A are *//*          unchanged. See Further Details. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= max(1,N). *//*  TAU     (output) COMPLEX*16 array, dimension (NB) *//*          The scalar factors of the elementary reflectors. See Further *//*          Details. *//*  T       (output) COMPLEX*16 array, dimension (LDT,NB) *//*          The upper triangular matrix T. *//*  LDT     (input) INTEGER *//*          The leading dimension of the array T.  LDT >= NB. *//*  Y       (output) COMPLEX*16 array, dimension (LDY,NB) *//*          The n-by-nb matrix Y. *//*  LDY     (input) INTEGER *//*          The leading dimension of the array Y. LDY >= N. *//*  Further Details *//*  =============== *//*  The matrix Q is represented as a product of nb elementary reflectors *//*     Q = H(1) H(2) . . . H(nb). *///.........这里部分代码省略.........

/* Subroutine */ int zlavsy_(char *uplo, char *trans, char *diag, integer *n, 	integer *nrhs, doublecomplex *a, integer *lda, integer *ipiv, 	doublecomplex *b, integer *ldb, integer *info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;    doublecomplex z__1, z__2, z__3;    /* Local variables */    static integer j, k;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zgeru_(integer *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *)	    , zswap_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    static doublecomplex t1, t2, d11, d12, d21, d22;    static integer kp;    extern /* Subroutine */ int xerbla_(char *, integer *);    static logical nounit;#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)]/*  -- 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       =======          ZLAVSY  performs one of the matrix-vector operations             x := A*x  or  x := A'*x,          where x is an N element vector and  A is one of the factors          from the symmetric factorization computed by ZSYTRF.          ZSYTRF produces a factorization of the form               U * D * U'      or     L * D * L' ,          where U (or L) is a product of permutation and unit upper (lower)          triangular matrices, U' (or L') is the transpose of          U (or L), and D is symmetric and block diagonal with 1 x 1 and          2 x 2 diagonal blocks.  The multipliers for the transformations          and the upper or lower triangular parts of the diagonal blocks          are stored in the leading upper or lower triangle of the 2-D          array A.          If TRANS = 'N' or 'n', ZLAVSY multiplies either by U or U * D          (or L or L * D).          If TRANS = 'T' or 't', ZLAVSY multiplies either by U' or D * U'          (or L' or D * L' ).       Arguments       ==========       UPLO   - CHARACTER*1                On entry, UPLO specifies whether the triangular matrix                stored in A is upper or lower triangular.                   UPLO = 'U' or 'u'   The matrix is upper triangular.                   UPLO = 'L' or 'l'   The matrix is lower triangular.                Unchanged on exit.       TRANS  - CHARACTER*1                On entry, TRANS specifies the operation to be performed as                follows:                   TRANS = 'N' or 'n'   x := A*x.                   TRANS = 'T' or 't'   x := A'*x.                Unchanged on exit.       DIAG   - CHARACTER*1                On entry, DIAG specifies whether the diagonal blocks are                assumed to be unit matrices:                   DIAG = 'U' or 'u'   Diagonal blocks are unit matrices.                   DIAG = 'N' or 'n'   Diagonal blocks are non-unit.                Unchanged on exit.       N      - INTEGER                On entry, N specifies the order of the matrix A.                N must be at least zero.                Unchanged on exit.       NRHS   - INTEGER                On entry, NRHS specifies the number of right hand sides,                i.e., the number of vectors x to be multiplied by A.                NRHS must be at least zero.                Unchanged on exit.       A      - COMPLEX*16 array, dimension( LDA, N )                On entry, A contains a block diagonal matrix and the                multipliers of the transformations used to obtain it,                stored as a 2-D triangular matrix.                Unchanged on exit.   //.........这里部分代码省略.........

/*! comple*_zgbmatrix operator */inline _zgbmatrix operator*(const comple& d, const _zgbmatrix& mat){VERBOSE_REPORT;  zscal_((mat.kl+mat.ku+1)*mat.n, d, mat.array, 1);  return mat;}

//.........这里部分代码省略.........	if (*ku > 0 && *m < *n) {/*           Annihilate a(m,m+1) by applying plane rotations from the *//*           right */	    i__1 = *ku + (*m + 1) * ab_dim1;	    rb.r = ab[i__1].r, rb.i = ab[i__1].i;	    for (i__ = *m; i__ >= 1; --i__) {		zlartg_(&ab[*ku + 1 + i__ * ab_dim1], &rb, &rc, &rs, &ra);		i__1 = *ku + 1 + i__ * ab_dim1;		ab[i__1].r = ra.r, ab[i__1].i = ra.i;		if (i__ > 1) {		    d_cnjg(&z__3, &rs);		    z__2.r = -z__3.r, z__2.i = -z__3.i;		    i__1 = *ku + i__ * ab_dim1;		    z__1.r = z__2.r * ab[i__1].r - z__2.i * ab[i__1].i, 			    z__1.i = z__2.r * ab[i__1].i + z__2.i * ab[i__1]			    .r;		    rb.r = z__1.r, rb.i = z__1.i;		    i__1 = *ku + i__ * ab_dim1;		    i__2 = *ku + i__ * ab_dim1;		    z__1.r = rc * ab[i__2].r, z__1.i = rc * ab[i__2].i;		    ab[i__1].r = z__1.r, ab[i__1].i = z__1.i;		}		if (wantpt) {		    d_cnjg(&z__1, &rs);		    zrot_(n, &pt[i__ + pt_dim1], ldpt, &pt[*m + 1 + pt_dim1], 			    ldpt, &rc, &z__1);		}	    }	}    }/*     Make diagonal and superdiagonal elements real, storing them in D *//*     and E */    i__1 = *ku + 1 + ab_dim1;    t.r = ab[i__1].r, t.i = ab[i__1].i;    i__1 = minmn;    for (i__ = 1; i__ <= i__1; ++i__) {	abst = z_abs(&t);	d__[i__] = abst;	if (abst != 0.) {	    z__1.r = t.r / abst, z__1.i = t.i / abst;	    t.r = z__1.r, t.i = z__1.i;	} else {	    t.r = 1., t.i = 0.;	}	if (wantq) {	    zscal_(m, &t, &q[i__ * q_dim1 + 1], &c__1);	}	if (wantc) {	    d_cnjg(&z__1, &t);	    zscal_(ncc, &z__1, &c__[i__ + c_dim1], ldc);	}	if (i__ < minmn) {	    if (*ku == 0 && *kl == 0) {		e[i__] = 0.;		i__2 = (i__ + 1) * ab_dim1 + 1;		t.r = ab[i__2].r, t.i = ab[i__2].i;	    } else {		if (*ku == 0) {		    i__2 = i__ * ab_dim1 + 2;		    d_cnjg(&z__2, &t);		    z__1.r = ab[i__2].r * z__2.r - ab[i__2].i * z__2.i, 			    z__1.i = ab[i__2].r * z__2.i + ab[i__2].i * 			    z__2.r;		    t.r = z__1.r, t.i = z__1.i;		} else {		    i__2 = *ku + (i__ + 1) * ab_dim1;		    d_cnjg(&z__2, &t);		    z__1.r = ab[i__2].r * z__2.r - ab[i__2].i * z__2.i, 			    z__1.i = ab[i__2].r * z__2.i + ab[i__2].i * 			    z__2.r;		    t.r = z__1.r, t.i = z__1.i;		}		abst = z_abs(&t);		e[i__] = abst;		if (abst != 0.) {		    z__1.r = t.r / abst, z__1.i = t.i / abst;		    t.r = z__1.r, t.i = z__1.i;		} else {		    t.r = 1., t.i = 0.;		}		if (wantpt) {		    zscal_(n, &t, &pt[i__ + 1 + pt_dim1], ldpt);		}		i__2 = *ku + 1 + (i__ + 1) * ab_dim1;		d_cnjg(&z__2, &t);		z__1.r = ab[i__2].r * z__2.r - ab[i__2].i * z__2.i, z__1.i = 			ab[i__2].r * z__2.i + ab[i__2].i * z__2.r;		t.r = z__1.r, t.i = z__1.i;	    }	}    }    return 0;/*     End of ZGBBRD */} /* zgbbrd_ */

//.........这里部分代码省略.........       ( a   h   a   a   a )          ( a   h   a   a   a )          ( a   h   a   a   a )          ( h   h   a   a   a )          ( v1  h   a   a   a )          ( v1  v2  a   a   a )          ( v1  v2  a   a   a )       where a denotes an element of the original matrix A, h denotes a       modified element of the upper Hessenberg matrix H, and vi denotes an       element of the vector defining H(i).       =====================================================================          Quick return if possible          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 a_dim1, a_offset, t_dim1, t_offset, y_dim1, y_offset, i__1, i__2, 	    i__3;    doublecomplex z__1;    /* Local variables */    static integer i;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zaxpy_(integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *), ztrmv_(char *, char *, 	    char *, integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    static doublecomplex ei;    extern /* Subroutine */ int zlarfg_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *), zlacgv_(integer *, 	    doublecomplex *, integer *);#define TAU(I) tau[(I)-1]#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]#define T(I,J) t[(I)-1 + ((J)-1)* ( *ldt)]#define Y(I,J) y[(I)-1 + ((J)-1)* ( *ldy)]    if (*n <= 1) {	return 0;    }    i__1 = *nb;    for (i = 1; i <= *nb; ++i) {	if (i > 1) {/*           Update A(1:n,i)                Compute i-th column of A - Y * V' */

 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 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;//.........这里部分代码省略.........

/* Subroutine */ int zung2r_(integer *m, integer *n, integer *k, 	doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *	work, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublecomplex z__1;    /* Local variables */    integer i__, j, l;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), zlarf_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *), xerbla_(char *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZUNG2R generates an m by n complex matrix Q with orthonormal columns, *//*  which is defined as the first n columns of a product of k elementary *//*  reflectors of order m *//*        Q  =  H(1) H(2) . . . H(k) *//*  as returned by ZGEQRF. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix Q. M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix Q. M >= N >= 0. *//*  K       (input) INTEGER *//*          The number of elementary reflectors whose product defines the *//*          matrix Q. N >= K >= 0. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the i-th column must contain the vector which *//*          defines the elementary reflector H(i), for i = 1,2,...,k, as *//*          returned by ZGEQRF in the first k columns of its array *//*          argument A. *//*          On exit, the m by n matrix Q. *//*  LDA     (input) INTEGER *//*          The first dimension of the array A. LDA >= max(1,M). *//*  TAU     (input) COMPLEX*16 array, dimension (K) *//*          TAU(i) must contain the scalar factor of the elementary *//*          reflector H(i), as returned by ZGEQRF. *//*  WORK    (workspace) COMPLEX*16 array, dimension (N) *//*  INFO    (output) INTEGER *//*          = 0: successful exit *//*          < 0: if INFO = -i, the i-th argument has an illegal value *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input arguments */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --tau;    --work;    /* Function Body */    *info = 0;    if (*m < 0) {	*info = -1;    } else if (*n < 0 || *n > *m) {	*info = -2;    } else if (*k < 0 || *k > *n) {	*info = -3;    } else if (*lda < max(1,*m)) {	*info = -5;//.........这里部分代码省略.........

/* Subroutine */ int ztptri_(char *uplo, char *diag, integer *n, 	doublecomplex *ap, integer *info){    /* System generated locals */    integer i__1, i__2;    doublecomplex z__1;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer j, jc, jj;    doublecomplex ajj;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    logical upper;    extern /* Subroutine */ int ztpmv_(char *, char *, char *, integer *, 	    doublecomplex *, doublecomplex *, integer *), xerbla_(char *, integer *);    integer jclast;    logical nounit;/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZTPTRI computes the inverse of a complex upper or lower triangular *//*  matrix A stored in packed format. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          = 'U':  A is upper triangular; *//*          = 'L':  A is lower triangular. *//*  DIAG    (input) CHARACTER*1 *//*          = 'N':  A is non-unit triangular; *//*          = 'U':  A is unit triangular. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) *//*          On entry, the upper or lower triangular matrix A, stored *//*          columnwise in a linear array.  The j-th column of A is stored *//*          in the array AP as follows: *//*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; *//*          if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n. *//*          See below for further details. *//*          On exit, the (triangular) inverse of the original matrix, in *//*          the same packed storage format. *//*  INFO    (output) INTEGER *//*          = 0:  successful exit *//*          < 0:  if INFO = -i, the i-th argument had an illegal value *//*          > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular *//*                matrix is singular and its inverse can not be computed. *//*  Further Details *//*  =============== *//*  A triangular matrix A can be transferred to packed storage using one *//*  of the following program segments: *//*  UPLO = 'U':                      UPLO = 'L': *//*        JC = 1                           JC = 1 *//*        DO 2 J = 1, N                    DO 2 J = 1, N *//*           DO 1 I = 1, J                    DO 1 I = J, N *//*              AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J) *//*      1    CONTINUE                    1    CONTINUE *//*           JC = JC + J                      JC = JC + N - J + 1 *//*      2 CONTINUE                       2 CONTINUE *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input parameters. */    /* Parameter adjustments */    --ap;//.........这里部分代码省略.........

//.........这里部分代码省略........./*     Compute the product L*L', overwriting L. */    } else {	if (*rank < *n) {	    i__1 = *n;	    for (j = *rank + 1; j <= i__1; ++j) {		i__2 = *n;		for (i__ = j; i__ <= i__2; ++i__) {		    i__3 = i__ + j * afac_dim1;		    afac[i__3].r = 0., afac[i__3].i = 0.;/* L140: */		}/* L150: */	    }	}	for (k = *n; k >= 1; --k) {/*           Add a multiple of column K of the factor L to each of *//*           columns K+1 through N. */	    if (k + 1 <= *n) {		i__1 = *n - k;		zher_("Lower", &i__1, &c_b20, &afac[k + 1 + k * afac_dim1], &			c__1, &afac[k + 1 + (k + 1) * afac_dim1], ldafac);	    }/*           Scale column K by the diagonal element. */	    i__1 = k + k * afac_dim1;	    tc.r = afac[i__1].r, tc.i = afac[i__1].i;	    i__1 = *n - k + 1;	    zscal_(&i__1, &tc, &afac[k + k * afac_dim1], &c__1);/* L160: */	}    }/*        Form P*L*L'*P' or P*U'*U*P' */    if (lsame_(uplo, "U")) {	i__1 = *n;	for (j = 1; j <= i__1; ++j) {	    i__2 = *n;	    for (i__ = 1; i__ <= i__2; ++i__) {		if (piv[i__] <= piv[j]) {		    if (i__ <= j) {			i__3 = piv[i__] + piv[j] * perm_dim1;			i__4 = i__ + j * afac_dim1;			perm[i__3].r = afac[i__4].r, perm[i__3].i = afac[i__4]				.i;		    } else {			i__3 = piv[i__] + piv[j] * perm_dim1;			d_cnjg(&z__1, &afac[j + i__ * afac_dim1]);			perm[i__3].r = z__1.r, perm[i__3].i = z__1.i;		    }		}/* L170: */	    }/* L180: */	}    } else {

void zscal( int n, doublecomplex* alpha, doublecomplex *x, int incx){    zscal_(&n, alpha, x, &incx);}

/*! _zrovector*comple operator */inline _zrovector operator*(const _zrovector& vec, const comple& d){VERBOSE_REPORT;  zscal_(vec.l, d, vec.array, 1);    return vec;}

/* Subroutine */ int zsptrf_(char *uplo, integer *n, doublecomplex *ap, 	integer *ipiv, integer *info){    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5, i__6;    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 z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, k;    doublecomplex t, r1, d11, d12, d21, d22;    integer kc, kk, kp;    doublecomplex wk;    integer kx, knc, kpc, npp;    doublecomplex wkm1, wkp1;    integer imax, jmax;    extern /* Subroutine */ int zspr_(char *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *);    doublereal alpha;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    integer kstep;    logical upper;    extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    doublereal absakk;    extern /* Subroutine */ int xerbla_(char *, integer *);    doublereal colmax;    extern integer izamax_(integer *, doublecomplex *, integer *);    doublereal rowmax;/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZSPTRF computes the factorization of a complex symmetric matrix A *//*  stored in packed format using the Bunch-Kaufman diagonal pivoting *//*  method: *//*     A = U*D*U**T  or  A = L*D*L**T *//*  where U (or L) is a product of permutation and unit upper (lower) *//*  triangular matrices, and D is symmetric and block diagonal with *//*  1-by-1 and 2-by-2 diagonal blocks. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          = 'U':  Upper triangle of A is stored; *//*          = 'L':  Lower triangle of A is stored. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) *//*          On entry, the upper or lower triangle of the symmetric matrix *//*          A, packed columnwise in a linear array.  The j-th column of A *//*          is stored in the array AP as follows: *//*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; *//*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. *//*          On exit, the block diagonal matrix D and the multipliers used *//*          to obtain the factor U or L, stored as a packed triangular *//*          matrix overwriting A (see below for further details). *//*  IPIV    (output) INTEGER array, dimension (N) *//*          Details of the interchanges and the block structure of D. *//*          If IPIV(k) > 0, then rows and columns k and IPIV(k) were *//*          interchanged and D(k,k) is a 1-by-1 diagonal block. *//*          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and *//*          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) *//*          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) = *//*          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were *//*          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. *//*  INFO    (output) INTEGER *//*          = 0: successful exit *//*          < 0: if INFO = -i, the i-th argument had an illegal value *//*          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization *//*               has been completed, but the block diagonal matrix D is *//*               exactly singular, and division by zero will occur if it *//*               is used to solve a system of equations. *//*  Further Details *//*  =============== *///.........这里部分代码省略.........

/* Subroutine */ int zgetf2_(integer *m, integer *n, doublecomplex *a, 	integer *lda, integer *ipiv, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublecomplex z__1;    /* Builtin functions */    double z_abs(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j, jp;    doublereal sfmin;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), zgeru_(integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), zswap_(integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGETF2 computes an LU factorization of a general m-by-n matrix A *//*  using partial pivoting with row interchanges. *//*  The factorization has the form *//*     A = P * L * U *//*  where P is a permutation matrix, L is lower triangular with unit *//*  diagonal elements (lower trapezoidal if m > n), and U is upper *//*  triangular (upper trapezoidal if m < n). *//*  This is the right-looking Level 2 BLAS version of the algorithm. *//*  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. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the m by n matrix to be factored. *//*          On exit, the factors L and U from the factorization *//*          A = P*L*U; the unit diagonal elements of L are not stored. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= max(1,M). *//*  IPIV    (output) INTEGER array, dimension (min(M,N)) *//*          The pivot indices; for 1 <= i <= min(M,N), row i of the *//*          matrix was interchanged with row IPIV(i). *//*  INFO    (output) INTEGER *//*          = 0: successful exit *//*          < 0: if INFO = -k, the k-th argument had an illegal value *//*          > 0: if INFO = k, U(k,k) is exactly zero. The factorization *//*               has been completed, but the factor U is exactly *//*               singular, and division by zero will occur if it is used *//*               to solve a system of equations. *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input parameters. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --ipiv;    /* Function Body */    *info = 0;//.........这里部分代码省略.........

/* Subroutine */ int zungr2_(integer *m, integer *n, integer *k, 	doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *	work, integer *info){/*  -- LAPACK 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       =======       ZUNGR2 generates an m by n complex matrix Q with orthonormal rows,       which is defined as the last m rows of a product of k elementary       reflectors of order n             Q  =  H(1)' H(2)' . . . H(k)'       as returned by ZGERQF.       Arguments       =========       M       (input) INTEGER               The number of rows of the matrix Q. M >= 0.       N       (input) INTEGER               The number of columns of the matrix Q. N >= M.       K       (input) INTEGER               The number of elementary reflectors whose product defines the               matrix Q. M >= K >= 0.       A       (input/output) COMPLEX*16 array, dimension (LDA,N)               On entry, the (m-k+i)-th row must contain the vector which               defines the elementary reflector H(i), for i = 1,2,...,k, as               returned by ZGERQF in the last k rows of its array argument               A.               On exit, the m-by-n matrix Q.       LDA     (input) INTEGER               The first dimension of the array A. LDA >= max(1,M).       TAU     (input) COMPLEX*16 array, dimension (K)               TAU(i) must contain the scalar factor of the elementary               reflector H(i), as returned by ZGERQF.       WORK    (workspace) COMPLEX*16 array, dimension (M)       INFO    (output) INTEGER               = 0: successful exit               < 0: if INFO = -i, the i-th argument has an illegal value       =====================================================================          Test the input arguments          Parameter adjustments */    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    doublecomplex z__1, z__2;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static integer i__, j, l;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), zlarf_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *);    static integer ii;    extern /* Subroutine */ int xerbla_(char *, integer *), zlacgv_(	    integer *, doublecomplex *, integer *);#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;    --tau;    --work;    /* Function Body */    *info = 0;    if (*m < 0) {	*info = -1;    } else if (*n < *m) {	*info = -2;    } else if (*k < 0 || *k > *m) {	*info = -3;    } else if (*lda < max(1,*m)) {	*info = -5;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("ZUNGR2", &i__1);	return 0;    }//.........这里部分代码省略.........

/* 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   //.........这里部分代码省略.........

/* Subroutine */ int ztptri_(char *uplo, char *diag, integer *n, 	doublecomplex *ap, integer *info){/*  -- LAPACK 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       =======       ZTPTRI computes the inverse of a complex upper or lower triangular       matrix A stored in packed format.       Arguments       =========       UPLO    (input) CHARACTER*1               = 'U':  A is upper triangular;               = 'L':  A is lower triangular.       DIAG    (input) CHARACTER*1               = 'N':  A is non-unit triangular;               = 'U':  A is unit triangular.       N       (input) INTEGER               The order of the matrix A.  N >= 0.       AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)               On entry, the upper or lower triangular matrix A, stored               columnwise in a linear array.  The j-th column of A is stored               in the array AP as follows:               if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;               if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.               See below for further details.               On exit, the (triangular) inverse of the original matrix, in               the same packed storage format.       INFO    (output) INTEGER               = 0:  successful exit               < 0:  if INFO = -i, the i-th argument had an illegal value               > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular                     matrix is singular and its inverse can not be computed.       Further Details       ===============       A triangular matrix A can be transferred to packed storage using one       of the following program segments:       UPLO = 'U':                      UPLO = 'L':             JC = 1                           JC = 1             DO 2 J = 1, N                    DO 2 J = 1, N                DO 1 I = 1, J                    DO 1 I = J, N                   AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J)           1    CONTINUE                    1    CONTINUE                JC = JC + J                      JC = JC + N - J + 1           2 CONTINUE                       2 CONTINUE       =====================================================================          Test the input parameters.          Parameter adjustments          Function Body */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer i__1, i__2;    doublecomplex z__1;    /* Builtin functions */    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);    /* Local variables */    static integer j;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static logical upper;    extern /* Subroutine */ int ztpmv_(char *, char *, char *, integer *, 	    doublecomplex *, doublecomplex *, integer *);    static integer jc, jj;    extern /* Subroutine */ int xerbla_(char *, integer *);    static integer jclast;    static logical nounit;    static doublecomplex ajj;//.........这里部分代码省略.........