--- rpl/lapack/lapack/dbdsvdx.f 2016/08/27 15:34:19 1.3 +++ rpl/lapack/lapack/dbdsvdx.f 2017/06/17 10:53:46 1.4 @@ -2,23 +2,23 @@ * * =========== 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 DBDSVDX + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> +*> Download DBDSVDX + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> *> [TXT] -*> \endhtmlonly +*> \endhtmlonly * * Definition: * =========== * -* SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, +* SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, * $ NS, S, Z, LDZ, WORK, IWORK, INFO ) * * .. Scalar Arguments .. @@ -28,45 +28,45 @@ * .. * .. Array Arguments .. * INTEGER IWORK( * ) -* DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), +* DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), * Z( LDZ, * ) * .. -* +* *> \par Purpose: * ============= *> *> \verbatim *> *> DBDSVDX computes the singular value decomposition (SVD) of a real -*> N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT, -*> where S is a diagonal matrix with non-negative diagonal elements -*> (the singular values of B), and U and VT are orthogonal matrices +*> N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT, +*> where S is a diagonal matrix with non-negative diagonal elements +*> (the singular values of B), and U and VT are orthogonal matrices *> of left and right singular vectors, respectively. *> -*> Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ] -*> and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the -*> singular value decompositon of B through the eigenvalues and +*> Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ] +*> and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the +*> singular value decompositon of B through the eigenvalues and *> eigenvectors of the N*2-by-N*2 tridiagonal matrix -*> -*> | 0 d_1 | -*> | d_1 0 e_1 | -*> TGK = | e_1 0 d_2 | -*> | d_2 . . | +*> +*> | 0 d_1 | +*> | d_1 0 e_1 | +*> TGK = | e_1 0 d_2 | +*> | d_2 . . | *> | . . . | *> -*> If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then -*> (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) / -*> sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and -*> P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ]. -*> -*> Given a TGK matrix, one can either a) compute -s,-v and change signs -*> so that the singular values (and corresponding vectors) are already in -*> descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder -*> the values (and corresponding vectors). DBDSVDX implements a) by -*> calling DSTEVX (bisection plus inverse iteration, to be replaced -*> with a version of the Multiple Relative Robust Representation -*> algorithm. (See P. Willems and B. Lang, A framework for the MR^3 -*> algorithm: theory and implementation, SIAM J. Sci. Comput., +*> If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then +*> (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) / +*> sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and +*> P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ]. +*> +*> Given a TGK matrix, one can either a) compute -s,-v and change signs +*> so that the singular values (and corresponding vectors) are already in +*> descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder +*> the values (and corresponding vectors). DBDSVDX implements a) by +*> calling DSTEVX (bisection plus inverse iteration, to be replaced +*> with a version of the Multiple Relative Robust Representation +*> algorithm. (See P. Willems and B. Lang, A framework for the MR^3 +*> algorithm: theory and implementation, SIAM J. Sci. Comput., *> 35:740-766, 2013.) *> \endverbatim * @@ -101,13 +101,13 @@ *> N is INTEGER *> The order of the bidiagonal matrix. N >= 0. *> \endverbatim -*> +*> *> \param[in] D *> \verbatim *> D is DOUBLE PRECISION array, dimension (N) *> The n diagonal elements of the bidiagonal matrix B. *> \endverbatim -*> +*> *> \param[in] E *> \verbatim *> E is DOUBLE PRECISION array, dimension (max(1,N-1)) @@ -167,14 +167,14 @@ *> \verbatim *> Z is DOUBLE PRECISION array, dimension (2*N,K) ) *> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z -*> contain the singular vectors of the matrix B corresponding to +*> contain the singular vectors of the matrix B corresponding to *> the selected singular values, with U in rows 1 to N and V *> in rows N+1 to N*2, i.e. -*> Z = [ U ] +*> Z = [ U ] *> [ V ] -*> If JOBZ = 'N', then Z is not referenced. -*> Note: The user must ensure that at least K = NS+1 columns are -*> supplied in the array Z; if RANGE = 'V', the exact value of +*> If JOBZ = 'N', then Z is not referenced. +*> Note: The user must ensure that at least K = NS+1 columns are +*> supplied in the array Z; if RANGE = 'V', the exact value of *> NS is not known in advance and an upper bound must be used. *> \endverbatim *> @@ -194,7 +194,7 @@ *> \verbatim *> IWORK is INTEGER array, dimension (12*N) *> If JOBZ = 'V', then if INFO = 0, the first NS elements of -*> IWORK are zero. If INFO > 0, then IWORK contains the indices +*> IWORK are zero. If INFO > 0, then IWORK contains the indices *> of the eigenvectors that failed to converge in DSTEVX. *> \endverbatim *> @@ -213,24 +213,24 @@ * 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 doubleOTHEReigen * -* ===================================================================== - SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, +* ===================================================================== + SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, $ NS, S, Z, LDZ, WORK, IWORK, INFO) * -* -- LAPACK driver routine (version 3.6.1) -- +* -- LAPACK driver 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..-- -* November 2016 -* +* December 2016 +* * .. Scalar Arguments .. CHARACTER JOBZ, RANGE, UPLO INTEGER IL, INFO, IU, LDZ, N, NS @@ -238,28 +238,28 @@ * .. * .. Array Arguments .. INTEGER IWORK( * ) - DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), + DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), $ Z( LDZ, * ) * .. * * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH - PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0, + DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0, $ HNDRD = 100.0D0, MEIGTH = -0.1250D0 ) DOUBLE PRECISION FUDGE PARAMETER ( FUDGE = 2.0D0 ) * .. -* .. Local Scalars .. +* .. Local Scalars .. CHARACTER RNGVX - LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ - INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR, - $ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV, - $ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K, + LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ + INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR, + $ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV, + $ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K, $ NTGK, NRU, NRV, NSL DOUBLE PRECISION ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX, - $ SMIN, SQRT2, THRESH, TOL, ULP, + $ SMIN, SQRT2, THRESH, TOL, ULP, $ VLTGK, VUTGK, ZJTJI * .. * .. External Functions .. @@ -274,7 +274,7 @@ * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, SIGN, SQRT * .. -* .. Executable Statements .. +* .. Executable Statements .. * * Test the input parameters. * @@ -321,7 +321,7 @@ * NS = 0 IF( N.EQ.0 ) RETURN -* +* IF( N.EQ.1 ) THEN IF( ALLSV .OR. INDSV ) THEN NS = 1 @@ -339,17 +339,17 @@ RETURN END IF * - ABSTOL = 2*DLAMCH( 'Safe Minimum' ) + ABSTOL = 2*DLAMCH( 'Safe Minimum' ) ULP = DLAMCH( 'Precision' ) EPS = DLAMCH( 'Epsilon' ) SQRT2 = SQRT( 2.0D0 ) ORTOL = SQRT( ULP ) -* +* * Criterion for splitting is taken from DBDSQR when singular -* values are computed to relative accuracy TOL. (See J. Demmel and -* W. Kahan, Accurate singular values of bidiagonal matrices, SIAM +* values are computed to relative accuracy TOL. (See J. Demmel and +* W. Kahan, Accurate singular values of bidiagonal matrices, SIAM * J. Sci. and Stat. Comput., 11:873–912, 1990.) -* +* TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS * * Compute approximate maximum, minimum singular values. @@ -390,19 +390,19 @@ IIWORK = IIFAIL + N*2 * * Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode. -* VL,VU or IL,IU are redefined to conform to implementation a) +* VL,VU or IL,IU are redefined to conform to implementation a) * described in the leading comments. * ILTGK = 0 - IUTGK = 0 + IUTGK = 0 VLTGK = ZERO VUTGK = ZERO * IF( ALLSV ) THEN * -* All singular values will be found. We aim at -s (see +* All singular values will be found. We aim at -s (see * leading comments) with RNGVX = 'I'. IL and IU are set -* later (as ILTGK and IUTGK) according to the dimension +* later (as ILTGK and IUTGK) according to the dimension * of the active submatrix. * RNGVX = 'I' @@ -419,7 +419,7 @@ WORK( IDTGK:IDTGK+2*N-1 ) = ZERO CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) - CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ), + CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ), $ VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S, $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), $ IWORK( IIFAIL ), INFO ) @@ -430,20 +430,20 @@ END IF ELSE IF( INDSV ) THEN * -* Find the IL-th through the IU-th singular values. We aim -* at -s (see leading comments) and indices are mapped into +* Find the IL-th through the IU-th singular values. We aim +* at -s (see leading comments) and indices are mapped into * values, therefore mimicking DSTEBZ, where * * GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN * GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN * ILTGK = IL - IUTGK = IU + IUTGK = IU RNGVX = 'V' WORK( IDTGK:IDTGK+2*N-1 ) = ZERO CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) - CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), + CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), $ VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S, $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), $ IWORK( IIFAIL ), INFO ) @@ -451,7 +451,7 @@ WORK( IDTGK:IDTGK+2*N-1 ) = ZERO CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) - CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), + CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), $ VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S, $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), $ IWORK( IIFAIL ), INFO ) @@ -459,12 +459,12 @@ VUTGK = MIN( VUTGK, ZERO ) * * If VLTGK=VUTGK, DSTEVX returns an error message, -* so if needed we change VUTGK slightly. +* so if needed we change VUTGK slightly. * IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL * IF( WANTZ ) CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ) - END IF + END IF * * Initialize variables and pointers for S, Z, and WORK. * @@ -483,7 +483,7 @@ IROWU = 2 IROWV = 1 SPLIT = .FALSE. - SVEQ0 = .FALSE. + SVEQ0 = .FALSE. * * Form the tridiagonal TGK matrix. * @@ -494,15 +494,15 @@ CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) * * -* Check for splits in two levels, outer level +* Check for splits in two levels, outer level * in E and inner level in D. * - DO IEPTR = 2, N*2, 2 - IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN + DO IEPTR = 2, N*2, 2 + IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN * * Split in E (this piece of B is square) or bottom * of the (input bidiagonal) matrix. -* +* ISPLT = IDBEG IDEND = IEPTR - 1 DO IDPTR = IDBEG, IDEND, 2 @@ -519,22 +519,22 @@ IF( IDBEG.EQ.IDEND) THEN NRU = 1 NRV = 1 - END IF + END IF ELSE IF( IDPTR.EQ.IDEND ) THEN * * D=0 at the bottom. * SVEQ0 = .TRUE. - NRU = (IDEND-ISPLT)/2 + 1 - NRV = NRU + NRU = (IDEND-ISPLT)/2 + 1 + NRV = NRU IF( ISPLT.NE.IDBEG ) THEN NRU = NRU + 1 - END IF + END IF ELSE IF( ISPLT.EQ.IDBEG ) THEN * * Split: top rectangular submatrix. -* +* NRU = (IDPTR-IDBEG)/2 NRV = NRU + 1 ELSE @@ -542,7 +542,7 @@ * Split: middle square submatrix. * NRU = (IDPTR-ISPLT)/2 + 1 - NRV = NRU + NRV = NRU END IF END IF ELSE IF( IDPTR.EQ.IDEND ) THEN @@ -560,7 +560,7 @@ * Split: bottom rectangular submatrix. * NRV = (IDEND-ISPLT)/2 + 1 - NRU = NRV + 1 + NRU = NRV + 1 END IF END IF * @@ -568,32 +568,32 @@ * IF( NTGK.GT.0 ) THEN * -* Compute eigenvalues/vectors of the active -* submatrix according to RANGE: +* Compute eigenvalues/vectors of the active +* submatrix according to RANGE: * if RANGE='A' (ALLSV) then RNGVX = 'I' * if RANGE='V' (VALSV) then RNGVX = 'V' * if RANGE='I' (INDSV) then RNGVX = 'V' * ILTGK = 1 - IUTGK = NTGK / 2 + IUTGK = NTGK / 2 IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN - IF( SVEQ0 .OR. - $ SMIN.LT.EPS .OR. + IF( SVEQ0 .OR. + $ SMIN.LT.EPS .OR. $ MOD(NTGK,2).GT.0 ) THEN * Special case: eigenvalue equal to zero or very * small, additional eigenvector is needed. IUTGK = IUTGK + 1 - END IF + END IF END IF * -* Workspace needed by DSTEVX: -* WORK( ITEMP: ): 2*5*NTGK +* Workspace needed by DSTEVX: +* WORK( ITEMP: ): 2*5*NTGK * IWORK( 1: ): 2*6*NTGK * - CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ), - $ WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK, - $ ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ), - $ Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ), + CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ), + $ WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK, + $ ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ), + $ Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ), $ IWORK( IIWORK ), IWORK( IIFAIL ), $ INFO ) IF( INFO.NE.0 ) THEN @@ -601,7 +601,7 @@ RETURN END IF EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) ) -* +* IF( NSL.GT.0 .AND. WANTZ ) THEN * * Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:), @@ -615,22 +615,22 @@ IF( NSL.GT.1 .AND. $ VUTGK.EQ.ZERO .AND. $ MOD(NTGK,2).EQ.0 .AND. - $ EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN + $ EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN * * D=0 at the top or bottom of the active submatrix: -* one eigenvalue is equal to zero; concatenate the -* eigenvectors corresponding to the two smallest +* one eigenvalue is equal to zero; concatenate the +* eigenvectors corresponding to the two smallest * eigenvalues. * Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) = $ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) + $ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) - Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) = - $ ZERO + Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) = + $ ZERO * IF( IUTGK*2.GT.NTGK ) THEN * Eigenvalue equal to zero or very small. * NSL = NSL - 1 -* END IF +* END IF END IF * DO I = 0, MIN( NSL-1, NRU-1 ) @@ -639,20 +639,20 @@ INFO = N*2 + 1 RETURN END IF - CALL DSCAL( NRU, ONE/NRMU, + CALL DSCAL( NRU, ONE/NRMU, $ Z( IROWU,ICOLZ+I ), 2 ) IF( NRMU.NE.ONE .AND. $ ABS( NRMU-ORTOL )*SQRT2.GT.ONE ) $ THEN DO J = 0, I-1 - ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ), + ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ), $ 2, Z( IROWU, ICOLZ+I ), 2 ) - CALL DAXPY( NRU, ZJTJI, + CALL DAXPY( NRU, ZJTJI, $ Z( IROWU, ICOLZ+J ), 2, $ Z( IROWU, ICOLZ+I ), 2 ) END DO NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 ) - CALL DSCAL( NRU, ONE/NRMU, + CALL DSCAL( NRU, ONE/NRMU, $ Z( IROWU,ICOLZ+I ), 2 ) END IF END DO @@ -662,7 +662,7 @@ INFO = N*2 + 1 RETURN END IF - CALL DSCAL( NRV, -ONE/NRMV, + CALL DSCAL( NRV, -ONE/NRMV, $ Z( IROWV,ICOLZ+I ), 2 ) IF( NRMV.NE.ONE .AND. $ ABS( NRMV-ORTOL )*SQRT2.GT.ONE ) @@ -670,12 +670,12 @@ DO J = 0, I-1 ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ), $ 2, Z( IROWV, ICOLZ+I ), 2 ) - CALL DAXPY( NRU, ZJTJI, + CALL DAXPY( NRU, ZJTJI, $ Z( IROWV, ICOLZ+J ), 2, $ Z( IROWV, ICOLZ+I ), 2 ) END DO NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 ) - CALL DSCAL( NRV, ONE/NRMV, + CALL DSCAL( NRV, ONE/NRMV, $ Z( IROWV,ICOLZ+I ), 2 ) END IF END DO @@ -684,18 +684,18 @@ $ MOD(NTGK,2).GT.0 ) THEN * * D=0 in the middle of the active submatrix (one -* eigenvalue is equal to zero): save the corresponding +* eigenvalue is equal to zero): save the corresponding * eigenvector for later use (when bottom of the * active submatrix is reached). * SPLIT = .TRUE. - Z( IROWZ:IROWZ+NTGK-1,N+1 ) = + Z( IROWZ:IROWZ+NTGK-1,N+1 ) = $ Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) - Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) = - $ ZERO - END IF + Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) = + $ ZERO + END IF END IF !** WANTZ **! -* +* NSL = MIN( NSL, NRU ) SVEQ0 = .FALSE. * @@ -706,17 +706,17 @@ END DO * * Update pointers for TGK, S and Z. -* +* ISBEG = ISBEG + NSL IROWZ = IROWZ + NTGK ICOLZ = ICOLZ + NSL IROWU = IROWZ - IROWV = IROWZ + 1 + IROWV = IROWZ + 1 ISPLT = IDPTR + 1 NS = NS + NSL NRU = 0 - NRV = 0 - END IF !** NTGK.GT.0 **! + NRV = 0 + END IF !** NTGK.GT.0 **! IF( IROWZ.LT.N*2 .AND. WANTZ ) THEN Z( 1:IROWZ-1, ICOLZ ) = ZERO END IF @@ -726,7 +726,7 @@ * Bring back eigenvector corresponding * to eigenvalue equal to zero. * - Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) = + Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) = $ Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) + $ Z( IDBEG:IDEND-NTGK+1,N+1 ) Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0 @@ -735,7 +735,7 @@ IROWU = IROWU + 1 IDBEG = IEPTR + 1 SVEQ0 = .FALSE. - SPLIT = .FALSE. + SPLIT = .FALSE. END IF !** Check for split in E **! END DO !** IEPTR loop **! * @@ -757,7 +757,7 @@ IF( WANTZ ) CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 ) END IF END DO -* +* * If RANGE=I, check for singular values/vectors to be discarded. * IF( INDSV ) THEN @@ -767,7 +767,7 @@ IF( WANTZ ) Z( 1:N*2,K+1:NS ) = ZERO NS = K END IF - END IF + END IF * * Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ). * If B is a lower diagonal, swap U and V.