From 5768f537085f3a27a3157df778af85b41c201491 Mon Sep 17 00:00:00 2001 From: julie Date: Sun, 8 Nov 2015 23:42:08 +0000 Subject: First commit for Zlatko Drmac Contribution - Fixing z precisions issues + modif sent by Zlatko --- SRC/cgesvj.f | 2 +- SRC/cgsvj0.f | 4 +- SRC/cgsvj1.f | 7 +-- SRC/zgesvj.f | 166 +++++++++++++++++++++++++++++------------------------------ SRC/zgsvj0.f | 110 +++++++++++++++++++-------------------- SRC/zgsvj1.f | 74 +++++++++++++------------- 6 files changed, 179 insertions(+), 184 deletions(-) diff --git a/SRC/cgesvj.f b/SRC/cgesvj.f index f7115d31..69d77048 100644 --- a/SRC/cgesvj.f +++ b/SRC/cgesvj.f @@ -385,7 +385,7 @@ * .. External Subroutines .. * .. * from BLAS - EXTERNAL CCOPY, CSROT, CSSCAL, CSWAP + EXTERNAL CCOPY, CROT, CSSCAL, CSWAP * from LAPACK EXTERNAL CLASCL, CLASET, CLASSQ, XERBLA EXTERNAL CGSVJ0, CGSVJ1 diff --git a/SRC/cgsvj0.f b/SRC/cgsvj0.f index 690bed01..79ffde62 100644 --- a/SRC/cgsvj0.f +++ b/SRC/cgsvj0.f @@ -255,7 +255,7 @@ * .. * .. * .. Intrinsic Functions .. - INTRINSIC ABS, AMAX1, CONJG, FLOAT, MIN0, MAX0, SIGN, SQRT + INTRINSIC ABS, AMAX1, CONJG, FLOAT, MIN0, SIGN, SQRT * .. * .. External Functions .. REAL SCNRM2 @@ -268,7 +268,7 @@ * .. External Subroutines .. * .. * from BLAS - EXTERNAL CCOPY, CSROT, CSSCAL, CSWAP + EXTERNAL CCOPY, CROT, CSSCAL, CSWAP * from LAPACK EXTERNAL CLASCL, CLASSQ, XERBLA * .. diff --git a/SRC/cgsvj1.f b/SRC/cgsvj1.f index 5e73b373..f4b1fc15 100644 --- a/SRC/cgsvj1.f +++ b/SRC/cgsvj1.f @@ -268,11 +268,9 @@ $ p, PSKIPPED, q, ROWSKIP, SWBAND LOGICAL APPLV, ROTOK, RSVEC * .. -* .. Local Arrays .. - REAL FASTR( 5 ) * .. * .. Intrinsic Functions .. - INTRINSIC ABS, AMAX1, FLOAT, MIN0, SIGN, SQRT + INTRINSIC ABS, AMAX1, CONJG, FLOAT, MIN0, SIGN, SQRT * .. * .. External Functions .. REAL SCNRM2 @@ -283,7 +281,7 @@ * .. * .. External Subroutines .. * .. from BLAS - EXTERNAL CCOPY, CSROT, CSSCAL, CSWAP + EXTERNAL CCOPY, CROT, CSSCAL, CSWAP * .. from LAPACK EXTERNAL CLASCL, CLASSQ, XERBLA * .. @@ -346,7 +344,6 @@ * EMPTSW = N1*( N-N1 ) NOTROT = 0 - FASTR( 1 ) = ZERO * * .. Row-cyclic pivot strategy with de Rijk's pivoting .. * diff --git a/SRC/zgesvj.f b/SRC/zgesvj.f index cd2e02f4..930a3530 100644 --- a/SRC/zgesvj.f +++ b/SRC/zgesvj.f @@ -360,10 +360,10 @@ * ===================================================================== * * .. Local Parameters .. - DOUBLE PRECISION ZERO, HALF, ONE - PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0) - COMPLEX*16 CZERO, CONE - PARAMETER ( CZERO = (0.0E0, 0.0E0), CONE = (1.0E0, 0.0E0) ) + DOUBLE PRECISION ZERO, HALF, ONE + PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = (0.0D0, 0.0D0), CONE = (1.0D0, 0.0D0) ) INTEGER NSWEEP PARAMETER ( NSWEEP = 30 ) * .. @@ -381,8 +381,8 @@ * .. * .. * .. Intrinsic Functions .. - INTRINSIC ABS, AMAX1, AMIN1, DCONJG, DBLE, MIN0, MAX0, - $ SIGN, SQRT + INTRINSIC ABS, DMAX1, DMIN1, DCONJG, DFLOAT, MIN0, MAX0, + $ DSIGN, DSQRT * .. * .. External Functions .. * .. @@ -390,8 +390,8 @@ DOUBLE PRECISION DZNRM2 COMPLEX*16 ZDOTC EXTERNAL ZDOTC, DZNRM2 - INTEGER IZAMAX - EXTERNAL IZAMAX + INTEGER IDAMAX + EXTERNAL IDAMAX * from LAPACK DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH @@ -401,7 +401,7 @@ * .. External Subroutines .. * .. * from BLAS - EXTERNAL ZCOPY, ZDROT, ZDSCAL, ZSWAP + EXTERNAL ZCOPY, ZROT, ZDSCAL, ZSWAP * from LAPACK EXTERNAL ZLASCL, ZLASET, ZLASSQ, XERBLA EXTERNAL ZGSVJ0, ZGSVJ1 @@ -467,29 +467,29 @@ ELSE * ... default IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN - CTOL = SQRT( DBLE( M ) ) + CTOL = DSQRT( DFLOAT( M ) ) ELSE - CTOL = DBLE( M ) + CTOL = DFLOAT( M ) END IF END IF * ... and the machine dependent parameters are *[!] (Make sure that DLAMCH() works properly on the target machine.) * EPSLN = DLAMCH( 'Epsilon' ) - ROOTEPS = SQRT( EPSLN ) + ROOTEPS = DSQRT( EPSLN ) SFMIN = DLAMCH( 'SafeMinimum' ) - ROOTSFMIN = SQRT( SFMIN ) + ROOTSFMIN = DSQRT( SFMIN ) SMALL = SFMIN / EPSLN BIG = DLAMCH( 'Overflow' ) * BIG = ONE / SFMIN ROOTBIG = ONE / ROOTSFMIN - LARGE = BIG / SQRT( DBLE( M*N ) ) + LARGE = BIG / DSQRT( DFLOAT( M*N ) ) BIGTHETA = ONE / ROOTEPS * TOL = CTOL*EPSLN - ROOTTOL = SQRT( TOL ) + ROOTTOL = DSQRT( TOL ) * - IF( DBLE( M )*EPSLN.GE.ONE ) THEN + IF( DFLOAT( M )*EPSLN.GE.ONE ) THEN INFO = -4 CALL XERBLA( 'ZGESVJ', -INFO ) RETURN @@ -514,7 +514,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 / SQRT( DBLE( M )*DBLE( N ) ) + SKL = ONE / DSQRT( DFLOAT( M )*DFLOAT( N ) ) NOSCALE = .TRUE. GOSCALE = .TRUE. * @@ -529,7 +529,7 @@ CALL XERBLA( 'ZGESVJ', -INFO ) RETURN END IF - AAQQ = SQRT( AAQQ ) + AAQQ = DSQRT( AAQQ ) IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN SVA( p ) = AAPP*AAQQ ELSE @@ -554,7 +554,7 @@ CALL XERBLA( 'ZGESVJ', -INFO ) RETURN END IF - AAQQ = SQRT( AAQQ ) + AAQQ = DSQRT( AAQQ ) IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN SVA( p ) = AAPP*AAQQ ELSE @@ -579,7 +579,7 @@ CALL XERBLA( 'ZGESVJ', -INFO ) RETURN END IF - AAQQ = SQRT( AAQQ ) + AAQQ = DSQRT( AAQQ ) IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN SVA( p ) = AAPP*AAQQ ELSE @@ -604,8 +604,8 @@ AAPP = ZERO AAQQ = BIG DO 4781 p = 1, N - IF( SVA( p ).NE.ZERO )AAQQ = AMIN1( AAQQ, SVA( p ) ) - AAPP = AMAX1( AAPP, SVA( p ) ) + IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) ) + AAPP = DMAX1( AAPP, SVA( p ) ) 4781 CONTINUE * * #:) Quick return for zero matrix @@ -642,23 +642,23 @@ * Protect small singular values from underflow, and try to * avoid underflows/overflows in computing Jacobi rotations. * - SN = SQRT( SFMIN / EPSLN ) - TEMP1 = SQRT( BIG / DBLE( N ) ) + SN = DSQRT( SFMIN / EPSLN ) + TEMP1 = DSQRT( BIG / DFLOAT( N ) ) IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. $ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN - TEMP1 = AMIN1( BIG, TEMP1 / AAPP ) + TEMP1 = DMIN1( BIG, TEMP1 / AAPP ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN - TEMP1 = AMIN1( SN / AAQQ, BIG / ( AAPP*SQRT( DBLE( N ) ) ) ) + TEMP1 = DMIN1( SN / AAQQ, BIG / (AAPP*DSQRT( DFLOAT(N)) ) ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN - TEMP1 = AMAX1( SN / AAQQ, TEMP1 / AAPP ) + TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN - TEMP1 = AMIN1( SN / AAQQ, BIG / ( SQRT( DBLE( N ) )*AAPP ) ) + TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DFLOAT( N ) )*AAPP ) ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE @@ -668,7 +668,7 @@ * Scale, if necessary * IF( TEMP1.NE.ONE ) THEN - CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR ) + CALL ZLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR ) END IF SKL = TEMP1*SKL IF( SKL.NE.ONE ) THEN @@ -690,7 +690,7 @@ SWBAND = 3 *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective * if ZGESVJ is used as a computational routine in the preconditioned -* Jacobi SVD algorithm CGEJSV. For sweeps i=1:SWBAND the procedure +* Jacobi SVD algorithm ZGEJSV. For sweeps i=1:SWBAND the procedure * works on pivots inside a band-like region around the diagonal. * The boundaries are determined dynamically, based on the number of * pivots above a threshold. @@ -754,7 +754,7 @@ $ CWORK( N2+1 ), SVA( N2+1 ), MVL, $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, $ CWORK( N+1 ), LWORK-N, IERR ) -* + CALL ZGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA, $ CWORK( N4+1 ), SVA( N4+1 ), MVL, $ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1, @@ -824,7 +824,7 @@ * * .. de Rijk's pivoting * - q = IZAMAX( N-p+1, SVA( p ), 1 ) + p - 1 + 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, @@ -843,12 +843,12 @@ * norm computation. *[!] Caveat: * Unfortunately, some BLAS implementations compute DZNRM2(M,A(1,p),1) -* as SQRT(S=ZDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to +* as SQRT(S=CDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to * overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to * underflow for ||A(:,p)||_2 < SQRT(underflow_threshold). * Hence, DZNRM2 cannot be trusted, not even in the case when * the true norm is far from the under(over)flow boundaries. -* If properly implemented DZNRM2 is available, the IF-THEN-ELSE-END IF +* 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. @@ -858,7 +858,7 @@ TEMP1 = ZERO AAPP = ONE CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) - SVA( p ) = TEMP1*SQRT( AAPP ) + SVA( p ) = TEMP1*DSQRT( AAPP ) END IF AAPP = SVA( p ) ELSE @@ -908,7 +908,7 @@ OMPQ = AAPQ / ABS(AAPQ) * AAPQ = AAPQ * DCONJG( CWORK(p) ) * CWORK(q) AAPQ1 = -ABS(AAPQ) - MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 ) + MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 ) * * TO rotate or NOT to rotate, THAT is the question ... * @@ -934,39 +934,39 @@ T = HALF / THETA CS = ONE - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), 1, $ CS, DCONJG(OMPQ)*T ) IF ( RSVEC ) THEN - CALL CROT( MVL, V(1,p), 1, + CALL ZROT( MVL, V(1,p), 1, $ V(1,q), 1, CS, DCONJG(OMPQ)*T ) END IF - SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, + SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*SQRT( AMAX1( ZERO, + AAPP = AAPP*DSQRT( DMAX1( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) - MXSINJ = AMAX1( MXSINJ, ABS( T ) ) + MXSINJ = DMAX1( MXSINJ, ABS( T ) ) * ELSE * * .. choose correct signum for THETA and rotate * - THSIGN = -SIGN( ONE, AAPQ1 ) + THSIGN = -DSIGN( ONE, AAPQ1 ) T = ONE / ( THETA+THSIGN* - $ SQRT( ONE+THETA*THETA ) ) - CS = SQRT( ONE / ( ONE+T*T ) ) + $ DSQRT( ONE+THETA*THETA ) ) + CS = DSQRT( ONE / ( ONE+T*T ) ) SN = T*CS * - MXSINJ = AMAX1( MXSINJ, ABS( SN ) ) - SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, + MXSINJ = DMAX1( MXSINJ, ABS( SN ) ) + SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*SQRT( AMAX1( ZERO, + AAPP = AAPP*DSQRT( DMAX1( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), 1, $ CS, DCONJG(OMPQ)*SN ) IF ( RSVEC ) THEN - CALL CROT( MVL, V(1,p), 1, + CALL ZROT( MVL, V(1,p), 1, $ V(1,q), 1, CS, DCONJG(OMPQ)*SN ) END IF END IF @@ -981,13 +981,13 @@ $ IERR ) CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, M, $ 1, A( 1, q ), LDA, IERR ) - CALL CAXPY( M, -AAPQ, CWORK(N+1), 1, + CALL ZAXPY( M, -AAPQ, CWORK(N+1), 1, $ A( 1, q ), 1 ) CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M, $ 1, A( 1, q ), LDA, IERR ) - SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, + SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) - MXSINJ = AMAX1( MXSINJ, SFMIN ) + MXSINJ = DMAX1( MXSINJ, SFMIN ) END IF * END IF ROTOK THEN ... ELSE * @@ -1004,7 +1004,7 @@ AAQQ = ONE CALL ZLASSQ( M, A( 1, q ), 1, T, $ AAQQ ) - SVA( q ) = T*SQRT( AAQQ ) + SVA( q ) = T*DSQRT( AAQQ ) END IF END IF IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN @@ -1016,7 +1016,7 @@ AAPP = ONE CALL ZLASSQ( M, A( 1, p ), 1, T, $ AAPP ) - AAPP = T*SQRT( AAPP ) + AAPP = T*DSQRT( AAPP ) END IF SVA( p ) = AAPP END IF @@ -1129,7 +1129,7 @@ OMPQ = AAPQ / ABS(AAPQ) * AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q) AAPQ1 = -ABS(AAPQ) - MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 ) + MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 ) * * TO rotate or NOT to rotate, THAT is the question ... * @@ -1149,37 +1149,37 @@ IF( ABS( THETA ).GT.BIGTHETA ) THEN T = HALF / THETA CS = ONE - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), 1, $ CS, DCONJG(OMPQ)*T ) IF( RSVEC ) THEN - CALL CROT( MVL, V(1,p), 1, + CALL ZROT( MVL, V(1,p), 1, $ V(1,q), 1, CS, DCONJG(OMPQ)*T ) END IF - SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, + SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*SQRT( AMAX1( ZERO, + AAPP = AAPP*DSQRT( DMAX1( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) - MXSINJ = AMAX1( MXSINJ, ABS( T ) ) + MXSINJ = DMAX1( MXSINJ, ABS( T ) ) ELSE * * .. choose correct signum for THETA and rotate * - THSIGN = -SIGN( ONE, AAPQ1 ) + THSIGN = -DSIGN( ONE, AAPQ1 ) IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN T = ONE / ( THETA+THSIGN* - $ SQRT( ONE+THETA*THETA ) ) - CS = SQRT( ONE / ( ONE+T*T ) ) + $ DSQRT( ONE+THETA*THETA ) ) + CS = DSQRT( ONE / ( ONE+T*T ) ) SN = T*CS - MXSINJ = AMAX1( MXSINJ, ABS( SN ) ) - SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, + MXSINJ = DMAX1( MXSINJ, ABS( SN ) ) + SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, $ ONE+T*APOAQ*AAPQ1 ) ) - AAPP = AAPP*SQRT( AMAX1( ZERO, + AAPP = AAPP*DSQRT( DMAX1( ZERO, $ ONE-T*AQOAP*AAPQ1 ) ) * - CALL CROT( M, A(1,p), 1, A(1,q), 1, + CALL ZROT( M, A(1,p), 1, A(1,q), 1, $ CS, DCONJG(OMPQ)*SN ) IF( RSVEC ) THEN - CALL CROT( MVL, V(1,p), 1, + CALL ZROT( MVL, V(1,p), 1, $ V(1,q), 1, CS, DCONJG(OMPQ)*SN ) END IF END IF @@ -1196,14 +1196,14 @@ CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, $ M, 1, A( 1, q ), LDA, $ IERR ) - CALL CAXPY( M, -AAPQ, CWORK(N+1), + CALL ZAXPY( M, -AAPQ, CWORK(N+1), $ 1, A( 1, q ), 1 ) CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, $ M, 1, A( 1, q ), LDA, $ IERR ) - SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, + SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) - MXSINJ = AMAX1( MXSINJ, SFMIN ) + MXSINJ = DMAX1( MXSINJ, SFMIN ) ELSE CALL ZCOPY( M, A( 1, q ), 1, $ CWORK(N+1), 1 ) @@ -1213,14 +1213,14 @@ CALL ZLASCL( 'G', 0, 0, AAPP, ONE, $ M, 1, A( 1, p ), LDA, $ IERR ) - CALL CAXPY( M, -DCONJG(AAPQ), + CALL ZAXPY( M, -DCONJG(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*SQRT( AMAX1( ZERO, + SVA( p ) = AAPP*DSQRT( DMAX1( ZERO, $ ONE-AAPQ1*AAPQ1 ) ) - MXSINJ = AMAX1( MXSINJ, SFMIN ) + MXSINJ = DMAX1( MXSINJ, SFMIN ) END IF END IF * END IF ROTOK THEN ... ELSE @@ -1237,7 +1237,7 @@ AAQQ = ONE CALL ZLASSQ( M, A( 1, q ), 1, T, $ AAQQ ) - SVA( q ) = T*SQRT( AAQQ ) + SVA( q ) = T*DSQRT( AAQQ ) END IF END IF IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN @@ -1249,7 +1249,7 @@ AAPP = ONE CALL ZLASSQ( M, A( 1, p ), 1, T, $ AAPP ) - AAPP = T*SQRT( AAPP ) + AAPP = T*DSQRT( AAPP ) END IF SVA( p ) = AAPP END IF @@ -1314,7 +1314,7 @@ T = ZERO AAPP = ONE CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP ) - SVA( N ) = T*SQRT( AAPP ) + SVA( N ) = T*DSQRT( AAPP ) END IF * * Additional steering devices @@ -1322,8 +1322,8 @@ IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. $ ( ISWROT.LE.N ) ) )SWBAND = i * - IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( DBLE( N ) )* - $ TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN + IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DFLOAT( N ) )* + $ TOL ) .AND. ( DFLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN GO TO 1994 END IF * @@ -1350,7 +1350,7 @@ N2 = 0 N4 = 0 DO 5991 p = 1, N - 1 - q = IZAMAX( N-p+1, SVA( p ), 1 ) + p - 1 + q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1 IF( p.NE.q ) THEN TEMP1 = SVA( p ) SVA( p ) = SVA( q ) @@ -1400,15 +1400,15 @@ * then some of the singular values may overflow or underflow and * the spectrum is given in this factored representation. * - RWORK( 2 ) = DBLE( N4 ) + RWORK( 2 ) = DFLOAT( N4 ) * N4 is the number of computed nonzero singular values of A. * - RWORK( 3 ) = DBLE( N2 ) + RWORK( 3 ) = DFLOAT( N2 ) * N2 is the number of singular values of A greater than SFMIN. * If N2