diff options
Diffstat (limited to 'SRC/zgesvj.f')
| -rw-r--r-- | SRC/zgesvj.f | 225 |
1 files changed, 120 insertions, 105 deletions
diff --git a/SRC/zgesvj.f b/SRC/zgesvj.f index 8ac489ba..2a43f1f3 100644 --- a/SRC/zgesvj.f +++ b/SRC/zgesvj.f @@ -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 @@ -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 *> @@ -272,7 +278,7 @@ * *> \date June 2016 * -*> \ingroup doubleGEcomputational +*> \ingroup complex16GEcomputational * *> \par Further Details: * ===================== @@ -298,14 +304,15 @@ *> 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: @@ -329,8 +336,8 @@ *> Department of Mathematics, University of Zagreb, 2008, 2015. *> \endverbatim * -*> \par Bugs, examples and comments: -* ================================= +*> \par Bugs, examples and comments: +* ================================= *> *> \verbatim *> =========================== @@ -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, $ 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,9 +443,9 @@ 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 + ELSE IF( ( LRWORK.LT.MAX( N, 6 ) ) .AND. ( .NOT.LQUERY ) ) THEN INFO = -15 ELSE INFO = 0 @@ -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 ) ) + 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 @@ -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,11 +827,11 @@ * 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 * @@ -858,7 +869,7 @@ 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 ) * @@ -893,7 +904,7 @@ ROTOK = AAPP.LE.( AAQQ / SMALL ) IF( AAPP.GT.( SMALL / AAQQ ) ) THEN AAPQ = ( ZDOTC( M, A( 1, p ), 1, - $ A( 1, q ), 1 ) / AAQQ ) / AAPP + $ A( 1, q ), 1 ) / AAPP ) / AAQQ ELSE CALL ZCOPY( M, A( 1, q ), 1, $ CWORK(N+1), 1 ) @@ -905,13 +916,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) * * .. rotate *[RTD] ROTATED = ROTATED + ONE @@ -924,8 +937,7 @@ * IF( ROTOK ) THEN * - OMPQ = AAPQ / ABS(AAPQ) - AQOAP = AAQQ / AAPP + AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1 * @@ -935,39 +947,39 @@ 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 ) + $ 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 ) 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 ) + $ V(1,q), 1, CS, CONJG(OMPQ)*SN ) END IF END IF CWORK(p) = -CWORK(q) * OMPQ @@ -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 @@ -1114,7 +1126,8 @@ END IF IF( AAPP.GT.( SMALL / AAQQ ) ) THEN AAPQ = ( ZDOTC( M, A( 1, p ), 1, - $ A( 1, q ), 1 ) / AAQQ ) / AAPP + $ 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 @@ -1150,37 +1164,37 @@ 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 ) + $ 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 ) + $ V(1,q), 1, CS, CONJG(OMPQ)*SN ) END IF END IF CWORK(p) = -CWORK(q) * OMPQ @@ -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 * @@ -1390,7 +1405,7 @@ $ .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 |