diff --git a/Makefile b/Makefile index 3aa02ef459..a739f7fbab 100644 --- a/Makefile +++ b/Makefile @@ -406,7 +406,7 @@ 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 test_cgges3_scaling test_sgges3_scaling test_dgges3_scaling test_zgges3_scaling + $(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 test_cgesvj_emptsw test_dgesvj_emptsw test_sgesvj_emptsw test_zgesvj_emptsw $(MAKE) -j 1 -C $(NETLIB_LAPACK_DIR)/TESTING/LIN xlintstc xlintstd xlintstds xlintstrfd xlintstrfz xlintsts xlintstz xlintstzc xlintstrfs xlintstrfc diff --git a/lapack-netlib/SRC/cgesvj.c b/lapack-netlib/SRC/cgesvj.c index 0e00dc5689..94576db291 100644 --- a/lapack-netlib/SRC/cgesvj.c +++ b/lapack-netlib/SRC/cgesvj.c @@ -1087,7 +1087,7 @@ f"> */ /* Row-cyclic Jacobi SVD algorithm with column pivoting */ - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; i__1 = *n; for (q = 1; q <= i__1; ++q) { diff --git a/lapack-netlib/SRC/cgesvj.f b/lapack-netlib/SRC/cgesvj.f index cae7a244f4..f1fca7d3ab 100644 --- a/lapack-netlib/SRC/cgesvj.f +++ b/lapack-netlib/SRC/cgesvj.f @@ -714,7 +714,7 @@ SUBROUTINE CGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, * * Row-cyclic Jacobi SVD algorithm with column pivoting * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 DO 1868 q = 1, N diff --git a/lapack-netlib/SRC/cgsvj0.c b/lapack-netlib/SRC/cgsvj0.c index d390bce7d6..7e8730f8f1 100644 --- a/lapack-netlib/SRC/cgsvj0.c +++ b/lapack-netlib/SRC/cgsvj0.c @@ -638,7 +638,7 @@ f"> */ roottol = sqrt(*tol); - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; diff --git a/lapack-netlib/SRC/cgsvj0.f b/lapack-netlib/SRC/cgsvj0.f index 7426ecdf14..d017fd6970 100644 --- a/lapack-netlib/SRC/cgsvj0.f +++ b/lapack-netlib/SRC/cgsvj0.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CGSVJ0 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -193,7 +191,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexOTHERcomputational +*> \ingroup gsvj0 * *> \par Further Details: * ===================== @@ -321,7 +319,7 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * * .. Row-cyclic Jacobi SVD algorithm with column pivoting .. * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 * * .. Row-cyclic pivot strategy with de Rijk's pivoting .. @@ -499,7 +497,8 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, T = HALF / THETA CS = ONE - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL CROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*T ) IF ( RSVEC ) THEN CALL CROT( MVL, V(1,p), 1, @@ -528,7 +527,8 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL CROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*SN ) IF ( RSVEC ) THEN CALL CROT( MVL, V(1,p), 1, @@ -541,14 +541,17 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * .. have to use modified Gram-Schmidt like transformation CALL CCOPY( M, A( 1, p ), 1, $ WORK, 1 ) - CALL CLASCL( 'G', 0, 0, AAPP, ONE, M, + CALL CLASCL( 'G', 0, 0, AAPP, ONE, + $ M, $ 1, WORK, LDA, $ IERR ) - CALL CLASCL( 'G', 0, 0, AAQQ, ONE, M, + CALL CLASCL( 'G', 0, 0, AAQQ, ONE, + $ M, $ 1, A( 1, q ), LDA, IERR ) CALL CAXPY( M, -AAPQ, WORK, 1, $ A( 1, q ), 1 ) - CALL CLASCL( 'G', 0, 0, ONE, AAQQ, M, + CALL CLASCL( 'G', 0, 0, ONE, AAQQ, + $ M, $ 1, A( 1, q ), LDA, IERR ) SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) @@ -563,7 +566,8 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = SCNRM2( M, A( 1, q ), 1 ) + SVA( q ) = SCNRM2( M, A( 1, q ), + $ 1 ) ELSE T = ZERO AAQQ = ONE @@ -715,7 +719,8 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, IF( ABS( THETA ).GT.BIGTHETA ) THEN T = HALF / THETA CS = ONE - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL CROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*T ) IF( RSVEC ) THEN CALL CROT( MVL, V(1,p), 1, @@ -742,7 +747,8 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL CROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*SN ) IF( RSVEC ) THEN CALL CROT( MVL, V(1,p), 1, @@ -756,15 +762,18 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, IF( AAPP.GT.AAQQ ) THEN CALL CCOPY( M, A( 1, p ), 1, $ WORK, 1 ) - CALL CLASCL( 'G', 0, 0, AAPP, ONE, + CALL CLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, WORK,LDA, $ IERR ) - CALL CLASCL( 'G', 0, 0, AAQQ, ONE, + CALL CLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, A( 1, q ), LDA, $ IERR ) CALL CAXPY( M, -AAPQ, WORK, $ 1, A( 1, q ), 1 ) - CALL CLASCL( 'G', 0, 0, ONE, AAQQ, + CALL CLASCL( 'G', 0, 0, ONE, + $ AAQQ, $ M, 1, A( 1, q ), LDA, $ IERR ) SVA( q ) = AAQQ*SQRT( MAX( ZERO, @@ -773,15 +782,18 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ELSE CALL CCOPY( M, A( 1, q ), 1, $ WORK, 1 ) - CALL CLASCL( 'G', 0, 0, AAQQ, ONE, + CALL CLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, WORK,LDA, $ IERR ) - CALL CLASCL( 'G', 0, 0, AAPP, ONE, + CALL CLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, A( 1, p ), LDA, $ IERR ) CALL CAXPY( M, -CONJG(AAPQ), $ WORK, 1, A( 1, p ), 1 ) - CALL CLASCL( 'G', 0, 0, ONE, AAPP, + CALL CLASCL( 'G', 0, 0, ONE, + $ AAPP, $ M, 1, A( 1, p ), LDA, $ IERR ) SVA( p ) = AAPP*SQRT( MAX( ZERO, @@ -797,7 +809,8 @@ SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = SCNRM2( M, A( 1, q ), 1) + SVA( q ) = SCNRM2( M, A( 1, q ), + $ 1) ELSE T = ZERO AAQQ = ONE diff --git a/lapack-netlib/SRC/dgesvj.c b/lapack-netlib/SRC/dgesvj.c index ec7420bcb3..3f1af4fb54 100644 --- a/lapack-netlib/SRC/dgesvj.c +++ b/lapack-netlib/SRC/dgesvj.c @@ -992,7 +992,7 @@ f"> */ /* Row-cyclic Jacobi SVD algorithm with column pivoting */ - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; fastr[0] = 0.; diff --git a/lapack-netlib/SRC/dgesvj.f b/lapack-netlib/SRC/dgesvj.f index ef17833baa..fc2e8c35d3 100644 --- a/lapack-netlib/SRC/dgesvj.f +++ b/lapack-netlib/SRC/dgesvj.f @@ -695,7 +695,7 @@ SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, * * Row-cyclic Jacobi SVD algorithm with column pivoting * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 FASTR( 1 ) = ZERO * diff --git a/lapack-netlib/SRC/dgsvj0.c b/lapack-netlib/SRC/dgsvj0.c index 36e9869806..132d98cbb4 100644 --- a/lapack-netlib/SRC/dgsvj0.c +++ b/lapack-netlib/SRC/dgsvj0.c @@ -574,7 +574,7 @@ f"> */ /* -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#- */ - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; fastr[0] = 0.; diff --git a/lapack-netlib/SRC/dgsvj0.f b/lapack-netlib/SRC/dgsvj0.f index 1008aa8c6e..3c3b17672e 100644 --- a/lapack-netlib/SRC/dgsvj0.f +++ b/lapack-netlib/SRC/dgsvj0.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DGSVJ0 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -193,7 +191,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleOTHERcomputational +*> \ingroup gsvj0 * *> \par Further Details: * ===================== @@ -215,6 +213,7 @@ * ===================================================================== SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -259,7 +258,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, EXTERNAL IDAMAX, LSAME, DDOT, DNRM2 * .. * .. External Subroutines .. - EXTERNAL DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP, + EXTERNAL DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, + $ DSWAP, $ XERBLA * .. * .. Executable Statements .. @@ -314,7 +314,7 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * * -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#- * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 FASTR( 1 ) = ZERO * @@ -424,12 +424,15 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, IF( AAQQ.GE.ONE ) THEN ROTOK = ( SMALL*AAPP ).LE.AAQQ IF( AAPP.LT.( BIG / AAQQ ) ) THEN - AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( DDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL DCOPY( M, A( 1, p ), 1, WORK, 1 ) - CALL DLASCL( 'G', 0, 0, AAPP, D( p ), + CALL DCOPY( M, A( 1, p ), 1, WORK, + $ 1 ) + CALL DLASCL( 'G', 0, 0, AAPP, + $ D( p ), $ M, 1, WORK, LDA, IERR ) AAPQ = DDOT( M, WORK, 1, A( 1, q ), $ 1 )*D( q ) / AAQQ @@ -437,12 +440,15 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ELSE ROTOK = AAPP.LE.( AAQQ / SMALL ) IF( AAPP.GT.( SMALL / AAQQ ) ) THEN - AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( DDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL DCOPY( M, A( 1, q ), 1, WORK, 1 ) - CALL DLASCL( 'G', 0, 0, AAQQ, D( q ), + CALL DCOPY( M, A( 1, q ), 1, WORK, + $ 1 ) + CALL DLASCL( 'G', 0, 0, AAQQ, + $ D( q ), $ M, 1, WORK, LDA, IERR ) AAPQ = DDOT( M, WORK, 1, A( 1, p ), $ 1 )*D( p ) / AAPP @@ -511,7 +517,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, FASTR( 4 ) = -T*AQOAP D( p ) = D( p )*CS D( q ) = D( q )*CS - CALL DROTM( M, A( 1, p ), 1, + CALL DROTM( M, A( 1, p ), + $ 1, $ A( 1, q ), 1, $ FASTR ) IF( RSVEC )CALL DROTM( MVL, @@ -527,7 +534,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, D( p ) = D( p )*CS D( q ) = D( q ) / CS IF( RSVEC ) THEN - CALL DAXPY( MVL, -T*AQOAP, + CALL DAXPY( MVL, + $ -T*AQOAP, $ V( 1, q ), 1, $ V( 1, p ), 1 ) CALL DAXPY( MVL, @@ -541,13 +549,15 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL DAXPY( M, T*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) - CALL DAXPY( M, -CS*SN*AQOAP, + CALL DAXPY( M, + $ -CS*SN*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) D( p ) = D( p ) / CS D( q ) = D( q )*CS IF( RSVEC ) THEN - CALL DAXPY( MVL, T*APOAQ, + CALL DAXPY( MVL, + $ T*APOAQ, $ V( 1, p ), 1, $ V( 1, q ), 1 ) CALL DAXPY( MVL, @@ -560,7 +570,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL DAXPY( M, -T*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) - CALL DAXPY( M, CS*SN*APOAQ, + CALL DAXPY( M, + $ CS*SN*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) D( p ) = D( p )*CS @@ -601,15 +612,19 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * ELSE * .. have to use modified Gram-Schmidt like transformation - CALL DCOPY( M, A( 1, p ), 1, WORK, 1 ) - CALL DLASCL( 'G', 0, 0, AAPP, ONE, M, + CALL DCOPY( M, A( 1, p ), 1, WORK, + $ 1 ) + CALL DLASCL( 'G', 0, 0, AAPP, ONE, + $ M, $ 1, WORK, LDA, IERR ) - CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M, + CALL DLASCL( 'G', 0, 0, AAQQ, ONE, + $ M, $ 1, A( 1, q ), LDA, IERR ) TEMP1 = -AAPQ*D( p ) / D( q ) CALL DAXPY( M, TEMP1, WORK, 1, $ A( 1, q ), 1 ) - CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M, + CALL DLASCL( 'G', 0, 0, ONE, AAQQ, + $ M, $ 1, A( 1, q ), LDA, IERR ) SVA( q ) = AAQQ*DSQRT( MAX( ZERO, $ ONE-AAPQ*AAPQ ) ) @@ -623,7 +638,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = DNRM2( M, A( 1, q ), 1 )* + SVA( q ) = DNRM2( M, A( 1, q ), + $ 1 )* $ D( q ) ELSE T = ZERO @@ -724,12 +740,15 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ROTOK = ( SMALL*AAQQ ).LE.AAPP END IF IF( AAPP.LT.( BIG / AAQQ ) ) THEN - AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( DDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL DCOPY( M, A( 1, p ), 1, WORK, 1 ) - CALL DLASCL( 'G', 0, 0, AAPP, D( p ), + CALL DCOPY( M, A( 1, p ), 1, WORK, + $ 1 ) + CALL DLASCL( 'G', 0, 0, AAPP, + $ D( p ), $ M, 1, WORK, LDA, IERR ) AAPQ = DDOT( M, WORK, 1, A( 1, q ), $ 1 )*D( q ) / AAQQ @@ -741,12 +760,15 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ROTOK = AAQQ.LE.( AAPP / SMALL ) END IF IF( AAPP.GT.( SMALL / AAQQ ) ) THEN - AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( DDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL DCOPY( M, A( 1, q ), 1, WORK, 1 ) - CALL DLASCL( 'G', 0, 0, AAQQ, D( q ), + CALL DCOPY( M, A( 1, q ), 1, WORK, + $ 1 ) + CALL DLASCL( 'G', 0, 0, AAQQ, + $ D( q ), $ M, 1, WORK, LDA, IERR ) AAPQ = DDOT( M, WORK, 1, A( 1, p ), $ 1 )*D( p ) / AAPP @@ -810,7 +832,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, FASTR( 4 ) = -T*AQOAP D( p ) = D( p )*CS D( q ) = D( q )*CS - CALL DROTM( M, A( 1, p ), 1, + CALL DROTM( M, A( 1, p ), + $ 1, $ A( 1, q ), 1, $ FASTR ) IF( RSVEC )CALL DROTM( MVL, @@ -824,7 +847,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ A( 1, p ), 1, $ A( 1, q ), 1 ) IF( RSVEC ) THEN - CALL DAXPY( MVL, -T*AQOAP, + CALL DAXPY( MVL, + $ -T*AQOAP, $ V( 1, q ), 1, $ V( 1, p ), 1 ) CALL DAXPY( MVL, @@ -840,11 +864,13 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL DAXPY( M, T*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) - CALL DAXPY( M, -CS*SN*AQOAP, + CALL DAXPY( M, + $ -CS*SN*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) IF( RSVEC ) THEN - CALL DAXPY( MVL, T*APOAQ, + CALL DAXPY( MVL, + $ T*APOAQ, $ V( 1, p ), 1, $ V( 1, q ), 1 ) CALL DAXPY( MVL, @@ -859,7 +885,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL DAXPY( M, -T*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) - CALL DAXPY( M, CS*SN*APOAQ, + CALL DAXPY( M, + $ CS*SN*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) D( p ) = D( p )*CS @@ -900,34 +927,42 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * ELSE IF( AAPP.GT.AAQQ ) THEN - CALL DCOPY( M, A( 1, p ), 1, WORK, + CALL DCOPY( M, A( 1, p ), 1, + $ WORK, $ 1 ) - CALL DLASCL( 'G', 0, 0, AAPP, ONE, + CALL DLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, WORK, LDA, IERR ) - CALL DLASCL( 'G', 0, 0, AAQQ, ONE, + CALL DLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, A( 1, q ), LDA, $ IERR ) TEMP1 = -AAPQ*D( p ) / D( q ) CALL DAXPY( M, TEMP1, WORK, 1, $ A( 1, q ), 1 ) - CALL DLASCL( 'G', 0, 0, ONE, AAQQ, + CALL DLASCL( 'G', 0, 0, ONE, + $ AAQQ, $ M, 1, A( 1, q ), LDA, $ IERR ) SVA( q ) = AAQQ*DSQRT( MAX( ZERO, $ ONE-AAPQ*AAPQ ) ) MXSINJ = MAX( MXSINJ, SFMIN ) ELSE - CALL DCOPY( M, A( 1, q ), 1, WORK, + CALL DCOPY( M, A( 1, q ), 1, + $ WORK, $ 1 ) - CALL DLASCL( 'G', 0, 0, AAQQ, ONE, + CALL DLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, WORK, LDA, IERR ) - CALL DLASCL( 'G', 0, 0, AAPP, ONE, + CALL DLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, A( 1, p ), LDA, $ IERR ) TEMP1 = -AAPQ*D( q ) / D( p ) CALL DAXPY( M, TEMP1, WORK, 1, $ A( 1, p ), 1 ) - CALL DLASCL( 'G', 0, 0, ONE, AAPP, + CALL DLASCL( 'G', 0, 0, ONE, + $ AAPP, $ M, 1, A( 1, p ), LDA, $ IERR ) SVA( p ) = AAPP*DSQRT( MAX( ZERO, @@ -943,7 +978,8 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = DNRM2( M, A( 1, q ), 1 )* + SVA( q ) = DNRM2( M, A( 1, q ), + $ 1 )* $ D( q ) ELSE T = ZERO diff --git a/lapack-netlib/SRC/sgesvj.c b/lapack-netlib/SRC/sgesvj.c index 55adf6a97d..b05348e5f7 100644 --- a/lapack-netlib/SRC/sgesvj.c +++ b/lapack-netlib/SRC/sgesvj.c @@ -1214,7 +1214,7 @@ f"> */ /* Row-cyclic Jacobi SVD algorithm with column pivoting */ - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; fastr[0] = 0.f; diff --git a/lapack-netlib/SRC/sgesvj.f b/lapack-netlib/SRC/sgesvj.f index 25eec31df2..7f7d263932 100644 --- a/lapack-netlib/SRC/sgesvj.f +++ b/lapack-netlib/SRC/sgesvj.f @@ -681,7 +681,7 @@ SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, * * Row-cyclic Jacobi SVD algorithm with column pivoting * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 FASTR( 1 ) = ZERO * diff --git a/lapack-netlib/SRC/sgsvj0.c b/lapack-netlib/SRC/sgsvj0.c index 747b9eff95..c5dcc4f69f 100644 --- a/lapack-netlib/SRC/sgsvj0.c +++ b/lapack-netlib/SRC/sgsvj0.c @@ -813,7 +813,7 @@ f"> */ roottol = sqrt(*tol); - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; fastr[0] = 0.f; diff --git a/lapack-netlib/SRC/sgsvj0.f b/lapack-netlib/SRC/sgsvj0.f index 7e5fddeb0b..923f33374a 100644 --- a/lapack-netlib/SRC/sgsvj0.f +++ b/lapack-netlib/SRC/sgsvj0.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SGSVJ0 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -193,7 +191,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup realOTHERcomputational +*> \ingroup gsvj0 * *> \par Further Details: * ===================== @@ -215,6 +213,7 @@ * ===================================================================== SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -259,7 +258,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, EXTERNAL ISAMAX, LSAME, SDOT, SNRM2 * .. * .. External Subroutines .. - EXTERNAL SAXPY, SCOPY, SLASCL, SLASSQ, SROTM, SSWAP, + EXTERNAL SAXPY, SCOPY, SLASCL, SLASSQ, SROTM, + $ SSWAP, $ XERBLA * .. * .. Executable Statements .. @@ -314,7 +314,7 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * * .. Row-cyclic Jacobi SVD algorithm with column pivoting .. * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 FASTR( 1 ) = ZERO * @@ -424,12 +424,15 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, IF( AAQQ.GE.ONE ) THEN ROTOK = ( SMALL*AAPP ).LE.AAQQ IF( AAPP.LT.( BIG / AAQQ ) ) THEN - AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( SDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL SCOPY( M, A( 1, p ), 1, WORK, 1 ) - CALL SLASCL( 'G', 0, 0, AAPP, D( p ), + CALL SCOPY( M, A( 1, p ), 1, WORK, + $ 1 ) + CALL SLASCL( 'G', 0, 0, AAPP, + $ D( p ), $ M, 1, WORK, LDA, IERR ) AAPQ = SDOT( M, WORK, 1, A( 1, q ), $ 1 )*D( q ) / AAQQ @@ -437,12 +440,15 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ELSE ROTOK = AAPP.LE.( AAQQ / SMALL ) IF( AAPP.GT.( SMALL / AAQQ ) ) THEN - AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( SDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL SCOPY( M, A( 1, q ), 1, WORK, 1 ) - CALL SLASCL( 'G', 0, 0, AAQQ, D( q ), + CALL SCOPY( M, A( 1, q ), 1, WORK, + $ 1 ) + CALL SLASCL( 'G', 0, 0, AAQQ, + $ D( q ), $ M, 1, WORK, LDA, IERR ) AAPQ = SDOT( M, WORK, 1, A( 1, p ), $ 1 )*D( p ) / AAPP @@ -511,7 +517,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, FASTR( 4 ) = -T*AQOAP D( p ) = D( p )*CS D( q ) = D( q )*CS - CALL SROTM( M, A( 1, p ), 1, + CALL SROTM( M, A( 1, p ), + $ 1, $ A( 1, q ), 1, $ FASTR ) IF( RSVEC )CALL SROTM( MVL, @@ -527,7 +534,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, D( p ) = D( p )*CS D( q ) = D( q ) / CS IF( RSVEC ) THEN - CALL SAXPY( MVL, -T*AQOAP, + CALL SAXPY( MVL, + $ -T*AQOAP, $ V( 1, q ), 1, $ V( 1, p ), 1 ) CALL SAXPY( MVL, @@ -541,13 +549,15 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL SAXPY( M, T*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) - CALL SAXPY( M, -CS*SN*AQOAP, + CALL SAXPY( M, + $ -CS*SN*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) D( p ) = D( p ) / CS D( q ) = D( q )*CS IF( RSVEC ) THEN - CALL SAXPY( MVL, T*APOAQ, + CALL SAXPY( MVL, + $ T*APOAQ, $ V( 1, p ), 1, $ V( 1, q ), 1 ) CALL SAXPY( MVL, @@ -560,7 +570,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL SAXPY( M, -T*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) - CALL SAXPY( M, CS*SN*APOAQ, + CALL SAXPY( M, + $ CS*SN*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) D( p ) = D( p )*CS @@ -601,15 +612,19 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * ELSE * .. have to use modified Gram-Schmidt like transformation - CALL SCOPY( M, A( 1, p ), 1, WORK, 1 ) - CALL SLASCL( 'G', 0, 0, AAPP, ONE, M, + CALL SCOPY( M, A( 1, p ), 1, WORK, + $ 1 ) + CALL SLASCL( 'G', 0, 0, AAPP, ONE, + $ M, $ 1, WORK, LDA, IERR ) - CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M, + CALL SLASCL( 'G', 0, 0, AAQQ, ONE, + $ M, $ 1, A( 1, q ), LDA, IERR ) TEMP1 = -AAPQ*D( p ) / D( q ) CALL SAXPY( M, TEMP1, WORK, 1, $ A( 1, q ), 1 ) - CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M, + CALL SLASCL( 'G', 0, 0, ONE, AAQQ, + $ M, $ 1, A( 1, q ), LDA, IERR ) SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE-AAPQ*AAPQ ) ) @@ -623,7 +638,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = SNRM2( M, A( 1, q ), 1 )* + SVA( q ) = SNRM2( M, A( 1, q ), + $ 1 )* $ D( q ) ELSE T = ZERO @@ -724,12 +740,15 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ROTOK = ( SMALL*AAQQ ).LE.AAPP END IF IF( AAPP.LT.( BIG / AAQQ ) ) THEN - AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( SDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL SCOPY( M, A( 1, p ), 1, WORK, 1 ) - CALL SLASCL( 'G', 0, 0, AAPP, D( p ), + CALL SCOPY( M, A( 1, p ), 1, WORK, + $ 1 ) + CALL SLASCL( 'G', 0, 0, AAPP, + $ D( p ), $ M, 1, WORK, LDA, IERR ) AAPQ = SDOT( M, WORK, 1, A( 1, q ), $ 1 )*D( q ) / AAQQ @@ -741,12 +760,15 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ROTOK = AAQQ.LE.( AAPP / SMALL ) END IF IF( AAPP.GT.( SMALL / AAQQ ) ) THEN - AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, + AAPQ = ( SDOT( M, A( 1, p ), 1, + $ A( 1, $ q ), 1 )*D( p )*D( q ) / AAQQ ) $ / AAPP ELSE - CALL SCOPY( M, A( 1, q ), 1, WORK, 1 ) - CALL SLASCL( 'G', 0, 0, AAQQ, D( q ), + CALL SCOPY( M, A( 1, q ), 1, WORK, + $ 1 ) + CALL SLASCL( 'G', 0, 0, AAQQ, + $ D( q ), $ M, 1, WORK, LDA, IERR ) AAPQ = SDOT( M, WORK, 1, A( 1, p ), $ 1 )*D( p ) / AAPP @@ -810,7 +832,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, FASTR( 4 ) = -T*AQOAP D( p ) = D( p )*CS D( q ) = D( q )*CS - CALL SROTM( M, A( 1, p ), 1, + CALL SROTM( M, A( 1, p ), + $ 1, $ A( 1, q ), 1, $ FASTR ) IF( RSVEC )CALL SROTM( MVL, @@ -824,7 +847,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ A( 1, p ), 1, $ A( 1, q ), 1 ) IF( RSVEC ) THEN - CALL SAXPY( MVL, -T*AQOAP, + CALL SAXPY( MVL, + $ -T*AQOAP, $ V( 1, q ), 1, $ V( 1, p ), 1 ) CALL SAXPY( MVL, @@ -840,11 +864,13 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL SAXPY( M, T*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) - CALL SAXPY( M, -CS*SN*AQOAP, + CALL SAXPY( M, + $ -CS*SN*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) IF( RSVEC ) THEN - CALL SAXPY( MVL, T*APOAQ, + CALL SAXPY( MVL, + $ T*APOAQ, $ V( 1, p ), 1, $ V( 1, q ), 1 ) CALL SAXPY( MVL, @@ -859,7 +885,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, CALL SAXPY( M, -T*AQOAP, $ A( 1, q ), 1, $ A( 1, p ), 1 ) - CALL SAXPY( M, CS*SN*APOAQ, + CALL SAXPY( M, + $ CS*SN*APOAQ, $ A( 1, p ), 1, $ A( 1, q ), 1 ) D( p ) = D( p )*CS @@ -900,34 +927,42 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * ELSE IF( AAPP.GT.AAQQ ) THEN - CALL SCOPY( M, A( 1, p ), 1, WORK, + CALL SCOPY( M, A( 1, p ), 1, + $ WORK, $ 1 ) - CALL SLASCL( 'G', 0, 0, AAPP, ONE, + CALL SLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, WORK, LDA, IERR ) - CALL SLASCL( 'G', 0, 0, AAQQ, ONE, + CALL SLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, A( 1, q ), LDA, $ IERR ) TEMP1 = -AAPQ*D( p ) / D( q ) CALL SAXPY( M, TEMP1, WORK, 1, $ A( 1, q ), 1 ) - CALL SLASCL( 'G', 0, 0, ONE, AAQQ, + CALL SLASCL( 'G', 0, 0, ONE, + $ AAQQ, $ M, 1, A( 1, q ), LDA, $ IERR ) SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE-AAPQ*AAPQ ) ) MXSINJ = MAX( MXSINJ, SFMIN ) ELSE - CALL SCOPY( M, A( 1, q ), 1, WORK, + CALL SCOPY( M, A( 1, q ), 1, + $ WORK, $ 1 ) - CALL SLASCL( 'G', 0, 0, AAQQ, ONE, + CALL SLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, WORK, LDA, IERR ) - CALL SLASCL( 'G', 0, 0, AAPP, ONE, + CALL SLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, A( 1, p ), LDA, $ IERR ) TEMP1 = -AAPQ*D( q ) / D( p ) CALL SAXPY( M, TEMP1, WORK, 1, $ A( 1, p ), 1 ) - CALL SLASCL( 'G', 0, 0, ONE, AAPP, + CALL SLASCL( 'G', 0, 0, ONE, + $ AAPP, $ M, 1, A( 1, p ), LDA, $ IERR ) SVA( p ) = AAPP*SQRT( MAX( ZERO, @@ -943,7 +978,8 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = SNRM2( M, A( 1, q ), 1 )* + SVA( q ) = SNRM2( M, A( 1, q ), + $ 1 )* $ D( q ) ELSE T = ZERO diff --git a/lapack-netlib/SRC/zgesvj.c b/lapack-netlib/SRC/zgesvj.c index dba1be2a36..9e068bfda5 100644 --- a/lapack-netlib/SRC/zgesvj.c +++ b/lapack-netlib/SRC/zgesvj.c @@ -1087,7 +1087,7 @@ f"> */ /* Row-cyclic Jacobi SVD algorithm with column pivoting */ - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; i__1 = *n; for (q = 1; q <= i__1; ++q) { diff --git a/lapack-netlib/SRC/zgesvj.f b/lapack-netlib/SRC/zgesvj.f index 6c69f8cf06..c13754ced7 100644 --- a/lapack-netlib/SRC/zgesvj.f +++ b/lapack-netlib/SRC/zgesvj.f @@ -715,7 +715,7 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, * * Row-cyclic Jacobi SVD algorithm with column pivoting * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 DO 1868 q = 1, N diff --git a/lapack-netlib/SRC/zgsvj0.c b/lapack-netlib/SRC/zgsvj0.c index e79fc8b5ba..7d68bae62b 100644 --- a/lapack-netlib/SRC/zgsvj0.c +++ b/lapack-netlib/SRC/zgsvj0.c @@ -649,7 +649,7 @@ f"> */ roottol = sqrt(*tol); - emptsw = *n * (*n - 1) / 2; + emptsw = (*n / 2) * (*n - 1 + *n % 2); notrot = 0; diff --git a/lapack-netlib/SRC/zgsvj0.f b/lapack-netlib/SRC/zgsvj0.f index 075c243c29..401d2ac243 100644 --- a/lapack-netlib/SRC/zgsvj0.f +++ b/lapack-netlib/SRC/zgsvj0.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZGSVJ0 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -193,7 +191,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHERcomputational +*> \ingroup gsvj0 *> *> \par Further Details: * ===================== @@ -321,7 +319,7 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * * .. Row-cyclic Jacobi SVD algorithm with column pivoting .. * - EMPTSW = ( N*( N-1 ) ) / 2 + EMPTSW = ( N / 2 )*( N-1 + MOD( N, 2 ) ) NOTROT = 0 * * .. Row-cyclic pivot strategy with de Rijk's pivoting .. @@ -499,7 +497,8 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, T = HALF / THETA CS = ONE - CALL ZROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*T ) IF ( RSVEC ) THEN CALL ZROT( MVL, V(1,p), 1, @@ -528,7 +527,8 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * - CALL ZROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*SN ) IF ( RSVEC ) THEN CALL ZROT( MVL, V(1,p), 1, @@ -541,14 +541,17 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, * .. have to use modified Gram-Schmidt like transformation CALL ZCOPY( M, A( 1, p ), 1, $ WORK, 1 ) - CALL ZLASCL( 'G', 0, 0, AAPP, ONE, M, + CALL ZLASCL( 'G', 0, 0, AAPP, ONE, + $ M, $ 1, WORK, LDA, $ IERR ) - CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, M, + CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, + $ M, $ 1, A( 1, q ), LDA, IERR ) CALL ZAXPY( M, -AAPQ, WORK, 1, $ A( 1, q ), 1 ) - CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M, + CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, + $ M, $ 1, A( 1, q ), LDA, IERR ) SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) @@ -563,7 +566,8 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = DZNRM2( M, A( 1, q ), 1 ) + SVA( q ) = DZNRM2( M, A( 1, q ), + $ 1 ) ELSE T = ZERO AAQQ = ONE @@ -715,7 +719,8 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, IF( ABS( THETA ).GT.BIGTHETA ) THEN T = HALF / THETA CS = ONE - CALL ZROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*T ) IF( RSVEC ) THEN CALL ZROT( MVL, V(1,p), 1, @@ -742,7 +747,8 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * - CALL ZROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), + $ 1, $ CS, CONJG(OMPQ)*SN ) IF( RSVEC ) THEN CALL ZROT( MVL, V(1,p), 1, @@ -756,15 +762,18 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, IF( AAPP.GT.AAQQ ) THEN CALL ZCOPY( M, A( 1, p ), 1, $ WORK, 1 ) - CALL ZLASCL( 'G', 0, 0, AAPP, ONE, + CALL ZLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, WORK,LDA, $ IERR ) - CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, + CALL ZLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, A( 1, q ), LDA, $ IERR ) CALL ZAXPY( M, -AAPQ, WORK, $ 1, A( 1, q ), 1 ) - CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, + CALL ZLASCL( 'G', 0, 0, ONE, + $ AAQQ, $ M, 1, A( 1, q ), LDA, $ IERR ) SVA( q ) = AAQQ*SQRT( MAX( ZERO, @@ -773,15 +782,18 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, ELSE CALL ZCOPY( M, A( 1, q ), 1, $ WORK, 1 ) - CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, + CALL ZLASCL( 'G', 0, 0, AAQQ, + $ ONE, $ M, 1, WORK,LDA, $ IERR ) - CALL ZLASCL( 'G', 0, 0, AAPP, ONE, + CALL ZLASCL( 'G', 0, 0, AAPP, + $ ONE, $ M, 1, A( 1, p ), LDA, $ IERR ) CALL ZAXPY( M, -CONJG(AAPQ), $ WORK, 1, A( 1, p ), 1 ) - CALL ZLASCL( 'G', 0, 0, ONE, AAPP, + CALL ZLASCL( 'G', 0, 0, ONE, + $ AAPP, $ M, 1, A( 1, p ), LDA, $ IERR ) SVA( p ) = AAPP*SQRT( MAX( ZERO, @@ -797,7 +809,8 @@ SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. $ ( AAQQ.GT.ROOTSFMIN ) ) THEN - SVA( q ) = DZNRM2( M, A( 1, q ), 1) + SVA( q ) = DZNRM2( M, A( 1, q ), + $ 1) ELSE T = ZERO AAQQ = ONE diff --git a/lapack-netlib/TESTING/CMakeLists.txt b/lapack-netlib/TESTING/CMakeLists.txt index f9c0a965f3..3ea59f6fc3 100644 --- a/lapack-netlib/TESTING/CMakeLists.txt +++ b/lapack-netlib/TESTING/CMakeLists.txt @@ -58,6 +58,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) +add_test(NAME LAPACK-test_sgesvj_emptsw COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_sgesvj_emptsw) endif() if(BUILD_DOUBLE) @@ -93,6 +94,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) +add_test(NAME LAPACK-test_dgesvj_emptsw COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_dgesvj_emptsw) endif() if(BUILD_COMPLEX) @@ -126,6 +128,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) +add_test(NAME LAPACK-test_cgesvj_emptsw COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_cgesvj_emptsw) endif() if(BUILD_COMPLEX16) @@ -161,6 +164,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) +add_test(NAME LAPACK-test_zgesvj_emptsw COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EIG/test_zgesvj_emptsw) endif() diff --git a/lapack-netlib/TESTING/EIG/CMakeLists.txt b/lapack-netlib/TESTING/EIG/CMakeLists.txt index c203161689..784323fa76 100644 --- a/lapack-netlib/TESTING/EIG/CMakeLists.txt +++ b/lapack-netlib/TESTING/EIG/CMakeLists.txt @@ -113,12 +113,14 @@ 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) +add_eig_executable(test_sgesvj_emptsw test_sgesvj_emptsw.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) +add_eig_executable(test_cgesvj_emptsw test_sgesvj_emptsw.f90) endif() @@ -126,6 +128,7 @@ 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) +add_eig_executable(test_dgesvj_emptsw test_sgesvj_emptsw.f90) endif() if(BUILD_COMPLEX16) @@ -133,4 +136,5 @@ 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) +add_eig_executable(test_zgesvj_emptsw test_sgesvj_emptsw.f90) endif() diff --git a/lapack-netlib/TESTING/EIG/Makefile b/lapack-netlib/TESTING/EIG/Makefile index 73661f1e46..ffb5493ae9 100644 --- a/lapack-netlib/TESTING/EIG/Makefile +++ b/lapack-netlib/TESTING/EIG/Makefile @@ -171,6 +171,18 @@ test_cgges3_scaling: test_cgges3_scaling.f90 ../$(LAPACKLIB) test_zgges3_scaling: test_zgges3_scaling.f90 ../$(LAPACKLIB) $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +test_sgesvj_emptsw: test_sgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + +test_dgesvj_emptsw: test_dgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + +test_cgesvj_emptsw: test_cgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + +test_zgesvj_emptsw: test_zgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ + $(SDMDEIGTST): $(FRC) $(CDMDEIGTST): $(FRC) $(DDMDEIGTST): $(FRC) @@ -191,7 +203,7 @@ clean: cleanobj cleanexe cleanobj: rm -f *.o cleanexe: - rm -f xeigtst* xdmdeigtst* + rm -f xeigtst* xdmdeigtst* test_?gges3_scaling test_?gesvj_emptsw schkee.o: schkee.F $(FC) $(FFLAGS_DRV) -c -o $@ $< diff --git a/lapack-netlib/TESTING/EIG/test_cgesvj_emptsw.f90 b/lapack-netlib/TESTING/EIG/test_cgesvj_emptsw.f90 new file mode 100644 index 0000000000..3f6426fbbf --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_cgesvj_emptsw.f90 @@ -0,0 +1,117 @@ +! Coverage driver for the EMPTSW integer-overflow fix in CGESVJ/CGSVJ0 +! (issue #1312: EMPTSW = N*(N-1)/2 overflows 32-bit arithmetic for +! N >= 46342, so it is computed in a factored form that never +! overflows). +! +! CGSVJ0 is only reached through CGESVJ's triangular (JOBA='L'/'U') +! blocked path for N > 64, which the xeigtstc SVD tests (JOBA='G') +! never exercise, even with the N=66 pair in TESTING/svd.in. This +! driver therefore calls CGESVJ with N > 64 lower- and upper-triangular +! matrices (even and odd N), which takes the blocked CGESVJ -> CGSVJ0 +! path and covers the fixed lines in both files. It validates the +! live factorization directly (residual and orthogonality), so there +! is no replica of the EMPTSW expression here to go stale. +program test_cgesvj_emptsw + use, intrinsic :: iso_fortran_env, only: sp => real32 + implicit none + integer, parameter :: neven = 66, nodd = 65 + integer, parameter :: lw = 3*neven, lrw = 2*neven + complex(sp) :: a_even(neven,neven), a0_even(neven,neven) + complex(sp) :: v_even(neven,neven) + complex(sp) :: a_odd(nodd,nodd), a0_odd(nodd,nodd) + complex(sp) :: v_odd(nodd,nodd) + complex(sp) :: cwork(lw) + real(sp) :: sva_even(neven), sva_odd(nodd), rwork(lrw) + real(sp), external :: slamch + external :: cgesvj + integer :: i, j, info + integer :: failures + + failures = 0 + + do j = 1, neven + do i = 1, neven + if (i == j) then + a_even(i,j) = cmplx(2.0_sp+0.01_sp*real(i,sp), 0.02_sp*real(i,sp), sp) + else if (i > j) then + a_even(i,j) = cmplx(0.1_sp/real(abs(i-j)+1,sp), & + 0.05_sp/real(abs(i-j)+1,sp), sp) + else + a_even(i,j) = cmplx(0.0_sp, 0.0_sp, sp) + end if + end do + end do + a0_even = a_even + call cgesvj('L','U','V',neven,neven,a_even,neven,sva_even,neven, & + v_even,neven,cwork,lw,rwork,lrw,info) + call check_svd('even', a0_even, a_even, v_even, sva_even, info, & + failures) + + do j = 1, nodd + do i = 1, nodd + if (i == j) then + a_odd(i,j) = cmplx(2.0_sp+0.01_sp*real(i,sp), 0.02_sp*real(i,sp), sp) + else if (i < j) then + a_odd(i,j) = cmplx(0.1_sp/real(abs(i-j)+1,sp), & + 0.05_sp/real(abs(i-j)+1,sp), sp) + else + a_odd(i,j) = cmplx(0.0_sp, 0.0_sp, sp) + end if + end do + end do + a0_odd = a_odd + call cgesvj('U','U','V',nodd,nodd,a_odd,nodd,sva_odd,nodd, & + v_odd,nodd,cwork,lw,rwork,lrw,info) + call check_svd('odd ', a0_odd, a_odd, v_odd, sva_odd, info, failures) + + print '(a,i0)', 'CGESVJ EMPTSW failures = ', failures + if (failures /= 0) stop 1 + +contains + + subroutine check_svd(tag, a0, u, v, sva, info, failures) + character(*), intent(in) :: tag + complex(sp), intent(in) :: a0(:,:), u(:,:), v(:,:) + real(sp), intent(in) :: sva(:) + integer, intent(in) :: info + integer, intent(inout) :: failures + integer :: n, i, j + complex(sp), allocatable :: us(:,:), res(:,:), eye(:,:) + real(sp) :: anorm, rerr, oerr_u, oerr_v, tol, eps + + n = size(sva) + eps = slamch('Precision') + tol = real(100*n, sp)*eps + if (info /= 0) then + print '(a,a,a,i0)', 'CGESVJ EMPTSW (', tag, ') INFO = ', info + failures = failures + 1 + return + end if + do i = 1, n + if (sva(i) <= 0.0_sp) failures = failures + 1 + if (i > 1) then + if (sva(i-1) < sva(i)) failures = failures + 1 + end if + end do + allocate(us(n,n), res(n,n), eye(n,n)) + anorm = maxval(abs(a0)) + do j = 1, n + do i = 1, n + us(i,j) = u(i,j)*sva(j) + end do + end do + res = a0 - matmul(us, conjg(transpose(v))) + rerr = maxval(abs(res)) / anorm + eye = (0.0_sp, 0.0_sp) + do i = 1, n + eye(i,i) = (1.0_sp, 0.0_sp) + end do + oerr_u = maxval(abs(matmul(conjg(transpose(u)), u) - eye)) + oerr_v = maxval(abs(matmul(conjg(transpose(v)), v) - eye)) + print '(a,a,a,3es12.4)', 'CGESVJ EMPTSW (', tag, ') res/orthU/orthV = ', & + rerr, oerr_u, oerr_v + if (.not. (rerr <= tol .and. oerr_u <= tol .and. oerr_v <= tol)) & + failures = failures + 1 + end subroutine + +end program diff --git a/lapack-netlib/TESTING/EIG/test_dgesvj_emptsw.f90 b/lapack-netlib/TESTING/EIG/test_dgesvj_emptsw.f90 new file mode 100644 index 0000000000..1254369a5f --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_dgesvj_emptsw.f90 @@ -0,0 +1,113 @@ +! Coverage driver for the EMPTSW integer-overflow fix in DGESVJ/DGSVJ0 +! (issue #1312: EMPTSW = N*(N-1)/2 overflows 32-bit arithmetic for +! N >= 46342, so it is computed in a factored form that never +! overflows). +! +! DGSVJ0 is only reached through DGESVJ's triangular (JOBA='L'/'U') +! blocked path for N > 64, which the xeigtstd SVD tests (JOBA='G') +! never exercise, even with the N=66 pair in TESTING/svd.in. This +! driver therefore calls DGESVJ with N > 64 lower- and upper-triangular +! matrices (even and odd N), which takes the blocked DGESVJ -> DGSVJ0 +! path and covers the fixed lines in both files. It validates the +! live factorization directly (residual and orthogonality), so there +! is no replica of the EMPTSW expression here to go stale. +program test_dgesvj_emptsw + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + integer, parameter :: neven = 66, nodd = 65, lw = 3*neven + real(dp) :: a_even(neven,neven), a0_even(neven,neven) + real(dp) :: v_even(neven,neven), sva_even(neven) + real(dp) :: a_odd(nodd,nodd), a0_odd(nodd,nodd) + real(dp) :: v_odd(nodd,nodd), sva_odd(nodd) + real(dp) :: work(lw) + real(dp), external :: dlamch + external :: dgesvj + integer :: i, j, info + integer :: failures + + failures = 0 + + do j = 1, neven + do i = 1, neven + if (i == j) then + a_even(i,j) = 2.0_dp + 0.01_dp*real(i,dp) + else if (i > j) then + a_even(i,j) = 0.1_dp/real(abs(i-j)+1,dp) + else + a_even(i,j) = 0.0_dp + end if + end do + end do + a0_even = a_even + call dgesvj('L','U','V',neven,neven,a_even,neven,sva_even,neven, & + v_even,neven,work,lw,info) + call check_svd('even', a0_even, a_even, v_even, sva_even, info, & + failures) + + do j = 1, nodd + do i = 1, nodd + if (i == j) then + a_odd(i,j) = 2.0_dp + 0.01_dp*real(i,dp) + else if (i < j) then + a_odd(i,j) = 0.1_dp/real(abs(i-j)+1,dp) + else + a_odd(i,j) = 0.0_dp + end if + end do + end do + a0_odd = a_odd + call dgesvj('U','U','V',nodd,nodd,a_odd,nodd,sva_odd,nodd, & + v_odd,nodd,work,lw,info) + call check_svd('odd ', a0_odd, a_odd, v_odd, sva_odd, info, failures) + + print '(a,i0)', 'DGESVJ EMPTSW failures = ', failures + if (failures /= 0) stop 1 + +contains + + subroutine check_svd(tag, a0, u, v, sva, info, failures) + character(*), intent(in) :: tag + real(dp), intent(in) :: a0(:,:), u(:,:), v(:,:) + real(dp), intent(in) :: sva(:) + integer, intent(in) :: info + integer, intent(inout) :: failures + integer :: n, i, j + real(dp), allocatable :: us(:,:), res(:,:), eye(:,:) + real(dp) :: anorm, rerr, oerr_u, oerr_v, tol, eps + + n = size(sva) + eps = dlamch('Precision') + tol = real(100*n, dp)*eps + if (info /= 0) then + print '(a,a,a,i0)', 'DGESVJ EMPTSW (', tag, ') INFO = ', info + failures = failures + 1 + return + end if + do i = 1, n + if (sva(i) <= 0.0_dp) failures = failures + 1 + if (i > 1) then + if (sva(i-1) < sva(i)) failures = failures + 1 + end if + end do + allocate(us(n,n), res(n,n), eye(n,n)) + anorm = maxval(abs(a0)) + do j = 1, n + do i = 1, n + us(i,j) = u(i,j)*sva(j) + end do + end do + res = a0 - matmul(us, transpose(v)) + rerr = maxval(abs(res)) / anorm + eye = 0.0_dp + do i = 1, n + eye(i,i) = 1.0_dp + end do + oerr_u = maxval(abs(matmul(transpose(u), u) - eye)) + oerr_v = maxval(abs(matmul(transpose(v), v) - eye)) + print '(a,a,a,3es12.4)', 'DGESVJ EMPTSW (', tag, ') res/orthU/orthV = ', & + rerr, oerr_u, oerr_v + if (.not. (rerr <= tol .and. oerr_u <= tol .and. oerr_v <= tol)) & + failures = failures + 1 + end subroutine + +end program diff --git a/lapack-netlib/TESTING/EIG/test_sgesvj_emptsw.f90 b/lapack-netlib/TESTING/EIG/test_sgesvj_emptsw.f90 new file mode 100644 index 0000000000..ee2ea8459a --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_sgesvj_emptsw.f90 @@ -0,0 +1,113 @@ +! Coverage driver for the EMPTSW integer-overflow fix in SGESVJ/SGSVJ0 +! (issue #1312: EMPTSW = N*(N-1)/2 overflows 32-bit arithmetic for +! N >= 46342, so it is computed in a factored form that never +! overflows). +! +! SGSVJ0 is only reached through SGESVJ's triangular (JOBA='L'/'U') +! blocked path for N > 64, which the xeigtsts SVD tests (JOBA='G') +! never exercise, even with the N=66 pair in TESTING/svd.in. This +! driver therefore calls SGESVJ with N > 64 lower- and upper-triangular +! matrices (even and odd N), which takes the blocked SGESVJ -> SGSVJ0 +! path and covers the fixed lines in both files. It validates the +! live factorization directly (residual and orthogonality), so there +! is no replica of the EMPTSW expression here to go stale. +program test_sgesvj_emptsw + use, intrinsic :: iso_fortran_env, only: sp => real32 + implicit none + integer, parameter :: neven = 66, nodd = 65, lw = 3*neven + real(sp) :: a_even(neven,neven), a0_even(neven,neven) + real(sp) :: v_even(neven,neven), sva_even(neven) + real(sp) :: a_odd(nodd,nodd), a0_odd(nodd,nodd) + real(sp) :: v_odd(nodd,nodd), sva_odd(nodd) + real(sp) :: work(lw) + real(sp), external :: slamch + external :: sgesvj + integer :: i, j, info + integer :: failures + + failures = 0 + + do j = 1, neven + do i = 1, neven + if (i == j) then + a_even(i,j) = 2.0_sp + 0.01_sp*real(i,sp) + else if (i > j) then + a_even(i,j) = 0.1_sp/real(abs(i-j)+1,sp) + else + a_even(i,j) = 0.0_sp + end if + end do + end do + a0_even = a_even + call sgesvj('L','U','V',neven,neven,a_even,neven,sva_even,neven, & + v_even,neven,work,lw,info) + call check_svd('even', a0_even, a_even, v_even, sva_even, info, & + failures) + + do j = 1, nodd + do i = 1, nodd + if (i == j) then + a_odd(i,j) = 2.0_sp + 0.01_sp*real(i,sp) + else if (i < j) then + a_odd(i,j) = 0.1_sp/real(abs(i-j)+1,sp) + else + a_odd(i,j) = 0.0_sp + end if + end do + end do + a0_odd = a_odd + call sgesvj('U','U','V',nodd,nodd,a_odd,nodd,sva_odd,nodd, & + v_odd,nodd,work,lw,info) + call check_svd('odd ', a0_odd, a_odd, v_odd, sva_odd, info, failures) + + print '(a,i0)', 'SGESVJ EMPTSW failures = ', failures + if (failures /= 0) stop 1 + +contains + + subroutine check_svd(tag, a0, u, v, sva, info, failures) + character(*), intent(in) :: tag + real(sp), intent(in) :: a0(:,:), u(:,:), v(:,:) + real(sp), intent(in) :: sva(:) + integer, intent(in) :: info + integer, intent(inout) :: failures + integer :: n, i, j + real(sp), allocatable :: us(:,:), res(:,:), eye(:,:) + real(sp) :: anorm, rerr, oerr_u, oerr_v, tol, eps + + n = size(sva) + eps = slamch('Precision') + tol = real(100*n, sp)*eps + if (info /= 0) then + print '(a,a,a,i0)', 'SGESVJ EMPTSW (', tag, ') INFO = ', info + failures = failures + 1 + return + end if + do i = 1, n + if (sva(i) <= 0.0_sp) failures = failures + 1 + if (i > 1) then + if (sva(i-1) < sva(i)) failures = failures + 1 + end if + end do + allocate(us(n,n), res(n,n), eye(n,n)) + anorm = maxval(abs(a0)) + do j = 1, n + do i = 1, n + us(i,j) = u(i,j)*sva(j) + end do + end do + res = a0 - matmul(us, transpose(v)) + rerr = maxval(abs(res)) / anorm + eye = 0.0_sp + do i = 1, n + eye(i,i) = 1.0_sp + end do + oerr_u = maxval(abs(matmul(transpose(u), u) - eye)) + oerr_v = maxval(abs(matmul(transpose(v), v) - eye)) + print '(a,a,a,3es12.4)', 'SGESVJ EMPTSW (', tag, ') res/orthU/orthV = ', & + rerr, oerr_u, oerr_v + if (.not. (rerr <= tol .and. oerr_u <= tol .and. oerr_v <= tol)) & + failures = failures + 1 + end subroutine + +end program diff --git a/lapack-netlib/TESTING/EIG/test_zgesvj_emptsw.f90 b/lapack-netlib/TESTING/EIG/test_zgesvj_emptsw.f90 new file mode 100644 index 0000000000..373415a9b9 --- /dev/null +++ b/lapack-netlib/TESTING/EIG/test_zgesvj_emptsw.f90 @@ -0,0 +1,117 @@ +! Coverage driver for the EMPTSW integer-overflow fix in ZGESVJ/ZGSVJ0 +! (issue #1312: EMPTSW = N*(N-1)/2 overflows 32-bit arithmetic for +! N >= 46342, so it is computed in a factored form that never +! overflows). +! +! ZGSVJ0 is only reached through ZGESVJ's triangular (JOBA='L'/'U') +! blocked path for N > 64, which the xeigtstz SVD tests (JOBA='G') +! never exercise, even with the N=66 pair in TESTING/svd.in. This +! driver therefore calls ZGESVJ with N > 64 lower- and upper-triangular +! matrices (even and odd N), which takes the blocked ZGESVJ -> ZGSVJ0 +! path and covers the fixed lines in both files. It validates the +! live factorization directly (residual and orthogonality), so there +! is no replica of the EMPTSW expression here to go stale. +program test_zgesvj_emptsw + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + integer, parameter :: neven = 66, nodd = 65 + integer, parameter :: lw = 3*neven, lrw = 2*neven + complex(dp) :: a_even(neven,neven), a0_even(neven,neven) + complex(dp) :: v_even(neven,neven) + complex(dp) :: a_odd(nodd,nodd), a0_odd(nodd,nodd) + complex(dp) :: v_odd(nodd,nodd) + complex(dp) :: cwork(lw) + real(dp) :: sva_even(neven), sva_odd(nodd), rwork(lrw) + real(dp), external :: dlamch + external :: zgesvj + integer :: i, j, info + integer :: failures + + failures = 0 + + do j = 1, neven + do i = 1, neven + if (i == j) then + a_even(i,j) = cmplx(2.0_dp+0.01_dp*real(i,dp), 0.02_dp*real(i,dp), dp) + else if (i > j) then + a_even(i,j) = cmplx(0.1_dp/real(abs(i-j)+1,dp), & + 0.05_dp/real(abs(i-j)+1,dp), dp) + else + a_even(i,j) = cmplx(0.0_dp, 0.0_dp, dp) + end if + end do + end do + a0_even = a_even + call zgesvj('L','U','V',neven,neven,a_even,neven,sva_even,neven, & + v_even,neven,cwork,lw,rwork,lrw,info) + call check_svd('even', a0_even, a_even, v_even, sva_even, info, & + failures) + + do j = 1, nodd + do i = 1, nodd + if (i == j) then + a_odd(i,j) = cmplx(2.0_dp+0.01_dp*real(i,dp), 0.02_dp*real(i,dp), dp) + else if (i < j) then + a_odd(i,j) = cmplx(0.1_dp/real(abs(i-j)+1,dp), & + 0.05_dp/real(abs(i-j)+1,dp), dp) + else + a_odd(i,j) = cmplx(0.0_dp, 0.0_dp, dp) + end if + end do + end do + a0_odd = a_odd + call zgesvj('U','U','V',nodd,nodd,a_odd,nodd,sva_odd,nodd, & + v_odd,nodd,cwork,lw,rwork,lrw,info) + call check_svd('odd ', a0_odd, a_odd, v_odd, sva_odd, info, failures) + + print '(a,i0)', 'ZGESVJ EMPTSW failures = ', failures + if (failures /= 0) stop 1 + +contains + + subroutine check_svd(tag, a0, u, v, sva, info, failures) + character(*), intent(in) :: tag + complex(dp), intent(in) :: a0(:,:), u(:,:), v(:,:) + real(dp), intent(in) :: sva(:) + integer, intent(in) :: info + integer, intent(inout) :: failures + integer :: n, i, j + complex(dp), allocatable :: us(:,:), res(:,:), eye(:,:) + real(dp) :: anorm, rerr, oerr_u, oerr_v, tol, eps + + n = size(sva) + eps = dlamch('Precision') + tol = real(100*n, dp)*eps + if (info /= 0) then + print '(a,a,a,i0)', 'ZGESVJ EMPTSW (', tag, ') INFO = ', info + failures = failures + 1 + return + end if + do i = 1, n + if (sva(i) <= 0.0_dp) failures = failures + 1 + if (i > 1) then + if (sva(i-1) < sva(i)) failures = failures + 1 + end if + end do + allocate(us(n,n), res(n,n), eye(n,n)) + anorm = maxval(abs(a0)) + do j = 1, n + do i = 1, n + us(i,j) = u(i,j)*sva(j) + end do + end do + res = a0 - matmul(us, conjg(transpose(v))) + rerr = maxval(abs(res)) / anorm + eye = (0.0_dp, 0.0_dp) + do i = 1, n + eye(i,i) = (1.0_dp, 0.0_dp) + end do + oerr_u = maxval(abs(matmul(conjg(transpose(u)), u) - eye)) + oerr_v = maxval(abs(matmul(conjg(transpose(v)), v) - eye)) + print '(a,a,a,3es12.4)', 'ZGESVJ EMPTSW (', tag, ') res/orthU/orthV = ', & + rerr, oerr_u, oerr_v + if (.not. (rerr <= tol .and. oerr_u <= tol .and. oerr_v <= tol)) & + failures = failures + 1 + end subroutine + +end program diff --git a/lapack-netlib/TESTING/Makefile b/lapack-netlib/TESTING/Makefile index 3963260ac0..d64582b217 100644 --- a/lapack-netlib/TESTING/Makefile +++ b/lapack-netlib/TESTING/Makefile @@ -1,639 +1,223 @@ -####################################################################### -# This makefile runs the test programs for the linear equation routines -# and the eigenvalue routines in LAPACK. The test output files -# are grouped as follows: -# -# SLINTST,SEIGTST -- Single precision real test routines -# CLINTST,CEIGTST -- Single precision complex test routines -# DLINTST,DEIGTST -- Double precision real test routines -# ZLINTST,ZEIGTST -- Double precision complex test routines -# -# Test programs can be executed for all or some of the four different -# precisions. Enter 'make' followed by one or more of the data types -# desired. -# Some examples: +######################################################################## +# This is the makefile for the eigenvalue test program from LAPACK. +# The test files are organized as follows: +# +# AEIGTST -- Auxiliary test routines used in all precisions +# SCIGTST -- Auxiliary test routines used in single precision +# DZIGTST -- Auxiliary test routines used in double precision +# SEIGTST -- Single precision real test routines +# CEIGTST -- Single precision complex test routines +# DEIGTST -- Double precision real test routines +# ZEIGTST -- Double precision complex test routines +# +# Test programs can be generated for all or some of the four different +# precisions. Enter make followed by one or more of the data types +# desired. Some examples: # make single # make single complex # make single double complex complex16 # Alternatively, the command # make -# without any arguments runs all eight test programs. -# The executable files are called: -# xlintsts, xlintstd, xlintstc, and xlintstz for LIN -# xeigtsts, xeigtstd, xeigtstc, and xeigtstz for EIG -# -# To remove the output files after the tests have been run, enter -# make cleantest -# -# To re-run specific tests after a make, enter (for example): -# 'rm ssvd.out; make' or: -# 'make ssvd.out' or: -# 'touch svd.in; make' (to re-run the single precision SVD tests.) +# without any arguments creates all four test programs. +# The executable files are called +# xeigtsts, xeigtstd, xeigtstc, and xeigtstz # -# 'rm *svd.out; make' (to re-run all the SVD tests.) +# To remove the object files after the executable files have been +# created, enter +# make cleanobj +# On some systems, you can force the source files to be recompiled by +# entering (for example) +# make single FRC=FRC # -####################################################################### +######################################################################## -TOPSRCDIR = .. +TOPSRCDIR = ../.. include $(TOPSRCDIR)/make.inc +AEIGTST = \ + alahdg.o \ + alasum.o \ + alasvm.o \ + alareq.o \ + ilaenv.o \ + xerbla.o \ + xlaenv.o \ + chkxer.o + +SCIGTST = slafts.o slahd2.o slasum.o slatb9.o sstech.o sstect.o \ + ssvdch.o ssvdct.o ssxt1.o + +SEIGTST = schkee.o \ + sbdt01.o sbdt02.o sbdt03.o sbdt04.o sbdt05.o \ + schkbb.o schkbd.o schkbk.o schkbl.o schkec.o \ + schkgg.o schkgk.o schkgl.o schkhs.o schksb.o schkst.o schkst2stg.o schksb2stg.o \ + sckcsd.o sckglm.o sckgqr.o sckgsv.o scklse.o scsdts.o \ + sdrges.o sdrgev.o sdrges3.o sdrgev3.o sdrgsx.o sdrgvx.o \ + sdrvbd.o sdrvmr3.o sdrves.o sdrvev.o sdrvsg.o sdrvsg2stg.o \ + sdrvst.o sdrvst2stg.o sdrvsx.o sdrvvx.o \ + serrbd.o serrec.o serred.o serrgg.o serrhs.o serrst.o \ + sget02.o sget10.o sget22.o sget23.o sget24.o sget31.o \ + sget32.o sget33.o sget34.o sget35.o sget36.o \ + sget37.o sget38.o sget39.o sget40.o sget51.o sget52.o sget53.o \ + sget54.o sglmts.o sgqrts.o sgrqts.o sgsvts3.o \ + shst01.o slarfy.o slarhs.o slatm4.o slctes.o slctsx.o slsets.o sort01.o \ + sort03.o ssbt21.o ssgt01.o sslect.o sspt21.o sstt21.o \ + sstt22.o ssyl01.o ssyt21.o ssyt22.o + +SDMDEIGTST = schkdmd.o + +CEIGTST = cchkee.o \ + cbdt01.o cbdt02.o cbdt03.o cbdt05.o \ + cchkbb.o cchkbd.o cchkbk.o cchkbl.o cchkec.o \ + cchkgg.o cchkgk.o cchkgl.o cchkhb.o cchkhs.o cchkst.o cchkst2stg.o cchkhb2stg.o \ + cckcsd.o cckglm.o cckgqr.o cckgsv.o ccklse.o ccsdts.o \ + cdrges.o cdrgev.o cdrges3.o cdrgev3.o cdrgsx.o cdrgvx.o \ + cdrvbd.o cdrvmr3.o cdrves.o cdrvev.o cdrvsg.o cdrvsg2stg.o \ + cdrvst.o cdrvst2stg.o cdrvsx.o cdrvvx.o \ + cerrbd.o cerrec.o cerred.o cerrgg.o cerrhs.o cerrst.o \ + cget02.o cget10.o cget22.o cget23.o cget24.o \ + cget35.o cget36.o cget37.o cget38.o cget51.o cget52.o \ + cget54.o cglmts.o cgqrts.o cgrqts.o cgsvts3.o \ + chbt21.o chet21.o chet22.o chpt21.o chst01.o \ + clarfy.o clarhs.o clatm4.o clctes.o clctsx.o clsets.o csbmv.o \ + csgt01.o cslect.o csyl01.o\ + cstt21.o cstt22.o cunt01.o cunt03.o + +CDMDEIGTST = cchkdmd.o + +DZIGTST = dlafts.o dlahd2.o dlasum.o dlatb9.o dstech.o dstect.o \ + dsvdch.o dsvdct.o dsxt1.o + +DEIGTST = dchkee.o \ + dbdt01.o dbdt02.o dbdt03.o dbdt04.o dbdt05.o \ + dchkbb.o dchkbd.o dchkbk.o dchkbl.o dchkec.o \ + dchkgg.o dchkgk.o dchkgl.o dchkhs.o dchksb.o dchkst.o dchkst2stg.o dchksb2stg.o \ + dckcsd.o dckglm.o dckgqr.o dckgsv.o dcklse.o dcsdts.o \ + ddrges.o ddrgev.o ddrges3.o ddrgev3.o ddrgsx.o ddrgvx.o \ + ddrvbd.o ddrvmr3.o ddrves.o ddrvev.o ddrvsg.o ddrvsg2stg.o \ + ddrvst.o ddrvst2stg.o ddrvsx.o ddrvvx.o \ + derrbd.o derrec.o derred.o derrgg.o derrhs.o derrst.o \ + dget02.o dget10.o dget22.o dget23.o dget24.o dget31.o \ + dget32.o dget33.o dget34.o dget35.o dget36.o \ + dget37.o dget38.o dget39.o dget40.o dget51.o dget52.o dget53.o \ + dget54.o dglmts.o dgqrts.o dgrqts.o dgsvts3.o \ + dhst01.o dlarfy.o dlarhs.o dlatm4.o dlctes.o dlctsx.o dlsets.o dort01.o \ + dort03.o dsbt21.o dsgt01.o dslect.o dspt21.o dstt21.o \ + dstt22.o dsyl01.o dsyt21.o dsyt22.o + +DDMDEIGTST = dchkdmd.o + +ZEIGTST = zchkee.o \ + zbdt01.o zbdt02.o zbdt03.o zbdt05.o \ + zchkbb.o zchkbd.o zchkbk.o zchkbl.o zchkec.o \ + zchkgg.o zchkgk.o zchkgl.o zchkhb.o zchkhs.o zchkst.o zchkst2stg.o zchkhb2stg.o \ + zckcsd.o zckglm.o zckgqr.o zckgsv.o zcklse.o zcsdts.o \ + zdrges.o zdrgev.o zdrges3.o zdrgev3.o zdrgsx.o zdrgvx.o \ + zdrvbd.o zdrvmr3.o zdrves.o zdrvev.o zdrvsg.o zdrvsg2stg.o \ + zdrvst.o zdrvst2stg.o zdrvsx.o zdrvvx.o \ + zerrbd.o zerrec.o zerred.o zerrgg.o zerrhs.o zerrst.o \ + zget02.o zget10.o zget22.o zget23.o zget24.o \ + zget35.o zget36.o zget37.o zget38.o zget51.o zget52.o \ + zget54.o zglmts.o zgqrts.o zgrqts.o zgsvts3.o \ + zhbt21.o zhet21.o zhet22.o zhpt21.o zhst01.o \ + zlarfy.o zlarhs.o zlatm4.o zlctes.o zlctsx.o zlsets.o zsbmv.o \ + zsgt01.o zslect.o zsyl01.o\ + zstt21.o zstt22.o zunt01.o zunt03.o + +ZDMDEIGTST = zchkdmd.o + .PHONY: all -all: single complex double complex16 singleproto doubleproto complexproto complex16proto - -SEIGTST= snep.out \ - ssep.out \ - sse2.out \ - ssvd.out \ - sec.out \ - sed.out \ - sgg.out \ - sgd.out \ - ssb.out \ - ssg.out \ - sbal.out \ - sbak.out \ - sgbal.out \ - sgbak.out \ - sbb.out \ - sglm.out \ - sgqr.out \ - sgsv.out \ - scsd.out \ - slse.out - -SDMDEIGTST= sdmd.out - -CEIGTST= cnep.out \ - csep.out \ - cse2.out \ - csvd.out \ - cec.out \ - ced.out \ - cgg.out \ - cgd.out \ - csb.out \ - csg.out \ - cbal.out \ - cbak.out \ - cgbal.out \ - cgbak.out \ - cbb.out \ - cglm.out \ - cgqr.out \ - cgsv.out \ - ccsd.out \ - clse.out - -CDMDEIGTST= cdmd.out - -DEIGTST= dnep.out \ - dsep.out \ - dse2.out \ - dsvd.out \ - dec.out \ - ded.out \ - dgg.out \ - dgd.out \ - dsb.out \ - dsg.out \ - dbal.out \ - dbak.out \ - dgbal.out \ - dgbak.out \ - dbb.out \ - dglm.out \ - dgqr.out \ - dgsv.out \ - dcsd.out \ - dlse.out - -DDMDEIGTST= ddmd.out - -ZEIGTST= znep.out \ - zsep.out \ - zse2.out \ - zsvd.out \ - zec.out \ - zed.out \ - zgg.out \ - zgd.out \ - zsb.out \ - zsg.out \ - zbal.out \ - zbak.out \ - zgbal.out \ - zgbak.out \ - zbb.out \ - zglm.out \ - zgqr.out \ - zgsv.out \ - zcsd.out \ - zlse.out - -ZDMDEIGTST= zdmd.out - -SLINTST= stest.out - -SLINTSTPROTO= stest_rfp.out - -CLINTST= ctest.out - -CLINTSTPROTO= ctest_rfp.out - -DLINTST= dtest.out - -DLINTSTPROTO= dstest.out dtest_rfp.out - -ZLINTST= ztest.out - -ZLINTSTPROTO= zctest.out ztest_rfp.out +all: single complex double complex16 .PHONY: single complex double complex16 -single: $(SLINTST) $(SEIGTST) $(SDMDEIGTST) -complex: $(CLINTST) $(CEIGTST) $(CDMDEIGTST) -double: $(DLINTST) $(DEIGTST) $(DDMDEIGTST) -complex16: $(ZLINTST) $(ZEIGTST) $(ZDMDEIGTST) - -.PHONY: singleproto complexproto doubleproto complex16proto -singleproto: $(SLINTSTPROTO) -complexproto: $(CLINTSTPROTO) -doubleproto: $(DLINTSTPROTO) -complex16proto: $(ZLINTSTPROTO) - -# -# ======== SINGLE LIN TESTS =========================== - -stest.out: stest.in LIN/xlintsts - @echo Testing REAL LAPACK linear equation routines - ./LIN/xlintsts < stest.in > $@ 2>&1 -# -# ======== COMPLEX LIN TESTS ========================== - -ctest.out: ctest.in LIN/xlintstc - @echo Testing COMPLEX LAPACK linear equation routines - ./LIN/xlintstc < ctest.in > $@ 2>&1 -# -# ======== DOUBLE LIN TESTS =========================== - -dtest.out: dtest.in LIN/xlintstd - @echo Testing DOUBLE PRECISION LAPACK linear equation routines - ./LIN/xlintstd < dtest.in > $@ 2>&1 -# -# ======== COMPLEX16 LIN TESTS ======================== - -ztest.out: ztest.in LIN/xlintstz - @echo Testing COMPLEX16 LAPACK linear equation routines - ./LIN/xlintstz < ztest.in > $@ 2>&1 -# -# ======== SINGLE-DOUBLE PROTO LIN TESTS ============== - -dstest.out: dstest.in LIN/xlintstds - @echo Testing SINGLE-DOUBLE PRECISION LAPACK prototype linear equation routines - ./LIN/xlintstds < dstest.in > $@ 2>&1 -# -# ======== COMPLEX-COMPLEX16 LIN TESTS ======================== - -zctest.out: zctest.in LIN/xlintstzc - @echo Testing COMPLEX-COMPLEX16 LAPACK prototype linear equation routines - ./LIN/xlintstzc < zctest.in > $@ 2>&1 -# -# ======== SINGLE RFP LIN TESTS ======================== - -stest_rfp.out: stest_rfp.in LIN/xlintstrfs - @echo Testing REAL LAPACK RFP prototype linear equation routines - ./LIN/xlintstrfs < stest_rfp.in > $@ 2>&1 -# -# ======== COMPLEX16 RFP LIN TESTS ======================== - -dtest_rfp.out: dtest_rfp.in LIN/xlintstrfd - @echo Testing DOUBLE PRECISION LAPACK RFP prototype linear equation routines - ./LIN/xlintstrfd < dtest_rfp.in > $@ 2>&1 -# -# ======== COMPLEX16 RFP LIN TESTS ======================== - -ctest_rfp.out: ctest_rfp.in LIN/xlintstrfc - @echo Testing COMPLEX LAPACK RFP prototype linear equation routines - ./LIN/xlintstrfc < ctest_rfp.in > $@ 2>&1 -# -# ======== COMPLEX16 RFP LIN TESTS ======================== - -ztest_rfp.out: ztest_rfp.in LIN/xlintstrfz - @echo Testing COMPLEX16 LAPACK RFP prototype linear equation routines - ./LIN/xlintstrfz < ztest_rfp.in > $@ 2>&1 -# -# -# ======== SINGLE EIG TESTS =========================== -# - -snep.out: nep.in EIG/xeigtsts - @echo NEP: Testing Nonsymmetric Eigenvalue Problem routines - ./EIG/xeigtsts < nep.in > $@ 2>&1 - -ssep.out: sep.in EIG/xeigtsts - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtsts < sep.in > $@ 2>&1 - -sse2.out: se2.in EIG/xeigtsts - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtsts < se2.in > $@ 2>&1 - -ssvd.out: svd.in EIG/xeigtsts - @echo SVD: Testing Singular Value Decomposition routines - ./EIG/xeigtsts < svd.in > $@ 2>&1 - -sec.out: sec.in EIG/xeigtsts - @echo SEC: Testing REAL Eigen Condition Routines - ./EIG/xeigtsts < sec.in > $@ 2>&1 - -sed.out: sed.in EIG/xeigtsts - @echo SEV: Testing REAL Nonsymmetric Eigenvalue Driver - ./EIG/xeigtsts < sed.in > $@ 2>&1 - -sgg.out: sgg.in EIG/xeigtsts - @echo SGG: Testing REAL Nonsymmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtsts < sgg.in > $@ 2>&1 - -sgd.out: sgd.in EIG/xeigtsts - @echo SGD: Testing REAL Nonsymmetric Generalized Eigenvalue Problem driver routines - ./EIG/xeigtsts < sgd.in > $@ 2>&1 - -ssb.out: ssb.in EIG/xeigtsts - @echo SSB: Testing REAL Symmetric Eigenvalue Problem routines - ./EIG/xeigtsts < ssb.in > $@ 2>&1 - -ssg.out: ssg.in EIG/xeigtsts - @echo SSG: Testing REAL Symmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtsts < ssg.in > $@ 2>&1 - -sbal.out: sbal.in EIG/xeigtsts - @echo SGEBAL: Testing the balancing of a REAL general matrix - ./EIG/xeigtsts < sbal.in > $@ 2>&1 - -sbak.out: sbak.in EIG/xeigtsts - @echo SGEBAK: Testing the back transformation of a REAL balanced matrix - ./EIG/xeigtsts < sbak.in > $@ 2>&1 - -sgbal.out: sgbal.in EIG/xeigtsts - @echo SGGBAL: Testing the balancing of a pair of REAL general matrices - ./EIG/xeigtsts < sgbal.in > $@ 2>&1 - -sgbak.out: sgbak.in EIG/xeigtsts - @echo SGGBAK: Testing the back transformation of a pair of REAL balanced matrices - ./EIG/xeigtsts < sgbak.in > $@ 2>&1 - -sbb.out: sbb.in EIG/xeigtsts - @echo SBB: Testing banded Singular Value Decomposition routines - ./EIG/xeigtsts < sbb.in > $@ 2>&1 - -sglm.out: glm.in EIG/xeigtsts - @echo GLM: Testing Generalized Linear Regression Model routines - ./EIG/xeigtsts < glm.in > $@ 2>&1 - -sgqr.out: gqr.in EIG/xeigtsts - @echo GQR: Testing Generalized QR and RQ factorization routines - ./EIG/xeigtsts < gqr.in > $@ 2>&1 - -sgsv.out: gsv.in EIG/xeigtsts - @echo GSV: Testing Generalized Singular Value Decomposition routines - ./EIG/xeigtsts < gsv.in > $@ 2>&1 - -scsd.out: csd.in EIG/xeigtsts - @echo CSD: Testing CS Decomposition routines - ./EIG/xeigtsts < csd.in > $@ 2>&1 - -slse.out: lse.in EIG/xeigtsts - @echo LSE: Testing Constrained Linear Least Squares routines - ./EIG/xeigtsts < lse.in > $@ 2>&1 - -sdmd.out: sdmd.in EIG/xdmdeigtsts - @echo DMD: Testing Dynamic Mode Decomposition routines - ./EIG/xdmdeigtsts < sdmd.in > $@ 2>&1 -# -# ======== COMPLEX EIG TESTS =========================== - -cnep.out: nep.in EIG/xeigtstc - @echo NEP: Testing Nonsymmetric Eigenvalue Problem routines - ./EIG/xeigtstc < nep.in > $@ 2>&1 - -csep.out: sep.in EIG/xeigtstc - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtstc < sep.in > $@ 2>&1 - -cse2.out: se2.in EIG/xeigtstc - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtstc < se2.in > $@ 2>&1 - -csvd.out: svd.in EIG/xeigtstc - @echo SVD: Testing Singular Value Decomposition routines - ./EIG/xeigtstc < svd.in > $@ 2>&1 - -cec.out: cec.in EIG/xeigtstc - @echo CEC: Testing COMPLEX Eigen Condition Routines - ./EIG/xeigtstc < cec.in > $@ 2>&1 - -ced.out: ced.in EIG/xeigtstc - @echo CES: Testing COMPLEX Nonsymmetric Schur Form Driver - ./EIG/xeigtstc < ced.in > $@ 2>&1 - -cgg.out: cgg.in EIG/xeigtstc - @echo CGG: Testing COMPLEX Nonsymmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtstc < cgg.in > $@ 2>&1 - -cgd.out: cgd.in EIG/xeigtstc - @echo CGD: Testing COMPLEX Nonsymmetric Generalized Eigenvalue Problem driver routines - ./EIG/xeigtstc < cgd.in > $@ 2>&1 - -csb.out: csb.in EIG/xeigtstc - @echo CHB: Testing Hermitian Eigenvalue Problem routines - ./EIG/xeigtstc < csb.in > $@ 2>&1 - -csg.out: csg.in EIG/xeigtstc - @echo CSG: Testing Symmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtstc < csg.in > $@ 2>&1 - -cbal.out: cbal.in EIG/xeigtstc - @echo CGEBAL: Testing the balancing of a COMPLEX general matrix - ./EIG/xeigtstc < cbal.in > $@ 2>&1 - -cbak.out: cbak.in EIG/xeigtstc - @echo CGEBAK: Testing the back transformation of a COMPLEX balanced matrix - ./EIG/xeigtstc < cbak.in > $@ 2>&1 - -cgbal.out: cgbal.in EIG/xeigtstc - @echo CGGBAL: Testing the balancing of a pair of COMPLEX general matrices - ./EIG/xeigtstc < cgbal.in > $@ 2>&1 - -cgbak.out: cgbak.in EIG/xeigtstc - @echo CGGBAK: Testing the back transformation of a pair of COMPLEX balanced matrices - ./EIG/xeigtstc < cgbak.in > $@ 2>&1 - -cbb.out: cbb.in EIG/xeigtstc - @echo CBB: Testing banded Singular Value Decomposition routines - ./EIG/xeigtstc < cbb.in > $@ 2>&1 - -cglm.out: glm.in EIG/xeigtstc - @echo GLM: Testing Generalized Linear Regression Model routines - ./EIG/xeigtstc < glm.in > $@ 2>&1 - -cgqr.out: gqr.in EIG/xeigtstc - @echo GQR: Testing Generalized QR and RQ factorization routines - ./EIG/xeigtstc < gqr.in > $@ 2>&1 - -cgsv.out: gsv.in EIG/xeigtstc - @echo GSV: Testing Generalized Singular Value Decomposition routines - ./EIG/xeigtstc < gsv.in > $@ 2>&1 - -ccsd.out: csd.in EIG/xeigtstc - @echo CSD: Testing CS Decomposition routines - ./EIG/xeigtstc < csd.in > $@ 2>&1 - -clse.out: lse.in EIG/xeigtstc - @echo LSE: Testing Constrained Linear Least Squares routines - ./EIG/xeigtstc < lse.in > $@ 2>&1 - -cdmd.out: cdmd.in EIG/xdmdeigtstc - @echo DMD: Testing Dynamic Mode Decomposition routines - ./EIG/xdmdeigtstc < cdmd.in > $@ 2>&1 -# -# ======== DOUBLE EIG TESTS =========================== - -dnep.out: nep.in EIG/xeigtstd - @echo NEP: Testing Nonsymmetric Eigenvalue Problem routines - ./EIG/xeigtstd < nep.in > $@ 2>&1 - -dsep.out: sep.in EIG/xeigtstd - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtstd < sep.in > $@ 2>&1 - -dse2.out: se2.in EIG/xeigtstd - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtstd < se2.in > $@ 2>&1 - -dsvd.out: svd.in EIG/xeigtstd - @echo SVD: Testing Singular Value Decomposition routines - ./EIG/xeigtstd < svd.in > $@ 2>&1 - -dec.out: dec.in EIG/xeigtstd - @echo DEC: Testing DOUBLE PRECISION Eigen Condition Routines - ./EIG/xeigtstd < dec.in > $@ 2>&1 - -ded.out: ded.in EIG/xeigtstd - @echo DEV: Testing DOUBLE PRECISION Nonsymmetric Eigenvalue Driver - ./EIG/xeigtstd < ded.in > $@ 2>&1 - -dgg.out: dgg.in EIG/xeigtstd - @echo DGG: Testing DOUBLE PRECISION Nonsymmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtstd < dgg.in > $@ 2>&1 - -dgd.out: dgd.in EIG/xeigtstd - @echo DGD: Testing DOUBLE PRECISION Nonsymmetric Generalized Eigenvalue Problem driver routines - ./EIG/xeigtstd < dgd.in > $@ 2>&1 - -dsb.out: dsb.in EIG/xeigtstd - @echo DSB: Testing DOUBLE PRECISION Symmetric Eigenvalue Problem routines - ./EIG/xeigtstd < dsb.in > $@ 2>&1 - -dsg.out: dsg.in EIG/xeigtstd - @echo DSG: Testing DOUBLE PRECISION Symmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtstd < dsg.in > $@ 2>&1 - -dbal.out: dbal.in EIG/xeigtstd - @echo DGEBAL: Testing the balancing of a DOUBLE PRECISION general matrix - ./EIG/xeigtstd < dbal.in > $@ 2>&1 - -dbak.out: dbak.in EIG/xeigtstd - @echo DGEBAK: Testing the back transformation of a DOUBLE PRECISION balanced matrix - ./EIG/xeigtstd < dbak.in > $@ 2>&1 - -dgbal.out: dgbal.in EIG/xeigtstd - @echo DGGBAL: Testing the balancing of a pair of DOUBLE PRECISION general matrices - ./EIG/xeigtstd < dgbal.in > $@ 2>&1 - -dgbak.out: dgbak.in EIG/xeigtstd - @echo DGGBAK: Testing the back transformation of a pair of DOUBLE PRECISION balanced matrices - ./EIG/xeigtstd < dgbak.in > $@ 2>&1 - -dbb.out: dbb.in EIG/xeigtstd - @echo DBB: Testing banded Singular Value Decomposition routines - ./EIG/xeigtstd < dbb.in > $@ 2>&1 - -dglm.out: glm.in EIG/xeigtstd - @echo GLM: Testing Generalized Linear Regression Model routines - ./EIG/xeigtstd < glm.in > $@ 2>&1 - -dgqr.out: gqr.in EIG/xeigtstd - @echo GQR: Testing Generalized QR and RQ factorization routines - ./EIG/xeigtstd < gqr.in > $@ 2>&1 - -dgsv.out: gsv.in EIG/xeigtstd - @echo GSV: Testing Generalized Singular Value Decomposition routines - ./EIG/xeigtstd < gsv.in > $@ 2>&1 - -dcsd.out: csd.in EIG/xeigtstd - @echo CSD: Testing CS Decomposition routines - ./EIG/xeigtstd < csd.in > $@ 2>&1 - -dlse.out: lse.in EIG/xeigtstd - @echo LSE: Testing Constrained Linear Least Squares routines - ./EIG/xeigtstd < lse.in > $@ 2>&1 - -ddmd.out: ddmd.in EIG/xdmdeigtstd - @echo DMD: Testing Dynamic Mode Decomposition routines - ./EIG/xdmdeigtstd < ddmd.in > $@ 2>&1 -# -# ======== COMPLEX16 EIG TESTS =========================== - -znep.out: nep.in EIG/xeigtstz - @echo NEP: Testing Nonsymmetric Eigenvalue Problem routines - ./EIG/xeigtstz < nep.in > $@ 2>&1 - -zsep.out: sep.in EIG/xeigtstz - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtstz < sep.in > $@ 2>&1 - -zse2.out: se2.in EIG/xeigtstz - @echo SEP: Testing Symmetric Eigenvalue Problem routines - ./EIG/xeigtstz < se2.in > $@ 2>&1 - -zsvd.out: svd.in EIG/xeigtstz - @echo SVD: Testing Singular Value Decomposition routines - ./EIG/xeigtstz < svd.in > $@ 2>&1 - -zec.out: zec.in EIG/xeigtstz - @echo ZEC: Testing COMPLEX16 Eigen Condition Routines - ./EIG/xeigtstz < zec.in > $@ 2>&1 - -zed.out: zed.in EIG/xeigtstz - @echo ZES: Testing COMPLEX16 Nonsymmetric Schur Form Driver - ./EIG/xeigtstz < zed.in > $@ 2>&1 - -zgg.out: zgg.in EIG/xeigtstz - @echo ZGG: Testing COMPLEX16 Nonsymmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtstz < zgg.in > $@ 2>&1 - -zgd.out: zgd.in EIG/xeigtstz - @echo ZGD: Testing COMPLEX16 Nonsymmetric Generalized Eigenvalue Problem driver routines - ./EIG/xeigtstz < zgd.in > $@ 2>&1 - -zsb.out: zsb.in EIG/xeigtstz - @echo ZHB: Testing Hermitian Eigenvalue Problem routines - ./EIG/xeigtstz < zsb.in > $@ 2>&1 - -zsg.out: zsg.in EIG/xeigtstz - @echo ZSG: Testing Symmetric Generalized Eigenvalue Problem routines - ./EIG/xeigtstz < zsg.in > $@ 2>&1 - -zbal.out: zbal.in EIG/xeigtstz - @echo ZGEBAL: Testing the balancing of a COMPLEX16 general matrix - ./EIG/xeigtstz < zbal.in > $@ 2>&1 - -zbak.out: zbak.in EIG/xeigtstz - @echo ZGEBAK: Testing the back transformation of a COMPLEX16 balanced matrix - ./EIG/xeigtstz < zbak.in > $@ 2>&1 - -zgbal.out: zgbal.in EIG/xeigtstz - @echo ZGGBAL: Testing the balancing of a pair of COMPLEX general matrices - ./EIG/xeigtstz < zgbal.in > $@ 2>&1 - -zgbak.out: zgbak.in EIG/xeigtstz - @echo ZGGBAK: Testing the back transformation of a pair of COMPLEX16 balanced matrices - ./EIG/xeigtstz < zgbak.in > $@ 2>&1 - -zbb.out: zbb.in EIG/xeigtstz - @echo ZBB: Testing banded Singular Value Decomposition routines - ./EIG/xeigtstz < zbb.in > $@ 2>&1 - -zglm.out: glm.in EIG/xeigtstz - @echo GLM: Testing Generalized Linear Regression Model routines - ./EIG/xeigtstz < glm.in > $@ 2>&1 - -zgqr.out: gqr.in EIG/xeigtstz - @echo GQR: Testing Generalized QR and RQ factorization routines - ./EIG/xeigtstz < gqr.in > $@ 2>&1 - -zgsv.out: gsv.in EIG/xeigtstz - @echo GSV: Testing Generalized Singular Value Decomposition routines - ./EIG/xeigtstz < gsv.in > $@ 2>&1 - -zcsd.out: csd.in EIG/xeigtstz - @echo CSD: Testing CS Decomposition routines - ./EIG/xeigtstz < csd.in > $@ 2>&1 - -zlse.out: lse.in EIG/xeigtstz - @echo LSE: Testing Constrained Linear Least Squares routines - ./EIG/xeigtstz < lse.in > $@ 2>&1 - -zdmd.out: zdmd.in EIG/xdmdeigtstz - @echo DMD: Testing Dynamic Mode Decomposition routines - ./EIG/xdmdeigtstz < zdmd.in > $@ 2>&1 -# ============================================================================== - -LIN/xlintsts: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintsts - -LIN/xlintstc: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstc - -LIN/xlintstd: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstd - -LIN/xlintstz: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstz +single: xeigtsts xdmdeigtsts test_sgges3_scaling +complex: xeigtstc xdmdeigtstc test_cgges3_scaling +double: xeigtstd xdmdeigtstd test_dgges3_scaling +complex16: xeigtstz xdmdeigtstz test_zgges3_scaling -LIN/xlintstrfs: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstrfs +xdmdeigtsts: $(SDMDEIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -LIN/xlintstrfc: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstrfc +xdmdeigtstc: $(CDMDEIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -LIN/xlintstrfd: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstrfd +xdmdeigtstd: $(DDMDEIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -LIN/xlintstrfz: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstrfz +xdmdeigtstz: $(ZDMDEIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -LIN/xlintstds: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstds +xeigtsts: $(SEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -LIN/xlintstzc: $(FRCLIN) $(FRC) - $(MAKE) -C LIN xlintstzc +xeigtstc: $(CEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xeigtsts: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xeigtsts +xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xdmdeigtsts: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xdmdeigtsts +xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xeigtstc: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xeigtstc +test_sgges3_scaling: test_sgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xdmdeigtstc: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xdmdeigtstc +test_dgges3_scaling: test_dgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xeigtstd: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xeigtstd +test_cgges3_scaling: test_cgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xdmdeigtstd: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xdmdeigtstd +test_zgges3_scaling: test_zgges3_scaling.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xeigtstz: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xeigtstz +test_sgesvj_emptsw: test_sgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -EIG/xdmdeigtstz: $(FRCEIG) $(FRC) - $(MAKE) -C EIG xdmdeigtstz +test_dgesvj_emptsw: test_dgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -.PHONY: clean cleantest -clean: cleantest -cleantest: - rm -f *.out core +test_cgesvj_emptsw: test_cgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -FRCLIN: - @FRCLIN=$(FRCLIN) +test_zgesvj_emptsw: test_zgesvj_emptsw.f90 ../$(LAPACKLIB) + $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ -FRCEIG: - @FRCEIG=$(FRCEIG) +$(SDMDEIGTST): $(FRC) +$(CDMDEIGTST): $(FRC) +$(DDMDEIGTST): $(FRC) +$(ZDMDEIGTST): $(FRC) +$(AEIGTST): $(FRC) +$(SCIGTST): $(FRC) +$(DZIGTST): $(FRC) +$(SEIGTST): $(FRC) +$(CEIGTST): $(FRC) +$(DEIGTST): $(FRC) +$(ZEIGTST): $(FRC) FRC: @FRC=$(FRC) -.NOTPARALLEL: +.PHONY: clean cleanobj cleanexe +clean: cleanobj cleanexe +cleanobj: + rm -f *.o +cleanexe: + rm -f xeigtst* xdmdeigtst* + +schkee.o: schkee.F + $(FC) $(FFLAGS_DRV) -c -o $@ $< +dchkee.o: dchkee.F + $(FC) $(FFLAGS_DRV) -c -o $@ $< +cchkee.o: cchkee.F + $(FC) $(FFLAGS_DRV) -c -o $@ $< +zchkee.o: zchkee.F + $(FC) $(FFLAGS_DRV) -c -o $@ $< +schkdmd.o: schkdmd.f90 + $(FC) $(FFLAGS_DRV) -c -o $@ $< +cchkdmd.o: cchkdmd.f90 + $(FC) $(FFLAGS_DRV) -c -o $@ $< +dchkdmd.o: dchkdmd.f90 + $(FC) $(FFLAGS_DRV) -c -o $@ $< +zchkdmd.o: zchkdmd.f90 + $(FC) $(FFLAGS_DRV) -c -o $@ $<