diff --git a/conda-recipe/bld.bat b/conda-recipe/bld.bat index 87b1481d..88bfedd3 100644 --- a/conda-recipe/bld.bat +++ b/conda-recipe/bld.bat @@ -1,4 +1,10 @@ -"%PYTHON%" setup.py install +:: Uncoment following two lines for local test build +cd %RECIPE_DIR% +cd .. + +"%PYTHON%" setup.py build --compiler=mingw32 +"%PYTHON%" setup.py install --skip-build + if errorlevel 1 exit 1 :: Add more build steps here, if they are necessary. diff --git a/conda-recipe/build.sh b/conda-recipe/build.sh new file mode 100644 index 00000000..bce90005 --- /dev/null +++ b/conda-recipe/build.sh @@ -0,0 +1,2 @@ +cd $RECIPE_DIR/.. +LDFLAGS="-shared" FFLAGS="-fPIC" $PYTHON setup.py install diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 74229a02..9644848a 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -5,15 +5,12 @@ package: build: number: 1 - script: - - cd $RECIPE_DIR/.. - - $PYTHON setup.py install - requirements: build: - python - - numpy + - numpy >=1.13.3 - lapack + - m2w64-gcc-fortran # [win] run: - python diff --git a/slycot/setup.py b/slycot/setup.py index f44f0b28..bcdb3a17 100644 --- a/slycot/setup.py +++ b/slycot/setup.py @@ -3,7 +3,7 @@ import glob import os import sys - +import sysconfig def configuration(parent_package='', top_path=None): from numpy.distutils.misc_util import Configuration @@ -21,14 +21,38 @@ def configuration(parent_package='', top_path=None): f2py_sources = ['src/_wrapper.pyf'] + pyver = sysconfig.get_config_var('VERSION') + if sys.platform == 'win32': - liblist = ['liblapack', 'libblas'] + liblist = [ + 'lapack', 'lapacke', 'blas', 'gfortran' + ] + extra_objects = [ + ] + ppath = os.sep.join(sys.executable.split(os.sep)[:-1]) + library_dirs = [r'\Library\lib', ] + library_dirs = [ppath + l for l in library_dirs] + extra_link_args = [ ] + extra_compile_args = [ ] else: - liblist = ['lapack'] + # this is needed on Py 3.x, and fails on Py 2.7 + try: + abiflags = sys.abiflags + except AttributeError: + abiflags = '' + liblist = ['lapack', 'blas', 'python'+pyver+abiflags] + extra_objects = [] + library_dirs = [] + extra_link_args = [] + extra_compile_args = [] config.add_extension( name='_wrapper', libraries=liblist, + extra_objects=extra_objects, + extra_link_args=extra_link_args, + library_dirs=library_dirs, + extra_compile_args=extra_compile_args, sources=fortran_sources + f2py_sources) config.make_config_py() # installs __config__.py diff --git a/slycot/src/AB08NX.f b/slycot/src/AB08NX.f index d67f6a19..ce7c9701 100644 --- a/slycot/src/AB08NX.f +++ b/slycot/src/AB08NX.f @@ -199,8 +199,8 @@ SUBROUTINE AB08NX( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD, INTEGER ILAENV EXTERNAL ILAENV C .. External Subroutines .. - EXTERNAL DLAPMT, DLARFG, DLASET, DLATZM, DORMQR, DORMRQ, - $ MB03OY, MB03PY, XERBLA + EXTERNAL DLAPMT, DLARFG, DLASET, SLCT_DLATZM, DORMQR, + $ DORMRQ, MB03OY, MB03PY, XERBLA C .. Intrinsic Functions .. INTRINSIC INT, MAX, MIN C .. Executable Statements .. @@ -308,7 +308,8 @@ SUBROUTINE AB08NX( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD, C DO 40 I1 = 1, SIGMA CALL DLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1, T ) - CALL DLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, T, + CALL SLCT_DLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), + $ 1, T, $ ABCD(IROW,I1+1), ABCD(IROW+1,I1+1), LDABCD, $ DWORK ) IROW = IROW + 1 diff --git a/slycot/src/AB13AD.f b/slycot/src/AB13AD.f index fb2b2018..b0889dd6 100644 --- a/slycot/src/AB13AD.f +++ b/slycot/src/AB13AD.f @@ -197,6 +197,7 @@ DOUBLE PRECISION FUNCTION AB13AD( DICO, EQUIL, N, M, P, ALPHA, A, C Nov. 1998, V. Sima, Research Institute for Informatics, Bucharest. C Dec. 1998, V. Sima, Katholieke Univ. Leuven, Leuven. C Oct. 2001, V. Sima, Research Institute for Informatics, Bucharest. +C Jun. 2017, RvP, made 1st error return value zero C C KEYWORDS C @@ -260,8 +261,10 @@ DOUBLE PRECISION FUNCTION AB13AD( DICO, EQUIL, N, M, P, ALPHA, A, C IF( INFO.NE.0 ) THEN C -C Error return. +C Error return. C +C + AB13AD = ZERO CALL XERBLA( 'AB13AD', -INFO ) RETURN END IF diff --git a/slycot/src/AB13AX.f b/slycot/src/AB13AX.f index 4053e2a7..ace56a1f 100644 --- a/slycot/src/AB13AX.f +++ b/slycot/src/AB13AX.f @@ -195,6 +195,7 @@ DOUBLE PRECISION FUNCTION AB13AX( DICO, N, M, P, A, LDA, B, LDB, C C Error return. C + AB13AX = ZERO CALL XERBLA( 'AB13AX', -INFO ) RETURN END IF diff --git a/slycot/src/AB13BD.f b/slycot/src/AB13BD.f index ac69fd7b..72a3b69b 100644 --- a/slycot/src/AB13BD.f +++ b/slycot/src/AB13BD.f @@ -292,6 +292,7 @@ DOUBLE PRECISION FUNCTION AB13BD( DICO, JOBN, N, M, P, A, LDA, C C Error return. C + AB13BD = ZERO CALL XERBLA( 'AB13BD', -INFO ) RETURN END IF @@ -300,6 +301,7 @@ DOUBLE PRECISION FUNCTION AB13BD( DICO, JOBN, N, M, P, A, LDA, C S2NORM = DLANGE( 'Frobenius', P, M, D, LDD, DWORK ) IF( .NOT.DISCR .AND. S2NORM.NE.ZERO ) THEN + AB13BD = ZERO INFO = 5 RETURN END IF diff --git a/slycot/src/AB13CD.f b/slycot/src/AB13CD.f index ec9fa255..fc9430f1 100644 --- a/slycot/src/AB13CD.f +++ b/slycot/src/AB13CD.f @@ -237,13 +237,17 @@ DOUBLE PRECISION FUNCTION AB13CD( N, M, NP, A, LDA, B, LDB, C, INFO = -17 END IF IF( INFO.NE.0 ) THEN + AB13CD = ZERO CALL XERBLA( 'AB13CD', -INFO ) RETURN END IF C C Quick return if possible. C - IF( M.EQ.0 .OR. NP.EQ.0 ) RETURN + IF( M.EQ.0 .OR. NP.EQ.0 ) THEN + AB13CD = ZERO + RETURN + END IF C C Workspace usage. C @@ -262,6 +266,7 @@ DOUBLE PRECISION FUNCTION AB13CD( N, M, NP, A, LDA, B, LDB, C, $ INFO2 ) IF( INFO2.GT.0 ) THEN INFO = 4 + AB13CD = ZERO RETURN END IF GAMMAL = DWORK( IW6+1 ) diff --git a/slycot/src/AB13DX.f b/slycot/src/AB13DX.f index 09362b7c..5de0bd78 100644 --- a/slycot/src/AB13DX.f +++ b/slycot/src/AB13DX.f @@ -324,6 +324,7 @@ DOUBLE PRECISION FUNCTION AB13DX( DICO, JOBE, JOBD, N, M, P, END IF C IF( INFO.NE.0 ) THEN + AB13DX = ZERO CALL XERBLA( 'AB13DX', -INFO ) RETURN END IF @@ -360,6 +361,7 @@ DOUBLE PRECISION FUNCTION AB13DX( DICO, JOBE, JOBD, N, M, P, $ LDWORK-IWRK+1, IERR ) IF( IERR.GT.0 ) THEN INFO = N + 1 + AB13DX = ZERO RETURN END IF AB13DX = DWORK( IS ) @@ -391,6 +393,7 @@ DOUBLE PRECISION FUNCTION AB13DX( DICO, JOBE, JOBD, N, M, P, INFO = IERR DWORK( 1 ) = ONE CWORK( 1 ) = ONE + AB13DX = ZERO RETURN END IF CALL MB02RD( 'No Transpose', N, M, A, LDA, IWORK, B, LDB, @@ -415,6 +418,7 @@ DOUBLE PRECISION FUNCTION AB13DX( DICO, JOBE, JOBD, N, M, P, END IF IF( IERR.GT.0 ) THEN INFO = N + 1 + AB13DX = ZERO RETURN END IF C @@ -515,6 +519,7 @@ DOUBLE PRECISION FUNCTION AB13DX( DICO, JOBE, JOBD, N, M, P, INFO = IERR DWORK( 1 ) = ONE CWORK( 1 ) = ICWK - 1 + AB13DX = ZERO RETURN END IF CALL MB02RZ( 'No Transpose', N, M, CWORK, N, IWORK, diff --git a/slycot/src/AB8NXZ.f b/slycot/src/AB8NXZ.f index 9ec0da56..a0976d23 100755 --- a/slycot/src/AB8NXZ.f +++ b/slycot/src/AB8NXZ.f @@ -197,7 +197,7 @@ SUBROUTINE AB8NXZ( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD, LOGICAL LQUERY INTEGER I1, IK, IROW, ITAU, IZ, JWORK, MM1, MNTAU, MNU, $ MPM, NB, NP, RANK, RO1, TAU, WRKOPT - COMPLEX*16 TC + COMPLEX*16 TC, TCCONJ C .. Local Arrays .. DOUBLE PRECISION SVAL(3) C .. External Functions .. @@ -205,7 +205,7 @@ SUBROUTINE AB8NXZ( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD, EXTERNAL ILAENV C .. External Subroutines .. EXTERNAL MB3OYZ, MB3PYZ, XERBLA, ZLAPMT, ZLARFG, ZLASET, - $ ZLATZM, ZUNMQR, ZUNMRQ + $ SLCT_ZLATZM, ZUNMQR, ZUNMRQ C .. Intrinsic Functions .. INTRINSIC DCONJG, INT, MAX, MIN C .. Executable Statements .. @@ -315,7 +315,9 @@ SUBROUTINE AB8NXZ( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD, DO 40 I1 = 1, SIGMA CALL ZLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1, $ TC ) - CALL ZLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, +C RvP 170623 slicot-specific ZLATZM + TCCONJ = DCONJG( TC ) + CALL SLCT_ZLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, $ DCONJG( TC ), ABCD(IROW,I1+1), $ ABCD(IROW+1,I1+1), LDABCD, ZWORK ) IROW = IROW + 1 diff --git a/slycot/src/AG08BY.f b/slycot/src/AG08BY.f index 7e980bf8..d6517dea 100644 --- a/slycot/src/AG08BY.f +++ b/slycot/src/AG08BY.f @@ -266,7 +266,8 @@ SUBROUTINE AG08BY( FIRST, N, M, P, SVLMAX, ABCD, LDABCD, E, LDE, EXTERNAL DLAMCH, DNRM2, IDAMAX, ILAENV C .. External Subroutines .. EXTERNAL DCOPY, DLAIC1, DLAPMT, DLARFG, DLARTG, DLASET, - $ DLATZM, DORMQR, DROT, DSWAP, MB03OY, XERBLA + $ SLCT_DLATZM, DORMQR, DROT, DSWAP, MB03OY, + $ XERBLA C .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT C .. Executable Statements .. @@ -398,7 +399,8 @@ SUBROUTINE AG08BY( FIRST, N, M, P, SVLMAX, ABCD, LDABCD, E, LDE, IROW = IROW + 1 CALL DLARFG( RO+1, ABCD(IROW,ICOL), ABCD(IROW+1,ICOL), 1, $ T ) - CALL DLATZM( 'L', RO+1, MNR-ICOL, ABCD(IROW+1,ICOL), 1, T, + CALL SLCT_DLATZM( 'L', RO+1, MNR-ICOL, ABCD(IROW+1,ICOL), + $ 1, T, $ ABCD(IROW,ICOL+1), ABCD(IROW+1,ICOL+1), $ LDABCD, DWORK ) CALL DCOPY( PR-ICOL, DUM, 0, ABCD(IROW+1,ICOL), 1 ) diff --git a/slycot/src/AG8BYZ.f b/slycot/src/AG8BYZ.f index c2dc7d5e..4f08c12e 100755 --- a/slycot/src/AG8BYZ.f +++ b/slycot/src/AG8BYZ.f @@ -274,7 +274,8 @@ SUBROUTINE AG8BYZ( FIRST, N, M, P, SVLMAX, ABCD, LDABCD, E, LDE, EXTERNAL DLAMCH, DZNRM2, IDAMAX, ILAENV C .. External Subroutines .. EXTERNAL MB3OYZ, XERBLA, ZCOPY, ZLAIC1, ZLAPMT, ZLARFG, - $ ZLARTG, ZLASET, ZLATZM, ZROT, ZSWAP, ZUNMQR + $ ZLARTG, ZLASET, SLCT_ZLATZM, ZROT, ZSWAP, + $ ZUNMQR C .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCONJG, INT, MAX, MIN, SQRT C .. Executable Statements .. @@ -406,7 +407,8 @@ SUBROUTINE AG8BYZ( FIRST, N, M, P, SVLMAX, ABCD, LDABCD, E, LDE, IROW = IROW + 1 CALL ZLARFG( RO+1, ABCD(IROW,ICOL), ABCD(IROW+1,ICOL), 1, $ TC ) - CALL ZLATZM( 'L', RO+1, MNR-ICOL, ABCD(IROW+1,ICOL), 1, +C RvP, replaced by slicot replacement for obsolete lapack routine + CALL SLCT_ZLATZM( 'L', RO+1, MNR-ICOL, ABCD(IROW+1,ICOL), 1, $ DCONJG( TC ), ABCD(IROW,ICOL+1), $ ABCD(IROW+1,ICOL+1), LDABCD, ZWORK ) CALL ZCOPY( PR-ICOL, DUM, 0, ABCD(IROW+1,ICOL), 1 ) diff --git a/slycot/src/MB03ND.f b/slycot/src/MB03ND.f index c681c2e5..a6cbbd52 100644 --- a/slycot/src/MB03ND.f +++ b/slycot/src/MB03ND.f @@ -182,6 +182,7 @@ INTEGER FUNCTION MB03ND( N, THETA, Q2, E2, PIVMIN, INFO ) C C Error return. C + MB03ND = ZERO CALL XERBLA( 'MB03ND', -INFO ) RETURN END IF diff --git a/slycot/src/MB03NY.f b/slycot/src/MB03NY.f index a6efae58..e507ea29 100644 --- a/slycot/src/MB03NY.f +++ b/slycot/src/MB03NY.f @@ -155,6 +155,7 @@ DOUBLE PRECISION FUNCTION MB03NY( N, OMEGA, A, LDA, S, DWORK, C C Error return. C + MB03NY = ZERO CALL XERBLA( 'MB03NY', -INFO ) RETURN END IF @@ -177,6 +178,7 @@ DOUBLE PRECISION FUNCTION MB03NY( N, OMEGA, A, LDA, S, DWORK, $ 1, DUMMY, 1, DWORK, LDWORK, INFO ) IF ( INFO.NE.0 ) THEN INFO = 2 + MB03NY = ZERO RETURN END IF ELSE @@ -196,6 +198,7 @@ DOUBLE PRECISION FUNCTION MB03NY( N, OMEGA, A, LDA, S, DWORK, $ DWORK, INFO ) IF ( INFO.NE.0 ) THEN INFO = 2 + MB03NY = ZERO RETURN END IF CWORK(1) = CWORK(N*N+1) + DBLE( N*N ) * CONE diff --git a/slycot/src/SB01BY.f b/slycot/src/SB01BY.f index 58b48013..592161f2 100644 --- a/slycot/src/SB01BY.f +++ b/slycot/src/SB01BY.f @@ -122,7 +122,7 @@ SUBROUTINE SB01BY( N, M, S, P, A, B, F, TOL, DWORK, INFO ) DOUBLE PRECISION DLAMC3, DLAMCH EXTERNAL DLAMC3, DLAMCH C .. External Subroutines .. - EXTERNAL DLANV2, DLARFG, DLASET, DLASV2, DLATZM, DROT + EXTERNAL DLANV2, DLARFG, DLASET, DLASV2, SLCT_DLATZM, DROT C .. Intrinsic Functions .. INTRINSIC ABS, MIN C .. Executable Statements .. @@ -148,7 +148,8 @@ SUBROUTINE SB01BY( N, M, S, P, A, B, F, TOL, DWORK, INFO ) F(1,1) = ( S - A(1,1) )/B1 IF( M.GT.1 ) THEN CALL DLASET( 'Full', M-1, 1, ZERO, ZERO, F(2,1), M ) - CALL DLATZM( 'Left', M, N, B(1,2), N, TAU1, F(1,1), F(2,1), + CALL SLCT_DLATZM( 'Left', M, N, B(1,2), N, TAU1, + $ F(1,1), F(2,1), $ M, DWORK ) END IF RETURN @@ -185,7 +186,8 @@ SUBROUTINE SB01BY( N, M, S, P, A, B, F, TOL, DWORK, INFO ) C and H2. C CALL DLARFG( M, B(1,1), B(1,2), N, TAU1 ) - CALL DLATZM( 'Right', N-1, M, B(1,2), N, TAU1, B(2,1), B(2,2), + CALL SLCT_DLATZM( 'Right', N-1, M, B(1,2), N, TAU1, + $ B(2,1), B(2,2), $ N, DWORK ) B1 = B(1,1) B21 = B(2,1) @@ -322,10 +324,10 @@ SUBROUTINE SB01BY( N, M, S, P, A, B, F, TOL, DWORK, INFO ) C Compute H1*H2*F. C IF( M.GT.2 ) - $ CALL DLATZM( 'Left', M-1, N, B(2,3), N, TAU2, F(2,1), F(3,1), - $ M, DWORK ) - CALL DLATZM( 'Left', M, N, B(1,2), N, TAU1, F(1,1), F(2,1), M, - $ DWORK ) + $ CALL SLCT_DLATZM( 'Left', M-1, N, B(2,3), N, TAU2, + $ F(2,1), F(3,1), M, DWORK ) + CALL SLCT_DLATZM( 'Left', M, N, B(1,2), N, TAU1, + $ F(1,1), F(2,1), M, DWORK ) C RETURN C *** Last line of SB01BY *** diff --git a/slycot/src/SB01FY.f b/slycot/src/SB01FY.f index 20a716ba..5cbb1ac6 100644 --- a/slycot/src/SB01FY.f +++ b/slycot/src/SB01FY.f @@ -128,8 +128,8 @@ SUBROUTINE SB01FY( DISCR, N, M, A, LDA, B, LDB, F, LDF, V, LDV, DOUBLE PRECISION DLAPY2, DLAPY3 EXTERNAL DLAPY2, DLAPY3 C .. External Subroutines .. - EXTERNAL DLARFG, DLASET, DLATZM, DROTG, DTRTRI, MA02AD, - $ MB04OX, SB03OY + EXTERNAL DLARFG, DLASET, SLCT_DLATZM, DROTG, DTRTRI, + $ MA02AD, MB04OX, SB03OY C .. Intrinsic Functions .. INTRINSIC ABS, SQRT C .. Executable Statements .. @@ -180,7 +180,7 @@ SUBROUTINE SB01FY( DISCR, N, M, A, LDA, B, LDB, F, LDF, V, LDV, C IF( M.GT.1 ) THEN CALL DLARFG( M, F(1,1), F(2,1), 1, TEMP ) - CALL DLATZM( 'Left', M, N-1, F(2,1), 1, TEMP, F(1,2), + CALL SLCT_DLATZM( 'Left', M, N-1, F(2,1), 1, TEMP, F(1,2), $ F(2,2), LDF, V ) END IF R11 = F(1,1) diff --git a/slycot/src/SB03MD.f b/slycot/src/SB03MD.f index 8097795d..f11b3ff4 100644 --- a/slycot/src/SB03MD.f +++ b/slycot/src/SB03MD.f @@ -1,4 +1,4 @@ - SUBROUTINE SB03MD( DICO, JOB, FACT, TRANA, N, C, LDC, A, LDA, U, + SUBROUTINE SB03MD( DICO, JOB, FACT, TRANA, N, C, LDC, A, LDA, U, $ LDU, SCALE, SEP, FERR, WR, WI, IWORK, DWORK, $ LDWORK, INFO ) C diff --git a/slycot/src/SB04OD.f b/slycot/src/SB04OD.f index 6a11ffa7..929b1d65 100644 --- a/slycot/src/SB04OD.f +++ b/slycot/src/SB04OD.f @@ -405,7 +405,7 @@ SUBROUTINE SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C, LOGICAL ILASCL, ILBSCL, ILDSCL, ILESCL, LJOB1, LJOB2, $ LJOBD, LJOBDF, LJOBF, LREDRA, LREDRB, LREDUA, $ LREDUB, LREDUC, LREDUR, LTRANN, SUFWRK - INTEGER I, IERR, IJOB, MINWRK, MN, WRKOPT + INTEGER I, IERR, IJOB, MINWRK, MN, WRKOPT, SDIM DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, DNRM, $ DNRMTO, ENRM, ENRMTO, SAFMAX, SAFMIN, SMLNUM C .. External Functions .. @@ -413,7 +413,7 @@ SUBROUTINE SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C, DOUBLE PRECISION DLAMCH, DLANGE EXTERNAL DLAMCH, DLANGE, LSAME C .. External Subroutines .. - EXTERNAL DCOPY, DGEGS, DGEMM, DGEMV, DLABAD, DLACPY, + EXTERNAL DCOPY, DGGES, DGEMM, DGEMV, DLABAD, DLACPY, $ DLASCL, DTGSYL, XERBLA C .. Intrinsic Functions .. INTRINSIC INT, MAX, SQRT @@ -563,9 +563,13 @@ SUBROUTINE SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C, C Workspace: need 7*M; C prefer 5*M + M*(NB+1). C - CALL DGEGS( 'Vectors left', 'Vectors right', M, A, LDA, D, - $ LDD, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q, - $ LDQ, DWORK(3*M+1), LDWORK-3*M, INFO ) +C CALL DGEGS( 'Vectors left', 'Vectors right', M, A, LDA, D, +C $ LDD, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q, +C $ LDQ, DWORK(3*M+1), LDWORK-3*M, INFO ) + CALL DGGES( 'Vectors left', 'Vectors right', 'N', 0, N, A, LDA, + $ D, LDD, SDIM, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q, + $ LDQ, DWORK(3*M+1), LDWORK-3*M, 0, INFO ) + C C Undo scaling C @@ -619,9 +623,12 @@ SUBROUTINE SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C, C Workspace: need 7*N; C prefer 5*N + N*(NB+1). C - CALL DGEGS( 'Vectors left', 'Vectors right', N, B, LDB, E, - $ LDE, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V, - $ LDV, DWORK(3*N+1), LDWORK-3*N, INFO ) +C CALL DGEGS( 'Vectors left', 'Vectors right', N, B, LDB, E, +C $ LDE, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V, +C $ LDV, DWORK(3*N+1), LDWORK-3*N, INFO ) + CALL DGGES( 'Vectors left', 'Vectors right', 'N', 0, N, B, LDB, E, + $ LDE, SDIM, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V, + $ LDV, DWORK(3*N+1), LDWORK-3*N, 0, INFO ) C C Undo scaling C diff --git a/slycot/src/SB06ND.f b/slycot/src/SB06ND.f index 3ea98637..a3774439 100644 --- a/slycot/src/SB06ND.f +++ b/slycot/src/SB06ND.f @@ -183,8 +183,8 @@ SUBROUTINE SB06ND( N, M, KMAX, A, LDA, B, LDB, KSTAIR, U, LDU, F, INTEGER J, J0, JCUR, JKCUR, JMKCUR, KCUR, KK, KMIN, $ KSTEP, MKCUR, NCONT C .. External Subroutines .. - EXTERNAL DCOPY, DGEMM, DLACPY, DLARFG, DLASET, DLATZM, - $ DTRSM, XERBLA + EXTERNAL DCOPY, DGEMM, DLACPY, DLARFG, DLASET, + $ SLCT_DLATZM, DTRSM, XERBLA C .. Executable Statements .. C INFO = 0 @@ -250,11 +250,11 @@ SUBROUTINE SB06ND( N, M, KMAX, A, LDA, B, LDB, KSTAIR, U, LDU, F, C C Backmultiply A and U with Ukk. C - CALL DLATZM( 'Right', JCUR-1, KCUR+1, F(1,JCUR), 1, + CALL SLCT_DLATZM( 'Right', JCUR-1, KCUR+1, F(1,JCUR), 1, $ DWORK(JCUR), A(1,JCUR), A(1,JMKCUR), LDA, $ DWORK ) C - CALL DLATZM( 'Right', N, KCUR+1, F(1,JCUR), 1, + CALL SLCT_DLATZM( 'Right', N, KCUR+1, F(1,JCUR), 1, $ DWORK(JCUR), U(1,JCUR), U(1,JMKCUR), LDU, $ DWORK(N+1) ) JCUR = JCUR - 1 @@ -291,7 +291,8 @@ SUBROUTINE SB06ND( N, M, KMAX, A, LDA, B, LDB, KSTAIR, U, LDU, F, JCUR = JKCUR - KCUR C DO 60 J = 1, KCUR - CALL DLATZM( 'Left', KCUR+1, N-JCUR+1, F(1,JKCUR), 1, + CALL SLCT_DLATZM( 'Left', KCUR+1, N-JCUR+1, + $ F(1,JKCUR), 1, $ DWORK(JKCUR), A(JKCUR,JCUR), $ A(JCUR,JCUR), LDA, DWORK(N+1) ) JCUR = JCUR - 1 @@ -306,7 +307,7 @@ SUBROUTINE SB06ND( N, M, KMAX, A, LDA, B, LDB, KSTAIR, U, LDU, F, JKCUR = JCUR + KCUR C DO 100 J = M, M - KCUR + 1, -1 - CALL DLATZM( 'Left', KCUR+1, M-J+1, F(1,JKCUR), 1, + CALL SLCT_DLATZM( 'Left', KCUR+1, M-J+1, F(1,JKCUR), 1, $ DWORK(JKCUR), B(JKCUR,J), B(JCUR,J), LDB, $ DWORK(N+1) ) JCUR = JCUR - 1 diff --git a/slycot/src/SG03AD.f b/slycot/src/SG03AD.f index a08e218c..92b5eaf1 100644 --- a/slycot/src/SG03AD.f +++ b/slycot/src/SG03AD.f @@ -197,9 +197,9 @@ SUBROUTINE SG03AD( DICO, JOB, FACT, TRANS, UPLO, N, A, LDA, E, C JOB FACT | LDWORK C -------------------+------------------- C 'X' 'F' | MAX(1,N) -C 'X' 'N' | MAX(1,4*N) +C 'X' 'N' | MAX(1,4*N,8*N+16) C 'B', 'S' 'F' | MAX(1,2*N**2) -C 'B', 'S' 'N' | MAX(1,2*N**2,4*N) +C 'B', 'S' 'N' | MAX(1,2*N**2,4*N,8*N+16) C C For optimum performance, LDWORK should be larger. C @@ -214,7 +214,7 @@ SUBROUTINE SG03AD( DICO, JOB, FACT, TRANS, UPLO, N, A, LDA, E, C quasitriangular form; C = 2: FACT = 'N' and the pencil A - lambda * E cannot be C reduced to generalized Schur form: LAPACK routine -C DGEGS has failed to converge; +C DGGES has failed to converge; C = 3: DICO = 'D' and the pencil A - lambda * E has a C pair of reciprocal eigenvalues. That is, lambda_i = C 1/lambda_j for some i and j, where lambda_i and @@ -392,7 +392,7 @@ SUBROUTINE SG03AD( DICO, JOB, FACT, TRANS, UPLO, N, A, LDA, E, LOGICAL LSAME EXTERNAL DLAMCH, DNRM2, LSAME C .. External Subroutines .. - EXTERNAL DCOPY, DGEGS, DLACON, MB01RD, MB01RW, SG03AX, + EXTERNAL DCOPY, DGGES, DLACON, MB01RD, MB01RW, SG03AX, $ SG03AY, XERBLA C .. Intrinsic Functions .. INTRINSIC DBLE, INT, MAX, MIN @@ -443,13 +443,13 @@ SUBROUTINE SG03AD( DICO, JOB, FACT, TRANS, UPLO, N, A, LDA, E, IF ( ISFACT ) THEN MINWRK = MAX( N, 1 ) ELSE - MINWRK = MAX( 4*N, 1 ) + MINWRK = MAX( 8*N+16, 1 ) END IF ELSE IF ( ISFACT ) THEN MINWRK = MAX( 2*N*N, 1 ) ELSE - MINWRK = MAX( 2*N*N, 4*N, 1 ) + MINWRK = MAX( 2*N*N, 8*N+16, 1 ) END IF END IF IF ( MINWRK .GT. LDWORK ) THEN @@ -492,9 +492,10 @@ SUBROUTINE SG03AD( DICO, JOB, FACT, TRANS, UPLO, N, A, LDA, E, C C ( Workspace: >= MAX(1,4*N) ) C - CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR, + CALL DGGES( 'Vectors', 'Vectors', 'N', 0, N, A, LDA, E, LDE, + $ SDIM, ALPHAR, $ ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK, - $ INFO1 ) + $ 0, INFO1 ) IF ( INFO1 .NE. 0 ) THEN INFO = 2 RETURN diff --git a/slycot/src/SG03BD.f b/slycot/src/SG03BD.f index 6bcd7400..36579be1 100644 --- a/slycot/src/SG03BD.f +++ b/slycot/src/SG03BD.f @@ -461,9 +461,10 @@ SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q, C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*), $ BETA(*), DWORK(*), E(LDE,*), Q(LDQ,*), Z(LDZ,*) + LOGICAL BWORK C .. Local Scalars .. DOUBLE PRECISION S1, S2, SAFMIN, WI, WR1, WR2 - INTEGER I, INFO1, MINMN, MINWRK, OPTWRK + INTEGER I, INFO1, MINMN, MINWRK, OPTWRK, SDIM LOGICAL ISDISC, ISFACT, ISTRAN C .. Local Arrays .. DOUBLE PRECISION E1(2,2) @@ -472,7 +473,7 @@ SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q, LOGICAL LSAME EXTERNAL DLAMCH, DLAPY2, LSAME C .. External Subroutines .. - EXTERNAL DCOPY, DGEGS, DGEMM, DGEMV, DGEQRF, DGERQF, + EXTERNAL DCOPY, DGGES, DGEMM, DGEMV, DGEQRF, DGERQF, $ DLACPY, DLAG2, DLASET, DSCAL, DTRMM, SG03BU, $ SG03BV, XERBLA C .. Intrinsic Functions .. @@ -559,9 +560,12 @@ SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q, C C ( Workspace: >= MAX(1,4*N) ) C - CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR, - $ ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK, - $ INFO1 ) +C CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR, +C $ ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK, +C $ INFO1 ) + CALL DGGES( 'Vectors', 'Vectors', 'N', 0, N, A, LDA, + $ E, LDE, SDIM, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, + $ DWORK, LDWORK, 0, INFO) IF ( INFO1 .NE. 0 ) THEN INFO = 4 RETURN diff --git a/slycot/src/SLCT_DLATZM.f b/slycot/src/SLCT_DLATZM.f new file mode 100644 index 00000000..06f140c2 --- /dev/null +++ b/slycot/src/SLCT_DLATZM.f @@ -0,0 +1,223 @@ +*> \brief \b DLATZM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLATZM + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE SLCT_DLATZM( SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK ) +* +* .. Scalar Arguments .. +* CHARACTER SIDE +* INTEGER INCV, LDC, M, N +* DOUBLE PRECISION TAU +* .. +* .. Array Arguments .. +* DOUBLE PRECISION C1( LDC, * ), C2( LDC, * ), V( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> This routine is deprecated and has been replaced by routine DORMRZ. +*> For compatibility, this copy re-named and included in SLYCOT +*> +*> DLATZM applies a Householder matrix generated by DTZRQF to a matrix. +*> +*> Let P = I - tau*u*u**T, u = ( 1 ), +*> ( v ) +*> where v is an (m-1) vector if SIDE = 'L', or a (n-1) vector if +*> SIDE = 'R'. +*> +*> If SIDE equals 'L', let +*> C = [ C1 ] 1 +*> [ C2 ] m-1 +*> n +*> Then C is overwritten by P*C. +*> +*> If SIDE equals 'R', let +*> C = [ C1, C2 ] m +*> 1 n-1 +*> Then C is overwritten by C*P. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': form P * C +*> = 'R': form C * P +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension +*> (1 + (M-1)*abs(INCV)) if SIDE = 'L' +*> (1 + (N-1)*abs(INCV)) if SIDE = 'R' +*> The vector v in the representation of P. V is not used +*> if TAU = 0. +*> \endverbatim +*> +*> \param[in] INCV +*> \verbatim +*> INCV is INTEGER +*> The increment between elements of v. INCV <> 0 +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is DOUBLE PRECISION +*> The value tau in the representation of P. +*> \endverbatim +*> +*> \param[in,out] C1 +*> \verbatim +*> C1 is DOUBLE PRECISION array, dimension +*> (LDC,N) if SIDE = 'L' +*> (M,1) if SIDE = 'R' +*> On entry, the n-vector C1 if SIDE = 'L', or the m-vector C1 +*> if SIDE = 'R'. +*> +*> On exit, the first row of P*C if SIDE = 'L', or the first +*> column of C*P if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in,out] C2 +*> \verbatim +*> C2 is DOUBLE PRECISION array, dimension +*> (LDC, N) if SIDE = 'L' +*> (LDC, N-1) if SIDE = 'R' +*> On entry, the (m - 1) x n matrix C2 if SIDE = 'L', or the +*> m x (n - 1) matrix C2 if SIDE = 'R'. +*> +*> On exit, rows 2:m of P*C if SIDE = 'L', or columns 2:m of C*P +*> if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the arrays C1 and C2. LDC >= (1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension +*> (N) if SIDE = 'L' +*> (M) if SIDE = 'R' +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleOTHERcomputational +* +* ===================================================================== + SUBROUTINE SLCT_DLATZM( SIDE, M, N, V, INCV, TAU, + & C1, C2, LDC, WORK ) +* +* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER SIDE + INTEGER INCV, LDC, M, N + DOUBLE PRECISION TAU +* .. +* .. Array Arguments .. + DOUBLE PRECISION C1( LDC, * ), C2( LDC, * ), V( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. External Subroutines .. + EXTERNAL DAXPY, DCOPY, DGEMV, DGER +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Intrinsic Functions .. + INTRINSIC MIN +* .. +* .. Executable Statements .. +* + IF( ( MIN( M, N ).EQ.0 ) .OR. ( TAU.EQ.ZERO ) ) + $ RETURN +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* w := (C1 + v**T * C2)**T +* + CALL DCOPY( N, C1, LDC, WORK, 1 ) + CALL DGEMV( 'Transpose', M-1, N, ONE, C2, LDC, V, INCV, ONE, + $ WORK, 1 ) +* +* [ C1 ] := [ C1 ] - tau* [ 1 ] * w**T +* [ C2 ] [ C2 ] [ v ] +* + CALL DAXPY( N, -TAU, WORK, 1, C1, LDC ) + CALL DGER( M-1, N, -TAU, V, INCV, WORK, 1, C2, LDC ) +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* w := C1 + C2 * v +* + CALL DCOPY( M, C1, 1, WORK, 1 ) + CALL DGEMV( 'No transpose', M, N-1, ONE, C2, LDC, V, INCV, ONE, + $ WORK, 1 ) +* +* [ C1, C2 ] := [ C1, C2 ] - tau* w * [ 1 , v**T] +* + CALL DAXPY( M, -TAU, WORK, 1, C1, 1 ) + CALL DGER( M, N-1, -TAU, WORK, 1, V, INCV, C2, LDC ) + END IF +* + RETURN +* +* End of DLATZM +* + END diff --git a/slycot/src/SLCT_ZLATZM.f b/slycot/src/SLCT_ZLATZM.f new file mode 100644 index 00000000..5768ed69 --- /dev/null +++ b/slycot/src/SLCT_ZLATZM.f @@ -0,0 +1,226 @@ +*> \brief \b ZLATZM +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLATZM + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLATZM( SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK ) +* +* .. Scalar Arguments .. +* CHARACTER SIDE +* INTEGER INCV, LDC, M, N +* COMPLEX*16 TAU +* .. +* .. Array Arguments .. +* COMPLEX*16 C1( LDC, * ), C2( LDC, * ), V( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> This routine is deprecated and has been replaced by routine ZUNMRZ. +*> +*> ZLATZM applies a Householder matrix generated by ZTZRQF to a matrix. +*> +*> Let P = I - tau*u*u**H, u = ( 1 ), +*> ( v ) +*> where v is an (m-1) vector if SIDE = 'L', or a (n-1) vector if +*> SIDE = 'R'. +*> +*> If SIDE equals 'L', let +*> C = [ C1 ] 1 +*> [ C2 ] m-1 +*> n +*> Then C is overwritten by P*C. +*> +*> If SIDE equals 'R', let +*> C = [ C1, C2 ] m +*> 1 n-1 +*> Then C is overwritten by C*P. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': form P * C +*> = 'R': form C * P +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is COMPLEX*16 array, dimension +*> (1 + (M-1)*abs(INCV)) if SIDE = 'L' +*> (1 + (N-1)*abs(INCV)) if SIDE = 'R' +*> The vector v in the representation of P. V is not used +*> if TAU = 0. +*> \endverbatim +*> +*> \param[in] INCV +*> \verbatim +*> INCV is INTEGER +*> The increment between elements of v. INCV <> 0 +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is COMPLEX*16 +*> The value tau in the representation of P. +*> \endverbatim +*> +*> \param[in,out] C1 +*> \verbatim +*> C1 is COMPLEX*16 array, dimension +*> (LDC,N) if SIDE = 'L' +*> (M,1) if SIDE = 'R' +*> On entry, the n-vector C1 if SIDE = 'L', or the m-vector C1 +*> if SIDE = 'R'. +*> +*> On exit, the first row of P*C if SIDE = 'L', or the first +*> column of C*P if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in,out] C2 +*> \verbatim +*> C2 is COMPLEX*16 array, dimension +*> (LDC, N) if SIDE = 'L' +*> (LDC, N-1) if SIDE = 'R' +*> On entry, the (m - 1) x n matrix C2 if SIDE = 'L', or the +*> m x (n - 1) matrix C2 if SIDE = 'R'. +*> +*> On exit, rows 2:m of P*C if SIDE = 'L', or columns 2:m of C*P +*> if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the arrays C1 and C2. +*> LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension +*> (N) if SIDE = 'L' +*> (M) if SIDE = 'R' +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complex16OTHERcomputational +* +* ===================================================================== + SUBROUTINE SLCT_ZLATZM( SIDE, M, N, V, INCV, TAU, C1, C2, LDC, + $ WORK ) +* +* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER SIDE + INTEGER INCV, LDC, M, N + COMPLEX*16 TAU +* .. +* .. Array Arguments .. + COMPLEX*16 C1( LDC, * ), C2( LDC, * ), V( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE, ZERO + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), + $ ZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. External Subroutines .. + EXTERNAL ZAXPY, ZCOPY, ZGEMV, ZGERC, ZGERU, ZLACGV +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Intrinsic Functions .. + INTRINSIC MIN +* .. +* .. Executable Statements .. +* + IF( ( MIN( M, N ).EQ.0 ) .OR. ( TAU.EQ.ZERO ) ) + $ RETURN +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* w := ( C1 + v**H * C2 )**H +* + CALL ZCOPY( N, C1, LDC, WORK, 1 ) + CALL ZLACGV( N, WORK, 1 ) + CALL ZGEMV( 'Conjugate transpose', M-1, N, ONE, C2, LDC, V, + $ INCV, ONE, WORK, 1 ) +* +* [ C1 ] := [ C1 ] - tau* [ 1 ] * w**H +* [ C2 ] [ C2 ] [ v ] +* + CALL ZLACGV( N, WORK, 1 ) + CALL ZAXPY( N, -TAU, WORK, 1, C1, LDC ) + CALL ZGERU( M-1, N, -TAU, V, INCV, WORK, 1, C2, LDC ) +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* w := C1 + C2 * v +* + CALL ZCOPY( M, C1, 1, WORK, 1 ) + CALL ZGEMV( 'No transpose', M, N-1, ONE, C2, LDC, V, INCV, ONE, + $ WORK, 1 ) +* +* [ C1, C2 ] := [ C1, C2 ] - tau* w * [ 1 , v**H] +* + CALL ZAXPY( M, -TAU, WORK, 1, C1, 1 ) + CALL ZGERC( M, N-1, -TAU, WORK, 1, V, INCV, C2, LDC ) + END IF +* + RETURN +* +* End of ZLATZM +* + END diff --git a/slycot/src/TB01MD.f b/slycot/src/TB01MD.f index b63aacee..b446e719 100644 --- a/slycot/src/TB01MD.f +++ b/slycot/src/TB01MD.f @@ -188,7 +188,7 @@ SUBROUTINE TB01MD( JOBU, UPLO, N, M, A, LDA, B, LDB, U, LDU, LOGICAL LSAME EXTERNAL LSAME C .. External Subroutines .. - EXTERNAL DLARFG, DLASET, DLATZM, XERBLA + EXTERNAL DLARFG, DLASET, SLCT_DLATZM, XERBLA C .. Intrinsic Functions .. INTRINSIC MAX, MIN C .. Executable Statements .. @@ -262,16 +262,16 @@ SUBROUTINE TB01MD( JOBU, UPLO, N, M, A, LDA, B, LDB, U, LDU, C C Update A. C - CALL DLATZM( 'Left', NJ+1, N, B(PAR3,PAR1), 1, DZ, A(PAR2,1), - $ A(PAR3,1), LDA, DWORK ) - CALL DLATZM( 'Right', N, NJ+1, B(PAR3,PAR1), 1, DZ, A(1,PAR2), - $ A(1,PAR3), LDA, DWORK ) + CALL SLCT_DLATZM( 'Left', NJ+1, N, B(PAR3,PAR1), 1, DZ, + $ A(PAR2,1), A(PAR3,1), LDA, DWORK ) + CALL SLCT_DLATZM( 'Right', N, NJ+1, B(PAR3,PAR1), 1, DZ, + $ A(1,PAR2), A(1,PAR3), LDA, DWORK ) C IF ( LJOBA ) THEN C C Update U. C - CALL DLATZM( 'Right', N, NJ+1, B(PAR3,PAR1), 1, DZ, + CALL SLCT_DLATZM( 'Right', N, NJ+1, B(PAR3,PAR1), 1, DZ, $ U(1,PAR2), U(1,PAR3), LDU, DWORK ) END IF C @@ -279,7 +279,8 @@ SUBROUTINE TB01MD( JOBU, UPLO, N, M, A, LDA, B, LDB, U, LDU, C C Update B C - CALL DLATZM( 'Left', NJ+1, PAR4-PAR3+1, B(PAR3,PAR1), 1, DZ, + CALL SLCT_DLATZM( 'Left', NJ+1, PAR4-PAR3+1, B(PAR3,PAR1), + $ 1, DZ, $ B(PAR2,PAR3), B(PAR3,PAR3), LDB, DWORK ) END IF C @@ -314,16 +315,17 @@ SUBROUTINE TB01MD( JOBU, UPLO, N, M, A, LDA, B, LDB, U, LDU, C C Update A. C - CALL DLATZM( 'Left', NJ+1, PAR6-PAR5+1, A(PAR3,PAR1), 1, DZ, + CALL SLCT_DLATZM( 'Left', NJ+1, PAR6-PAR5+1, A(PAR3,PAR1), + $ 1, DZ, $ A(PAR2,PAR5), A(PAR3,PAR5), LDA, DWORK ) - CALL DLATZM( 'Right', N, NJ+1, A(PAR3,PAR1), 1, DZ, + CALL SLCT_DLATZM( 'Right', N, NJ+1, A(PAR3,PAR1), 1, DZ, $ A(1,PAR2), A(1,PAR3), LDA, DWORK ) C IF ( LJOBA ) THEN C C Update U. C - CALL DLATZM( 'Right', N, NJ+1, A(PAR3,PAR1), 1, DZ, + CALL SLCT_DLATZM( 'Right', N, NJ+1, A(PAR3,PAR1), 1, DZ, $ U(1,PAR2), U(1,PAR3), LDU, DWORK ) END IF C diff --git a/slycot/src/TB01ND.f b/slycot/src/TB01ND.f index cc93dd3a..4bd4fbad 100644 --- a/slycot/src/TB01ND.f +++ b/slycot/src/TB01ND.f @@ -195,7 +195,7 @@ SUBROUTINE TB01ND( JOBU, UPLO, N, P, A, LDA, C, LDC, U, LDU, LOGICAL LSAME EXTERNAL LSAME C .. External Subroutines .. - EXTERNAL DLARFG, DLASET, DLATZM, XERBLA + EXTERNAL DLARFG, DLASET, SLCT_DLATZM, XERBLA C .. Intrinsic Functions .. INTRINSIC MAX, MIN C .. Executable Statements .. @@ -269,16 +269,16 @@ SUBROUTINE TB01ND( JOBU, UPLO, N, P, A, LDA, C, LDC, U, LDU, C C Update A. C - CALL DLATZM( 'Left', NJ+1, N, C(PAR1,PAR3), LDC, DZ, A(PAR2,1), - $ A(PAR3,1), LDA, DWORK ) - CALL DLATZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ, + CALL SLCT_DLATZM( 'Left', NJ+1, N, C(PAR1,PAR3), LDC, DZ, + $ A(PAR2,1), A(PAR3,1), LDA, DWORK ) + CALL SLCT_DLATZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ, $ A(1,PAR2), A(1,PAR3), LDA, DWORK ) C IF ( LJOBA ) THEN C C Update U. C - CALL DLATZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ, + CALL SLCT_DLATZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ, $ U(1,PAR2), U(1,PAR3), LDU, DWORK ) END IF C @@ -286,7 +286,8 @@ SUBROUTINE TB01ND( JOBU, UPLO, N, P, A, LDA, C, LDC, U, LDU, C C Update C. C - CALL DLATZM( 'Right', PAR4-PAR3+1, NJ+1, C(PAR1,PAR3), LDC, + CALL SLCT_DLATZM( 'Right', PAR4-PAR3+1, NJ+1, + $ C(PAR1,PAR3), LDC, $ DZ, C(PAR3,PAR2), C(PAR3,PAR3), LDC, DWORK ) END IF C @@ -323,16 +324,18 @@ SUBROUTINE TB01ND( JOBU, UPLO, N, P, A, LDA, C, LDC, U, LDU, C C Update A. C - CALL DLATZM( 'Left', NJ+1, N, A(PAR1,PAR3), LDA, DZ, + CALL SLCT_DLATZM( 'Left', NJ+1, N, A(PAR1,PAR3), LDA, DZ, $ A(PAR2,1), A(PAR3,1), LDA, DWORK ) - CALL DLATZM( 'Right', PAR6-PAR5+1, NJ+1, A(PAR1,PAR3), LDA, + CALL SLCT_DLATZM( 'Right', PAR6-PAR5+1, NJ+1, A(PAR1,PAR3), + $ LDA, $ DZ, A(PAR5,PAR2), A(PAR5,PAR3), LDA, DWORK ) C IF ( LJOBA ) THEN C C Update U. C - CALL DLATZM( 'Right', N, NJ+1, A(PAR1,PAR3), LDA, DZ, + CALL SLCT_DLATZM( 'Right', N, NJ+1, A(PAR1,PAR3), + $ LDA, DZ, $ U(1,PAR2), U(1,PAR3), LDU, DWORK ) END IF C diff --git a/slycot/src/synthesis.pyf b/slycot/src/synthesis.pyf index bc95229a..2fdb1308 100644 --- a/slycot/src/synthesis.pyf +++ b/slycot/src/synthesis.pyf @@ -550,7 +550,7 @@ subroutine sg03ad(dico,job,fact,trans,uplo,n,a,lda,e,lde,q,ldq,z,ldz,x,ldx,scale double precision intent(out), dimension(n), depend(n) :: beta integer intent(hide,cache),dimension(n*n),depend(n) :: iwork double precision intent(hide,cache),dimension(ldwork) :: dwork - integer optional,check(ldwork>=max(1,n)),depend(n) :: ldwork=max(1,max(2*n*n,4*n)) + integer optional,check(ldwork>=max(1,max(2*n*n,8*n+16))),depend(n) :: ldwork=max(1,max(2*n*n,8*n+16)) integer intent(out) :: info end subroutine sg03ad subroutine sg02ad_g(dico,jobb,fact,uplo,jobl,scal,sort,acc,n,m,p,a,lda,e,lde,b,ldb,q,ldq,r,ldr,l,ldl,rcondu,x,ldx,alfar,alfai,beta,s,lds,t,ldt,u,ldu,tol,iwork,dwork,ldwork,bwork,iwarn,info) ! in SG02AD.f @@ -593,7 +593,7 @@ subroutine sg02ad_g(dico,jobb,fact,uplo,jobl,scal,sort,acc,n,m,p,a,lda,e,lde,b,l double precision intent(hide) :: tol = 0 ! Unused for jobb = 'G' integer intent(hide,cache), dimension(2*n), depend(n) :: iwork ! Since jobb = 'G' double precision intent(hide,cache), dimension(ldwork) :: dwork - integer optional,check(ldwork>=max(7*(2*n+1)+16,16*n)),depend(n) :: ldwork=max(7*(2*n+1)+16,16*n) + integer optional,check(ldwork>=max(7*(2*n+1)+16,max(16*n,8*n+16))),depend(n) :: ldwork=max(7*(2*n+1)+16,max(16*n,8*n+16)) logical intent(hide,cache), dimension(2*n), depend(n) :: bwork integer intent(out) :: iwarn integer intent(out) :: info @@ -638,7 +638,7 @@ subroutine sg02ad_bn(dico,jobb,fact,uplo,jobl,scal,sort,acc,n,m,p,a,lda,e,lde,b, double precision optional :: tol = -1 integer intent(hide,cache), dimension(max(m,2*n)), depend(n,m) :: iwork ! Since jobb = 'B' double precision intent(hide,cache), dimension(ldwork) :: dwork - integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n)))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n))) + integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(16*n,8*n+16))))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(16*n,8*n+16)))) logical intent(hide,cache), dimension(2*n), depend(n) :: bwork integer intent(out) :: iwarn integer intent(out) :: info @@ -683,7 +683,7 @@ subroutine sg02ad_bc(dico,jobb,fact,uplo,jobl,scal,sort,acc,n,m,p,a,lda,e,lde,b, double precision optional :: tol = -1 integer intent(hide,cache), dimension(max(m,2*n)), depend(n,m) :: iwork ! Since jobb = 'B' double precision intent(hide,cache), dimension(ldwork) :: dwork - integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n)))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n))) + integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(16*n,8*n+16))))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(16*n,8*n+16)))) logical intent(hide,cache), dimension(2*n), depend(n) :: bwork integer intent(out) :: iwarn integer intent(out) :: info @@ -728,7 +728,7 @@ subroutine sg02ad_bd(dico,jobb,fact,uplo,jobl,scal,sort,acc,n,m,p,a,lda,e,lde,b, double precision optional :: tol = -1 integer intent(hide,cache), dimension(max(m,2*n)), depend(n,m) :: iwork ! Since jobb = 'B' double precision intent(hide,cache), dimension(ldwork) :: dwork - integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n)))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n))) + integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(16*n,8*n+16))))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(16*n,8*n+16)))) logical intent(hide,cache), dimension(2*n), depend(n) :: bwork integer intent(out) :: iwarn integer intent(out) :: info @@ -773,7 +773,7 @@ subroutine sg02ad_bb(dico,jobb,fact,uplo,jobl,scal,sort,acc,n,m,p,a,lda,e,lde,b, double precision optional :: tol = -1 integer intent(hide,cache), dimension(max(m,2*n)), depend(n) :: iwork ! Since jobb = 'B' double precision intent(hide,cache), dimension(ldwork) :: dwork - integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n)))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,16*n))) + integer optional,check(ldwork>=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(16*n,8*n+16))))),depend(n,m) :: ldwork=max(3*m,max(2*n+m,max(7*(2*n+1)+16,max(8*n+16,16*n)))) logical intent(hide,cache), dimension(2*n), depend(n) :: bwork integer intent(out) :: iwarn integer intent(out) :: info diff --git a/slycot/synthesis.py b/slycot/synthesis.py index 0ea2e5d7..75daa347 100644 --- a/slycot/synthesis.py +++ b/slycot/synthesis.py @@ -2064,11 +2064,11 @@ def sg03ad(dico,job,fact,trans,uplo,N,A,E,Q,Z,X,ldwork=None): if job == 'X' and fact == 'F': ldwork = max(1,N) elif job == 'X' and fact == 'N': - ldwork = max(1,4*N) + ldwork = max(1,8*N+16) elif (job == 'B' or job == 'S') and fact == 'F': ldwork = max(1,2*N**2) elif (job == 'B' or job == 'S') and fact == 'N': - ldwork = max(1,2*N**2,4*N) + ldwork = max(1,2*N**2,8*N+16) out = _wrapper.sg03ad(dico,job,fact,trans,uplo,N,A,E,Q,Z,X,ldwork) diff --git a/slycot/tests/test_sg02ad.py b/slycot/tests/test_sg02ad.py new file mode 100644 index 00000000..7fe0bf65 --- /dev/null +++ b/slycot/tests/test_sg02ad.py @@ -0,0 +1,47 @@ +#!/usr/bin/env python +# +# test_sg02ad.py - test suite for ricatti equation solving +# RvP, 19 Jun 2017 +from __future__ import print_function + +import unittest +from slycot import synthesis +import numpy as np +from numpy import linalg + +from numpy.testing import assert_raises, assert_almost_equal + +class test_sg03ad(unittest.TestCase): + + def test_sg02ad_case1(self): + n = 3 + m = 1 + # from a discussion here: + # https://github.com/scipy/scipy/issues/2251 + A = np.matrix([[ 0.63399379, 0.54906824, 0.76253406], + [ 0.5404729 , 0.53745766, 0.08731853], + [ 0.27524045, 0.84922129, 0.4681622 ]]) + B = np.matrix([[ 0.96861695],[ 0.05532739],[ 0.78934047]]) + Q = np.matrix(np.eye(3)) + E = np.matrix(np.eye(3)) + R = np.matrix(np.ones((1,1), dtype=float)) + S = np.matrix([[-2.67522766, -5.39447418, 2.19128542], + [-1.94918951, -3.15480639, 5.24379117], + [ 4.29133973, 8.10585767, -5.88895897]]) + L = np.matrix(np.zeros((3,1))) + rcondu, X, alphar, alphai, beta, S, T, U, iwarn = \ + synthesis.sg02ad('D', 'B', 'N', 'U', 'Z', 'N', 'S', 'R', + n, m, 1, + A, E, B, Q, R, L) + assert_almost_equal( + A.T*X*A - E.T*X*E - + (L + A.T*X*B) * np.linalg.solve (R+B.T*X*B, (L+A.T*X*B).T) + + Q, + np.zeros((n,n))) + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + + +if __name__ == "__main__": + unittest.main() diff --git a/slycot/tests/test_sg03ad.py b/slycot/tests/test_sg03ad.py new file mode 100644 index 00000000..40f9684c --- /dev/null +++ b/slycot/tests/test_sg03ad.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python +# +# test_sg03ad.py - test suit for stability margin commands +# RvP, 15 Jun 2017 +from __future__ import print_function + +import unittest +from slycot import synthesis +import numpy as np + +from numpy.testing import assert_raises, assert_almost_equal + +# test cases from +# http://www.qucosa.de/fileadmin/data/qucosa/documents/4168/data/b002.pdf + +class test_sg03ad(unittest.TestCase): + + def test_sg03ad_a(self): + # Example 1 + n = 100 + Xref = np.ones((n,n)) + U = np.tril(Xref) + for t in range(0, 50, 10): + A = np.matrix( + 2.0**(-t) - np.eye(n) + np.diag(range(1,n+1)) + U.T) + #print(A) + #print(t, n) + E = np.matrix(np.eye(n) + 2**(-t)*U) + Y = A.T*Xref*E + E.T*Xref*A + Q = np.zeros((n,n)) + Z = np.zeros((n,n)) + A, E, Q, Z, X, scale, sep, ferr, alphar, alphai, beta = \ + synthesis.sg03ad('C', 'B', 'N', 'N', 'L', n, A, E, Q, Z, Y) + assert_almost_equal(Xref, X) + + def test_sg03ad_3(self): + n = 3 + A = np.matrix([[3.0, 1.0, 1.0], + [1.0, 3.0, 0.0], + [1.0, 0.0, 2.0]]) + E = np.matrix([[1.0, 3.0, 0.0], + [3.0, 2.0, 1.0], + [1.0, 0.0, 1.0]]) + Y = np.matrix([[64.0, 73.0, 28.0], + [73.0, 70.0, 25.0], + [28.0, 25.0, 18.0]]) + Xref = np.array([[-2.0000, -1.0000, 0.0000], + [-1.0000, -3.0000, -1.0000], + [0.0000, -1.0000, -3.0000]]) + Q = np.matrix(np.zeros((3,3))) + Z = np.matrix(np.zeros((3,3))) + A, E, Q, Z, X, scale, sep, ferr, alphar, alphai, beta = \ + synthesis.sg03ad('C', 'B', 'N', 'N', 'L', n, A, E, Q, Z, -Y) + #print(A, E, Q, Z, X, scale, sep) + assert_almost_equal(X, Xref) + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestConvert) + + +if __name__ == "__main__": + unittest.main()