--- rpl/lapack/lapack/zgesvj.f 2016/08/27 15:34:47 1.3 +++ rpl/lapack/lapack/zgesvj.f 2017/06/17 10:54:11 1.4 @@ -1,26 +1,26 @@ -*> \brief \b ZGESVJ +*> \brief ZGESVJ * * =========== DOCUMENTATION =========== * -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ * *> \htmlonly -*> Download ZGESVJ + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> +*> Download ZGESVJ + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> *> [TXT] -*> \endhtmlonly +*> \endhtmlonly * * Definition: * =========== * * SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, * LDV, CWORK, LWORK, RWORK, LRWORK, INFO ) -* +* * .. Scalar Arguments .. * INTEGER INFO, LDA, LDV, LWORK, LRWORK, M, MV, N * CHARACTER*1 JOBA, JOBU, JOBV @@ -29,7 +29,7 @@ * COMPLEX*16 A( LDA, * ), V( LDV, * ), CWORK( LWORK ) * DOUBLE PRECISION RWORK( LRWORK ), SVA( N ) * .. -* +* * *> \par Purpose: * ============= @@ -64,11 +64,11 @@ *> JOBU is CHARACTER*1 *> Specifies whether to compute the left singular vectors *> (columns of U): -*> = 'U': The left singular vectors corresponding to the nonzero +*> = 'U' or 'F': The left singular vectors corresponding to the nonzero *> singular values are computed and returned in the leading *> columns of A. See more details in the description of A. *> The default numerical orthogonality threshold is set to -*> approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). +*> approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=DLAMCH('E'). *> = 'C': Analogous to JOBU='U', except that user can control the *> level of numerical orthogonality of the computed left *> singular vectors. TOL can be set to TOL = CTOL*EPS, where @@ -88,10 +88,10 @@ *> JOBV is CHARACTER*1 *> Specifies whether to compute the right singular vectors, that *> is, the matrix V: -*> = 'V' : the matrix V is computed and returned in the array V +*> = 'V' or 'J': the matrix V is computed and returned in the array V *> = 'A' : the Jacobi rotations are applied to the MV-by-N *> array V. In other words, the right singular vector -*> matrix V is not computed explicitly, instead it is +*> matrix V is not computed explicitly; instead it is *> applied to an MV-by-N matrix initially stored in the *> first MV rows of V. *> = 'N' : the matrix V is not computed and the array V is not @@ -101,7 +101,7 @@ *> \param[in] M *> \verbatim *> M is INTEGER -*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. +*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. *> \endverbatim *> *> \param[in] N @@ -206,8 +206,11 @@ *> *> \param[in,out] CWORK *> \verbatim -*> CWORK is COMPLEX*16 array, dimension M+N. -*> Used as work space. +*> CWORK is COMPLEX*16 array, dimension max(1,LWORK). +*> Used as workspace. +*> If on entry LWORK .EQ. -1, then a workspace query is assumed and +*> no computation is done; CWORK(1) is set to the minial (and optimal) +*> length of CWORK. *> \endverbatim *> *> \param[in] LWORK @@ -218,7 +221,7 @@ *> *> \param[in,out] RWORK *> \verbatim -*> RWORK is DOUBLE PRECISION array, dimension max(6,M+N). +*> RWORK is DOUBLE PRECISION array, dimension max(6,LRWORK). *> On entry, *> If JOBU .EQ. 'C' : *> RWORK(1) = CTOL, where CTOL defines the threshold for convergence. @@ -244,11 +247,14 @@ *> RWORK(6) = the largest absolute value over all sines of the *> Jacobi rotation angles in the last sweep. It can be *> useful for a post festum analysis. +*> If on entry LRWORK .EQ. -1, then a workspace query is assumed and +*> no computation is done; RWORK(1) is set to the minial (and optimal) +*> length of RWORK. *> \endverbatim *> *> \param[in] LRWORK *> \verbatim -*> LRWORK is INTEGER +*> LRWORK is INTEGER *> Length of RWORK, LRWORK >= MAX(6,N). *> \endverbatim *> @@ -257,22 +263,22 @@ *> INFO is INTEGER *> = 0 : successful exit. *> < 0 : if INFO = -i, then the i-th argument had an illegal value -*> > 0 : ZGESVJ did not converge in the maximal allowed number -*> (NSWEEP=30) of sweeps. The output may still be useful. +*> > 0 : ZGESVJ did not converge in the maximal allowed number +*> (NSWEEP=30) of sweeps. The output may still be useful. *> See the description of RWORK. *> \endverbatim *> * Authors: * ======== * -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. * *> \date June 2016 * -*> \ingroup doubleGEcomputational +*> \ingroup complex16GEcomputational * *> \par Further Details: * ===================== @@ -291,29 +297,30 @@ *> procedure is achieved if used in an accelerated version of Drmac and *> Veselic [4,5], and it is the kernel routine in the SIGMA library [6]. *> Some tunning parameters (marked with [TP]) are available for the -*> implementer. +*> implementer. *> The computational range for the nonzero singular values is the machine *> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even *> denormalized singular values can be computed with the corresponding *> gradual loss of accurate digits. *> \endverbatim * -*> \par Contributors: +*> \par Contributor: * ================== *> *> \verbatim *> *> ============ *> -*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) +*> Zlatko Drmac (Zagreb, Croatia) +*> *> \endverbatim * *> \par References: * ================ *> *> [1] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the -*> singular value decomposition on a vector computer. -*> SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. +*> singular value decomposition on a vector computer. +*> SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. *> [2] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. *> [3] Z. Drmac: Implementation of Jacobi rotations for accurate singular *> value computation in floating point arithmetic. @@ -329,8 +336,8 @@ *> Department of Mathematics, University of Zagreb, 2008, 2015. *> \endverbatim * -*> \par Bugs, examples and comments: -* ================================= +*> \par Bugs, examples and comments: +* ================================= *> *> \verbatim *> =========================== @@ -339,15 +346,15 @@ *> \endverbatim *> * ===================================================================== - SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, + SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, $ LDV, CWORK, LWORK, RWORK, LRWORK, INFO ) * -* -- LAPACK computational routine (version 3.6.1) -- +* -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * June 2016 * - IMPLICIT NONE + IMPLICIT NONE * .. Scalar Arguments .. INTEGER INFO, LDA, LDV, LWORK, LRWORK, M, MV, N CHARACTER*1 JOBA, JOBU, JOBV @@ -369,20 +376,19 @@ * .. * .. Local Scalars .. COMPLEX*16 AAPQ, OMPQ - DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG, - $ BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ, - $ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, - $ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, THSIGN, TOL + DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG, + $ BIGTHETA, CS, CTOL, EPSLN, MXAAPQ, + $ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, + $ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, THSIGN, TOL INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, - $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, + $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, $ N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND - LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK, + LOGICAL APPLV, GOSCALE, LOWER, LQUERY, LSVEC, NOSCALE, ROTOK, $ RSVEC, UCTOL, UPPER * .. * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DMAX1, DMIN1, DCONJG, DBLE, MIN0, MAX0, - $ DSIGN, DSQRT + INTRINSIC ABS, MAX, MIN, CONJG, DBLE, SIGN, SQRT * .. * .. External Functions .. * .. @@ -410,13 +416,14 @@ * * Test the input arguments * - LSVEC = LSAME( JOBU, 'U' ) + LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' ) UCTOL = LSAME( JOBU, 'C' ) - RSVEC = LSAME( JOBV, 'V' ) + RSVEC = LSAME( JOBV, 'V' ) .OR. LSAME( JOBV, 'J' ) APPLV = LSAME( JOBV, 'A' ) UPPER = LSAME( JOBA, 'U' ) LOWER = LSAME( JOBA, 'L' ) * + LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 ) IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN INFO = -1 ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN @@ -436,10 +443,10 @@ INFO = -11 ELSE IF( UCTOL .AND. ( RWORK( 1 ).LE.ONE ) ) THEN INFO = -12 - ELSE IF( LWORK.LT.( M+N ) ) THEN + ELSE IF( ( LWORK.LT.( M+N ) ) .AND. ( .NOT.LQUERY ) ) THEN INFO = -13 - ELSE IF( LRWORK.LT.MAX0( N, 6 ) ) THEN - INFO = -15 + ELSE IF( ( LRWORK.LT.MAX( N, 6 ) ) .AND. ( .NOT.LQUERY ) ) THEN + INFO = -15 ELSE INFO = 0 END IF @@ -448,6 +455,10 @@ IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZGESVJ', -INFO ) RETURN + ELSE IF ( LQUERY ) THEN + CWORK(1) = M + N + RWORK(1) = MAX( N, 6 ) + RETURN END IF * * #:) Quick return for void matrix @@ -467,27 +478,27 @@ ELSE * ... default IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN - CTOL = DSQRT( DBLE( M ) ) + CTOL = SQRT( DBLE( M ) ) ELSE CTOL = DBLE( M ) END IF END IF * ... and the machine dependent parameters are -*[!] (Make sure that DLAMCH() works properly on the target machine.) +*[!] (Make sure that SLAMCH() works properly on the target machine.) * EPSLN = DLAMCH( 'Epsilon' ) - ROOTEPS = DSQRT( EPSLN ) + ROOTEPS = SQRT( EPSLN ) SFMIN = DLAMCH( 'SafeMinimum' ) - ROOTSFMIN = DSQRT( SFMIN ) + ROOTSFMIN = SQRT( SFMIN ) SMALL = SFMIN / EPSLN BIG = DLAMCH( 'Overflow' ) * BIG = ONE / SFMIN ROOTBIG = ONE / ROOTSFMIN - LARGE = BIG / DSQRT( DBLE( M*N ) ) +* LARGE = BIG / SQRT( DBLE( M*N ) ) BIGTHETA = ONE / ROOTEPS * TOL = CTOL*EPSLN - ROOTTOL = DSQRT( TOL ) + ROOTTOL = SQRT( TOL ) * IF( DBLE( M )*EPSLN.GE.ONE ) THEN INFO = -4 @@ -514,7 +525,7 @@ * SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries * in A are detected, the procedure returns with INFO=-6. * - SKL = ONE / DSQRT( DBLE( M )*DBLE( N ) ) + SKL = ONE / SQRT( DBLE( M )*DBLE( N ) ) NOSCALE = .TRUE. GOSCALE = .TRUE. * @@ -529,7 +540,7 @@ CALL XERBLA( 'ZGESVJ', -INFO ) RETURN END IF - AAQQ = DSQRT( AAQQ ) + AAQQ = SQRT( AAQQ ) IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN SVA( p ) = AAPP*AAQQ ELSE @@ -554,7 +565,7 @@ CALL XERBLA( 'ZGESVJ', -INFO ) RETURN END IF - AAQQ = DSQRT( AAQQ ) + AAQQ = SQRT( AAQQ ) IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN SVA( p ) = AAPP*AAQQ ELSE @@ -579,7 +590,7 @@ CALL XERBLA( 'ZGESVJ', -INFO ) RETURN END IF - AAQQ = DSQRT( AAQQ ) + AAQQ = SQRT( AAQQ ) IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN SVA( p ) = AAPP*AAQQ ELSE @@ -604,8 +615,8 @@ AAPP = ZERO AAQQ = BIG DO 4781 p = 1, N - IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) ) - AAPP = DMAX1( AAPP, SVA( p ) ) + IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) ) + AAPP = MAX( AAPP, SVA( p ) ) 4781 CONTINUE * * #:) Quick return for zero matrix @@ -642,23 +653,23 @@ * Protect small singular values from underflow, and try to * avoid underflows/overflows in computing Jacobi rotations. * - SN = DSQRT( SFMIN / EPSLN ) - TEMP1 = DSQRT( BIG / DBLE( N ) ) - IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. + SN = SQRT( SFMIN / EPSLN ) + TEMP1 = SQRT( BIG / DBLE( N ) ) + IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. $ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN - TEMP1 = DMIN1( BIG, TEMP1 / AAPP ) + TEMP1 = MIN( BIG, TEMP1 / AAPP ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN - TEMP1 = DMIN1( SN / AAQQ, BIG / (AAPP*DSQRT( DBLE(N)) ) ) + TEMP1 = MIN( SN / AAQQ, BIG / (AAPP*SQRT( DBLE(N)) ) ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN - TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP ) + TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN - TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) ) + TEMP1 = MIN( SN / AAQQ, BIG / ( SQRT( DBLE( N ) )*AAPP ) ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE @@ -680,10 +691,10 @@ * EMPTSW = ( N*( N-1 ) ) / 2 NOTROT = 0 - + DO 1868 q = 1, N CWORK( q ) = CONE - 1868 CONTINUE + 1868 CONTINUE * * * @@ -695,7 +706,7 @@ * The boundaries are determined dynamically, based on the number of * pivots above a threshold. * - KBL = MIN0( 8, N ) + KBL = MIN( 8, N ) *[TP] KBL is a tuning parameter that defines the tile size in the * tiling of the p-q loops of pivot pairs. In general, an optimal * value of KBL depends on the matrix dimensions and on the @@ -707,7 +718,7 @@ BLSKIP = KBL**2 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL. * - ROWSKIP = MIN0( 5, KBL ) + ROWSKIP = MIN( 5, KBL ) *[TP] ROWSKIP is a tuning parameter. * LKAHEAD = 1 @@ -718,7 +729,7 @@ * invokes cubic convergence. Big part of this cycle is done inside * canonical subspaces of dimensions less than M. * - IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN + IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN *[TP] The number of partition levels and the actual partition are * tuning parameters. N4 = N / 4 @@ -816,18 +827,18 @@ * igl = ( ibr-1 )*KBL + 1 * - DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr ) + DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr ) * igl = igl + ir1*KBL * - DO 2001 p = igl, MIN0( igl+KBL-1, N-1 ) + DO 2001 p = igl, MIN( igl+KBL-1, N-1 ) * * .. de Rijk's pivoting * q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1 IF( p.NE.q ) THEN CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) - IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1, + IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1, $ V( 1, q ), 1 ) TEMP1 = SVA( p ) SVA( p ) = SVA( q ) @@ -851,14 +862,14 @@ * If properly implemented SCNRM2 is available, the IF-THEN-ELSE-END IF * below should be replaced with "AAPP = DZNRM2( M, A(1,p), 1 )". * - IF( ( SVA( p ).LT.ROOTBIG ) .AND. + IF( ( SVA( p ).LT.ROOTBIG ) .AND. $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN SVA( p ) = DZNRM2( M, A( 1, p ), 1 ) ELSE TEMP1 = ZERO AAPP = ONE CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) - SVA( p ) = TEMP1*DSQRT( AAPP ) + SVA( p ) = TEMP1*SQRT( AAPP ) END IF AAPP = SVA( p ) ELSE @@ -869,7 +880,7 @@ * PSKIPPED = 0 * - DO 2002 q = p + 1, MIN0( igl+KBL-1, N ) + DO 2002 q = p + 1, MIN( igl+KBL-1, N ) * AAQQ = SVA( q ) * @@ -879,12 +890,12 @@ IF( AAQQ.GE.ONE ) THEN ROTOK = ( SMALL*AAPP ).LE.AAQQ IF( AAPP.LT.( BIG / AAQQ ) ) THEN - AAPQ = ( ZDOTC( M, A( 1, p ), 1, + AAPQ = ( ZDOTC( M, A( 1, p ), 1, $ A( 1, q ), 1 ) / AAQQ ) / AAPP ELSE - CALL ZCOPY( M, A( 1, p ), 1, + CALL ZCOPY( M, A( 1, p ), 1, $ CWORK(N+1), 1 ) - CALL ZLASCL( 'G', 0, 0, AAPP, ONE, + CALL ZLASCL( 'G', 0, 0, AAPP, ONE, $ M, 1, CWORK(N+1), LDA, IERR ) AAPQ = ZDOTC( M, CWORK(N+1), 1, $ A( 1, q ), 1 ) / AAQQ @@ -892,10 +903,10 @@ ELSE ROTOK = AAPP.LE.( AAQQ / SMALL ) IF( AAPP.GT.( SMALL / AAQQ ) ) THEN - AAPQ = ( ZDOTC( M, A( 1, p ), 1, - $ A( 1, q ), 1 ) / AAQQ ) / AAPP + AAPQ = ( ZDOTC( M, A( 1, p ), 1, + $ A( 1, q ), 1 ) / AAPP ) / AAQQ ELSE - CALL ZCOPY( M, A( 1, q ), 1, + CALL ZCOPY( M, A( 1, q ), 1, $ CWORK(N+1), 1 ) CALL ZLASCL( 'G', 0, 0, AAQQ, $ ONE, M, 1, @@ -905,13 +916,15 @@ END IF END IF * -* AAPQ = AAPQ * DCONJG( CWORK(p) ) * CWORK(q) - AAPQ1 = -ABS(AAPQ) - MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 ) + +* AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q) + AAPQ1 = -ABS(AAPQ) + MXAAPQ = MAX( MXAAPQ, -AAPQ1 ) * * TO rotate or NOT to rotate, THAT is the question ... * IF( ABS( AAPQ1 ).GT.TOL ) THEN + OMPQ = AAPQ / ABS(AAPQ) * * .. rotate *[RTD] ROTATED = ROTATED + ONE @@ -924,53 +937,52 @@ * IF( ROTOK ) THEN * - OMPQ = AAPQ / ABS(AAPQ) - AQOAP = AAQQ / AAPP + AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1 * IF( ABS( THETA ).GT.BIGTHETA ) THEN -* +* T = HALF / THETA CS = ONE CALL ZROT( M, A(1,p), 1, A(1,q), 1, - $ CS, DCONJG(OMPQ)*T ) + $ CS, CONJG(OMPQ)*T ) IF ( RSVEC ) THEN - CALL ZROT( MVL, V(1,p), 1, - $ V(1,q), 1, CS, DCONJG(OMPQ)*T ) + CALL ZROT( MVL, V(1,p), 1, + $ V(1,q), 1, CS, CONJG(OMPQ)*T ) END IF - - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + + SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*DSQRT( DMAX1( ZERO, + AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) - MXSINJ = DMAX1( MXSINJ, ABS( T ) ) + MXSINJ = MAX( MXSINJ, ABS( T ) ) * ELSE * * .. choose correct signum for THETA and rotate * - THSIGN = -DSIGN( ONE, AAPQ1 ) - T = ONE / ( THETA+THSIGN* - $ DSQRT( ONE+THETA*THETA ) ) - CS = DSQRT( ONE / ( ONE+T*T ) ) + THSIGN = -SIGN( ONE, AAPQ1 ) + T = ONE / ( THETA+THSIGN* + $ SQRT( ONE+THETA*THETA ) ) + CS = SQRT( ONE / ( ONE+T*T ) ) SN = T*CS * - MXSINJ = DMAX1( MXSINJ, ABS( SN ) ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + MXSINJ = MAX( MXSINJ, ABS( SN ) ) + SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*DSQRT( DMAX1( ZERO, + AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * CALL ZROT( M, A(1,p), 1, A(1,q), 1, - $ CS, DCONJG(OMPQ)*SN ) + $ CS, CONJG(OMPQ)*SN ) IF ( RSVEC ) THEN - CALL ZROT( MVL, V(1,p), 1, - $ V(1,q), 1, CS, DCONJG(OMPQ)*SN ) - END IF - END IF - CWORK(p) = -CWORK(q) * OMPQ + CALL ZROT( MVL, V(1,p), 1, + $ V(1,q), 1, CS, CONJG(OMPQ)*SN ) + END IF + END IF + CWORK(p) = -CWORK(q) * OMPQ * ELSE * .. have to use modified Gram-Schmidt like transformation @@ -985,9 +997,9 @@ $ A( 1, q ), 1 ) CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M, $ 1, A( 1, q ), LDA, IERR ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) - MXSINJ = DMAX1( MXSINJ, SFMIN ) + MXSINJ = MAX( MXSINJ, SFMIN ) END IF * END IF ROTOK THEN ... ELSE * @@ -1004,7 +1016,7 @@ AAQQ = ONE CALL ZLASSQ( M, A( 1, q ), 1, T, $ AAQQ ) - SVA( q ) = T*DSQRT( AAQQ ) + SVA( q ) = T*SQRT( AAQQ ) END IF END IF IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN @@ -1016,7 +1028,7 @@ AAPP = ONE CALL ZLASSQ( M, A( 1, p ), 1, T, $ AAPP ) - AAPP = T*DSQRT( AAPP ) + AAPP = T*SQRT( AAPP ) END IF SVA( p ) = AAPP END IF @@ -1051,7 +1063,7 @@ ELSE SVA( p ) = AAPP IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) ) - $ NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p + $ NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p END IF * 2001 CONTINUE @@ -1071,14 +1083,14 @@ * doing the block at ( ibr, jbc ) * IJBLSK = 0 - DO 2100 p = igl, MIN0( igl+KBL-1, N ) + DO 2100 p = igl, MIN( igl+KBL-1, N ) * AAPP = SVA( p ) IF( AAPP.GT.ZERO ) THEN * PSKIPPED = 0 * - DO 2200 q = jgl, MIN0( jgl+KBL-1, N ) + DO 2200 q = jgl, MIN( jgl+KBL-1, N ) * AAQQ = SVA( q ) IF( AAQQ.GT.ZERO ) THEN @@ -1095,7 +1107,7 @@ ROTOK = ( SMALL*AAQQ ).LE.AAPP END IF IF( AAPP.LT.( BIG / AAQQ ) ) THEN - AAPQ = ( ZDOTC( M, A( 1, p ), 1, + AAPQ = ( ZDOTC( M, A( 1, p ), 1, $ A( 1, q ), 1 ) / AAQQ ) / AAPP ELSE CALL ZCOPY( M, A( 1, p ), 1, @@ -1113,8 +1125,9 @@ ROTOK = AAQQ.LE.( AAPP / SMALL ) END IF IF( AAPP.GT.( SMALL / AAQQ ) ) THEN - AAPQ = ( ZDOTC( M, A( 1, p ), 1, - $ A( 1, q ), 1 ) / AAQQ ) / AAPP + AAPQ = ( ZDOTC( M, A( 1, p ), 1, + $ A( 1, q ), 1 ) / MAX(AAQQ,AAPP) ) + $ / MIN(AAQQ,AAPP) ELSE CALL ZCOPY( M, A( 1, q ), 1, $ CWORK(N+1), 1 ) @@ -1126,13 +1139,15 @@ END IF END IF * -* AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q) + +* AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q) AAPQ1 = -ABS(AAPQ) - MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 ) + MXAAPQ = MAX( MXAAPQ, -AAPQ1 ) * * TO rotate or NOT to rotate, THAT is the question ... * IF( ABS( AAPQ1 ).GT.TOL ) THEN + OMPQ = AAPQ / ABS(AAPQ) NOTROT = 0 *[RTD] ROTATED = ROTATED + 1 PSKIPPED = 0 @@ -1140,7 +1155,6 @@ * IF( ROTOK ) THEN * - OMPQ = AAPQ / ABS(AAPQ) AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1 @@ -1148,42 +1162,42 @@ * IF( ABS( THETA ).GT.BIGTHETA ) THEN T = HALF / THETA - CS = ONE + CS = ONE CALL ZROT( M, A(1,p), 1, A(1,q), 1, - $ CS, DCONJG(OMPQ)*T ) + $ CS, CONJG(OMPQ)*T ) IF( RSVEC ) THEN - CALL ZROT( MVL, V(1,p), 1, - $ V(1,q), 1, CS, DCONJG(OMPQ)*T ) + CALL ZROT( MVL, V(1,p), 1, + $ V(1,q), 1, CS, CONJG(OMPQ)*T ) END IF - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*DSQRT( DMAX1( ZERO, + AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) - MXSINJ = DMAX1( MXSINJ, ABS( T ) ) + MXSINJ = MAX( MXSINJ, ABS( T ) ) ELSE * * .. choose correct signum for THETA and rotate * - THSIGN = -DSIGN( ONE, AAPQ1 ) + THSIGN = -SIGN( ONE, AAPQ1 ) IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN T = ONE / ( THETA+THSIGN* - $ DSQRT( ONE+THETA*THETA ) ) - CS = DSQRT( ONE / ( ONE+T*T ) ) + $ SQRT( ONE+THETA*THETA ) ) + CS = SQRT( ONE / ( ONE+T*T ) ) SN = T*CS - MXSINJ = DMAX1( MXSINJ, ABS( SN ) ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + MXSINJ = MAX( MXSINJ, ABS( SN ) ) + SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*DSQRT( DMAX1( ZERO, + AAPP = AAPP*SQRT( MAX( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * CALL ZROT( M, A(1,p), 1, A(1,q), 1, - $ CS, DCONJG(OMPQ)*SN ) + $ CS, CONJG(OMPQ)*SN ) IF( RSVEC ) THEN - CALL ZROT( MVL, V(1,p), 1, - $ V(1,q), 1, CS, DCONJG(OMPQ)*SN ) + CALL ZROT( MVL, V(1,p), 1, + $ V(1,q), 1, CS, CONJG(OMPQ)*SN ) END IF END IF - CWORK(p) = -CWORK(q) * OMPQ + CWORK(p) = -CWORK(q) * OMPQ * ELSE * .. have to use modified Gram-Schmidt like transformation @@ -1201,9 +1215,9 @@ CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, $ M, 1, A( 1, q ), LDA, $ IERR ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + SVA( q ) = AAQQ*SQRT( MAX( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) - MXSINJ = DMAX1( MXSINJ, SFMIN ) + MXSINJ = MAX( MXSINJ, SFMIN ) ELSE CALL ZCOPY( M, A( 1, q ), 1, $ CWORK(N+1), 1 ) @@ -1213,14 +1227,14 @@ CALL ZLASCL( 'G', 0, 0, AAPP, ONE, $ M, 1, A( 1, p ), LDA, $ IERR ) - CALL ZAXPY( M, -DCONJG(AAPQ), + CALL ZAXPY( M, -CONJG(AAPQ), $ CWORK(N+1), 1, A( 1, p ), 1 ) CALL ZLASCL( 'G', 0, 0, ONE, AAPP, $ M, 1, A( 1, p ), LDA, $ IERR ) - SVA( p ) = AAPP*DSQRT( DMAX1( ZERO, + SVA( p ) = AAPP*SQRT( MAX( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) - MXSINJ = DMAX1( MXSINJ, SFMIN ) + MXSINJ = MAX( MXSINJ, SFMIN ) END IF END IF * END IF ROTOK THEN ... ELSE @@ -1237,7 +1251,7 @@ AAQQ = ONE CALL ZLASSQ( M, A( 1, q ), 1, T, $ AAQQ ) - SVA( q ) = T*DSQRT( AAQQ ) + SVA( q ) = T*SQRT( AAQQ ) END IF END IF IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN @@ -1249,7 +1263,7 @@ AAPP = ONE CALL ZLASSQ( M, A( 1, p ), 1, T, $ AAPP ) - AAPP = T*DSQRT( AAPP ) + AAPP = T*SQRT( AAPP ) END IF SVA( p ) = AAPP END IF @@ -1288,7 +1302,7 @@ ELSE * IF( AAPP.EQ.ZERO )NOTROT = NOTROT + - $ MIN0( jgl+KBL-1, N ) - jgl + 1 + $ MIN( jgl+KBL-1, N ) - jgl + 1 IF( AAPP.LT.ZERO )NOTROT = 0 * END IF @@ -1299,7 +1313,7 @@ * end of the jbc-loop 2011 CONTINUE *2011 bailed out of the jbc-loop - DO 2012 p = igl, MIN0( igl+KBL-1, N ) + DO 2012 p = igl, MIN( igl+KBL-1, N ) SVA( p ) = ABS( SVA( p ) ) 2012 CONTINUE *** @@ -1314,7 +1328,7 @@ T = ZERO AAPP = ONE CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP ) - SVA( N ) = T*DSQRT( AAPP ) + SVA( N ) = T*SQRT( AAPP ) END IF * * Additional steering devices @@ -1322,7 +1336,7 @@ IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. $ ( ISWROT.LE.N ) ) )SWBAND = i * - IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )* + IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( DBLE( N ) )* $ TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN GO TO 1994 END IF @@ -1371,8 +1385,9 @@ * Normalize the left singular vectors. * IF( LSVEC .OR. UCTOL ) THEN - DO 1998 p = 1, N2 - CALL ZDSCAL( M, ONE / SVA( p ), A( 1, p ), 1 ) + DO 1998 p = 1, N4 +* CALL ZDSCAL( M, ONE / SVA( p ), A( 1, p ), 1 ) + CALL ZLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M, IERR ) 1998 CONTINUE END IF * @@ -1386,11 +1401,11 @@ END IF * * Undo scaling, if necessary (and possible). - IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL ) ) ) + IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL ) ) ) $ .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT. $ ( SFMIN / SKL ) ) ) ) THEN DO 2400 p = 1, N - SVA( P ) = SKL*SVA( P ) + SVA( p ) = SKL*SVA( p ) 2400 CONTINUE SKL = ONE END IF