int chgeqz_(char *job, char *compq, char *compz, int *n,
int *ilo, int *ihi, complex *h__, int *ldh, complex *t,
int *ldt, complex *alpha, complex *beta, complex *q, int *ldq,
complex *z__, int *ldz, complex *work, int *lwork, float *
rwork, int *info)
{
/* System generated locals */
int h_dim1, h_offset, q_dim1, q_offset, t_dim1, t_offset, z_dim1,
z_offset, i__1, i__2, i__3, i__4, i__5, i__6;
float r__1, r__2, r__3, r__4, r__5, r__6;
complex q__1, q__2, q__3, q__4, q__5, q__6;
/* Builtin functions */
double c_abs(complex *);
void r_cnjg(complex *, complex *);
double r_imag(complex *);
void c_div(complex *, complex *, complex *), pow_ci(complex *, complex *,
int *), c_sqrt(complex *, complex *);
/* Local variables */
float c__;
int j;
complex s, t1;
int jc, in;
complex u12;
int jr;
complex ad11, ad12, ad21, ad22;
int jch;
int ilq, ilz;
float ulp;
complex abi22;
float absb, atol, btol, temp;
extern int crot_(int *, complex *, int *,
complex *, int *, float *, complex *);
float temp2;
extern int cscal_(int *, complex *, complex *,
int *);
extern int lsame_(char *, char *);
complex ctemp;
int iiter, ilast, jiter;
float anorm, bnorm;
int maxit;
complex shift;
float tempr;
complex ctemp2, ctemp3;
int ilazr2;
float ascale, bscale;
complex signbc;
extern double slamch_(char *), clanhs_(char *, int *,
complex *, int *, float *);
extern int claset_(char *, int *, int *, complex
*, complex *, complex *, int *), clartg_(complex *,
complex *, float *, complex *, complex *);
float safmin;
extern int xerbla_(char *, int *);
complex eshift;
int ilschr;
int icompq, ilastm;
complex rtdisc;
int ischur;
int ilazro;
int icompz, ifirst, ifrstm, istart;
int lquery;
/* -- LAPACK routine (version 3.2) -- */
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CHGEQZ computes the eigenvalues of a complex matrix pair (H,T), */
/* where H is an upper Hessenberg matrix and T is upper triangular, */
/* using the single-shift QZ method. */
/* Matrix pairs of this type are produced by the reduction to */
/* generalized upper Hessenberg form of a complex matrix pair (A,B): */
/* A = Q1*H*Z1**H, B = Q1*T*Z1**H, */
/* as computed by CGGHRD. */
/* If JOB='S', then the Hessenberg-triangular pair (H,T) is */
/* also reduced to generalized Schur form, */
/* H = Q*S*Z**H, T = Q*P*Z**H, */
/* where Q and Z are unitary matrices and S and P are upper triangular. */
/* Optionally, the unitary matrix Q from the generalized Schur */
/* factorization may be postmultiplied into an input matrix Q1, and the */
/* unitary matrix Z may be postmultiplied into an input matrix Z1. */
/* If Q1 and Z1 are the unitary matrices from CGGHRD that reduced */
/* the matrix pair (A,B) to generalized Hessenberg form, then the output */
//.........这里部分代码省略.........
/* Subroutine */ int clabrd_(integer *m, integer *n, integer *nb, complex *a,
integer *lda, real *d__, real *e, complex *tauq, complex *taup,
complex *x, integer *ldx, complex *y, integer *ldy)
{
/* System generated locals */
integer a_dim1, a_offset, x_dim1, x_offset, y_dim1, y_offset, i__1, i__2,
i__3;
complex q__1;
/* Local variables */
integer i__;
complex alpha;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *), cgemv_(char *, integer *, integer *, complex *,
complex *, integer *, complex *, integer *, complex *, complex *,
integer *), clarfg_(integer *, complex *, complex *,
integer *, complex *), clacgv_(integer *, complex *, integer *);
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CLABRD reduces the first NB rows and columns of a complex general */
/* m by n matrix A to upper or lower real 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 CGEBRD */
/* 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 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) REAL array, dimension (NB) */
/* The diagonal elements of the first NB rows and columns of */
/* the reduced matrix. D(i) = A(i,i). */
/* E (output) REAL array, dimension (NB) */
/* The off-diagonal elements of the first NB rows and columns of */
/* the reduced matrix. */
/* TAUQ (output) COMPLEX array dimension (NB) */
/* The scalar factors of the elementary reflectors which */
/* represent the unitary matrix Q. See Further Details. */
/* TAUP (output) COMPLEX array, dimension (NB) */
/* The scalar factors of the elementary reflectors which */
/* represent the unitary matrix P. See Further Details. */
/* X (output) COMPLEX 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 array, dimension (LDY,NB) */
/* The n-by-nb matrix Y required to update the unreduced part */
/* of A. */
//.........这里部分代码省略.........
int cgbbrd_(char *vect, int *m, int *n, int *ncc,
int *kl, int *ku, complex *ab, int *ldab, float *d__,
float *e, complex *q, int *ldq, complex *pt, int *ldpt,
complex *c__, int *ldc, complex *work, float *rwork, int *info)
{
/* System generated locals */
int ab_dim1, ab_offset, c_dim1, c_offset, pt_dim1, pt_offset, q_dim1,
q_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7;
complex q__1, q__2, q__3;
/* Builtin functions */
void r_cnjg(complex *, complex *);
double c_abs(complex *);
/* Local variables */
int i__, j, l;
complex t;
int j1, j2, kb;
complex ra, rb;
float rc;
int kk, ml, nr, mu;
complex rs;
int kb1, ml0, mu0, klm, kun, nrt, klu1, inca;
float abst;
extern int crot_(int *, complex *, int *,
complex *, int *, float *, complex *), cscal_(int *,
complex *, complex *, int *);
extern int lsame_(char *, char *);
int wantb, wantc;
int minmn;
int wantq;
extern int claset_(char *, int *, int *, complex
*, complex *, complex *, int *), clartg_(complex *,
complex *, float *, complex *, complex *), xerbla_(char *, int
*), clargv_(int *, complex *, int *, complex *,
int *, float *, int *), clartv_(int *, complex *,
int *, complex *, int *, float *, complex *, int *);
int wantpt;
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CGBBRD reduces a complex general m-by-n band matrix A to float upper */
/* bidiagonal form B by a unitary transformation: Q' * A * P = B. */
/* The routine computes B, and optionally forms Q or P', or computes */
/* Q'*C for a given matrix C. */
/* Arguments */
/* ========= */
/* VECT (input) CHARACTER*1 */
/* Specifies whether or not the matrices Q and P' are to be */
/* formed. */
/* = 'N': do not form Q or P'; */
/* = 'Q': form Q only; */
/* = 'P': form P' only; */
/* = 'B': form both. */
/* 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. */
/* NCC (input) INTEGER */
/* The number of columns of the matrix C. NCC >= 0. */
/* KL (input) INTEGER */
/* The number of subdiagonals of the matrix A. KL >= 0. */
/* KU (input) INTEGER */
/* The number of superdiagonals of the matrix A. KU >= 0. */
/* AB (input/output) COMPLEX array, dimension (LDAB,N) */
/* On entry, the m-by-n band matrix A, stored in rows 1 to */
/* KL+KU+1. The j-th column of A is stored in the j-th column of */
/* the array AB as follows: */
/* AB(ku+1+i-j,j) = A(i,j) for MAX(1,j-ku)<=i<=MIN(m,j+kl). */
/* On exit, A is overwritten by values generated during the */
/* reduction. */
/* LDAB (input) INTEGER */
/* The leading dimension of the array A. LDAB >= KL+KU+1. */
/* D (output) REAL array, dimension (MIN(M,N)) */
/* The diagonal elements of the bidiagonal matrix B. */
/* E (output) REAL array, dimension (MIN(M,N)-1) */
/* The superdiagonal elements of the bidiagonal matrix B. */
//.........这里部分代码省略.........
/*< subroutine csvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) >*/
/* Subroutine */ int csvdc_(complex *x, integer *ldx, integer *n, integer *p,
complex *s, complex *e, complex *u, integer *ldu, complex *v, integer
*ldv, complex *work, integer *job, integer *info)
{
/* System generated locals */
integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
i__3, i__4;
real r__1, r__2, r__3, r__4;
complex q__1, q__2, q__3;
/* Builtin functions */
double r_imag(complex *), c_abs(complex *);
void c_div(complex *, complex *, complex *), r_cnjg(complex *, complex *);
double sqrt(doublereal);
/* Local variables */
real b, c__, f, g;
integer i__, j, k, l=0, m;
complex r__, t;
real t1, el;
integer kk;
real cs;
integer ll, mm, ls=0;
real sl;
integer lu;
real sm, sn;
integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
real emm1, smm1;
integer kase, jobu, iter;
real test;
integer nctp1, nrtp1;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *);
real scale;
extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer
*, complex *, integer *);
real shift;
extern /* Subroutine */ int cswap_(integer *, complex *, integer *,
complex *, integer *);
integer maxit;
extern /* Subroutine */ int caxpy_(integer *, complex *, complex *,
integer *, complex *, integer *), csrot_(integer *, complex *,
integer *, complex *, integer *, real *, real *);
logical wantu, wantv;
extern /* Subroutine */ int srotg_(real *, real *, real *, real *);
real ztest;
extern doublereal scnrm2_(integer *, complex *, integer *);
/*< integer ldx,n,p,ldu,ldv,job,info >*/
/*< complex x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1) >*/
/* csvdc is a subroutine to reduce a complex nxp matrix x by */
/* unitary transformations u and v to diagonal form. the */
/* diagonal elements s(i) are the singular values of x. the */
/* columns of u are the corresponding left singular vectors, */
/* and the columns of v the right singular vectors. */
/* on entry */
/* x complex(ldx,p), where ldx.ge.n. */
/* x contains the matrix whose singular value */
/* decomposition is to be computed. x is */
/* destroyed by csvdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* ldu integer. */
/* ldu is the leading dimension of the array u */
/* (see below). */
/* ldv integer. */
/* ldv is the leading dimension of the array v */
/* (see below). */
/* work complex(n). */
/* work is a scratch array. */
/* job integer. */
/* job controls the computation of the singular */
/* vectors. it has the decimal expansion ab */
/* with the following meaning */
/* a.eq.0 do not compute the left singular */
/* vectors. */
/* a.eq.1 return the n left singular vectors */
/* in u. */
/* a.ge.2 returns the first min(n,p) */
/* left singular vectors in u. */
/* b.eq.0 do not compute the right singular */
/* vectors. */
/* b.eq.1 return the right singular vectors */
//.........这里部分代码省略.........
开发者ID:151706061，项目名称:ITK，代码行数:101，代码来源:csvdc.c
示例8: r_imag
/* Subroutine */ int chseqr_(char *job, char *compz, integer *n, integer *ilo,
integer *ihi, complex *h__, integer *ldh, complex *w, complex *z__,
integer *ldz, complex *work, integer *lwork, integer *info)
{
/* System generated locals */
address a__1[2];
integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4[2],
i__5, i__6;
real r__1, r__2, r__3, r__4;
complex q__1;
char ch__1[2];
/* Builtin functions */
double r_imag(complex *);
void r_cnjg(complex *, complex *);
/* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
/* Local variables */
static integer maxb, ierr;
static real unfl;
static complex temp;
static real ovfl, opst;
static integer i__, j, k, l;
static complex s[225] /* was [15][15] */;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *);
static complex v[16];
extern logical lsame_(char *, char *);
extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
, complex *, integer *, complex *, integer *, complex *, complex *
, integer *), ccopy_(integer *, complex *, integer *,
complex *, integer *);
static integer itemp;
static real rtemp;
static integer i1, i2;
static logical initz, wantt, wantz;
static real rwork[1];
extern doublereal slapy2_(real *, real *);
static integer ii, nh;
extern /* Subroutine */ int slabad_(real *, real *), clarfg_(integer *,
complex *, complex *, integer *, complex *);
static integer nr, ns;
extern integer icamax_(integer *, complex *, integer *);
static integer nv;
extern doublereal slamch_(char *), clanhs_(char *, integer *,
complex *, integer *, real *);
extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer
*), clahqr_(logical *, logical *, integer *, integer *, integer *,
complex *, integer *, complex *, integer *, integer *, complex *,
integer *, integer *), clacpy_(char *, integer *, integer *,
complex *, integer *, complex *, integer *);
static complex vv[16];
extern /* Subroutine */ int claset_(char *, integer *, integer *, complex
*, complex *, complex *, integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
extern /* Subroutine */ int clarfx_(char *, integer *, integer *, complex
*, complex *, complex *, integer *, complex *), xerbla_(
char *, integer *);
static real smlnum;
static logical lquery;
static integer itn;
static complex tau;
static integer its;
static real ulp, tst1;
#define h___subscr(a_1,a_2) (a_2)*h_dim1 + a_1
#define h___ref(a_1,a_2) h__[h___subscr(a_1,a_2)]
#define s_subscr(a_1,a_2) (a_2)*15 + a_1 - 16
#define s_ref(a_1,a_2) s[s_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
Common block to return operation count.
Purpose
=======
CHSEQR computes the eigenvalues of a complex upper Hessenberg
matrix H, and, optionally, the matrices T and Z from the Schur
decomposition H = Z T Z**H, where T is an upper triangular matrix
(the Schur form), and Z is the unitary matrix of Schur vectors.
Optionally Z may be postmultiplied into an input unitary matrix Q,
so that this routine can give the Schur factorization of a matrix A
which has been reduced to the Hessenberg form H by the unitary
matrix Q: A = Q*H*Q**H = (QZ)*T*(QZ)**H.
Arguments
=========
JOB (input) CHARACTER*1
= 'E': compute eigenvalues only;
//.........这里部分代码省略.........
/* Subroutine */ int cgetf2_(integer *m, integer *n, complex *a, integer *lda,
integer *ipiv, integer *info)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
complex q__1;
/* Builtin functions */
double c_abs(complex *);
void c_div(complex *, complex *, complex *);
/* Local variables */
integer i__, j, jp;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *), cgeru_(integer *, integer *, complex *, complex *,
integer *, complex *, integer *, complex *, integer *);
real sfmin;
extern /* Subroutine */ int cswap_(integer *, complex *, integer *,
complex *, integer *);
extern integer icamax_(integer *, complex *, integer *);
extern doublereal slamch_(char *);
extern /* Subroutine */ int xerbla_(char *, integer *);
/* -- LAPACK routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CGETF2 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 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 ctrti2_(char *uplo, char *diag, integer *n, complex *a,
integer *lda, 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
=======
CTRTI2 computes the inverse of a complex upper or lower triangular
matrix.
This is the Level 2 BLAS version of the algorithm.
Arguments
=========
UPLO (input) CHARACTER*1
Specifies whether the matrix A is upper or lower triangular.
= 'U': Upper triangular
= 'L': Lower triangular
DIAG (input) CHARACTER*1
Specifies whether or not the matrix A is unit triangular.
= 'N': Non-unit triangular
= 'U': Unit triangular
N (input) INTEGER
The order of the matrix A. N >= 0.
A (input/output) COMPLEX array, dimension (LDA,N)
On entry, the triangular matrix A. If UPLO = 'U', the
leading n by n upper triangular part of the array A contains
the upper triangular matrix, and the strictly lower
triangular part of A is not referenced. If UPLO = 'L', the
leading n by n lower triangular part of the array A contains
the lower triangular matrix, and the strictly upper
triangular part of A is not referenced. If DIAG = 'U', the
diagonal elements of A are also not referenced and are
assumed to be 1.
On exit, the (triangular) inverse of the original matrix, in
the same storage format.
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 complex c_b1 = {1.f,0.f};
static integer c__1 = 1;
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
complex q__1;
/* Builtin functions */
void c_div(complex *, complex *, complex *);
/* Local variables */
static integer j;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *);
extern logical lsame_(char *, char *);
static logical upper;
extern /* Subroutine */ int ctrmv_(char *, char *, char *, integer *,
complex *, integer *, complex *, integer *), xerbla_(char *, integer *);
static logical nounit;
static complex ajj;
#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");
nounit = lsame_(diag, "N");
if (! upper && ! lsame_(uplo, "L")) {
*info = -1;
} else if (! nounit && ! lsame_(diag, "U")) {
*info = -2;
} else if (*n < 0) {
*info = -3;
} else if (*lda < max(1,*n)) {
*info = -5;
}
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论