#include #include #include #include #include #ifdef complex #undef complex #endif #ifdef I #undef I #endif #if defined(_WIN64) typedef long long BLASLONG; typedef unsigned long long BLASULONG; #else typedef long BLASLONG; typedef unsigned long BLASULONG; #endif #ifdef LAPACK_ILP64 typedef BLASLONG blasint; #if defined(_WIN64) #define blasabs(x) llabs(x) #else #define blasabs(x) labs(x) #endif #else typedef int blasint; #define blasabs(x) abs(x) #endif typedef blasint integer; typedef unsigned int uinteger; typedef char *address; typedef short int shortint; typedef float real; typedef double doublereal; typedef struct { real r, i; } complex; typedef struct { doublereal r, i; } doublecomplex; #ifdef _MSC_VER static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;} static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;} static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;} static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;} #else static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;} static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;} static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;} static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;} #endif #define pCf(z) (*_pCf(z)) #define pCd(z) (*_pCd(z)) typedef blasint logical; typedef char logical1; typedef char integer1; #define TRUE_ (1) #define FALSE_ (0) /* Extern is for use with -E */ #ifndef Extern #define Extern extern #endif /* I/O stuff */ typedef int flag; typedef int ftnlen; typedef int ftnint; /*external read, write*/ typedef struct { flag cierr; ftnint ciunit; flag ciend; char *cifmt; ftnint cirec; } cilist; /*internal read, write*/ typedef struct { flag icierr; char *iciunit; flag iciend; char *icifmt; ftnint icirlen; ftnint icirnum; } icilist; /*open*/ typedef struct { flag oerr; ftnint ounit; char *ofnm; ftnlen ofnmlen; char *osta; char *oacc; char *ofm; ftnint orl; char *oblnk; } olist; /*close*/ typedef struct { flag cerr; ftnint cunit; char *csta; } cllist; /*rewind, backspace, endfile*/ typedef struct { flag aerr; ftnint aunit; } alist; /* inquire */ typedef struct { flag inerr; ftnint inunit; char *infile; ftnlen infilen; ftnint *inex; /*parameters in standard's order*/ ftnint *inopen; ftnint *innum; ftnint *innamed; char *inname; ftnlen innamlen; char *inacc; ftnlen inacclen; char *inseq; ftnlen inseqlen; char *indir; ftnlen indirlen; char *infmt; ftnlen infmtlen; char *inform; ftnint informlen; char *inunf; ftnlen inunflen; ftnint *inrecl; ftnint *innrec; char *inblank; ftnlen inblanklen; } inlist; #define VOID void union Multitype { /* for multiple entry points */ integer1 g; shortint h; integer i; /* longint j; */ real r; doublereal d; complex c; doublecomplex z; }; typedef union Multitype Multitype; struct Vardesc { /* for Namelist */ char *name; char *addr; ftnlen *dims; int type; }; typedef struct Vardesc Vardesc; struct Namelist { char *name; Vardesc **vars; int nvars; }; typedef struct Namelist Namelist; #define abs(x) ((x) >= 0 ? (x) : -(x)) #define dabs(x) (fabs(x)) #define f2cmin(a,b) ((a) <= (b) ? (a) : (b)) #define f2cmax(a,b) ((a) >= (b) ? (a) : (b)) #define dmin(a,b) (f2cmin(a,b)) #define dmax(a,b) (f2cmax(a,b)) #define bit_test(a,b) ((a) >> (b) & 1) #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b))) #define bit_set(a,b) ((a) | ((uinteger)1 << (b))) #define abort_() { sig_die("Fortran abort routine called", 1); } #define c_abs(z) (cabsf(Cf(z))) #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); } #ifdef _MSC_VER #define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);} #define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);} #else #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} #endif #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));} #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));} #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));} //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));} #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));} #define d_abs(x) (fabs(*(x))) #define d_acos(x) (acos(*(x))) #define d_asin(x) (asin(*(x))) #define d_atan(x) (atan(*(x))) #define d_atn2(x, y) (atan2(*(x),*(y))) #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); } #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); } #define d_cos(x) (cos(*(x))) #define d_cosh(x) (cosh(*(x))) #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 ) #define d_exp(x) (exp(*(x))) #define d_imag(z) (cimag(Cd(z))) #define r_imag(z) (cimagf(Cf(z))) #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x))) #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x))) #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) ) #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) ) #define d_log(x) (log(*(x))) #define d_mod(x, y) (fmod(*(x), *(y))) #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x))) #define d_nint(x) u_nint(*(x)) #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a))) #define d_sign(a,b) u_sign(*(a),*(b)) #define r_sign(a,b) u_sign(*(a),*(b)) #define d_sin(x) (sin(*(x))) #define d_sinh(x) (sinh(*(x))) #define d_sqrt(x) (sqrt(*(x))) #define d_tan(x) (tan(*(x))) #define d_tanh(x) (tanh(*(x))) #define i_abs(x) abs(*(x)) #define i_dnnt(x) ((integer)u_nint(*(x))) #define i_len(s, n) (n) #define i_nint(x) ((integer)u_nint(*(x))) #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b))) #define pow_dd(ap, bp) ( pow(*(ap), *(bp))) #define pow_si(B,E) spow_ui(*(B),*(E)) #define pow_ri(B,E) spow_ui(*(B),*(E)) #define pow_di(B,E) dpow_ui(*(B),*(E)) #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));} #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));} #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));} #define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; } #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d)))) #define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; } #define sig_die(s, kill) { exit(1); } #define s_stop(s, n) {exit(0);} static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n"; #define z_abs(z) (cabs(Cd(z))) #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));} #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));} #define myexit_() break; #define mycycle_() continue; #define myceiling_(w) {ceil(w)} #define myhuge_(w) {HUGE_VAL} //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n) /* procedure parameter types for -A and -C++ */ #ifdef __cplusplus typedef logical (*L_fp)(...); #else typedef logical (*L_fp)(); #endif static float spow_ui(float x, integer n) { float pow=1.0; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x = 1/x; for(u = n; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return pow; } static double dpow_ui(double x, integer n) { double pow=1.0; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x = 1/x; for(u = n; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return pow; } #ifdef _MSC_VER static _Fcomplex cpow_ui(complex x, integer n) { complex pow={1.0,0.0}; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i; for(u = n; ; ) { if(u & 01) pow.r *= x.r, pow.i *= x.i; if(u >>= 1) x.r *= x.r, x.i *= x.i; else break; } } _Fcomplex p={pow.r, pow.i}; return p; } #else static _Complex float cpow_ui(_Complex float x, integer n) { _Complex float pow=1.0; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x = 1/x; for(u = n; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return pow; } #endif #ifdef _MSC_VER static _Dcomplex zpow_ui(_Dcomplex x, integer n) { _Dcomplex pow={1.0,0.0}; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1]; for(u = n; ; ) { if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1]; if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1]; else break; } } _Dcomplex p = {pow._Val[0], pow._Val[1]}; return p; } #else static _Complex double zpow_ui(_Complex double x, integer n) { _Complex double pow=1.0; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x = 1/x; for(u = n; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return pow; } #endif static integer pow_ii(integer x, integer n) { integer pow; unsigned long int u; if (n <= 0) { if (n == 0 || x == 1) pow = 1; else if (x != -1) pow = x == 0 ? 1/x : 0; else n = -n; } if ((n > 0) || !(n == 0 || x == 1 || x != -1)) { u = n; for(pow = 1; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return pow; } static integer dmaxloc_(double *w, integer s, integer e, integer *n) { double m; integer i, mi; for(m=w[s-1], mi=s, i=s+1; i<=e; i++) if (w[i-1]>m) mi=i ,m=w[i-1]; return mi-s+1; } static integer smaxloc_(float *w, integer s, integer e, integer *n) { float m; integer i, mi; for(m=w[s-1], mi=s, i=s+1; i<=e; i++) if (w[i-1]>m) mi=i ,m=w[i-1]; return mi-s+1; } static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) { integer n = *n_, incx = *incx_, incy = *incy_, i; #ifdef _MSC_VER _Fcomplex zdotc = {0.0, 0.0}; if (incx == 1 && incy == 1) { for (i=0;i= 0 */ /* The state space dimension (the row dimension of X, Y). */ /* ..... */ /* N (input) INTEGER, 0 <= N <= M */ /* The number of data snapshot pairs */ /* (the number of columns of X and Y). */ /* ..... */ /* X (input/output) REAL(KIND=WP) M-by-N array */ /* > On entry, X contains the data snapshot matrix X. It is */ /* assumed that the column norms of X are in the range of */ /* the normalized floating point numbers. */ /* < On exit, the leading K columns of X contain a POD basis, */ /* i.e. the leading K left singular vectors of the input */ /* data matrix X, U(:,1:K). All N columns of X contain all */ /* left singular vectors of the input matrix X. */ /* See the descriptions of K, Z and W. */ /* ..... */ /* LDX (input) INTEGER, LDX >= M */ /* The leading dimension of the array X. */ /* ..... */ /* Y (input/workspace/output) REAL(KIND=WP) M-by-N array */ /* > On entry, Y contains the data snapshot matrix Y */ /* < On exit, */ /* If JOBR == 'R', the leading K columns of Y contain */ /* the residual vectors for the computed Ritz pairs. */ /* See the description of RES. */ /* If JOBR == 'N', Y contains the original input data, */ /* scaled according to the value of JOBS. */ /* ..... */ /* LDY (input) INTEGER , LDY >= M */ /* The leading dimension of the array Y. */ /* ..... */ /* NRNK (input) INTEGER */ /* Determines the mode how to compute the numerical rank, */ /* i.e. how to truncate small singular values of the input */ /* matrix X. On input, if */ /* NRNK = -1 :: i-th singular value sigma(i) is truncated */ /* if sigma(i) <= TOL*sigma(1). */ /* This option is recommended. */ /* NRNK = -2 :: i-th singular value sigma(i) is truncated */ /* if sigma(i) <= TOL*sigma(i-1) */ /* This option is included for R&D purposes. */ /* It requires highly accurate SVD, which */ /* may not be feasible. */ /* The numerical rank can be enforced by using positive */ /* value of NRNK as follows: */ /* 0 < NRNK <= N :: at most NRNK largest singular values */ /* will be used. If the number of the computed nonzero */ /* singular values is less than NRNK, then only those */ /* nonzero values will be used and the actually used */ /* dimension is less than NRNK. The actual number of */ /* the nonzero singular values is returned in the variable */ /* K. See the descriptions of TOL and K. */ /* ..... */ /* TOL (input) REAL(KIND=WP), 0 <= TOL < 1 */ /* The tolerance for truncating small singular values. */ /* See the description of NRNK. */ /* ..... */ /* K (output) INTEGER, 0 <= K <= N */ /* The dimension of the POD basis for the data snapshot */ /* matrix X and the number of the computed Ritz pairs. */ /* The value of K is determined according to the rule set */ /* by the parameters NRNK and TOL. */ /* See the descriptions of NRNK and TOL. */ /* ..... */ /* REIG (output) REAL(KIND=WP) N-by-1 array */ /* The leading K (K<=N) entries of REIG contain */ /* the real parts of the computed eigenvalues */ /* REIG(1:K) + sqrt(-1)*IMEIG(1:K). */ /* See the descriptions of K, IMEIG, and Z. */ /* ..... */ /* IMEIG (output) REAL(KIND=WP) N-by-1 array */ /* The leading K (K<=N) entries of IMEIG contain */ /* the imaginary parts of the computed eigenvalues */ /* REIG(1:K) + sqrt(-1)*IMEIG(1:K). */ /* The eigenvalues are determined as follows: */ /* If IMEIG(i) == 0, then the corresponding eigenvalue is */ /* real, LAMBDA(i) = REIG(i). */ /* If IMEIG(i)>0, then the corresponding complex */ /* conjugate pair of eigenvalues reads */ /* LAMBDA(i) = REIG(i) + sqrt(-1)*IMAG(i) */ /* LAMBDA(i+1) = REIG(i) - sqrt(-1)*IMAG(i) */ /* That is, complex conjugate pairs have consecutive */ /* indices (i,i+1), with the positive imaginary part */ /* listed first. */ /* See the descriptions of K, REIG, and Z. */ /* ..... */ /* Z (workspace/output) REAL(KIND=WP) M-by-N array */ /* If JOBZ =='V' then */ /* Z contains real Ritz vectors as follows: */ /* If IMEIG(i)=0, then Z(:,i) is an eigenvector of */ /* the i-th Ritz value; ||Z(:,i)||_2=1. */ /* If IMEIG(i) > 0 (and IMEIG(i+1) < 0) then */ /* [Z(:,i) Z(:,i+1)] span an invariant subspace and */ /* the Ritz values extracted from this subspace are */ /* REIG(i) + sqrt(-1)*IMEIG(i) and */ /* REIG(i) - sqrt(-1)*IMEIG(i). */ /* The corresponding eigenvectors are */ /* Z(:,i) + sqrt(-1)*Z(:,i+1) and */ /* Z(:,i) - sqrt(-1)*Z(:,i+1), respectively. */ /* || Z(:,i:i+1)||_F = 1. */ /* If JOBZ == 'F', then the above descriptions hold for */ /* the columns of X(:,1:K)*W(1:K,1:K), where the columns */ /* of W(1:k,1:K) are the computed eigenvectors of the */ /* K-by-K Rayleigh quotient. The columns of W(1:K,1:K) */ /* are similarly structured: If IMEIG(i) == 0 then */ /* X(:,1:K)*W(:,i) is an eigenvector, and if IMEIG(i)>0 */ /* then X(:,1:K)*W(:,i)+sqrt(-1)*X(:,1:K)*W(:,i+1) and */ /* X(:,1:K)*W(:,i)-sqrt(-1)*X(:,1:K)*W(:,i+1) */ /* are the eigenvectors of LAMBDA(i), LAMBDA(i+1). */ /* See the descriptions of REIG, IMEIG, X and W. */ /* ..... */ /* LDZ (input) INTEGER , LDZ >= M */ /* The leading dimension of the array Z. */ /* ..... */ /* RES (output) REAL(KIND=WP) N-by-1 array */ /* RES(1:K) contains the residuals for the K computed */ /* Ritz pairs. */ /* If LAMBDA(i) is real, then */ /* RES(i) = || A * Z(:,i) - LAMBDA(i)*Z(:,i))||_2. */ /* If [LAMBDA(i), LAMBDA(i+1)] is a complex conjugate pair */ /* then */ /* RES(i)=RES(i+1) = || A * Z(:,i:i+1) - Z(:,i:i+1) *B||_F */ /* where B = [ real(LAMBDA(i)) imag(LAMBDA(i)) ] */ /* [-imag(LAMBDA(i)) real(LAMBDA(i)) ]. */ /* It holds that */ /* RES(i) = || A*ZC(:,i) - LAMBDA(i) *ZC(:,i) ||_2 */ /* RES(i+1) = || A*ZC(:,i+1) - LAMBDA(i+1)*ZC(:,i+1) ||_2 */ /* where ZC(:,i) = Z(:,i) + sqrt(-1)*Z(:,i+1) */ /* ZC(:,i+1) = Z(:,i) - sqrt(-1)*Z(:,i+1) */ /* See the description of REIG, IMEIG and Z. */ /* ..... */ /* B (output) REAL(KIND=WP) M-by-N array. */ /* IF JOBF =='R', B(1:M,1:K) contains A*U(:,1:K), and can */ /* be used for computing the refined vectors; see further */ /* details in the provided references. */ /* If JOBF == 'E', B(1:M,1;K) contains */ /* A*U(:,1:K)*W(1:K,1:K), which are the vectors from the */ /* Exact DMD, up to scaling by the inverse eigenvalues. */ /* If JOBF =='N', then B is not referenced. */ /* See the descriptions of X, W, K. */ /* ..... */ /* LDB (input) INTEGER, LDB >= M */ /* The leading dimension of the array B. */ /* ..... */ /* W (workspace/output) REAL(KIND=WP) N-by-N array */ /* On exit, W(1:K,1:K) contains the K computed */ /* eigenvectors of the matrix Rayleigh quotient (real and */ /* imaginary parts for each complex conjugate pair of the */ /* eigenvalues). The Ritz vectors (returned in Z) are the */ /* product of X (containing a POD basis for the input */ /* matrix X) and W. See the descriptions of K, S, X and Z. */ /* W is also used as a workspace to temporarily store the */ /* right singular vectors of X. */ /* ..... */ /* LDW (input) INTEGER, LDW >= N */ /* The leading dimension of the array W. */ /* ..... */ /* S (workspace/output) REAL(KIND=WP) N-by-N array */ /* The array S(1:K,1:K) is used for the matrix Rayleigh */ /* quotient. This content is overwritten during */ /* the eigenvalue decomposition by DGEEV. */ /* See the description of K. */ /* ..... */ /* LDS (input) INTEGER, LDS >= N */ /* The leading dimension of the array S. */ /* ..... */ /* WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array */ /* On exit, WORK(1:N) contains the singular values of */ /* X (for JOBS=='N') or column scaled X (JOBS=='S', 'C'). */ /* If WHTSVD==4, then WORK(N+1) and WORK(N+2) contain */ /* scaling factor WORK(N+2)/WORK(N+1) used to scale X */ /* and Y to avoid overflow in the SVD of X. */ /* This may be of interest if the scaling option is off */ /* and as many as possible smallest eigenvalues are */ /* desired to the highest feasible accuracy. */ /* If the call to DGEDMD is only workspace query, then */ /* WORK(1) contains the minimal workspace length and */ /* WORK(2) is the optimal workspace length. Hence, the */ /* leng of work is at least 2. */ /* See the description of LWORK. */ /* ..... */ /* LWORK (input) INTEGER */ /* The minimal length of the workspace vector WORK. */ /* LWORK is calculated as follows: */ /* If WHTSVD == 1 :: */ /* If JOBZ == 'V', then */ /* LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)). */ /* If JOBZ == 'N' then */ /* LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)). */ /* Here LWORK_SVD = MAX(1,3*N+M,5*N) is the minimal */ /* workspace length of DGESVD. */ /* If WHTSVD == 2 :: */ /* If JOBZ == 'V', then */ /* LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,4*N)) */ /* If JOBZ == 'N', then */ /* LWORK >= MAX(2, N + LWORK_SVD, N+MAX(1,3*N)) */ /* Here LWORK_SVD = MAX(M, 5*N*N+4*N)+3*N*N is the */ /* minimal workspace length of DGESDD. */ /* If WHTSVD == 3 :: */ /* If JOBZ == 'V', then */ /* LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N)) */ /* If JOBZ == 'N', then */ /* LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N)) */ /* Here LWORK_SVD = N+M+MAX(3*N+1, */ /* MAX(1,3*N+M,5*N),MAX(1,N)) */ /* is the minimal workspace length of DGESVDQ. */ /* If WHTSVD == 4 :: */ /* If JOBZ == 'V', then */ /* LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,4*N)) */ /* If JOBZ == 'N', then */ /* LWORK >= MAX(2, N+LWORK_SVD,N+MAX(1,3*N)) */ /* Here LWORK_SVD = MAX(7,2*M+N,6*N+2*N*N) is the */ /* minimal workspace length of DGEJSV. */ /* The above expressions are not simplified in order to */ /* make the usage of WORK more transparent, and for */ /* easier checking. In any case, LWORK >= 2. */ /* If on entry LWORK = -1, then a workspace query is */ /* assumed and the procedure only computes the minimal */ /* and the optimal workspace lengths for both WORK and */ /* IWORK. See the descriptions of WORK and IWORK. */ /* ..... */ /* IWORK (workspace/output) INTEGER LIWORK-by-1 array */ /* Workspace that is required only if WHTSVD equals */ /* 2 , 3 or 4. (See the description of WHTSVD). */ /* If on entry LWORK =-1 or LIWORK=-1, then the */ /* minimal length of IWORK is computed and returned in */ /* IWORK(1). See the description of LIWORK. */ /* ..... */ /* LIWORK (input) INTEGER */ /* The minimal length of the workspace vector IWORK. */ /* If WHTSVD == 1, then only IWORK(1) is used; LIWORK >=1 */ /* If WHTSVD == 2, then LIWORK >= MAX(1,8*MIN(M,N)) */ /* If WHTSVD == 3, then LIWORK >= MAX(1,M+N-1) */ /* If WHTSVD == 4, then LIWORK >= MAX(3,M+3*N) */ /* If on entry LIWORK = -1, then a workspace query is */ /* assumed and the procedure only computes the minimal */ /* and the optimal workspace lengths for both WORK and */ /* IWORK. See the descriptions of WORK and IWORK. */ /* ..... */ /* INFO (output) INTEGER */ /* -i < 0 :: On entry, the i-th argument had an */ /* illegal value */ /* = 0 :: Successful return. */ /* = 1 :: Void input. Quick exit (M=0 or N=0). */ /* = 2 :: The SVD computation of X did not converge. */ /* Suggestion: Check the input data and/or */ /* repeat with different WHTSVD. */ /* = 3 :: The computation of the eigenvalues did not */ /* converge. */ /* = 4 :: If data scaling was requested on input and */ /* the procedure found inconsistency in the data */ /* such that for some column index i, */ /* X(:,i) = 0 but Y(:,i) /= 0, then Y(:,i) is set */ /* to zero if JOBS=='C'. The computation proceeds */ /* with original or modified data and warning */ /* flag is set with INFO=4. */ /* ............................................................. */ /* ............................................................. */ /* Parameters */ /* ~~~~~~~~~~ */ /* Local scalars */ /* ~~~~~~~~~~~~~ */ /* Local arrays */ /* ~~~~~~~~~~~~ */ /* External functions (BLAS and LAPACK) */ /* ~~~~~~~~~~~~~~~~~ */ /* External subroutines (BLAS and LAPACK) */ /* ~~~~~~~~~~~~~~~~~~~~ */ /* Intrinsic functions */ /* ~~~~~~~~~~~~~~~~~~~ */ /* ............................................................ */ /* Parameter adjustments */ x_dim1 = *ldx; x_offset = 1 + x_dim1 * 1; x -= x_offset; y_dim1 = *ldy; y_offset = 1 + y_dim1 * 1; y -= y_offset; --reig; --imeig; z_dim1 = *ldz; z_offset = 1 + z_dim1 * 1; z__ -= z_offset; --res; b_dim1 = *ldb; b_offset = 1 + b_dim1 * 1; b -= b_offset; w_dim1 = *ldw; w_offset = 1 + w_dim1 * 1; w -= w_offset; s_dim1 = *lds; s_offset = 1 + s_dim1 * 1; s -= s_offset; --work; --iwork; /* Function Body */ one = 1.f; zero = 0.f; /* Test the input arguments */ wntres = lsame_(jobr, "R"); sccolx = lsame_(jobs, "S") || lsame_(jobs, "C"); sccoly = lsame_(jobs, "Y"); wntvec = lsame_(jobz, "V"); wntref = lsame_(jobf, "R"); wntex = lsame_(jobf, "E"); *info = 0; lquery = *lwork == -1 || *liwork == -1; if (! (sccolx || sccoly || lsame_(jobs, "N"))) { *info = -1; } else if (! (wntvec || lsame_(jobz, "N") || lsame_( jobz, "F"))) { *info = -2; } else if (! (wntres || lsame_(jobr, "N")) || wntres && ! wntvec) { *info = -3; } else if (! (wntref || wntex || lsame_(jobf, "N"))) { *info = -4; } else if (! (*whtsvd == 1 || *whtsvd == 2 || *whtsvd == 3 || *whtsvd == 4)) { *info = -5; } else if (*m < 0) { *info = -6; } else if (*n < 0 || *n > *m) { *info = -7; } else if (*ldx < *m) { *info = -9; } else if (*ldy < *m) { *info = -11; } else if (! (*nrnk == -2 || *nrnk == -1 || *nrnk >= 1 && *nrnk <= *n)) { *info = -12; } else if (*tol < zero || *tol >= one) { *info = -13; } else if (*ldz < *m) { *info = -18; } else if ((wntref || wntex) && *ldb < *m) { *info = -21; } else if (*ldw < *n) { *info = -23; } else if (*lds < *n) { *info = -25; } if (*info == 0) { /* Compute the minimal and the optimal workspace */ /* requirements. Simulate running the code and */ /* determine minimal and optimal sizes of the */ /* workspace at any moment of the run. */ if (*n == 0) { /* Quick return. All output except K is void. */ /* INFO=1 signals the void input. */ /* In case of a workspace query, the default */ /* minimal workspace lengths are returned. */ if (lquery) { iwork[1] = 1; work[1] = 2.; work[2] = 2.; } else { *k = 0; } *info = 1; return 0; } mlwork = f2cmax(2,*n); olwork = f2cmax(2,*n); iminwr = 1; /* SELECT CASE ( WHTSVD ) */ if (*whtsvd == 1) { /* The following is specified as the minimal */ /* length of WORK in the definition of DGESVD: */ /* MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) */ /* Computing MAX */ i__1 = 1, i__2 = f2cmin(*m,*n) * 3 + f2cmax(*m,*n), i__1 = f2cmax(i__1, i__2), i__2 = f2cmin(*m,*n) * 5; mwrsvd = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = mlwork, i__2 = *n + mwrsvd; mlwork = f2cmax(i__1,i__2); if (lquery) { dgesvd_("O", "S", m, n, &x[x_offset], ldx, &work[1], &b[ b_offset], ldb, &w[w_offset], ldw, rdummy, &c_n1, & info1); /* Computing MAX */ i__1 = mwrsvd, i__2 = (integer) rdummy[0]; lwrsvd = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = olwork, i__2 = *n + lwrsvd; olwork = f2cmax(i__1,i__2); } } else if (*whtsvd == 2) { /* The following is specified as the minimal */ /* length of WORK in the definition of DGESDD: */ /* MWRSDD = 3*MIN(M,N)*MIN(M,N) + */ /* MAX( MAX(M,N),5*MIN(M,N)*MIN(M,N)+4*MIN(M,N) ) */ /* IMINWR = 8*MIN(M,N) */ /* Computing MAX */ i__1 = f2cmax(*m,*n), i__2 = f2cmin(*m,*n) * 5 * f2cmin(*m,*n) + (f2cmin(*m,* n) << 2); mwrsdd = f2cmin(*m,*n) * 3 * f2cmin(*m,*n) + f2cmax(i__1,i__2); /* Computing MAX */ i__1 = mlwork, i__2 = *n + mwrsdd; mlwork = f2cmax(i__1,i__2); iminwr = f2cmin(*m,*n) << 3; if (lquery) { dgesdd_("O", m, n, &x[x_offset], ldx, &work[1], &b[b_offset], ldb, &w[w_offset], ldw, rdummy, &c_n1, &iwork[1], & info1); /* Computing MAX */ i__1 = mwrsdd, i__2 = (integer) rdummy[0]; lwrsdd = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = olwork, i__2 = *n + lwrsdd; olwork = f2cmax(i__1,i__2); } } else if (*whtsvd == 3) { /* LWQP3 = 3*N+1 */ /* LWORQ = MAX(N, 1) */ /* MWRSVD = MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) */ /* MWRSVQ = N + MAX( LWQP3, MWRSVD, LWORQ ) + MAX(M,2) */ /* MLWORK = N + MWRSVQ */ /* IMINWR = M+N-1 */ dgesvdq_("H", "P", "N", "R", "R", m, n, &x[x_offset], ldx, &work[ 1], &z__[z_offset], ldz, &w[w_offset], ldw, &numrnk, & iwork[1], liwork, rdummy, &c_n1, rdummy2, &c_n1, &info1); iminwr = iwork[1]; mwrsvq = (integer) rdummy[1]; /* Computing MAX */ i__1 = mlwork, i__2 = *n + mwrsvq + (integer) rdummy2[0]; mlwork = f2cmax(i__1,i__2); if (lquery) { /* Computing MAX */ i__1 = mwrsvq, i__2 = (integer) rdummy[0]; lwrsvq = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = olwork, i__2 = *n + lwrsvq + (integer) rdummy2[0]; olwork = f2cmax(i__1,i__2); } } else if (*whtsvd == 4) { *(unsigned char *)jsvopt = 'J'; /* MWRSVJ = MAX( 7, 2*M+N, 6*N+2*N*N ) ! for JSVOPT='V' */ /* Computing MAX */ i__1 = 7, i__2 = (*m << 1) + *n, i__1 = f2cmax(i__1,i__2), i__2 = (* n << 2) + *n * *n, i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + *n * *n + 6; mwrsvj = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = mlwork, i__2 = *n + mwrsvj; mlwork = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = 3, i__2 = *m + *n * 3; iminwr = f2cmax(i__1,i__2); if (lquery) { /* Computing MAX */ i__1 = olwork, i__2 = *n + mwrsvj; olwork = f2cmax(i__1,i__2); } /* END SELECT */ } if (wntvec || wntex || lsame_(jobz, "F")) { *(unsigned char *)jobzl = 'V'; } else { *(unsigned char *)jobzl = 'N'; } /* Workspace calculation to the DGEEV call */ if (lsame_(jobzl, "V")) { /* Computing MAX */ i__1 = 1, i__2 = *n << 2; mwrkev = f2cmax(i__1,i__2); } else { /* Computing MAX */ i__1 = 1, i__2 = *n * 3; mwrkev = f2cmax(i__1,i__2); } /* Computing MAX */ i__1 = mlwork, i__2 = *n + mwrkev; mlwork = f2cmax(i__1,i__2); if (lquery) { dgeev_("N", jobzl, n, &s[s_offset], lds, &reig[1], &imeig[1], &w[ w_offset], ldw, &w[w_offset], ldw, rdummy, &c_n1, &info1); /* Computing MAX */ i__1 = mwrkev, i__2 = (integer) rdummy[0]; lwrkev = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = olwork, i__2 = *n + lwrkev; olwork = f2cmax(i__1,i__2); } if (*liwork < iminwr && ! lquery) { *info = -29; } if (*lwork < mlwork && ! lquery) { *info = -27; } } if (*info != 0) { i__1 = -(*info); xerbla_("DGEDMD", &i__1); return 0; } else if (lquery) { /* Return minimal and optimal workspace sizes */ iwork[1] = iminwr; work[1] = (doublereal) mlwork; work[2] = (doublereal) olwork; return 0; } /* ............................................................ */ ofl = dlamch_("O"); small = dlamch_("S"); badxy = FALSE_; /* <1> Optional scaling of the snapshots (columns of X, Y) */ /* ========================================================== */ if (sccolx) { /* The columns of X will be normalized. */ /* To prevent overflows, the column norms of X are */ /* carefully computed using DLASSQ. */ *k = 0; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* WORK(i) = DNRM2( M, X(1,i), 1 ) */ scale = zero; dlassq_(m, &x[i__ * x_dim1 + 1], &c__1, &scale, &ssum); if (disnan_(&scale) || disnan_(&ssum)) { *k = 0; *info = -8; i__2 = -(*info); xerbla_("DGEDMD", &i__2); } if (scale != zero && ssum != zero) { rootsc = sqrt(ssum); if (scale >= ofl / rootsc) { /* Norm of X(:,i) overflows. First, X(:,i) */ /* is scaled by */ /* ( ONE / ROOTSC ) / SCALE = 1/||X(:,i)||_2. */ /* Next, the norm of X(:,i) is stored without */ /* overflow as WORK(i) = - SCALE * (ROOTSC/M), */ /* the minus sign indicating the 1/M factor. */ /* Scaling is performed without overflow, and */ /* underflow may occur in the smallest entries */ /* of X(:,i). The relative backward and forward */ /* errors are small in the ell_2 norm. */ d__1 = one / rootsc; dlascl_("G", &c__0, &c__0, &scale, &d__1, m, &c__1, &x[ i__ * x_dim1 + 1], m, &info2); work[i__] = -scale * (rootsc / (doublereal) (*m)); } else { /* X(:,i) will be scaled to unit 2-norm */ work[i__] = scale * rootsc; dlascl_("G", &c__0, &c__0, &work[i__], &one, m, &c__1, &x[ i__ * x_dim1 + 1], m, &info2); /* X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC */ /* LAPACK */ } } else { work[i__] = zero; ++(*k); } } if (*k == *n) { /* All columns of X are zero. Return error code -8. */ /* (the 8th input variable had an illegal value) */ *k = 0; *info = -8; i__1 = -(*info); xerbla_("DGEDMD", &i__1); return 0; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* Now, apply the same scaling to the columns of Y. */ if (work[i__] > zero) { d__1 = one / work[i__]; dscal_(m, &d__1, &y[i__ * y_dim1 + 1], &c__1); /* Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC */ /* BLAS CALL */ } else if (work[i__] < zero) { d__1 = -work[i__]; d__2 = one / (doublereal) (*m); dlascl_("G", &c__0, &c__0, &d__1, &d__2, m, &c__1, &y[i__ * y_dim1 + 1], m, &info2); /* LAPACK CAL */ } else if (y[idamax_(m, &y[i__ * y_dim1 + 1], &c__1) + i__ * y_dim1] != zero) { /* X(:,i) is zero vector. For consistency, */ /* Y(:,i) should also be zero. If Y(:,i) is not */ /* zero, then the data might be inconsistent or */ /* corrupted. If JOBS == 'C', Y(:,i) is set to */ /* zero and a warning flag is raised. */ /* The computation continues but the */ /* situation will be reported in the output. */ badxy = TRUE_; if (lsame_(jobs, "C")) { dscal_(m, &zero, &y[i__ * y_dim1 + 1], &c__1); } /* BLAS CALL */ } } } if (sccoly) { /* The columns of Y will be normalized. */ /* To prevent overflows, the column norms of Y are */ /* carefully computed using DLASSQ. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* WORK(i) = DNRM2( M, Y(1,i), 1 ) */ scale = zero; dlassq_(m, &y[i__ * y_dim1 + 1], &c__1, &scale, &ssum); if (disnan_(&scale) || disnan_(&ssum)) { *k = 0; *info = -10; i__2 = -(*info); xerbla_("DGEDMD", &i__2); } if (scale != zero && ssum != zero) { rootsc = sqrt(ssum); if (scale >= ofl / rootsc) { /* Norm of Y(:,i) overflows. First, Y(:,i) */ /* is scaled by */ /* ( ONE / ROOTSC ) / SCALE = 1/||Y(:,i)||_2. */ /* Next, the norm of Y(:,i) is stored without */ /* overflow as WORK(i) = - SCALE * (ROOTSC/M), */ /* the minus sign indicating the 1/M factor. */ /* Scaling is performed without overflow, and */ /* underflow may occur in the smallest entries */ /* of Y(:,i). The relative backward and forward */ /* errors are small in the ell_2 norm. */ d__1 = one / rootsc; dlascl_("G", &c__0, &c__0, &scale, &d__1, m, &c__1, &y[ i__ * y_dim1 + 1], m, &info2); work[i__] = -scale * (rootsc / (doublereal) (*m)); } else { /* X(:,i) will be scaled to unit 2-norm */ work[i__] = scale * rootsc; dlascl_("G", &c__0, &c__0, &work[i__], &one, m, &c__1, &y[ i__ * y_dim1 + 1], m, &info2); /* Y(1:M,i) = (ONE/WORK(i)) * Y(1:M,i) ! INTRINSIC */ /* LAPACK */ } } else { work[i__] = zero; } } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* Now, apply the same scaling to the columns of X. */ if (work[i__] > zero) { d__1 = one / work[i__]; dscal_(m, &d__1, &x[i__ * x_dim1 + 1], &c__1); /* X(1:M,i) = (ONE/WORK(i)) * X(1:M,i) ! INTRINSIC */ /* BLAS CALL */ } else if (work[i__] < zero) { d__1 = -work[i__]; d__2 = one / (doublereal) (*m); dlascl_("G", &c__0, &c__0, &d__1, &d__2, m, &c__1, &x[i__ * x_dim1 + 1], m, &info2); /* LAPACK CAL */ } else if (x[idamax_(m, &x[i__ * x_dim1 + 1], &c__1) + i__ * x_dim1] != zero) { /* Y(:,i) is zero vector. If X(:,i) is not */ /* zero, then a warning flag is raised. */ /* The computation continues but the */ /* situation will be reported in the output. */ badxy = TRUE_; } } } /* <2> SVD of the data snapshot matrix X. */ /* ===================================== */ /* The left singular vectors are stored in the array X. */ /* The right singular vectors are in the array W. */ /* The array W will later on contain the eigenvectors */ /* of a Rayleigh quotient. */ numrnk = *n; /* SELECT CASE ( WHTSVD ) */ if (*whtsvd == 1) { i__1 = *lwork - *n; dgesvd_("O", "S", m, n, &x[x_offset], ldx, &work[1], &b[b_offset], ldb, &w[w_offset], ldw, &work[*n + 1], &i__1, &info1); /* LAPACK CAL */ *(unsigned char *)t_or_n__ = 'T'; } else if (*whtsvd == 2) { i__1 = *lwork - *n; dgesdd_("O", m, n, &x[x_offset], ldx, &work[1], &b[b_offset], ldb, &w[ w_offset], ldw, &work[*n + 1], &i__1, &iwork[1], &info1); /* LAPACK CAL */ *(unsigned char *)t_or_n__ = 'T'; } else if (*whtsvd == 3) { i__1 = *lwork - *n - f2cmax(2,*m); i__2 = f2cmax(2,*m); dgesvdq_("H", "P", "N", "R", "R", m, n, &x[x_offset], ldx, &work[1], & z__[z_offset], ldz, &w[w_offset], ldw, &numrnk, &iwork[1], liwork, &work[*n + f2cmax(2,*m) + 1], &i__1, &work[*n + 1], & i__2, &info1); /* L */ dlacpy_("A", m, &numrnk, &z__[z_offset], ldz, &x[x_offset], ldx); /* LAPACK C */ *(unsigned char *)t_or_n__ = 'T'; } else if (*whtsvd == 4) { i__1 = *lwork - *n; dgejsv_("F", "U", jsvopt, "N", "N", "P", m, n, &x[x_offset], ldx, & work[1], &z__[z_offset], ldz, &w[w_offset], ldw, &work[*n + 1] , &i__1, &iwork[1], &info1); /* LAPACK CALL */ dlacpy_("A", m, n, &z__[z_offset], ldz, &x[x_offset], ldx); /* LAPACK CALL */ *(unsigned char *)t_or_n__ = 'N'; xscl1 = work[*n + 1]; xscl2 = work[*n + 2]; if (xscl1 != xscl2) { /* This is an exceptional situation. If the */ /* data matrices are not scaled and the */ /* largest singular value of X overflows. */ /* In that case DGEJSV can return the SVD */ /* in scaled form. The scaling factor can be used */ /* to rescale the data (X and Y). */ dlascl_("G", &c__0, &c__0, &xscl1, &xscl2, m, n, &y[y_offset], ldy, &info2); } /* END SELECT */ } if (info1 > 0) { /* The SVD selected subroutine did not converge. */ /* Return with an error code. */ *info = 2; return 0; } if (work[1] == zero) { /* The largest computed singular value of (scaled) */ /* X is zero. Return error code -8 */ /* (the 8th input variable had an illegal value). */ *k = 0; *info = -8; i__1 = -(*info); xerbla_("DGEDMD", &i__1); return 0; } /* <3> Determine the numerical rank of the data */ /* snapshots matrix X. This depends on the */ /* parameters NRNK and TOL. */ /* SELECT CASE ( NRNK ) */ if (*nrnk == -1) { *k = 1; i__1 = numrnk; for (i__ = 2; i__ <= i__1; ++i__) { if (work[i__] <= work[1] * *tol || work[i__] <= small) { myexit_(); } ++(*k); } } else if (*nrnk == -2) { *k = 1; i__1 = numrnk - 1; for (i__ = 1; i__ <= i__1; ++i__) { if (work[i__ + 1] <= work[i__] * *tol || work[i__] <= small) { myexit_(); } ++(*k); } } else { *k = 1; i__1 = *nrnk; for (i__ = 2; i__ <= i__1; ++i__) { if (work[i__] <= small) { myexit_(); } ++(*k); } /* END SELECT */ } /* Now, U = X(1:M,1:K) is the SVD/POD basis for the */ /* snapshot data in the input matrix X. */ /* <4> Compute the Rayleigh quotient S = U^T * A * U. */ /* Depending on the requested outputs, the computation */ /* is organized to compute additional auxiliary */ /* matrices (for the residuals and refinements). */ /* In all formulas below, we need V_k*Sigma_k^(-1) */ /* where either V_k is in W(1:N,1:K), or V_k^T is in */ /* W(1:K,1:N). Here Sigma_k=diag(WORK(1:K)). */ if (lsame_(t_or_n__, "N")) { i__1 = *k; for (i__ = 1; i__ <= i__1; ++i__) { d__1 = one / work[i__]; dscal_(n, &d__1, &w[i__ * w_dim1 + 1], &c__1); /* W(1:N,i) = (ONE/WORK(i)) * W(1:N,i) ! INTRINSIC */ /* BLAS CALL */ } } else { /* This non-unit stride access is due to the fact */ /* that DGESVD, DGESVDQ and DGESDD return the */ /* transposed matrix of the right singular vectors. */ /* DO i = 1, K */ /* CALL DSCAL( N, ONE/WORK(i), W(i,1), LDW ) ! BLAS CALL */ /* ! W(i,1:N) = (ONE/WORK(i)) * W(i,1:N) ! INTRINSIC */ /* END DO */ i__1 = *k; for (i__ = 1; i__ <= i__1; ++i__) { work[*n + i__] = one / work[i__]; } i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *k; for (i__ = 1; i__ <= i__2; ++i__) { w[i__ + j * w_dim1] = work[*n + i__] * w[i__ + j * w_dim1]; } } } if (wntref) { /* Need A*U(:,1:K)=Y*V_k*inv(diag(WORK(1:K))) */ /* for computing the refined Ritz vectors */ /* (optionally, outside DGEDMD). */ dgemm_("N", t_or_n__, m, k, n, &one, &y[y_offset], ldy, &w[w_offset], ldw, &zero, &z__[z_offset], ldz); /* Z(1:M,1:K)=MATMUL(Y(1:M,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRI */ /* Z(1:M,1:K)=MATMUL(Y(1:M,1:N),W(1:N,1:K)) ! INTRI */ /* At this point Z contains */ /* A * U(:,1:K) = Y * V_k * Sigma_k^(-1), and */ /* this is needed for computing the residuals. */ /* This matrix is returned in the array B and */ /* it can be used to compute refined Ritz vectors. */ /* BLAS */ dlacpy_("A", m, k, &z__[z_offset], ldz, &b[b_offset], ldb); /* B(1:M,1:K) = Z(1:M,1:K) ! INTRINSIC */ /* BLAS CALL */ dgemm_("T", "N", k, k, m, &one, &x[x_offset], ldx, &z__[z_offset], ldz, &zero, &s[s_offset], lds); /* S(1:K,1:K) = MATMUL(TANSPOSE(X(1:M,1:K)),Z(1:M,1:K)) ! INTRI */ /* At this point S = U^T * A * U is the Rayleigh quotient. */ /* BLAS */ } else { /* A * U(:,1:K) is not explicitly needed and the */ /* computation is organized differently. The Rayleigh */ /* quotient is computed more efficiently. */ dgemm_("T", "N", k, n, m, &one, &x[x_offset], ldx, &y[y_offset], ldy, &zero, &z__[z_offset], ldz); /* Z(1:K,1:N) = MATMUL( TRANSPOSE(X(1:M,1:K)), Y(1:M,1:N) ) ! IN */ /* In the two DGEMM calls here, can use K for LDZ. */ /* B */ dgemm_("N", t_or_n__, k, k, n, &one, &z__[z_offset], ldz, &w[w_offset] , ldw, &zero, &s[s_offset], lds); /* S(1:K,1:K) = MATMUL(Z(1:K,1:N),TRANSPOSE(W(1:K,1:N))) ! INTRIN */ /* S(1:K,1:K) = MATMUL(Z(1:K,1:N),(W(1:N,1:K))) ! INTRIN */ /* At this point S = U^T * A * U is the Rayleigh quotient. */ /* If the residuals are requested, save scaled V_k into Z. */ /* Recall that V_k or V_k^T is stored in W. */ /* BLAS */ if (wntres || wntex) { if (lsame_(t_or_n__, "N")) { dlacpy_("A", n, k, &w[w_offset], ldw, &z__[z_offset], ldz); } else { dlacpy_("A", k, n, &w[w_offset], ldw, &z__[z_offset], ldz); } } } /* <5> Compute the Ritz values and (if requested) the */ /* right eigenvectors of the Rayleigh quotient. */ i__1 = *lwork - *n; dgeev_("N", jobzl, k, &s[s_offset], lds, &reig[1], &imeig[1], &w[w_offset] , ldw, &w[w_offset], ldw, &work[*n + 1], &i__1, &info1); /* W(1:K,1:K) contains the eigenvectors of the Rayleigh */ /* quotient. Even in the case of complex spectrum, all */ /* computation is done in real arithmetic. REIG and */ /* IMEIG are the real and the imaginary parts of the */ /* eigenvalues, so that the spectrum is given as */ /* REIG(:) + sqrt(-1)*IMEIG(:). Complex conjugate pairs */ /* are listed at consecutive positions. For such a */ /* complex conjugate pair of the eigenvalues, the */ /* corresponding eigenvectors are also a complex */ /* conjugate pair with the real and imaginary parts */ /* stored column-wise in W at the corresponding */ /* consecutive column indices. See the description of Z. */ /* Also, see the description of DGEEV. */ /* LAPACK C */ if (info1 > 0) { /* DGEEV failed to compute the eigenvalues and */ /* eigenvectors of the Rayleigh quotient. */ *info = 3; return 0; } /* <6> Compute the eigenvectors (if requested) and, */ /* the residuals (if requested). */ if (wntvec || wntex) { if (wntres) { if (wntref) { /* Here, if the refinement is requested, we have */ /* A*U(:,1:K) already computed and stored in Z. */ /* For the residuals, need Y = A * U(:,1;K) * W. */ dgemm_("N", "N", m, k, k, &one, &z__[z_offset], ldz, &w[ w_offset], ldw, &zero, &y[y_offset], ldy); /* Y(1:M,1:K) = Z(1:M,1:K) * W(1:K,1:K) ! INTRINSIC */ /* This frees Z; Y contains A * U(:,1:K) * W. */ /* BLAS CALL */ } else { /* Compute S = V_k * Sigma_k^(-1) * W, where */ /* V_k * Sigma_k^(-1) is stored in Z */ dgemm_(t_or_n__, "N", n, k, k, &one, &z__[z_offset], ldz, &w[ w_offset], ldw, &zero, &s[s_offset], lds); /* Then, compute Z = Y * S = */ /* = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) = */ /* = A * U(:,1:K) * W(1:K,1:K) */ dgemm_("N", "N", m, k, n, &one, &y[y_offset], ldy, &s[ s_offset], lds, &zero, &z__[z_offset], ldz); /* Save a copy of Z into Y and free Z for holding */ /* the Ritz vectors. */ dlacpy_("A", m, k, &z__[z_offset], ldz, &y[y_offset], ldy); if (wntex) { dlacpy_("A", m, k, &z__[z_offset], ldz, &b[b_offset], ldb); } } } else if (wntex) { /* Compute S = V_k * Sigma_k^(-1) * W, where */ /* V_k * Sigma_k^(-1) is stored in Z */ dgemm_(t_or_n__, "N", n, k, k, &one, &z__[z_offset], ldz, &w[ w_offset], ldw, &zero, &s[s_offset], lds); /* Then, compute Z = Y * S = */ /* = Y * V_k * Sigma_k^(-1) * W(1:K,1:K) = */ /* = A * U(:,1:K) * W(1:K,1:K) */ dgemm_("N", "N", m, k, n, &one, &y[y_offset], ldy, &s[s_offset], lds, &zero, &b[b_offset], ldb); /* The above call replaces the following two calls */ /* that were used in the developing-testing phase. */ /* CALL DGEMM( 'N', 'N', M, K, N, ONE, Y, LDY, S, & */ /* LDS, ZERO, Z, LDZ) */ /* Save a copy of Z into B and free Z for holding */ /* the Ritz vectors. */ /* CALL DLACPY( 'A', M, K, Z, LDZ, B, LDB ) */ } /* Compute the real form of the Ritz vectors */ if (wntvec) { dgemm_("N", "N", m, k, k, &one, &x[x_offset], ldx, &w[w_offset], ldw, &zero, &z__[z_offset], ldz); } /* Z(1:M,1:K) = MATMUL(X(1:M,1:K), W(1:K,1:K)) ! INTRINSIC */ /* BLAS CALL */ if (wntres) { i__ = 1; while(i__ <= *k) { if (imeig[i__] == zero) { /* have a real eigenvalue with real eigenvector */ d__1 = -reig[i__]; daxpy_(m, &d__1, &z__[i__ * z_dim1 + 1], &c__1, &y[i__ * y_dim1 + 1], &c__1); /* Y(1:M,i) = Y(1:M,i) - REIG(i) * Z(1:M,i) ! */ res[i__] = dnrm2_(m, &y[i__ * y_dim1 + 1], &c__1); ++i__; } else { /* Have a complex conjugate pair */ /* REIG(i) +- sqrt(-1)*IMEIG(i). */ /* Since all computation is done in real */ /* arithmetic, the formula for the residual */ /* is recast for real representation of the */ /* complex conjugate eigenpair. See the */ /* description of RES. */ ab[0] = reig[i__]; ab[1] = -imeig[i__]; ab[2] = imeig[i__]; ab[3] = reig[i__]; d__1 = -one; dgemm_("N", "N", m, &c__2, &c__2, &d__1, &z__[i__ * z_dim1 + 1], ldz, ab, &c__2, &one, &y[i__ * y_dim1 + 1], ldy); /* Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB ! INT */ /* BL */ res[i__] = dlange_("F", m, &c__2, &y[i__ * y_dim1 + 1], ldy, &work[*n + 1]); /* LA */ res[i__ + 1] = res[i__]; i__ += 2; } } } } if (*whtsvd == 4) { work[*n + 1] = xscl1; work[*n + 2] = xscl2; } /* Successful exit. */ if (! badxy) { *info = 0; } else { /* A warning on possible data inconsistency. */ /* This should be a rare event. */ *info = 4; } /* ............................................................ */ return 0; /* ...... */ } /* dgedmd_ */