#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 number of rows of F). */ /* ..... */ /* N (input) INTEGER, 0 <= N <= M */ /* The number of data snapshots from a single trajectory, */ /* taken at equidistant discrete times. This is the */ /* number of columns of F. */ /* ..... */ /* F (input/output) COMPLEX(KIND=WP) M-by-N array */ /* > On entry, */ /* the columns of F are the sequence of data snapshots */ /* from a single trajectory, taken at equidistant discrete */ /* times. It is assumed that the column norms of F are */ /* in the range of the normalized floating point numbers. */ /* < On exit, */ /* If JOBQ == 'Q', the array F contains the orthogonal */ /* matrix/factor of the QR factorization of the initial */ /* data snapshots matrix F. See the description of JOBQ. */ /* If JOBQ == 'N', the entries in F strictly below the main */ /* diagonal contain, column-wise, the information on the */ /* Householder vectors, as returned by CGEQRF. The */ /* remaining information to restore the orthogonal matrix */ /* of the initial QR factorization is stored in ZWORK(1:MIN(M,N)). */ /* See the description of ZWORK. */ /* ..... */ /* LDF (input) INTEGER, LDF >= M */ /* The leading dimension of the array F. */ /* ..... */ /* X (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N-1) array */ /* X is used as workspace to hold representations of the */ /* leading N-1 snapshots in the orthonormal basis computed */ /* in the QR factorization of F. */ /* On exit, the leading K columns of X contain the leading */ /* K left singular vectors of the above described content */ /* of X. To lift them to the space of the left singular */ /* vectors U(:,1:K) of the input data, pre-multiply with the */ /* Q factor from the initial QR factorization. */ /* See the descriptions of F, K, V and Z. */ /* ..... */ /* LDX (input) INTEGER, LDX >= N */ /* The leading dimension of the array X. */ /* ..... */ /* Y (workspace/output) COMPLEX(KIND=WP) MIN(M,N)-by-(N) array */ /* Y is used as workspace to hold representations of the */ /* trailing N-1 snapshots in the orthonormal basis computed */ /* in the QR factorization of F. */ /* On exit, */ /* If JOBT == 'R', Y contains the MIN(M,N)-by-N upper */ /* triangular factor from the QR factorization of the data */ /* snapshot matrix F. */ /* ..... */ /* LDY (input) INTEGER , LDY >= N */ /* 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-1 :: 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 description of 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 SVD/POD basis for the leading N-1 */ /* data snapshots (columns of F) 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. */ /* ..... */ /* EIGS (output) COMPLEX(KIND=WP) (N-1)-by-1 array */ /* The leading K (K<=N-1) entries of EIGS contain */ /* the computed eigenvalues (Ritz values). */ /* See the descriptions of K, and Z. */ /* ..... */ /* Z (workspace/output) COMPLEX(KIND=WP) M-by-(N-1) array */ /* If JOBZ =='V' then Z contains the Ritz vectors. Z(:,i) */ /* is an eigenvector of the i-th Ritz value; ||Z(:,i)||_2=1. */ /* If JOBZ == 'F', then the Z(:,i)'s are given implicitly as */ /* Z*V, where Z contains orthonormal matrix (the product of */ /* Q from the initial QR factorization and the SVD/POD_basis */ /* returned by CGEDMD in X) and the second factor (the */ /* eigenvectors of the Rayleigh quotient) is in the array V, */ /* as returned by CGEDMD. That is, X(:,1:K)*V(:,i) */ /* is an eigenvector corresponding to EIGS(i). The columns */ /* of V(1:K,1:K) are the computed eigenvectors of the */ /* K-by-K Rayleigh quotient. */ /* See the descriptions of EIGS, X and V. */ /* ..... */ /* LDZ (input) INTEGER , LDZ >= M */ /* The leading dimension of the array Z. */ /* ..... */ /* RES (output) REAL(KIND=WP) (N-1)-by-1 array */ /* RES(1:K) contains the residuals for the K computed */ /* Ritz pairs, */ /* RES(i) = || A * Z(:,i) - EIGS(i)*Z(:,i))||_2. */ /* See the description of EIGS and Z. */ /* ..... */ /* B (output) COMPLEX(KIND=WP) MIN(M,N)-by-(N-1) array. */ /* IF JOBF =='R', B(1:N,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:N,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. */ /* In both cases, the content of B can be lifted to the */ /* original dimension of the input data by pre-multiplying */ /* with the Q factor from the initial QR factorization. */ /* Here A denotes a compression of the underlying operator. */ /* See the descriptions of F and X. */ /* If JOBF =='N', then B is not referenced. */ /* ..... */ /* LDB (input) INTEGER, LDB >= MIN(M,N) */ /* The leading dimension of the array B. */ /* ..... */ /* V (workspace/output) COMPLEX(KIND=WP) (N-1)-by-(N-1) array */ /* On exit, V(1:K,1:K) V contains the K eigenvectors of */ /* the Rayleigh quotient. The Ritz vectors */ /* (returned in Z) are the product of Q from the initial QR */ /* factorization (see the description of F) X (see the */ /* description of X) and V. */ /* ..... */ /* LDV (input) INTEGER, LDV >= N-1 */ /* The leading dimension of the array V. */ /* ..... */ /* S (output) COMPLEX(KIND=WP) (N-1)-by-(N-1) array */ /* The array S(1:K,1:K) is used for the matrix Rayleigh */ /* quotient. This content is overwritten during */ /* the eigenvalue decomposition by CGEEV. */ /* See the description of K. */ /* ..... */ /* LDS (input) INTEGER, LDS >= N-1 */ /* The leading dimension of the array S. */ /* ..... */ /* ZWORK (workspace/output) COMPLEX(KIND=WP) LWORK-by-1 array */ /* On exit, */ /* ZWORK(1:MIN(M,N)) contains the scalar factors of the */ /* elementary reflectors as returned by CGEQRF of the */ /* M-by-N input matrix F. */ /* If the call to CGEDMDQ is only workspace query, then */ /* ZWORK(1) contains the minimal complex workspace length and */ /* ZWORK(2) is the optimal complex workspace length. */ /* Hence, the length of work is at least 2. */ /* See the description of LZWORK. */ /* ..... */ /* LZWORK (input) INTEGER */ /* The minimal length of the workspace vector ZWORK. */ /* LZWORK is calculated as follows: */ /* Let MLWQR = N (minimal workspace for CGEQRF[M,N]) */ /* MLWDMD = minimal workspace for CGEDMD (see the */ /* description of LWORK in CGEDMD) */ /* MLWMQR = N (minimal workspace for */ /* ZUNMQR['L','N',M,N,N]) */ /* MLWGQR = N (minimal workspace for ZUNGQR[M,N,N]) */ /* MINMN = MIN(M,N) */ /* Then */ /* LZWORK = MAX(2, MIN(M,N)+MLWQR, MINMN+MLWDMD) */ /* is further updated as follows: */ /* if JOBZ == 'V' or JOBZ == 'F' THEN */ /* LZWORK = MAX( LZWORK, MINMN+MLWMQR ) */ /* if JOBQ == 'Q' THEN */ /* LZWORK = MAX( ZLWORK, MINMN+MLWGQR) */ /* ..... */ /* WORK (workspace/output) REAL(KIND=WP) LWORK-by-1 array */ /* On exit, */ /* WORK(1:N-1) contains the singular values of */ /* the input submatrix F(1:M,1:N-1). */ /* If the call to CGEDMDQ is only workspace query, then */ /* WORK(1) contains the minimal workspace length and */ /* WORK(2) is the optimal workspace length. hence, the */ /* length of work is at least 2. */ /* See the description of LWORK. */ /* ..... */ /* LWORK (input) INTEGER */ /* The minimal length of the workspace vector WORK. */ /* LWORK is the same as in CGEDMD, because in CGEDMDQ */ /* only CGEDMD requires real workspace for snapshots */ /* of dimensions MIN(M,N)-by-(N-1). */ /* 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 */ /* Let M1=MIN(M,N), N1=N-1. Then */ /* 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 */ /* ~~~~~~~~~~ */ /* COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP ) */ /* Local scalars */ /* ~~~~~~~~~~~~~ */ /* External functions (BLAS and LAPACK) */ /* ~~~~~~~~~~~~~~~~~ */ /* External subroutines (BLAS and LAPACK) */ /* ~~~~~~~~~~~~~~~~~~~~ */ /* External subroutines */ /* ~~~~~~~~~~~~~~~~~~~~ */ /* Intrinsic functions */ /* ~~~~~~~~~~~~~~~~~~~ */ /* .......................................................... */ /* Parameter adjustments */ f_dim1 = *ldf; f_offset = 1 + f_dim1 * 1; f -= f_offset; x_dim1 = *ldx; x_offset = 1 + x_dim1 * 1; x -= x_offset; y_dim1 = *ldy; y_offset = 1 + y_dim1 * 1; y -= y_offset; --eigs; 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; v_dim1 = *ldv; v_offset = 1 + v_dim1 * 1; v -= v_offset; s_dim1 = *lds; s_offset = 1 + s_dim1 * 1; s -= s_offset; --zwork; --work; --iwork; /* Function Body */ one = 1.f; zero = 0.f; zzero.r = 0.f, zzero.i = 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"); wntvcf = lsame_(jobz, "F"); wntvcq = lsame_(jobz, "Q"); wntref = lsame_(jobf, "R"); wntex = lsame_(jobf, "E"); wantq = lsame_(jobq, "Q"); wnttrf = lsame_(jobt, "R"); minmn = f2cmin(*m,*n); *info = 0; lquery = *lwork == -1 || *liwork == -1; if (! (sccolx || sccoly || lsame_(jobs, "N"))) { *info = -1; } else if (! (wntvec || wntvcf || wntvcq || lsame_(jobz, "N"))) { *info = -2; } else if (! (wntres || lsame_(jobr, "N")) || wntres && lsame_(jobz, "N")) { *info = -3; } else if (! (wantq || lsame_(jobq, "N"))) { *info = -4; } else if (! (wnttrf || lsame_(jobt, "N"))) { *info = -5; } else if (! (wntref || wntex || lsame_(jobf, "N"))) { *info = -6; } else if (! (*whtsvd == 1 || *whtsvd == 2 || *whtsvd == 3 || *whtsvd == 4)) { *info = -7; } else if (*m < 0) { *info = -8; } else if (*n < 0 || *n > *m + 1) { *info = -9; } else if (*ldf < *m) { *info = -11; } else if (*ldx < minmn) { *info = -13; } else if (*ldy < minmn) { *info = -15; } else if (! (*nrnk == -2 || *nrnk == -1 || *nrnk >= 1 && *nrnk <= *n)) { *info = -16; } else if (*tol < zero || *tol >= one) { *info = -17; } else if (*ldz < *m) { *info = -21; } else if ((wntref || wntex) && *ldb < minmn) { *info = -24; } else if (*ldv < *n - 1) { *info = -26; } else if (*lds < *n - 1) { *info = -28; } if (wntvec || wntvcf || wntvcq) { *(unsigned char *)jobvl = 'V'; } else { *(unsigned char *)jobvl = 'N'; } 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 || *n == 1) { /* All output except K is void. INFO=1 signals */ /* the void input. In case of a workspace query, */ /* the minimal workspace lengths are returned. */ if (lquery) { iwork[1] = 1; work[1] = 2.f; work[2] = 2.f; } else { *k = 0; } *info = 1; return 0; } mlrwrk = 2; mlwork = 2; olwork = 2; iminwr = 1; mlwqr = f2cmax(1,*n); /* Minimal workspace length for CGEQRF. */ /* Computing MAX */ i__1 = mlwork, i__2 = minmn + mlwqr; mlwork = f2cmax(i__1,i__2); if (lquery) { cgeqrf_(m, n, &f[f_offset], ldf, &zwork[1], &zwork[1], &c_n1, & info1); olwqr = (integer) zwork[1].r; /* Computing MAX */ i__1 = olwork, i__2 = minmn + olwqr; olwork = f2cmax(i__1,i__2); } i__1 = *n - 1; cgedmd_(jobs, jobvl, jobr, jobf, whtsvd, &minmn, &i__1, &x[x_offset], ldx, &y[y_offset], ldy, nrnk, tol, k, &eigs[1], &z__[z_offset] , ldz, &res[1], &b[b_offset], ldb, &v[v_offset], ldv, &s[ s_offset], lds, &zwork[1], lzwork, &work[1], &c_n1, &iwork[1], liwork, &info1); mlwdmd = (integer) zwork[1].r; /* Computing MAX */ i__1 = mlwork, i__2 = minmn + mlwdmd; mlwork = f2cmax(i__1,i__2); /* Computing MAX */ i__1 = mlrwrk, i__2 = (integer) work[1]; mlrwrk = f2cmax(i__1,i__2); iminwr = f2cmax(iminwr,iwork[1]); if (lquery) { olwdmd = (integer) zwork[2].r; /* Computing MAX */ i__1 = olwork, i__2 = minmn + olwdmd; olwork = f2cmax(i__1,i__2); } if (wntvec || wntvcf) { mlwmqr = f2cmax(1,*n); /* Computing MAX */ i__1 = mlwork, i__2 = minmn + mlwmqr; mlwork = f2cmax(i__1,i__2); if (lquery) { cunmqr_("L", "N", m, n, &minmn, &f[f_offset], ldf, &zwork[1], &z__[z_offset], ldz, &zwork[1], &c_n1, &info1); olwmqr = (integer) zwork[1].r; /* Computing MAX */ i__1 = olwork, i__2 = minmn + olwmqr; olwork = f2cmax(i__1,i__2); } } if (wantq) { mlwgqr = f2cmax(1,*n); /* Computing MAX */ i__1 = mlwork, i__2 = minmn + mlwgqr; mlwork = f2cmax(i__1,i__2); if (lquery) { cungqr_(m, &minmn, &minmn, &f[f_offset], ldf, &zwork[1], & zwork[1], &c_n1, &info1); olwgqr = (integer) zwork[1].r; /* Computing MAX */ i__1 = olwork, i__2 = minmn + olwgqr; olwork = f2cmax(i__1,i__2); } } if (*liwork < iminwr && ! lquery) { *info = -34; } if (*lwork < mlrwrk && ! lquery) { *info = -32; } if (*lzwork < mlwork && ! lquery) { *info = -30; } } if (*info != 0) { i__1 = -(*info); xerbla_("CGEDMDQ", &i__1); return 0; } else if (lquery) { /* Return minimal and optimal workspace sizes */ iwork[1] = iminwr; zwork[1].r = (real) mlwork, zwork[1].i = 0.f; zwork[2].r = (real) olwork, zwork[2].i = 0.f; work[1] = (real) mlrwrk; work[2] = (real) mlrwrk; return 0; } /* ..... */ /* Initial QR factorization that is used to represent the */ /* snapshots as elements of lower dimensional subspace. */ /* For large scale computation with M >>N , at this place */ /* one can use an out of core QRF. */ i__1 = *lzwork - minmn; cgeqrf_(m, n, &f[f_offset], ldf, &zwork[1], &zwork[minmn + 1], &i__1, & info1); /* Define X and Y as the snapshots representations in the */ /* orthogonal basis computed in the QR factorization. */ /* X corresponds to the leading N-1 and Y to the trailing */ /* N-1 snapshots. */ i__1 = *n - 1; claset_("L", &minmn, &i__1, &zzero, &zzero, &x[x_offset], ldx); i__1 = *n - 1; clacpy_("U", &minmn, &i__1, &f[f_offset], ldf, &x[x_offset], ldx); i__1 = *n - 1; clacpy_("A", &minmn, &i__1, &f[(f_dim1 << 1) + 1], ldf, &y[y_offset], ldy); if (*m >= 3) { i__1 = minmn - 2; i__2 = *n - 2; claset_("L", &i__1, &i__2, &zzero, &zzero, &y[y_dim1 + 3], ldy); } /* Compute the DMD of the projected snapshot pairs (X,Y) */ i__1 = *n - 1; i__2 = *lzwork - minmn; cgedmd_(jobs, jobvl, jobr, jobf, whtsvd, &minmn, &i__1, &x[x_offset], ldx, &y[y_offset], ldy, nrnk, tol, k, &eigs[1], &z__[z_offset], ldz, & res[1], &b[b_offset], ldb, &v[v_offset], ldv, &s[s_offset], lds, & zwork[minmn + 1], &i__2, &work[1], lwork, &iwork[1], liwork, & info1); if (info1 == 2 || info1 == 3) { /* Return with error code. See CGEDMD for details. */ *info = info1; return 0; } else { *info = info1; } /* The Ritz vectors (Koopman modes) can be explicitly */ /* formed or returned in factored form. */ if (wntvec) { /* Compute the eigenvectors explicitly. */ if (*m > minmn) { i__1 = *m - minmn; claset_("A", &i__1, k, &zzero, &zzero, &z__[minmn + 1 + z_dim1], ldz); } i__1 = *lzwork - minmn; cunmqr_("L", "N", m, k, &minmn, &f[f_offset], ldf, &zwork[1], &z__[ z_offset], ldz, &zwork[minmn + 1], &i__1, &info1); } else if (wntvcf) { /* Return the Ritz vectors (eigenvectors) in factored */ /* form Z*V, where Z contains orthonormal matrix (the */ /* product of Q from the initial QR factorization and */ /* the SVD/POD_basis returned by CGEDMD in X) and the */ /* second factor (the eigenvectors of the Rayleigh */ /* quotient) is in the array V, as returned by CGEDMD. */ clacpy_("A", n, k, &x[x_offset], ldx, &z__[z_offset], ldz); if (*m > *n) { i__1 = *m - *n; claset_("A", &i__1, k, &zzero, &zzero, &z__[*n + 1 + z_dim1], ldz); } i__1 = *lzwork - minmn; cunmqr_("L", "N", m, k, &minmn, &f[f_offset], ldf, &zwork[1], &z__[ z_offset], ldz, &zwork[minmn + 1], &i__1, &info1); } /* Some optional output variables: */ /* The upper triangular factor R in the initial QR */ /* factorization is optionally returned in the array Y. */ /* This is useful if this call to CGEDMDQ is to be */ /* followed by a streaming DMD that is implemented in a */ /* QR compressed form. */ if (wnttrf) { /* Return the upper triangular R in Y */ claset_("A", &minmn, n, &zzero, &zzero, &y[y_offset], ldy); clacpy_("U", &minmn, n, &f[f_offset], ldf, &y[y_offset], ldy); } /* The orthonormal/unitary factor Q in the initial QR */ /* factorization is optionally returned in the array F. */ /* Same as with the triangular factor above, this is */ /* useful in a streaming DMD. */ if (wantq) { /* Q overwrites F */ i__1 = *lzwork - minmn; cungqr_(m, &minmn, &minmn, &f[f_offset], ldf, &zwork[1], &zwork[minmn + 1], &i__1, &info1); } return 0; } /* cgedmdq_ */