diff --git a/Makefile b/Makefile index ae1da5e9fe..3aa02ef459 100644 --- a/Makefile +++ b/Makefile @@ -406,8 +406,10 @@ endif lapack-test : (cd $(NETLIB_LAPACK_DIR)/TESTING && rm -f x* *.out) - $(MAKE) -j 1 -C $(NETLIB_LAPACK_DIR)/TESTING/EIG xeigtstc xeigtstd xeigtsts xeigtstz + $(MAKE) -j 1 -C $(NETLIB_LAPACK_DIR)/TESTING/EIG xeigtstc xeigtstd xeigtsts xeigtstz test_cgges3_scaling test_sgges3_scaling test_dgges3_scaling test_zgges3_scaling + $(MAKE) -j 1 -C $(NETLIB_LAPACK_DIR)/TESTING/LIN xlintstc xlintstd xlintstds xlintstrfd xlintstrfz xlintsts xlintstz xlintstzc xlintstrfs xlintstrfc + ifneq ($(CROSS), 1) ( cd $(NETLIB_LAPACK_DIR)/INSTALL; $(MAKE) all; ./testlsame; ./testslamch; ./testdlamch; \ ./testsecond; ./testdsecnd; ./testieee; ./testversion ) diff --git a/lapack-netlib/SRC/claqz0.c b/lapack-netlib/SRC/claqz0.c index 7d73d5f3ce..48f66ab890 100644 --- a/lapack-netlib/SRC/claqz0.c +++ b/lapack-netlib/SRC/claqz0.c @@ -252,271 +252,6 @@ typedef struct Namelist Namelist; //#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++ */ - -#define F2C_proc_par_types 1 -#ifdef __cplusplus -typedef logical (*L_fp)(...); -#else -typedef logical (*L_fp)(); -#endif -#if 0 -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= istart2) { if (c_abs(&b[k + k * b_dim1]) < btol) { + i__2 = k + k * b_dim1; + b[i__2].r = 0.f, b[i__2].i = 0.f; /* A diagonal element of B is negligible, move it */ /* to the top and deflate it */ i__2 = istart2 + 1; diff --git a/lapack-netlib/SRC/claqz0.f b/lapack-netlib/SRC/claqz0.f index 3433b4c8a5..680f787c7b 100644 --- a/lapack-netlib/SRC/claqz0.f +++ b/lapack-netlib/SRC/claqz0.f @@ -534,6 +534,7 @@ RECURSIVE SUBROUTINE CLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, DO WHILE ( K.GE.ISTART2 ) IF( ABS( B( K, K ) ) .LT. BTOL ) THEN + B( K, K ) = CZERO * A diagonal element of B is negligible, move it * to the top and deflate it diff --git a/lapack-netlib/SRC/dlaqz0.c b/lapack-netlib/SRC/dlaqz0.c index e001a77d6b..db5d65e499 100644 --- a/lapack-netlib/SRC/dlaqz0.c +++ b/lapack-netlib/SRC/dlaqz0.c @@ -226,21 +226,6 @@ typedef struct Namelist Namelist; */ - -/* -- translated by f2c (version 20200916). - You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip -*/ - - - /* Table of constant values */ static integer c__12 = 12; @@ -603,6 +588,7 @@ static integer c__1 = 1; k = istop; while(k >= istart2) { if ((d__1 = b[k + k * b_dim1], abs(d__1)) < btol) { + b[k + k * b_dim1] = 0.; /* A diagonal element of B is negligible, move it */ /* to the top and deflate it */ i__2 = istart2 + 1; diff --git a/lapack-netlib/SRC/dlaqz0.f b/lapack-netlib/SRC/dlaqz0.f index f91ac447b2..5b62839768 100644 --- a/lapack-netlib/SRC/dlaqz0.f +++ b/lapack-netlib/SRC/dlaqz0.f @@ -569,6 +569,7 @@ RECURSIVE SUBROUTINE DLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, DO WHILE ( K.GE.ISTART2 ) IF( ABS( B( K, K ) ) .LT. BTOL ) THEN + B( K, K ) = ZERO * A diagonal element of B is negligible, move it * to the top and deflate it diff --git a/lapack-netlib/SRC/slaqz0.c b/lapack-netlib/SRC/slaqz0.c index 7ed8c38eb1..b4fdefcd67 100644 --- a/lapack-netlib/SRC/slaqz0.c +++ b/lapack-netlib/SRC/slaqz0.c @@ -220,271 +220,6 @@ typedef struct Namelist Namelist; //#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++ */ - -#define F2C_proc_par_types 1 -#ifdef __cplusplus -typedef logical (*L_fp)(...); -#else -typedef logical (*L_fp)(); -#endif -#if 0 -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= istart2) { if ((r__1 = b[k + k * b_dim1], abs(r__1)) < btol) { + b[k + k * b_dim1] = 0.f; /* A diagonal element of B is negligible, move it */ /* to the top and deflate it */ i__2 = istart2 + 1; diff --git a/lapack-netlib/SRC/slaqz0.f b/lapack-netlib/SRC/slaqz0.f index 504bf488eb..17cd8bb884 100644 --- a/lapack-netlib/SRC/slaqz0.f +++ b/lapack-netlib/SRC/slaqz0.f @@ -567,6 +567,7 @@ RECURSIVE SUBROUTINE SLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, DO WHILE ( K.GE.ISTART2 ) IF( ABS( B( K, K ) ) .LT. BTOL ) THEN + B( K, K ) = ZERO * A diagonal element of B is negligible, move it * to the top and deflate it diff --git a/lapack-netlib/SRC/zlaqz0.c b/lapack-netlib/SRC/zlaqz0.c index b8a9f994f7..3081e59478 100644 --- a/lapack-netlib/SRC/zlaqz0.c +++ b/lapack-netlib/SRC/zlaqz0.c @@ -252,271 +252,6 @@ typedef struct Namelist Namelist; //#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++ */ - -#define F2C_proc_par_types 1 -#ifdef __cplusplus -typedef logical (*L_fp)(...); -#else -typedef logical (*L_fp)(); -#endif -#if 0 -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= istart2) { if (z_abs(&b[k + k * b_dim1]) < btol) { + i__2 = k + k * b_dim1; + b[i__2].r = 0., b[i__2].i = 0.; /* A diagonal element of B is negligible, move it */ /* to the top and deflate it */ i__2 = istart2 + 1; diff --git a/lapack-netlib/SRC/zlaqz0.f b/lapack-netlib/SRC/zlaqz0.f index c5f58294dc..beaea54f09 100644 --- a/lapack-netlib/SRC/zlaqz0.f +++ b/lapack-netlib/SRC/zlaqz0.f @@ -1,280 +1,3 @@ -*> \brief \b ZLAQZ0 -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> Download ZLAQZ0 + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -* -* Definition: -* =========== -* -* SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, -* $ LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, -* $ INFO ) -* IMPLICIT NONE -* -* Arguments -* CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ -* INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, -* $ REC -* INTEGER, INTENT( OUT ) :: INFO -* COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, -* $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * ), WORK( * ) -* DOUBLE PRECISION, INTENT( OUT ) :: RWORK( * ) -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> ZLAQZ0 computes the eigenvalues of a real matrix pair (H,T), -*> where H is an upper Hessenberg matrix and T is upper triangular, -*> using the double-shift QZ method. -*> Matrix pairs of this type are produced by the reduction to -*> generalized upper Hessenberg form of a real matrix pair (A,B): -*> -*> A = Q1*H*Z1**H, B = Q1*T*Z1**H, -*> -*> as computed by ZGGHRD. -*> -*> 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, P and S are an upper triangular -*> matrices. -*> -*> 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 ZGGHRD that reduced -*> the matrix pair (A,B) to generalized upper Hessenberg form, then the -*> output matrices Q1*Q and Z1*Z are the unitary factors from the -*> generalized Schur factorization of (A,B): -*> -*> A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H. -*> -*> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently, -*> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is -*> complex and beta real. -*> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the -*> generalized nonsymmetric eigenvalue problem (GNEP) -*> A*x = lambda*B*x -*> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the -*> alternate form of the GNEP -*> mu*A*y = B*y. -*> Eigenvalues can be read directly from the generalized Schur -*> form: -*> alpha = S(i,i), beta = P(i,i). -*> -*> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix -*> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), -*> pp. 241--256. -*> -*> Ref: B. Kagstrom, D. Kressner, "Multishift Variants of the QZ -*> Algorithm with Aggressive Early Deflation", SIAM J. Numer. -*> Anal., 29(2006), pp. 199--227. -*> -*> Ref: T. Steel, D. Camps, K. Meerbergen, R. Vandebril "A multishift, -*> multipole rational QZ method with aggressive early deflation" -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] WANTS -*> \verbatim -*> WANTS is CHARACTER*1 -*> = 'E': Compute eigenvalues only; -*> = 'S': Compute eigenvalues and the Schur form. -*> \endverbatim -*> -*> \param[in] WANTQ -*> \verbatim -*> WANTQ is CHARACTER*1 -*> = 'N': Left Schur vectors (Q) are not computed; -*> = 'I': Q is initialized to the unit matrix and the matrix Q -*> of left Schur vectors of (A,B) is returned; -*> = 'V': Q must contain an unitary matrix Q1 on entry and -*> the product Q1*Q is returned. -*> \endverbatim -*> -*> \param[in] WANTZ -*> \verbatim -*> WANTZ is CHARACTER*1 -*> = 'N': Right Schur vectors (Z) are not computed; -*> = 'I': Z is initialized to the unit matrix and the matrix Z -*> of right Schur vectors of (A,B) is returned; -*> = 'V': Z must contain an unitary matrix Z1 on entry and -*> the product Z1*Z is returned. -*> \endverbatim -*> -*> \param[in] N -*> \verbatim -*> N is INTEGER -*> The order of the matrices A, B, Q, and Z. N >= 0. -*> \endverbatim -*> -*> \param[in] ILO -*> \verbatim -*> ILO is INTEGER -*> \endverbatim -*> -*> \param[in] IHI -*> \verbatim -*> IHI is INTEGER -*> ILO and IHI mark the rows and columns of A which are in -*> Hessenberg form. It is assumed that A is already upper -*> triangular in rows and columns 1:ILO-1 and IHI+1:N. -*> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. -*> \endverbatim -*> -*> \param[in,out] A -*> \verbatim -*> A is COMPLEX*16 array, dimension (LDA, N) -*> On entry, the N-by-N upper Hessenberg matrix A. -*> On exit, if JOB = 'S', A contains the upper triangular -*> matrix S from the generalized Schur factorization. -*> If JOB = 'E', the diagonal blocks of A match those of S, but -*> the rest of A is unspecified. -*> \endverbatim -*> -*> \param[in] LDA -*> \verbatim -*> LDA is INTEGER -*> The leading dimension of the array A. LDA >= max( 1, N ). -*> \endverbatim -*> -*> \param[in,out] B -*> \verbatim -*> B is COMPLEX*16 array, dimension (LDB, N) -*> On entry, the N-by-N upper triangular matrix B. -*> On exit, if JOB = 'S', B contains the upper triangular -*> matrix P from the generalized Schur factorization; -*> If JOB = 'E', the diagonal blocks of B match those of P, but -*> the rest of B is unspecified. -*> \endverbatim -*> -*> \param[in] LDB -*> \verbatim -*> LDB is INTEGER -*> The leading dimension of the array B. LDB >= max( 1, N ). -*> \endverbatim -*> -*> \param[out] ALPHA -*> \verbatim -*> ALPHA is COMPLEX*16 array, dimension (N) -*> Each scalar alpha defining an eigenvalue -*> of GNEP. -*> \endverbatim -*> -*> \param[out] BETA -*> \verbatim -*> BETA is COMPLEX*16 array, dimension (N) -*> The scalars beta that define the eigenvalues of GNEP. -*> Together, the quantities alpha = ALPHA(j) and -*> beta = BETA(j) represent the j-th eigenvalue of the matrix -*> pair (A,B), in one of the forms lambda = alpha/beta or -*> mu = beta/alpha. Since either lambda or mu may overflow, -*> they should not, in general, be computed. -*> \endverbatim -*> -*> \param[in,out] Q -*> \verbatim -*> Q is COMPLEX*16 array, dimension (LDQ, N) -*> On entry, if COMPQ = 'V', the unitary matrix Q1 used in -*> the reduction of (A,B) to generalized Hessenberg form. -*> On exit, if COMPQ = 'I', the unitary matrix of left Schur -*> vectors of (A,B), and if COMPQ = 'V', the unitary matrix -*> of left Schur vectors of (A,B). -*> Not referenced if COMPQ = 'N'. -*> \endverbatim -*> -*> \param[in] LDQ -*> \verbatim -*> LDQ is INTEGER -*> The leading dimension of the array Q. LDQ >= 1. -*> If COMPQ='V' or 'I', then LDQ >= N. -*> \endverbatim -*> -*> \param[in,out] Z -*> \verbatim -*> Z is COMPLEX*16 array, dimension (LDZ, N) -*> On entry, if COMPZ = 'V', the unitary matrix Z1 used in -*> the reduction of (A,B) to generalized Hessenberg form. -*> On exit, if COMPZ = 'I', the unitary matrix of -*> right Schur vectors of (H,T), and if COMPZ = 'V', the -*> unitary matrix of right Schur vectors of (A,B). -*> Not referenced if COMPZ = 'N'. -*> \endverbatim -*> -*> \param[in] LDZ -*> \verbatim -*> LDZ is INTEGER -*> The leading dimension of the array Z. LDZ >= 1. -*> If COMPZ='V' or 'I', then LDZ >= N. -*> \endverbatim -*> -*> \param[out] WORK -*> \verbatim -*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) -*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. -*> \endverbatim -*> -*> \param[out] RWORK -*> \verbatim -*> RWORK is DOUBLE PRECISION array, dimension (N) -*> \endverbatim -*> -*> \param[in] LWORK -*> \verbatim -*> LWORK is INTEGER -*> The dimension of the array WORK. LWORK >= max(1,N). -*> -*> If LWORK = -1, then a workspace query is assumed; the routine -*> only calculates the optimal size of the WORK array, returns -*> this value as the first entry of the WORK array, and no error -*> message related to LWORK is issued by XERBLA. -*> \endverbatim -*> -*> \param[in] REC -*> \verbatim -*> REC is INTEGER -*> REC indicates the current recursion level. Should be set -*> to 0 on first call. -*> \endverbatim -*> -*> \param[out] INFO -*> \verbatim -*> INFO is INTEGER -*> = 0: successful exit -*> < 0: if INFO = -i, the i-th argument had an illegal value -*> = 1,...,N: the QZ iteration did not converge. (A,B) is not -*> in Schur form, but ALPHA(i) and -*> BETA(i), i=INFO+1,...,N should be correct. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Thijs Steel, KU Leuven -* -*> \date May 2020 -* -*> \ingroup laqz0 -*> -* ===================================================================== RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, $ A, $ LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, @@ -283,39 +6,42 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, IMPLICIT NONE * Arguments - CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ - INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, + CHARACTER WANTS, WANTQ, WANTZ + INTEGER N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, $ REC - INTEGER, INTENT( OUT ) :: INFO - COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, + INTEGER INFO + COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * ), WORK( * ) - DOUBLE PRECISION, INTENT( OUT ) :: RWORK( * ) + DOUBLE PRECISION RWORK( * ) * Parameters COMPLEX*16 CZERO, CONE PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), CONE = ( 1.0D+0, $ 0.0D+0 ) ) - DOUBLE PRECISION :: ZERO, ONE, HALF + DOUBLE PRECISION ZERO, ONE, HALF PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 ) * Local scalars - DOUBLE PRECISION :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR, + DOUBLE PRECISION SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR, $ BNORM, BTOL - COMPLEX*16 :: ESHIFT, S1, TEMP - INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS, + COMPLEX*16 ESHIFT, S1, TEMP + INTEGER ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS, $ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED, $ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM, $ ISTOPM, IWANTS, IWANTQ, IWANTZ, NORM_INFO, AED_INFO, $ NWR, NBR, NSR, ITEMP1, ITEMP2, RCOST - LOGICAL :: ILSCHUR, ILQ, ILZ - CHARACTER :: JBCMPZ*3 + LOGICAL ILSCHUR, ILQ, ILZ + CHARACTER JBCMPZ*3 * External Functions - EXTERNAL :: XERBLA, ZHGEQZ, ZLAQZ2, ZLAQZ3, ZLASET, + EXTERNAL XERBLA, ZHGEQZ, ZLAQZ2, ZLAQZ3, ZLASET, $ ZLARTG, ZROT - DOUBLE PRECISION, EXTERNAL :: DLAMCH, ZLANHS - LOGICAL, EXTERNAL :: LSAME - INTEGER, EXTERNAL :: ILAENV + DOUBLE PRECISION DLAMCH, ZLANHS + EXTERNAL DLAMCH, ZLANHS + LOGICAL LSAME + EXTERNAL LSAME + INTEGER ILAENV + EXTERNAL ILAENV * * Decode wantS,wantQ,wantZ @@ -535,6 +261,7 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, DO WHILE ( K.GE.ISTART2 ) IF( ABS( B( K, K ) ) .LT. BTOL ) THEN + B( K, K ) = CZERO * A diagonal element of B is negligible, move it * to the top and deflate it @@ -570,7 +297,7 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, $ K2 ), LDB, C1, S1 ) IF( ILQ ) THEN CALL ZROT( N, Q( 1, K2 ), 1, Q( 1, K2+1 ), 1, - $ C1, DCONJG( S1 ) ) + $ C1, CONJG( S1 ) ) END IF END IF @@ -590,7 +317,7 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, $ ISTART2+1 ), LDB, C1, S1 ) IF( ILQ ) THEN CALL ZROT( N, Q( 1, ISTART2 ), 1, Q( 1, - $ ISTART2+1 ), 1, C1, DCONJG( S1 ) ) + $ ISTART2+1 ), 1, C1, CONJG( S1 ) ) END IF END IF diff --git a/lapack-netlib/TESTING/CMakeLists.txt b/lapack-netlib/TESTING/CMakeLists.txt index d4e6f970d7..f9c0a965f3 100644 --- a/lapack-netlib/TESTING/CMakeLists.txt +++ b/lapack-netlib/TESTING/CMakeLists.txt @@ -57,6 +57,7 @@ add_lapack_test(slse.out lse.in xeigtsts) # # ======== SINGLE DMD EIG TESTS =========================== add_lapack_test(sdmd.out sdmd.in xdmdeigtsts) +add_test(NAME LAPACK-test_sgges3_scaling COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_sgges3_scaling) endif() if(BUILD_DOUBLE) @@ -91,6 +92,7 @@ add_lapack_test(dlse.out lse.in xeigtstd) # # ======== DOUBLE DMD EIG TESTS =========================== add_lapack_test(ddmd.out ddmd.in xdmdeigtstd) +add_test(NAME LAPACK-test_dgges3_scaling COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_dgges3_scaling) endif() if(BUILD_COMPLEX) @@ -123,6 +125,7 @@ add_lapack_test(clse.out lse.in xeigtstc) # # ======== COMPLEX DMD EIG TESTS =========================== add_lapack_test(cdmd.out cdmd.in xdmdeigtstc) +add_test(NAME LAPACK-test_cgges3_scaling COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_cgges3_scaling) endif() if(BUILD_COMPLEX16) @@ -157,6 +160,7 @@ add_lapack_test(zlse.out lse.in xeigtstz) # # ======== COMPLEX16 DMD EIG TESTS =========================== add_lapack_test(zdmd.out zdmd.in xdmdeigtstz) +add_test(NAME LAPACK-test_zgges3_scaling COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_zgges3_scaling) endif() diff --git a/lapack-netlib/TESTING/EIG/CMakeLists.txt b/lapack-netlib/TESTING/EIG/CMakeLists.txt index 6b7b0b7693..c203161689 100644 --- a/lapack-netlib/TESTING/EIG/CMakeLists.txt +++ b/lapack-netlib/TESTING/EIG/CMakeLists.txt @@ -112,21 +112,25 @@ endmacro() if(BUILD_SINGLE) add_eig_executable(xeigtsts ${SEIGTST} ${SCIGTST} ${AEIGTST}) add_eig_executable(xdmdeigtsts ${SDMDEIGTST} ${AEIGTST}) +add_eig_executable(test_sgges3_scaling test_sgges3_scaling.f90) endif() if(BUILD_COMPLEX) add_eig_executable(xeigtstc ${CEIGTST} ${SCIGTST} ${AEIGTST}) add_eig_executable(xdmdeigtstc ${CDMDEIGTST} ${AEIGTST}) +add_eig_executable(test_cgges3_scaling test_cgges3_scaling.f90) endif() if(BUILD_DOUBLE) add_eig_executable(xeigtstd ${DEIGTST} ${DZIGTST} ${AEIGTST}) add_eig_executable(xdmdeigtstd ${DDMDEIGTST} ${AEIGTST}) +add_eig_executable(test_dgges3_scaling test_dgges3_scaling.f90) endif() if(BUILD_COMPLEX16) add_eig_executable(xeigtstz ${ZEIGTST} ${DZIGTST} ${AEIGTST}) add_eig_executable(xdmdeigtstz ${ZDMDEIGTST} ${AEIGTST}) #add_eig_executable(test_wq_zrwork test_wq_rwork.f) +add_eig_executable(test_zgges3_scaling test_zgges3_scaling.f90) endif() diff --git a/lapack-netlib/TESTING/EIG/Makefile b/lapack-netlib/TESTING/EIG/Makefile index f685c7b0e5..73661f1e46 100644 --- a/lapack-netlib/TESTING/EIG/Makefile +++ b/lapack-netlib/TESTING/EIG/Makefile @@ -130,10 +130,10 @@ ZDMDEIGTST = zchkdmd.o all: single complex double complex16 .PHONY: single complex double complex16 -single: xeigtsts xdmdeigtsts -complex: xeigtstc xdmdeigtstc -double: xeigtstd xdmdeigtstd -complex16: xeigtstz xdmdeigtstz +single: xeigtsts xdmdeigtsts test_sgges3_scaling +complex: xeigtstc xdmdeigtstc test_cgges3_scaling +double: xeigtstd xdmdeigtstd test_dgges3_scaling +complex16: xeigtstz xdmdeigtstz test_zgges3_scaling xdmdeigtsts: $(SDMDEIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ @@ -159,6 +159,18 @@ xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +test_sgges3_scaling: test_sgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + +test_dgges3_scaling: test_dgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + +test_cgges3_scaling: test_cgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + +test_zgges3_scaling: test_zgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + $(SDMDEIGTST): $(FRC) $(CDMDEIGTST): $(FRC) $(DDMDEIGTST): $(FRC) diff --git a/lapack-netlib/TESTING/EIG/test_cgges3_scaling.f90 b/lapack-netlib/TESTING/EIG/test_cgges3_scaling.f90 new file mode 100644 index 0000000000..819a1f2612 --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_cgges3_scaling.f90 @@ -0,0 +1,69 @@ +! Check that a negligible nonzero pivot in CLAQZ0 does not leave fill +! below B's diagonal. Such fill is not rescaled by CGGES3's final +! CLASCL('U',...) call and corrupts the decomposition of a small pencil. +program test_cgges3_scaling + use, intrinsic :: iso_fortran_env, only: sp => real32 + implicit none + ! Pad the active 3-by-3 block to enable multishift QZ with the default + ! IPARMQ crossover NMIN=75, without overriding the tuning routines. + integer, parameter :: n=80, lwork=200000 + complex(sp) :: a(n,n), b(n,n), a0(n,n), b0(n,n), q(n,n), z(n,n) + complex(sp) :: alpha(n), beta(n), residual(n,n) + complex(sp), allocatable :: work(:) + real(sp) :: rwork(8*n), scale, anorm, bnorm, aerr, berr, lower, tol + real(sp), external :: slamch + logical :: bwork(n), select_none + external :: select_none, cgges3 + integer :: i, j, info, sdim + + ! Keep the large workspace off the stack, including recursive builds. + allocate(work(lwork)) + scale=1.0e-30_sp + a=0; b=0 + do i=1,n + a(i,i)=1; b(i,i)=1 + end do + ! A = [0 -1 -1; 1e-10 0 -1; 0 1 -1] + ! B = [1 1 1; 0 (1+2i)*1e-8 1; 0 0 1] + a(1,1)=0; a(1,2)=-1; a(1,3)=-1 + a(2,1)=1.0e-10_sp; a(2,2)=0; a(2,3)=-1 + a(3,2)=1; a(3,3)=-1 + b(1,2)=1; b(1,3)=1; b(2,3)=1 + b(2,2)=cmplx(1.0e-8_sp,2.0e-8_sp,sp) + a=a*scale; b=b*scale + a0=a; b0=b + + call cgges3('V','V','N',select_none,n,a,n,b,n,sdim,alpha,beta, & + q,n,z,n,work,lwork,rwork,bwork,info) + if(info/=0) then + print '(a,i0)', 'CGGES3 returned INFO = ',info + stop 1 + end if + + ! Normalize before reconstruction to avoid underflow in the checks. + anorm=maxval(abs(a0)); bnorm=maxval(abs(b0)) + residual=a0/anorm-matmul(q,matmul(a/anorm,conjg(transpose(z)))) + aerr=maxval(abs(residual)) + residual=b0/bnorm-matmul(q,matmul(b/bnorm,conjg(transpose(z)))) + berr=maxval(abs(residual)) + lower=0 + do j=1,n + do i=j+1,n + lower=max(lower,abs(b(i,j))) + end do + end do + tol=100*n*slamch('Precision') + print '(a,es24.16e3)', 'A relative max-entry residual = ',aerr + print '(a,es24.16e3)', 'B relative max-entry residual = ',berr + print '(a,es24.16e3)', 'B strict lower triangle maximum = ',lower + deallocate(work) + ! Negated comparisons also reject NaNs. + if(.not.(aerr<=tol .and. berr<=tol .and. lower==0.0_sp)) stop 1 +end program + +logical function select_none(alpha,beta) + use, intrinsic :: iso_fortran_env, only: sp => real32 + implicit none + complex(sp), intent(in) :: alpha,beta + select_none=.false. +end function diff --git a/lapack-netlib/TESTING/EIG/test_dgges3_scaling.f90 b/lapack-netlib/TESTING/EIG/test_dgges3_scaling.f90 new file mode 100644 index 0000000000..22cfc08bdb --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_dgges3_scaling.f90 @@ -0,0 +1,60 @@ +! Check that a negligible nonzero pivot in DLAQZ0 does not leave fill +! below B's diagonal before DGGES3 restores the original scale. +program test_dgges3_scaling + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + integer, parameter :: n=80, lwork=200000 + real(dp) :: a(n,n), b(n,n), a0(n,n), b0(n,n), q(n,n), z(n,n) + real(dp) :: alphar(n), alphai(n), beta(n), residual(n,n) + real(dp), allocatable :: work(:) + real(dp) :: scale, anorm, bnorm, aerr, berr, lower, tol + real(dp), external :: dlamch + logical :: bwork(n), select_none + external :: select_none, dgges3 + integer :: i, j, info, sdim + + allocate(work(lwork)) + scale=1.0e-280_dp + a=0; b=0 + do i=1,n + a(i,i)=1; b(i,i)=1 + end do + a(1,1)=0; a(1,2)=-1; a(1,3)=-1 + a(2,1)=1.0e-20_dp; a(2,2)=0; a(2,3)=-1 + a(3,2)=1; a(3,3)=-1 + b(1,2)=1; b(1,3)=1; b(2,3)=1 + b(2,2)=1.0e-18_dp + a=a*scale; b=b*scale + a0=a; b0=b + + call dgges3('V','V','N',select_none,n,a,n,b,n,sdim,alphar,alphai, & + beta,q,n,z,n,work,lwork,bwork,info) + if(info/=0) then + print '(a,i0)', 'DGGES3 returned INFO = ',info + stop 1 + end if + anorm=maxval(abs(a0)); bnorm=maxval(abs(b0)) + residual=a0/anorm-matmul(q,matmul(a/anorm,transpose(z))) + aerr=maxval(abs(residual)) + residual=b0/bnorm-matmul(q,matmul(b/bnorm,transpose(z))) + berr=maxval(abs(residual)) + lower=0 + do j=1,n + do i=j+1,n + lower=max(lower,abs(b(i,j))) + end do + end do + tol=100*n*dlamch('Precision') + print '(a,es24.16e3)', 'A relative max-entry residual = ',aerr + print '(a,es24.16e3)', 'B relative max-entry residual = ',berr + print '(a,es24.16e3)', 'B strict lower triangle maximum = ',lower + deallocate(work) + if(.not.(aerr<=tol .and. berr<=tol .and. lower==0.0_dp)) stop 1 +end program + +logical function select_none(alphar,alphai,beta) + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + real(dp), intent(in) :: alphar,alphai,beta + select_none=.false. +end function diff --git a/lapack-netlib/TESTING/EIG/test_sgges3_scaling.f90 b/lapack-netlib/TESTING/EIG/test_sgges3_scaling.f90 new file mode 100644 index 0000000000..5ec358f99c --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_sgges3_scaling.f90 @@ -0,0 +1,60 @@ +! Check that a negligible nonzero pivot in SLAQZ0 does not leave fill +! below B's diagonal before SGGES3 restores the original scale. +program test_sgges3_scaling + use, intrinsic :: iso_fortran_env, only: sp => real32 + implicit none + integer, parameter :: n=80, lwork=200000 + real(sp) :: a(n,n), b(n,n), a0(n,n), b0(n,n), q(n,n), z(n,n) + real(sp) :: alphar(n), alphai(n), beta(n), residual(n,n) + real(sp), allocatable :: work(:) + real(sp) :: scale, anorm, bnorm, aerr, berr, lower, tol + real(sp), external :: slamch + logical :: bwork(n), select_none + external :: select_none, sgges3 + integer :: i, j, info, sdim + + allocate(work(lwork)) + scale=1.0e-30_sp + a=0; b=0 + do i=1,n + a(i,i)=1; b(i,i)=1 + end do + a(1,1)=0; a(1,2)=-1; a(1,3)=-1 + a(2,1)=1.0e-10_sp; a(2,2)=0; a(2,3)=-1 + a(3,2)=1; a(3,3)=-1 + b(1,2)=1; b(1,3)=1; b(2,3)=1 + b(2,2)=1.0e-8_sp + a=a*scale; b=b*scale + a0=a; b0=b + + call sgges3('V','V','N',select_none,n,a,n,b,n,sdim,alphar,alphai, & + beta,q,n,z,n,work,lwork,bwork,info) + if(info/=0) then + print '(a,i0)', 'SGGES3 returned INFO = ',info + stop 1 + end if + anorm=maxval(abs(a0)); bnorm=maxval(abs(b0)) + residual=a0/anorm-matmul(q,matmul(a/anorm,transpose(z))) + aerr=maxval(abs(residual)) + residual=b0/bnorm-matmul(q,matmul(b/bnorm,transpose(z))) + berr=maxval(abs(residual)) + lower=0 + do j=1,n + do i=j+1,n + lower=max(lower,abs(b(i,j))) + end do + end do + tol=100*n*slamch('Precision') + print '(a,es24.16e3)', 'A relative max-entry residual = ',aerr + print '(a,es24.16e3)', 'B relative max-entry residual = ',berr + print '(a,es24.16e3)', 'B strict lower triangle maximum = ',lower + deallocate(work) + if(.not.(aerr<=tol .and. berr<=tol .and. lower==0.0_sp)) stop 1 +end program + +logical function select_none(alphar,alphai,beta) + use, intrinsic :: iso_fortran_env, only: sp => real32 + implicit none + real(sp), intent(in) :: alphar,alphai,beta + select_none=.false. +end function diff --git a/lapack-netlib/TESTING/EIG/test_zgges3_scaling.f90 b/lapack-netlib/TESTING/EIG/test_zgges3_scaling.f90 new file mode 100644 index 0000000000..15d98630ae --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_zgges3_scaling.f90 @@ -0,0 +1,69 @@ +! Check that a negligible nonzero pivot in ZLAQZ0 does not leave fill +! below B's diagonal. Such fill is not rescaled by ZGGES3's final +! ZLASCL('U',...) call and corrupts the decomposition of a small pencil. +program test_zgges3_scaling + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + ! Pad the active 3-by-3 block to enable multishift QZ with the default + ! IPARMQ crossover NMIN=75, without overriding the tuning routines. + integer, parameter :: n=80, lwork=200000 + complex(dp) :: a(n,n), b(n,n), a0(n,n), b0(n,n), q(n,n), z(n,n) + complex(dp) :: alpha(n), beta(n), residual(n,n) + complex(dp), allocatable :: work(:) + real(dp) :: rwork(8*n), scale, anorm, bnorm, aerr, berr, lower, tol + real(dp), external :: dlamch + logical :: bwork(n), select_none + external :: select_none, zgges3 + integer :: i, j, info, sdim + + ! Keep the large workspace off the stack, including recursive builds. + allocate(work(lwork)) + scale=1.0e-280_dp + a=0; b=0 + do i=1,n + a(i,i)=1; b(i,i)=1 + end do + ! A = [0 -1 -1; 1e-20 0 -1; 0 1 -1] + ! B = [1 1 1; 0 (1+2i)*1e-18 1; 0 0 1] + a(1,1)=0; a(1,2)=-1; a(1,3)=-1 + a(2,1)=1.0e-20_dp; a(2,2)=0; a(2,3)=-1 + a(3,2)=1; a(3,3)=-1 + b(1,2)=1; b(1,3)=1; b(2,3)=1 + b(2,2)=cmplx(1.0e-18_dp,2.0e-18_dp,dp) + a=a*scale; b=b*scale + a0=a; b0=b + + call zgges3('V','V','N',select_none,n,a,n,b,n,sdim,alpha,beta, & + q,n,z,n,work,lwork,rwork,bwork,info) + if(info/=0) then + print '(a,i0)', 'ZGGES3 returned INFO = ',info + stop 1 + end if + + ! Normalize before reconstruction to avoid underflow in the checks. + anorm=maxval(abs(a0)); bnorm=maxval(abs(b0)) + residual=a0/anorm-matmul(q,matmul(a/anorm,conjg(transpose(z)))) + aerr=maxval(abs(residual)) + residual=b0/bnorm-matmul(q,matmul(b/bnorm,conjg(transpose(z)))) + berr=maxval(abs(residual)) + lower=0 + do j=1,n + do i=j+1,n + lower=max(lower,abs(b(i,j))) + end do + end do + tol=100*n*dlamch('Precision') + print '(a,es24.16e3)', 'A relative max-entry residual = ',aerr + print '(a,es24.16e3)', 'B relative max-entry residual = ',berr + print '(a,es24.16e3)', 'B strict lower triangle maximum = ',lower + deallocate(work) + ! Negated comparisons also reject NaNs. + if(.not.(aerr<=tol .and. berr<=tol .and. lower==0.0_dp)) stop 1 +end program + +logical function select_none(alpha,beta) + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + complex(dp), intent(in) :: alpha,beta + select_none=.false. +end function diff --git a/lapack-netlib/TESTING/cec.in b/lapack-netlib/TESTING/cec.in index 622f05f1f1..ded23c8f83 100644 --- a/lapack-netlib/TESTING/cec.in +++ b/lapack-netlib/TESTING/cec.in @@ -1,5 +1,6 @@ CEC Key indicating type of input 20.0E0 Threshold value for test ratios +T Put T to test the error exits 1 1 ( 2.0E0, 0.0E0) ( 2.0E0, 0.0E0) diff --git a/lapack-netlib/TESTING/dec.in b/lapack-netlib/TESTING/dec.in index 4a2aa52a13..a40fdc4ef2 100644 --- a/lapack-netlib/TESTING/dec.in +++ b/lapack-netlib/TESTING/dec.in @@ -1,5 +1,6 @@ DEC Key indicating type of input 20.0 Threshold value for test ratios +T Put T to test the error exits 8 2 7 1.0D+00 1.0D+00 1.1D+00 1.3D+00 2.0D+00 3.0D+00 -4.7D+00 3.3D+00 -1.0D+00 1.0D+00 3.7D+00 7.9D+00 4.0D+00 5.3D+00 3.3D+00 -9.0D-01 @@ -1102,4 +1103,4 @@ DEC Key indicating type of input 0.000000000000000D+00 1.000000000000000D+00 4.839061798610604D-01 6.207808731846947D-01 0.000000000000000D+00 0.000000000000000D+00 1.000000000000000D+00 0.000000000000000D+00 0.000000000000000D+00 0.000000000000000D+00 0.000000000000000D+00 1.000000000000000D+00 - 0 0 0 \ No newline at end of file + 0 0 0 diff --git a/lapack-netlib/TESTING/sec.in b/lapack-netlib/TESTING/sec.in index 0b2f4e68c7..c1a72b3378 100644 --- a/lapack-netlib/TESTING/sec.in +++ b/lapack-netlib/TESTING/sec.in @@ -1,5 +1,6 @@ SEC Key indicating type of input 20.0 Threshold value for test ratios +T Put T to test the error exits 8 2 7 1.0E+00 1.0E+00 1.1E+00 1.3E+00 2.0E+00 3.0E+00 -4.7E+00 3.3E+00 -1.0E+00 1.0E+00 3.7E+00 7.9E+00 4.0E+00 5.3E+00 3.3E+00 -9.0E-01 diff --git a/lapack-netlib/TESTING/zec.in b/lapack-netlib/TESTING/zec.in index 32a1b890c9..e72ab512ef 100644 --- a/lapack-netlib/TESTING/zec.in +++ b/lapack-netlib/TESTING/zec.in @@ -1,5 +1,6 @@ ZEC Key indicating type of input 20.0D0 Threshold value for test ratios +T Put T to test the error exits 1 1 ( 2.0D0, 0.0D0) ( 2.0D0, 0.0D0)