diff options
author | Julien Langou <julien.langou@ucdenver.edu> | 2016-12-19 11:27:35 +0100 |
---|---|---|
committer | Julien Langou <julien.langou@ucdenver.edu> | 2016-12-19 11:27:35 +0100 |
commit | ad5bc21cb50535d66d628a309d60128db96c8851 (patch) | |
tree | cc7b72b0795c8c64ebf18cf28c984c41cfbedc54 /SRC/zgsvj1.f | |
parent | 5f3f247a5876ae4d5c67a765ffe8a35ef7944211 (diff) | |
download | lapack-ad5bc21cb50535d66d628a309d60128db96c8851.tar.gz lapack-ad5bc21cb50535d66d628a309d60128db96c8851.tar.bz2 lapack-ad5bc21cb50535d66d628a309d60128db96c8851.zip |
contribution from Zlatko Drmac
Note: I still need to work on merging [C/Z]GEJSV, but there is much more work
on these two files. We will see when this can be done.
Diffstat (limited to 'SRC/zgsvj1.f')
-rw-r--r-- | SRC/zgsvj1.f | 85 |
1 files changed, 43 insertions, 42 deletions
diff --git a/SRC/zgsvj1.f b/SRC/zgsvj1.f index 89ce3d01..9c764c89 100644 --- a/SRC/zgsvj1.f +++ b/SRC/zgsvj1.f @@ -27,8 +27,8 @@ * CHARACTER*1 JOBV * .. * .. Array Arguments .. -* COMPLEX*16 A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK ) -* DOUBLE PRECISION SVA( N ) +* COMPLEX*16 A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK ) +* DOUBLE PRECISION SVA( N ) * .. * * @@ -227,10 +227,10 @@ * *> \ingroup complex16OTHERcomputational * -*> \par Contributors: +*> \par Contributor: * ================== *> -*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) +*> Zlatko Drmac (Zagreb, Croatia) * * ===================================================================== SUBROUTINE ZGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, @@ -261,7 +261,7 @@ * .. Local Scalars .. COMPLEX*16 AAPQ, OMPQ DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG, - $ BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG, + $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, $ ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T, $ TEMP1, THETA, THSIGN INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK, @@ -271,7 +271,7 @@ * .. * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DCONJG, DMAX1, DBLE, MIN0, DSIGN, DSQRT + INTRINSIC ABS, CONJG, MAX, DBLE, MIN, SIGN, SQRT * .. * .. External Functions .. DOUBLE PRECISION DZNRM2 @@ -330,14 +330,14 @@ END IF RSVEC = RSVEC .OR. APPLV - ROOTEPS = DSQRT( EPS ) - ROOTSFMIN = DSQRT( SFMIN ) + ROOTEPS = SQRT( EPS ) + ROOTSFMIN = SQRT( SFMIN ) SMALL = SFMIN / EPS BIG = ONE / SFMIN ROOTBIG = ONE / ROOTSFMIN - LARGE = BIG / DSQRT( DBLE( M*N ) ) +* LARGE = BIG / SQRT( DBLE( M*N ) ) BIGTHETA = ONE / ROOTEPS - ROOTTOL = DSQRT( TOL ) + ROOTTOL = SQRT( TOL ) * * .. Initialize the right singular vector matrix .. * @@ -348,7 +348,7 @@ * * .. Row-cyclic pivot strategy with de Rijk's pivoting .. * - KBL = MIN0( 8, N ) + KBL = MIN( 8, N ) NBLR = N1 / KBL IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1 @@ -359,7 +359,7 @@ BLSKIP = ( KBL**2 ) + 1 *[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. SWBAND = 0 *[TP] SWBAND is a tuning parameter. It is meaningful and effective @@ -409,14 +409,14 @@ * doing the block at ( ibr, jbc ) * IJBLSK = 0 - DO 2100 p = igl, MIN0( igl+KBL-1, N1 ) + DO 2100 p = igl, MIN( igl+KBL-1, N1 ) * 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 @@ -452,7 +452,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, $ WORK, 1 ) @@ -464,14 +465,14 @@ END IF END IF * - OMPQ = AAPQ / ABS(AAPQ) -* 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 @@ -488,37 +489,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 D(p) = -D(q) * OMPQ @@ -539,9 +540,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, $ WORK, 1 ) @@ -551,14 +552,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), $ WORK, 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 @@ -575,7 +576,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 @@ -587,7 +588,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 @@ -626,7 +627,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 @@ -637,7 +638,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 *** @@ -652,7 +653,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 @@ -660,7 +661,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 |