# To unpack, sh this file. # Please, read the file readme.doc first. cat > readme.doc <<'CUT HERE..........' Accurate Symmetric Eigensolvers ------------------------------- DSYEVJ and SSYEVJ ------------------- Ivan Slapnicar Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture R. Boskovica b.b 58000 Split, Croatia e-mail: islapnicar at uni-zg.ac.mail.yu December 4, 1992 1. Contents of the Package ======================= You have received this file, readme.doc, and the following FORTRAN files: dsyevj.f - contains the subroutine DSYEVJ dgjgt.f - contains the subroutine DGJGT djgjf.f - contains the subroutine DJGJF ssyevj.f - contains the subroutine SSYEVJ sgjgt.f - contains the subroutine SGJGT sjgjf.f - contains the subroutine SJGJF djac.f - contains the subroutine DJAC sjac.f - contains the subroutine SJAC blasj.f - contains BLAS-type subroutines (single and double precision) lamch.f - contains LAPACK auxiliary functions DLAMCH and SLAMCH which determine the machine precision gensym.f - contains the subroutine GENSYM auxil.f - contains auxiliary subroutines dtestj.f - contains the test program DTESTJ stestj.f - contains the test program STESTJ compar.f - contains the test program COMPAR 2. Description =========== Subroutine DSYEVJ computes the non-zero eigenvalues and the corresponding eigenvectors of a DOUBLE PRECISION symmetric matrix. DSYEVJ implements symmetric indefinite decomposition (subroutine DGJGT) followed by implicit (one-sided) Jacobi method (subroutine DJGJF). Subroutines DGJGT and DJGJF can also be used separately. DSYEVJ computes the eigenvalues of non-singular matrices with high relative accuracy. For more details and for references, see the comments of DSYEVJ, DGJGT and DJGJF. SSYEVJ is a SINGLE PRECISION version of DSYEVJ. It uses subroutines SGJGT and SJGJF. DJAC and SJAC are double and single precision implementations of the standard Jacobi algorithms. DTESTJ and STESTJ are short test programs used to test the subroutines DSYEVJ and SSYEVJ, respectively. COMPAR is a more sophisticated test program. It compares our algorithm with the QR and the standard Jacobi algorithm. More precisely, COMPAR compares eigensolutions obtained by DSYEVJ, SSYEVJ, DJAC, SJAC, and the LAPACK driver routines DSYEV and SSYEV which can be obtained from NETLIB. (DSYEV and SSYEV implement tridiagonalization followed by QR iteration.) In COMPAR, user can add one more routine or substitute some routines. For details see the comments of COMPAR. 3. Linking Instructions ==================== All files need to be linked with the level one BLAS. a) Subroutines ----------- DSYEVJ - link the files dsyevj.o, dgjgt.o, djgjf.o, blasj.o, lamch.o. SSYEVJ - link the files ssyevj.o, sgjgt.o, sjgjf.o, blasj.o, lamch.o. DGJGT - link the files dgjgt.o and blasj.o. SGJGT - link the files sgjgt.o and blasj.o. DJGJF - link the files djgjf.o, blasj.o and lamch.o. SJGJF - link the files sjgjf.o, blasj.o and lamch.o. DJAC - link the files djac.o, blasj.o and lamch.o. SJAC - link the files sjac.o, blasj.o and lamch.o. b) Test Programs ------------- DTESTJ - link the file dtestj.o with the file auxil.o, and with everything needed for DSYEVJ above. STESTJ - link the file stestj.o with the file auxil.o, and with everything needed for SSYEVJ above. COMPAR - link the file compar.o with the files dsyevj.o, dgjgt.o, djgjf.o, ssyevj.o, sgjgt.o, sjgjf.o, (*) djac.o, sjac.o, blasj.o, gensym.o, auxil.o, and LAPACK driver routines DSYEV and SSYEV which can be obtained from NETLIB. (DSYEV and SSYEV already contain DLAMCH and SLAMCH, respectively, so there is no need to link the file lamch.o) User can also make a library from all files in (*). 4. Remark ====== We do not claim that our programs are optimal and free of errors. We appreciate all comments, and especially the ones concerning detected errors and/or bugs, and test results on various machines (our tests were preformed mainly on an IBM RISC/6000). CUT HERE.......... cat > dsyevj.f <<'CUT HERE..........' SUBROUTINE DSYEVJ( N, RANK, LDA, H, EV, V, PIVOT, JOB ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992. * * .. Scalar Arguments INTEGER N, RANK, LDA CHARACTER JOB * .. * .. Array Arguments DOUBLE PRECISION H( LDA, * ), EV( * ), V( * ) INTEGER PIVOT( * ) * .. * * Purpose * ======= * * DSYEVJ computes the non-zero eigenvalues and the corresponding * eigenvectors of a real symmetric matrix H of order N. * * ! If H is non-singular and highly conditioned, but the measure * ! C( A, A ) is small, DSYEVJ computes all eigenvalues with high * ! relative accuracy. More precisely, the computed eigenvalues * ! have approximately LOG10( 1 / EPS ) - LOG10( C( A, A ) ) * ! accurate decimal digits, where EPS is the machine precision. * ! For details see [1]. * * If, on exit, RANK < N, the computed eigensolution may be inaccurate. * * * DSYEVJ first uses the subroutine DGJGT (symmetric indefinite * decomposition with complete pivoting) to decompose H as * * P * H * P' = G * J * G' * * where ' denotes the transpose, P is a permutation matrix, * G is a N x RANK real matrix with full column rank, * and J is direct sum of I( NPOS ) and ( - I( RANK - NPOS ) ). * Here I denotes an identity matrix, and NPOS is the number of * positive eigenvalues of H. * * On exit from DGJGT, RANK < N only if DGJGT encountered a * submatrix which was exactly zero. * * * DSYEVJ then uses BLASJ subroutine DIPERR to apply inverse * of the permutation stored in PIVOT to rows of the factor G, * so that H = G * J * G'. * * * DSYEVJ finally uses the subroutine DJGJF (implicit J-orthogonal * Jacobi on the pair G, J) to compute the RANK of the matrix H, * non-zero eigenvalues of H, and the corresponding eigenvectors. * * * This method for solving the symmetric eigenvalue problem was * originally proposed in [2]. The formal algorithms and the error * analysis are given in [1]. * * * References * ========== * * [1] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * [2] K. Veselic: * An Eigenreduction Algorithm for Definite Matrix Pairs, preprint, * FB Mathematik, Fernuniversitaet Hagen, Hagen 1986, 1989, 1992. * * * Arguments * ========= * * N (input) INTEGER * Dimension of the matrix H. N >= 0. * * RANK (output) INTEGER * Rank of the matrix H. N >= RANK >= 0. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * H (input/output) DOUBLE PRECISION array, dimension (LDA,N). * On entry a N x N real symmetric matrix. * On exit, the first RANK columns of H contain the normalized * eigenvectors which correspond to the non-zero eigenvalues * of H. Here Rank( H ) = RANK. * * EV (output) DOUBLE PRECISION array, dimension (N) * The first RANK entries contain the non-zero eigenvalues of H. * * V (workspace) DOUBLE PRECISION array, dimension (N) * Used by DGJGT and DJGJF. * * PIVOT (workspace) INTEGER array, dimension (N) * Used by DGJGT, DIPERR and DJGJF. * * JOB (input) CHARACTER * If JOB = 'f' or 'F', then DJGJF uses fast rotations. * (This is the default choice.) * If JOB = 's' or 'S', then DJGJF uses fast self-scaling * rotations. (These rotations should be used if * underflow/owerflow occur in keeping the diagonal of the * fast rotations, which is unlikely to happen.) * * Further Details * =============== * * DSYEVJ uses subroutines DGJGT, DIPERR and DJGJF, which further use * BLAS1 subroutines DSWAP, DCOPY, DSCAL and DROT, * BLAS1 functions IDAMAX and DDOT, * BLASJ subroutines DROTJJ, DROTF and DROTFS, and * LAPACK auxiliary function DLAMCH. * * BLASJ contains BLAS-type routines for use by J-orthogonal Jacobi * methods. DLAMCH determines the machine precision. * * ============================================================= * * .. * .. Local Scalars .. INTEGER INFO, NPOS * .. * .. External Subroutines .. EXTERNAL DGJGT, DIPERR, DJGJF * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 3 ELSE IF( ( JOB .NE. 'f') .AND. ( JOB .NE. 'F' ) .AND. $ ( JOB .NE. 's' ) .AND. ( JOB .NE. 'S' ) ) THEN INFO = 8 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to DSYEVJ parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( N .EQ. 0 ) RETURN * * Decompose the matrix H. On exit from DGJGT, factor G is * stored in the first RANK columns of H, where Rank ( H ) = RANK. * The matrix J is implicitly defined by RANK and NPOS. * CALL DGJGT( N, RANK, LDA, H, NPOS, PIVOT, V ) * * Apply inverse of the permutation stored in PIVOT to rows of * the factor G (stored in array H), so that H = G * J * G'. * CALL DIPERR( N, RANK, LDA, H, PIVOT ) * * Calculate the rank of H, non-zero eigenvalues of H, and the * corresponding eigenvectors. * CALL DJGJF( N, RANK, LDA, H, NPOS, EV, V, PIVOT, JOB ) RETURN * * End of DSYEVJ * END CUT HERE.......... cat > dgjgt.f <<'CUT HERE..........' SUBROUTINE DGJGT( N, RANK, LDA, H, NPOS, PIVOT, JJ ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 20, 1992 * * .. Scalar Arguments INTEGER N, RANK, LDA, NPOS * .. * .. Array Arguments DOUBLE PRECISION H( LDA, * ), JJ( * ) INTEGER PIVOT( * ) * .. * * Purpose * ======= * * Subroutine DGJGT (symmetric indefinite decomposition with complete * pivoting) decomposes a real symmetric matrix H( N,N ) as * * P * H * P' = G * J * G' * * where ' denotes the transpose, P is the permutation matrix * of complete pivoting, G is a N x RANK real matrix with full * column rank, rank( H ) = RANK, and J is direct sum of I( NPOS ) * and ( - I( RANK - NPOS ) ). Here I denotes an identity matrix, * and NPOS is the number of the positive eigenvalues of H. * * RANK < N only if DGJGT encountered a submatrix which was * exactly zero. * * The algorithm is a modification of the symmetric indefinite * decomposition from [1]. The algorithm and its error annalysis are * given in [2]. * * If the initial matrix H is positive definite, DGJGT reduces to the * Cholesky decomposition with complete pivoting. * * If we follow DGJGT by DJGJF (implicit J-orthogonal Jacobi on * the pair G, J), we obtain non-zero eigenvalues and the * corresponding eigenvectors of the initial matrix H. * * * References * ========== * * [1] J. R. Bunch, B. N. Parlett: * Direct Methods for Solving Symmetric Indefinite Systems of Linear * Equations, SIAM J. Num. Anal., Vol. 8, No. 4, (639-655) 1971, * * [2] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * * Arguments * ========= * * N (input) INTEGER * The order of the matrix H. N >= 0. * * RANK (output) INTEGER * Rank of the matrix H. N >= RANK >= 0. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * H (input/output) DOUBLE PRECISION array, dimension (LDA,N). * On entry real symmetric matrix of order N. * On exit, the first RANK columns of H contain the factor G. * * NPOS (output) INTEGER * number of positive eigenvalues of H. * (R and NPOS define implicitly the matrix J.) * * PIVOT (output) INTEGER array, dimension (N) * defines implicitly the pivot matrix P, such that, * in the MATLAB sense, G J G' = H( PIVOT, PIVOT) . * * JJ (workspace) DOUBLE PRECISION array, dimension (N) * contains the matrix J. * * Further Details * =============== * * DGJGT uses BLAS1 subroutines DSWAP, DSCAL and DROT, * BLAS1 integer function IDAMAX, and * BLASJ subroutine DROTJJ. * * BLASJ contains BLAS-type routines for use by J-orthogonal * Jacobi methods. * * ============================================================= * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER INFO, NN, NN1, NN2, ISTEP, ISTEP1, ISTEP2 INTEGER IDIAG, IDIAG1, IDIAG2, I, J, IOFF, JOFF DOUBLE PRECISION ALPHA, MAXDIA, MAXOFF, TEMP, A, B, C, S, T * .. * .. External Subroutines .. EXTERNAL DSWAP, DSCAL, DROT, DROTJJ * .. * .. External Functions .. INTEGER IDAMAX EXTERNAL IDAMAX * .. * .. Intrinsic Functions .. INTRINSIC MAX, DSQRT, DABS * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 3 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to DGJGT parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( N .EQ. 0 ) THEN RANK = 0 RETURN END IF * * Initialize JJ, PIVOT, NPOS and ALPHA. * DO 10 I = 1, N JJ( I ) = - ONE 10 PIVOT( I ) = I NPOS=0 ALPHA=( ONE + DSQRT(17.0D0) ) / 8.0D0 * * Start the main loop from N downto 1. * NN = N ISTEP = 1 20 CONTINUE * * Find the absolutely largest diagonal and off-diagonal elements * of the current block of H, and their indices. * IDIAG = ISTEP - 1 + IDAMAX( NN, H( ISTEP, ISTEP ), LDA + 1 ) MAXDIA = DABS( H( IDIAG, IDIAG ) ) MAXOFF = ZERO ISTEP1 = ISTEP + 1 DO 30 I = ISTEP1, N DO 30 J = ISTEP, I - 1 IF( DABS( H( I, J ) ) .LE. MAXOFF ) GO TO 35 MAXOFF = DABS( H( I, J ) ) IOFF = I JOFF = J 35 CONTINUE 30 CONTINUE * * Choose 1 by 1 or 2 by 2 pivot. * IF( MAXDIA. LT . ALPHA * MAXOFF ) GO TO 100 * * 1 by 1 pivot. * * If MAXDIA = 0, then the rest of the current matrix H is zero. * This means that the initial H is singular with rank(H) = ISTEP-1. * IF( MAXDIA .LE. ZERO) THEN RANK = ISTEP - 1 IF( RANK.LT.0) RANK = 0 GO TO 1000 END IF * * Bring the element H( IDIAG, IDIAG ) to the position ISTEP, ISTEP * and notify this in the vector PIVOT. * CALL DSWAP( NN, H( ISTEP, ISTEP ), 1, H( ISTEP, IDIAG ), 1 ) CALL DSWAP( N, H( ISTEP,1 ), LDA, H( IDIAG,1 ), LDA ) I = PIVOT( ISTEP ) PIVOT( ISTEP ) = PIVOT( IDIAG ) PIVOT( IDIAG ) = I * * Update JJ( ISTEP ), and NPOS when necessary. * TEMP = H( ISTEP, ISTEP ) IF( TEMP .GT. ZERO ) THEN JJ( ISTEP ) = ONE NPOS = NPOS + 1 ENDIF * * Update the elements H( ISTEP:N, ISTEP) of the ISTEP-th column of H. * TEMP = DSQRT( DABS( TEMP ) ) H( ISTEP, ISTEP ) = TEMP NN1 = NN - 1 CALL DSCAL( NN1, JJ( ISTEP ) / TEMP, H( ISTEP1, ISTEP ), 1 ) * * Set H( ISTEP, ISTEP1:N ) to zero. * CALL DSCAL( NN1, ZERO, H( ISTEP, ISTEP1 ), LDA ) DO 50 I = ISTEP1, N DO 50 J = ISTEP1, I H( I,J ) = H( I,J ) - H( I,ISTEP ) * H( J,ISTEP ) * JJ(ISTEP) H( J,I ) = H( I,J ) 50 CONTINUE * * Either go to the next step ( GO TO 20), or finish the * decomposition ( GO TO 1000). * NN = NN1 ISTEP = ISTEP1 IF( NN .GT. 0 ) GO TO 20 RANK = N GO TO 1000 * * 2 by 2 pivot. * * Bring the element H( JOFF, JOFF ) to the position ISTEP, ISTEP, * bring the element H( IOFF, IOFF ) to the position ISTEP1, ISTEP1, * and notify this in vector PIVOT. * 100 CONTINUE CALL DSWAP( NN, H( ISTEP, ISTEP ), 1, H( ISTEP, JOFF ), 1 ) CALL DSWAP( NN, H( ISTEP, ISTEP1 ), 1, H( ISTEP, IOFF ), 1 ) I = PIVOT( ISTEP ) PIVOT( ISTEP ) = PIVOT( JOFF ) PIVOT( JOFF ) = I CALL DSWAP( N, H( ISTEP, 1 ), LDA, H( JOFF, 1 ), LDA ) CALL DSWAP( N, H( ISTEP1, 1 ), LDA, H( IOFF, 1 ), LDA ) I = PIVOT( ISTEP1 ) PIVOT( ISTEP1 ) = PIVOT( IOFF ) PIVOT( IOFF ) = I * * Calculate the eigenvalue decomposition of the pivot submatrix, * and initialize ISTEP2 and NN2. * A = H( ISTEP, ISTEP ) B = H( ISTEP1, ISTEP1 ) TEMP = H( ISTEP1, ISTEP ) CALL DROTJJ( A, B, TEMP, C, S, T, -ONE ) A = A - T * TEMP B = B + T * TEMP ISTEP2 = ISTEP1 + 1 NN2 = NN - 2 * * Update JJ( ISTEP ), JJ( ISTEP1 ) and NPOS. * IF( A .GT. ZERO ) JJ( ISTEP ) = ONE IF( B .GT. ZERO ) JJ( ISTEP1 ) = ONE NPOS = NPOS + 1 * * Update the elements H( ISTEP:N, ISTEP:ISTEP1) of the columns * ISTEP and ISTEP1 of the current matrix H. * A = DSQRT( DABS( A ) ) B = DSQRT( DABS( B ) ) H( ISTEP, ISTEP ) = C * A H( ISTEP1, ISTEP) = -S * A H( ISTEP, ISTEP1) = S * B H( ISTEP1, ISTEP1 ) = C * B CALL DROT( NN2, H( ISTEP2,ISTEP ), 1, H( ISTEP2,ISTEP1 ), 1,C,-S) CALL DSCAL( NN2, JJ( ISTEP ) / A, H( ISTEP2, ISTEP ), 1) CALL DSCAL( NN2, JJ( ISTEP1 ) / B, H( ISTEP2, ISTEP1 ), 1) * * Set H( ISTEP : ISTEP1, ISTEP2 : N ) to zero. * CALL DSCAL( NN2, ZERO, H( ISTEP, ISTEP2 ), LDA ) CALL DSCAL( NN2, ZERO, H( ISTEP1, ISTEP2 ), LDA ) * * Update the remaining part of the matrix H, * H( ISTEP2 : N, ISTEP2 : N ). * DO 60 I = ISTEP2, N DO 60 J = ISTEP2, I H( I,J ) = H( I,J ) - H( I, ISTEP ) * H( J, ISTEP )*JJ(ISTEP) $ - H( I, ISTEP1 ) * H( J, ISTEP1 ) * JJ( ISTEP1 ) H( J,I ) = H( I,J ) 60 CONTINUE * * Go to the next step, if necessary. * NN = NN2 ISTEP = ISTEP2 IF( NN .GT. 0) GO TO 20 RANK = N * * Finally, permute the columns of H, in order to sort the matrix J. * 1000 CONTINUE J = NPOS + 1 DO 70 I = 1, NPOS IF( JJ( I ) .GT. ZERO ) GO TO 70 75 IF (JJ( J ) .GT. ZERO ) GO TO 80 J = J + 1 GO TO 75 80 CONTINUE CALL DSWAP( N, H( 1,I ), 1, H( 1,J ), 1 ) J = J + 1 70 CONTINUE RETURN * * End of DGJGT. * END CUT HERE.......... cat > djgjf.f <<'CUT HERE..........' SUBROUTINE DJGJF( N, M, LDA, G, NPOS, EV, V, JJ, JOB ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, NPOS CHARACTER JOB * .. * .. Array Arguments INTEGER JJ( * ) DOUBLE PRECISION G( LDA, * ), EV( * ), V( * ) * .. * * Purpose * ======= * * Subroutine DJGJF (implicit J-orthogonal Jacobi method on the * pair G, J) computes the rank, the non-zero eigenvalues, and the * corresponding eigenvectors of the implicitly defined N x N real * symmetric matrix * * H = G * J * G' , * * where ' denotes the transpose, G is a N x M real matrix, * and J is direct sum of I( NPOS ) and ( -I( M - NPOS ) ). Here * I denotes an identity matrix. * * * ! If the matrix H is highly conditioned, but the matrix * ! Scal( G ), where G = Scal( G ) * D, has small condition, then * ! SJGJF computes all eigenvalues with high relative accuracy. * ! Here D is diagonal positive definite, and the columns of * ! Scal( G ) have unit norms. More precisely, the computed * ! eigenvalues have approximately * ! LOG10( 1 / EPS ) - LOG10( Cond( B ) ) accurate decimal digits, * ! where EPS is the machine precision. For more details see [1]. * * * On entry, the parameter M is the number of columns of G. On exit, * M returns the rank of G, which is also the rank of H. If, on exit, * M is smaller than its starting value (deflation), then the initial * G was (almost) singular, and the results may be inaccurate. * * DJGJF implements the implicit (one-sided) J-orthogonal Jacobi method * for solving Det( G' * G - Lambda * J ) = 0, originally defined in * [2]. The formal algorithm and its error analysis are given in [1]. * The global and the quadratic convergence have been proved in [2] and * [3], respectively. * * The implicit method works only on the matrix G, thus obtaining the * eigenvectors automatically. The diagonal of the matrix G' * G, which * is needed to determine the transformation parameters, is kept and * updated in a separate vector. * * DJGJF performs hyperbolic plane rotation if the pivot indices I < J * satisfy ( I .LE. NPOS ) .AND. ( J .GT. NPOS ), and trigonometric * rotation otherwise. E.g. for N = M = 5 and NPOS = 3 the rotation * scheme looks like * * ( x T T H H ) * ( x T H H ) T = trigonometric rotation * ( x H H ) * ( x T ) H = hyperbolic rotation * ( x ) * * DJGJF uses either fast or fast self-scaling rotations. * * Treshold and convergence tests are both made after the relative * error analysis of [1]. Let * * ( A CC ) * ( CC B ) * * be the pivot submatrix. We rotate if * * SQRT( ABS (CC) ) * ------------------------ > N * EPS , * SQRT( A ) * SQRT ( B ) * * where EPS is the machine precision. The algorithm has converged * if M * ( M - 1 ) / 2 succesive rotations were not performed. * * * References * ========== * * [1] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * [2] K. Veselic: * An Eigenreduction Algorithm for Definite Matrix Pairs, preprint, * FB Mathematik, Fernuniversitaet Hagen, Hagen 1986, 1989, 1992. * * [3] Z. Drmac, V. Hari, On quadratic convergence bounds for the * J-symmetric Jacobi method, to appear in Numer. Math. * * * Arguments * ========= * * N (input) INTEGER * Number of rows of the matrix G. N >= 0. * * M (input/output) INTEGER * On entry, number of columns of the matrix G. N >= M >= 0. * On exit, M = rank( G ) = rank ( H ) = number of non-zero * eigenvalues of H. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * G (input/output) DOUBLE PRECISION array, dimension (LDA,M). * On entry a N x M real matrix. * On exit, the first M (M can change!) columns of G contain * the normalized eigenvectors corresponding to the non-zero * eigenvalues of the matrix H = G * J * G'. * * NPOS (input/output) INTEGER, M >= NPOS >= 0 * On entry, number of +1's in the matrix J. * On exit, number of +1's in the matrix J and number of the * positive eigenvalues of the matrix H. * (On entry and on exit, M and NPOS define implicitly the * matrix J.) * * EV (output) DOUBLE PRECISION array, dimension (M) * On exit, EV contains the non-zero eigenvalues of the matrix * H = G J G'. During the computation EV contains the diagonal * of the matrix G' * G. * * V (workspace) DOUBLE PRECISION array, dimension (M) * accumulates the diagonal of the fast rotations. * * JJ (workspace) INTEGER array, dimension (M) * JJ( I ) = 0 if the I-th column of G is non-zero, * JJ( I ) = 1 if the I-th column of G is zero. * The rotations which involve zero columns are skipped. * * JOB (input) CHARACTER * If JOB = 'f' or 'F', then DJGJ uses fast rotations. * If JOB = 's' or 'S', then DJGJ uses fast self-scaling * rotations. (These rotations should be used if * underflow/owerflow occur in the diagonal of fast * rotations, which is unlikely to happen.) * * * Further Details * =============== * * DJGJF uses BLAS1 subroutines DSCAL and DCOPY, * BLAS1 function DDOT, * BLASJ subroutines DROTJJ, DROTF and DROTFS, and * LAPACK auxiliary function DLAMCH. * * BLASJ contains BLAS-type routines for use by J-orthogonal JACOBI * methods. DLAMCH determines the machine precision. * * ============================================================= * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL NOTSS INTEGER INFO, RANK, NP, I, J, ISWEEP INTEGER NTEMP, NMAX, SINGUL, IDI, IDJ, IDIDJ DOUBLE PRECISION EPS, ONEEPS, A, B, CC, HH, HYP, C, S, T * .. * .. External Subroutines .. EXTERNAL DROTJJ, DROTF, DROTFS, DCOPY, DSCAL * .. * .. External Functions .. DOUBLE PRECISION DDOT, DLAMCH EXTERNAL DDOT, DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC MAX, DSQRT, DABS, DBLE, IDINT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( ( M .LT. 0 ) .OR. ( M .GT. N ) ) THEN INFO = 2 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 3 ELSE IF( ( NPOS. LT. 0 ) .OR. ( NPOS .GT. M ) ) THEN INFO = 5 ELSE IF( ( JOB .NE. 'f') .AND. ( JOB .NE. 'F' ) .AND. $ ( JOB .NE. 's' ) .AND. ( JOB .NE. 'S' ) ) THEN INFO = 9 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to DJGJF parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( ( N .EQ. 0 ) .OR. ( M .EQ. 0 ) ) RETURN IF( M .EQ. 1 ) THEN EV( 1 ) = DDOT( N, G( 1, 1 ), 1, G( 1, 1 ), 1) IF( EV( 1 ) .LE. ZERO ) THEN * * We have a zero matrix. * M = 0 ELSE * * Compute the only eigenvector and set the sign of the eigenvalue. * CALL DSCAL( N, ONE / DSQRT( EV( 1 ) ), G( 1,1 ), 1 ) IF( NPOS .EQ. 0 ) EV( 1 ) = - EV( 1 ) ENDIF RETURN ENDIF * * Calculate the treshold N * EPS, and ONEEPS = ONE - SQRT( EPS ). * EPS = DLAMCH('E') ONEEPS = ONE - DSQRT( EPS ) EPS = DBLE( N ) * EPS * * Set NMAX to size of the sweep, and set NTEMP to 0. * NMAX = M * ( M - 1 ) / 2 NTEMP = 0 * * Set parameter NOTSS to .FALSE. if self-scaling rotations are to * be used. * NOTSS = .TRUE. IF( ( JOB .EQ. 's' ) .OR. ( JOB .EQ. 'S' ) ) NOTSS = .FALSE. * * Set vector V where we accumulate the diagonal of fast rotations * to 1, and set vector JJ to zero (at the beginning all columns of * G are treated as non-zero columns). * DO 9 I = 1, M V(I) = ONE JJ( I ) = 0 9 CONTINUE * * Set RANK to M, NP to NPOS, and SINGUL to 0. * SINGUL = 0 means that no singularity is detected in the current * rotation. * SINGUL = 1 means that the updated diagonal may be inaccurate. * SINGUL = 2 means that the singularity is detected. * RANK = M NP = NPOS SINGUL = 0 * * Start the iterations. * DO 40 ISWEEP = 1, 30 * * Refresh the diagonal of the matrix G' * G. If EV( I ) = 0, then * set JJ( I ) = 1 and decrease RANK by 1. If I-th column * corresponds to a positive eigenvalue, decrease NP by 1. * DO 30 I = 1, M IF( JJ( I ) .EQ. 1 ) GO TO 30 EV(I) = DDOT( N, G( 1,I ), 1, G( 1,I ), 1) * V( I ) ** 2 IF( EV( I ) . LE. ZERO ) THEN JJ( I ) = 1 RANK = RANK - 1 IF( I .LE. NPOS ) NP = NP - 1 END IF 30 CONTINUE * * Start the main loop. * DO 50 I = 1, M - 1 DO 50 J = I + 1, M * * If I-th or J-th columns are zero, skip this rotation and test * for convergence. * IF( ( JJ( I ) + JJ( J ) ) .GT. 0 ) GO TO 60 * * Choose trigonometric or hyperbolic rotation. * HYP = - ONE IF( ( I .LE. NPOS ) .AND. ( J .GT. NPOS ) ) HYP = ONE * * Calculate the pivot submatrix. * A = EV( I ) B = EV( J ) CC = DDOT( N, G( 1,I ), 1, G( 1,J ), 1 ) * V( I ) * V( J ) * * Test whether the matrix Scal( G )' * Scal( G ) is positive * definite. This means that HH must be smaller than 1. If * HH >= 1, we compute the diagonal elements A and B from the * columns of G once more and repeat the test. If the test fails * again, G is singular and the computed eigensolution might be * inaccurate. * HH = DABS( CC ) / DSQRT( A ) / DSQRT( B ) IF( HH .GE. ONE ) THEN A = DDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) * V(I) ** 2 B = DDOT( N, G( 1,J ), 1, G( 1,J ), 1 ) * V(J) ** 2 HH = DABS( CC ) / DSQRT( A ) / DSQRT( B ) IF( HH .GE. ONE ) THEN * * G' * G is singular. Set SINGUL = 2. * SINGUL = 2 GO TO 55 END IF END IF * * If the following test is true, the update of the diagonal may be * inaccurate. In this case set SINGUL = 1, which means that the * diagonal elements are later calculted from the columns of G. * IF( HH .GE. ONEEPS ) SINGUL = 1 * * Test whether to perform this rotation. If the rotation is * not performed, set NTEMP = NTEMP + 1. If NTEMP = NMAX, * the algorithm has converged. * IF( HH .GT. EPS ) GO TO 55 60 NTEMP = NTEMP + 1 IF( NTEMP .GT. NMAX ) GO TO 100 GO TO 50 * * Calculate the rotation parameters and set NTEMP = 0. * 55 CALL DROTJJ( A, B, CC, C, S, T, HYP ) NTEMP = 0 * * If in the hyperbolic case, SROTJJ could not calculate the * transformation parameters, it returns C = - 1. In this case * I-th and J-th columns of G both correspond to zero * eigenvalues. Set JJ( I ) and JJ( J ) to 1, decrease RANK * by 2, decrease NP by 1, and go to the next rotation. On exit * NPOS = NP is the number of the positive eigenvalues of H. * IF( C .LT. ZERO ) THEN JJ( I ) = 1 JJ( J ) = 1 RANK = RANK - 2 NP = NP - 1 SINGUL = 0 GO TO 50 END IF * * Update the columns I and J of G, the elements of the * accumulated diagonal of fast rotations, V( I ) and V( J ). * DROTF performs fast J-orthogonal Jacobi rotation. * DROTFS performs fast self-scaling J-orthogonal Jacobi rotation. * * Self--scaling rotation, if used, pushes the element of V, * which is further away from 1, towards 1. For details see [1]. * IF ( NOTSS ) GO TO 44 IDI = 1 IDJ = 1 IDIDJ = 1 IF( V( I ) .LT. ONE ) IDI = - 1 IF( V( J ) .LT. ONE ) IDJ = - 1 IF( V( I ) .LE. V( J ) ) IDIDJ = - 1 IF( ( IDI + IDJ + 2 * IDINT( HYP ) ) * IDIDJ ) 43, 44, 42 * 42 CALL DROTFS ( N, G( 1,I ), 1, G( 1,J ), 1, T * V(I) / V(J), $ HYP * C * S * V( J ) / V( I ) ) V( I ) = V( I ) / C V( J ) = V( J ) * C GO TO 45 * 43 CALL DROTFS( N, G( 1,J ), 1, G( 1,I ), 1, $ HYP * T * V( J ) / V( I ), C * S * V( I ) / V( J ) ) V( I ) = V( I ) * C V( J ) = V( J ) / C GO TO 45 * 44 CALL DROTF( N, G( 1,I ), 1, G( 1,J ), 1, T * V( I ) / V( J ), $ HYP * T * V( J ) / V( I ) ) V( I ) = V( I ) * C V( J ) = V( J ) * C * * Both, hyperbolic and trigonometric rotations decrease the smaller * diagonal element. If the matrix is singular (SINGUL = 2), then * the column corresponding to the smaller diagonal element becomes * after rotation a zero-column. Set the corresponding element of JJ * to 1 and decrease RANK by 1. If this column corresponds to a * positive eigenvalue, decrease NP by 1. The non-zero diagonal * element is computed from the corresponding column of G. * 45 CONTINUE IF( SINGUL .EQ. 2 ) THEN IF( A .LE. B ) THEN JJ( I ) = 1 IF( I .LE. NPOS ) NP = NP - 1 EV( J ) = DDOT( N, G( 1,J ), 1, G( 1,J ), 1 ) * V(J) ** 2 ELSE JJ( J ) = 1 IF( J .LE. NPOS ) NP = NP - 1 EV( I ) = DDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) * V(I) ** 2 END IF RANK = RANK - 1 SINGUL = 0 GO TO 50 END IF * * Update the diagonal elements of G' * G, EV( I ) and EV( J ). * IF( SINGUL .EQ. 0 ) THEN EV( I ) = A + HYP * CC * T EV( J ) = B + CC * T ELSE * * This point is reached if SINGUL = 1. This means that the * updates of the diagonal elements might be inaccurate. The * diagonal elements are therefore computed from the I-th and * J-th column of G. Although no singularity was detected, one * of the diagonal elements may have become zero. In this case, * set the corresponding element of JJ to 1 and decrease RANK * by 1. If this column corresponds to a positive eigenvalue, * decrease NP by 1. * EV( I ) = DDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) * V(I) ** 2 IF( EV( I ) .LE. ZERO ) THEN JJ( I ) = 1 RANK = RANK - 1 IF( I .LE. NPOS ) NP = NP - 1 END IF EV( J ) = DDOT( N, G( 1,J ), 1, G( 1,J ), 1 ) * V(J) ** 2 IF( EV( J ) .LE. ZERO ) THEN JJ( J ) = 1 RANK = RANK - 1 IF( J .LE. NPOS ) NP = NP - 1 END IF SINGUL = 0 ENDIF * * End of the main loop. * 50 CONTINUE * * End of the ISWEEP loop. * 40 CONTINUE * * If this point is reached, the algorithm did not converge. * WRITE(*,*) ' DJGJF executed ', ISWEEP, ' sweeps without', $ ' convergence. ' 100 CONTINUE * * If no deflation occured, skip the next part. * IF( M .EQ. RANK ) GO TO 200 * * The deflation has occured. The results might be inaccurate! * Move NP non-zero columns of G to the first NP columns, * and do the same for the corresponding elements of the * vectors EV and V. * DO 110 I = 1, NP IF( JJ( I ) .EQ. 1 ) THEN J = I + 1 120 CONTINUE IF( JJ( J ) .EQ. 0 ) THEN CALL DCOPY( N, G( 1, J ), 1, G( 1, I ), 1 ) EV( I ) = EV( J ) V( I ) = V( J ) JJ( J ) = 1 GO TO 110 END IF J = J + 1 GO TO 120 END IF 110 CONTINUE * * Move the rest of the RANK - NP non-zero columns of G to * the columns NP + 1, ... , RANK, and do the same for the * corresponding elements of the vectors EV and V. * DO 140 I = NP + 1, RANK IF( JJ( I ) .EQ. 1 ) THEN J = I + 1 150 CONTINUE IF( JJ( J ) .EQ. 0 ) THEN CALL DCOPY( N, G( 1, J ), 1, G( 1, I ), 1 ) EV( I ) = EV( J ) V( I ) = V( J ) JJ( J ) = 1 GO TO 140 END IF J = J + 1 GO TO 150 END IF 140 CONTINUE * * Set M to rank of G, and NPOS to the number of positive * eigenvalues of H. * M = RANK NPOS = NP * * Calculate the eigenvalues and the corresponding eigenvectors. * Absolute values of the eigenvalues are the squares of the norms * of the columns of G. The corresponding eigenvectors are the * columns of G divided by their norms. * 200 CONTINUE DO 80 I = 1, M C = DDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) CALL DSCAL( N, ONE / DSQRT( C ), G( 1,I ), 1 ) EV( I ) = C * V( I ) ** 2 80 CONTINUE * * Change the signs for the negative eigenvalues. * DO 70 I = NPOS + 1, M EV( I ) = - EV( I ) 70 CONTINUE RETURN * * End of DJGJF * END CUT HERE.......... cat > ssyevj.f <<'CUT HERE..........' SUBROUTINE SSYEVJ( N, RANK, LDA, H, EV, V, PIVOT, JOB ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, RANK, LDA CHARACTER JOB * .. * .. Array Arguments REAL H( LDA, * ), EV( * ), V( * ) INTEGER PIVOT( * ) * .. * * Purpose * ======= * * SSYEVJ computes the non-zero eigenvalues and the corresponding * eigenvectors of a real symmetric matrix H of order N. * * ! If H is non-singular and highly conditioned, but the measure * ! C( A, A ) is small, SSYEVJ computes all eigenvalues with high * ! relative accuracy. More precisely, the computed eigenvalues * ! have approximately LOG10( 1 / EPS ) - LOG10( C( A, A ) ) * ! accurate decimal digits, where EPS is the machine precision. * ! For details see [1]. * * If, on exit, RANK < N, the computed eigensolution may be inaccurate. * * * SSYEVJ first uses the subroutine SGJGT (symmetric indefinite * decomposition with complete pivoting) to decompose H as * * P * H * P' = G * J * G' * * where ' denotes the transpose, P is a permutation matrix, * G is a N x RANK real matrix with full column rank, * and J is direct sum of I (NPOS) and ( -I (RANK-NPOS) ). * Here I denotes an identity matrix, and NPOS is the number of * positive eigenvalues of H. * * On exit from SGJGT, RANK < N only if SGJGT encountered a * submatrix which was exactly zero. * * * SSYEVJ then uses BLASJ subroutine SIPERR to apply inverse * of the permutation stored in PIVOT to rows of the factor G, * so that H = G * J G'. * * * SSYEVJ finally uses the subroutine SJGJF (implicit J-orthogonal * Jacobi on the pair G, J) to compute the RANK of the matrix H, * non-zero eigenvalues of H, and the corresponding eigenvectors. * * * This method for solving the symmetric eigenvalue problem was * originally proposed in [2]. The formal algorithms and the error * analysis are given in [1]. * * * References * ========== * * [1] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * [2] K. Veselic: * An Eigenreduction Algorithm for Definite Matrix Pairs, preprint, * FB Mathematik, Fernuniversitaet Hagen, Hagen 1986, 1989, 1992. * * * Arguments * ========= * * N (input) INTEGER * Dimension of the matrix H. N >= 0. * * RANK (output) INTEGER * Rank of the matrix H. N >= RANK >= 0. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * H (input/output) REAL array, dimension (LDA,N). * On entry a N x N real symmetric matrix. * On exit, the first RANK columns of H contain the normalized * eigenvectors which correspond to the non-zero eigenvalues * of H. Here Rank( H ) = RANK. * * EV (output) REAL array, dimension (N) * The first RANK entries contain the non-zero eigenvalues of H. * * V (workspace) REAL array, dimension (N) * Used by SGJGT and SJGJF. * * PIVOT (workspace) INTEGER array, dimension (N) * Used by SGJGT, SIPERR and SJGJF. * * JOB (input) CHARACTER * If JOB = 'f' or 'F', then SJGJF uses fast rotations. * (This is the default choice.) * If JOB = 's' or 'S', then SJGJF uses fast self-scaling * rotations. (These rotations should be used if * underflow/owerflow occur in keeping the diagonal of the * fast rotations, which is unlikely to happen.) * * Further Details * =============== * * SSYEVJ uses subroutines SGJGT, SIPERR and SJGJF, which further use * BLAS1 subroutines SSWAP, SCOPY, SSCAL and SROT, * BLAS1 functions ISAMAX and SDOT, * BLASJ subroutines SROTJJ, SROTF and SROTFS, and * LAPACK auxiliary function SLAMCH. * * BLASJ contains BLAS-type routines for use by J-orthogonal Jacobi * methods. SLAMCH determines the machine precision. * * ============================================================= * * .. * .. Local Scalars .. INTEGER INFO, NPOS * .. * .. External Subroutines .. EXTERNAL SGJGT, SIPERR, SJGJF * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 3 ELSE IF( ( JOB .NE. 'f') .AND. ( JOB .NE. 'F' ) .AND. $ ( JOB .NE. 's' ) .AND. ( JOB .NE. 'S' ) ) THEN INFO = 8 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to SSYEVJ parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( N .EQ. 0 ) RETURN * * Decompose the matrix H. On exit from DGJGT, factor G is * stored in the first RANK columns of H, where Rank ( H ) = RANK. * The matrix J is implicitly defined by RANK and NPOS. * CALL SGJGT( N, RANK, LDA, H, NPOS, PIVOT, V ) * * Apply inverse of the permutation stored in PIVOT to rows of * the factor G (stored in array H), so that H = G * J * G'. * CALL SIPERR( N, RANK, LDA, H, PIVOT ) * * Calculate the rank of H, non-zero eigenvelues of H, and the * corresponding eigenvectors. * CALL SJGJF( N, RANK, LDA, H, NPOS, EV, V, PIVOT, JOB ) RETURN * * End of SSYEVJ * END CUT HERE.......... cat > sgjgt.f <<'CUT HERE..........' SUBROUTINE SGJGT( N, RANK, LDA, H, NPOS, PIVOT, JJ ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 20, 1992 * * .. Scalar Arguments INTEGER N, RANK, LDA, NPOS * .. * .. Array Arguments REAL H( LDA, * ), JJ( * ) INTEGER PIVOT( * ) * .. * * Purpose * ======= * * Subroutine SGJGT (symmetric indefinite decomposition with complete * pivoting) decomposes a real symmetric matrix H( N,N ) as * * P * H * P' = G * J * G' * * where ' denotes the transpose, P is the permutation matrix * of complete pivoting, G is a N x RANK real matrix with full * column rank, rank( H ) = RANK, and J is direct sum of I( NPOS ) * and ( -I( RANK - NPOS ) ). Here I denotes an identity matrix, * and NPOS is the number of the positive eigenvalues of H. * * RANK < N only if DGJGT encountered a submatrix which was * exactly zero. * * The algorithm is a modification of the symmetric indefinite * decomposition from [1]. The algorithm and its error annalysis are * given in [2]. * * If the initial matrix H is positive definite, DGJGT reduces to the * Cholesky decomposition with complete pivoting. * * If we follow DGJGT by DJGJF (implicit J-orthogonal Jacobi on * the pair G, J), we obtain non-zero eigenvalues and the * corresponding eigenvectors of the initial matrix H. * * * References * ========== * * [1] J. R. Bunch, B. N. Parlett: * Direct Methods for Solving Symmetric Indefinite Systems of Linear * Equations, SIAM J. Num. Anal., Vol. 8, No. 4, (639-655) 1971, * * [2] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * * Arguments * ========= * * N (input) INTEGER * The order of the matrix H. N >= 0. * * RANK (output) INTEGER * Rank of the matrix H. N >= RANK >= 0. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * H (input/output) REAL array, dimension (LDA,N). * On entry real symmetric matrix of order N. * On exit, the first RANK columns of H contain the factor G. * * NPOS (output) INTEGER * number of positive eigenvalues of H. * (R and NPOS define implicitly the matrix J.) * * PIVOT (output) INTEGER array, dimension (N) * defines implicitly the pivot matrix P, such that, * in the MATLAB sense, G J G' = H( PIVOT, PIVOT) . * * JJ (workspace) REAL array, dimension (N) * contains the matrix J. * * Further Details * =============== * * SGJGT uses BLAS1 subroutines SSWAP, SSCAL and SROT, * BLAS1 integer function ISAMAX, and * BLASJ subroutine SROTJJ. * * BLASJ contains BLAS-type routines for use by J-orthogonal * Jacobi methods. * * ============================================================= * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER INFO, NN, NN1, NN2, ISTEP, ISTEP1, ISTEP2 INTEGER IDIAG, IDIAG1, IDIAG2, I, J, IOFF, JOFF REAL ALPHA, MAXDIA, MAXOFF, TEMP, A, B, C, S, T * .. * .. External Subroutines .. EXTERNAL SSWAP, SSCAL, SROT, SROTJJ * .. * .. External Functions .. INTEGER IDAMAX EXTERNAL IDAMAX * .. * .. Intrinsic Functions .. INTRINSIC MAX, SQRT, ABS * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 3 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to SGJGT parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( N .EQ. 0 ) THEN RANK = 0 RETURN END IF * * Initialize JJ, PIVOT, NPOS and ALPHA. * DO 10 I = 1, N JJ( I ) = - ONE 10 PIVOT( I ) = I NPOS=0 ALPHA=( ONE + SQRT(17.0E0) ) / 8.0E0 * * Start the main loop from N downto 1. * NN = N ISTEP = 1 20 CONTINUE * * Find the absolutely largest diagonal and off-diagonal elements * of the current block of H, and their indices. * IDIAG = ISTEP - 1 + ISAMAX( NN, H( ISTEP, ISTEP ), LDA + 1 ) MAXDIA = ABS( H( IDIAG, IDIAG ) ) MAXOFF = ZERO ISTEP1 = ISTEP + 1 DO 30 I = ISTEP1, N DO 30 J = ISTEP, I - 1 IF( ABS( H( I, J ) ) .LE. MAXOFF ) GO TO 35 MAXOFF = ABS( H( I, J ) ) IOFF = I JOFF = J 35 CONTINUE 30 CONTINUE * * Choose 1 by 1 or 2 by 2 pivot. * IF( MAXDIA. LT . ALPHA * MAXOFF ) GO TO 100 * * 1 by 1 pivot. * * If MAXDIA = 0, then the rest of the current matrix H is zero. * This means that the initial H is singular with rank(H) = ISTEP-1. * IF( MAXDIA .LE. ZERO) THEN RANK = ISTEP - 1 IF( RANK.LT.0) RANK = 0 GO TO 1000 END IF * * Bring the element H( IDIAG, IDIAG ) to the position ISTEP, ISTEP * and notify this in the vector PIVOT. * CALL SSWAP( NN, H( ISTEP, ISTEP ), 1, H( ISTEP, IDIAG ), 1 ) CALL SSWAP( N, H( ISTEP,1 ), LDA, H( IDIAG,1 ), LDA ) I = PIVOT( ISTEP ) PIVOT( ISTEP ) = PIVOT( IDIAG ) PIVOT( IDIAG ) = I * * Update JJ( ISTEP ), and NPOS when necessary. * TEMP = H( ISTEP, ISTEP ) IF( TEMP .GT. ZERO ) THEN JJ( ISTEP ) = ONE NPOS = NPOS + 1 ENDIF * * Update the elements H( ISTEP:N, ISTEP) of the ISTEP-th column of H. * TEMP = SQRT( ABS( TEMP ) ) H( ISTEP, ISTEP ) = TEMP NN1 = NN - 1 CALL SSCAL( NN1, JJ( ISTEP ) / TEMP, H( ISTEP1, ISTEP ), 1 ) * * Set H( ISTEP, ISTEP1:N ) to zero. * CALL SSCAL( NN1, ZERO, H( ISTEP, ISTEP1 ), LDA ) DO 50 I = ISTEP1, N DO 50 J = ISTEP1, I H( I,J ) = H( I,J ) - H( I,ISTEP ) * H( J,ISTEP ) * JJ(ISTEP) H( J,I ) = H( I,J ) 50 CONTINUE * * Either go to the next step ( GO TO 20), or finish the * decomposition ( GO TO 1000). * NN = NN1 ISTEP = ISTEP1 IF( NN .GT. 0 ) GO TO 20 RANK = N GO TO 1000 * * 2 by 2 pivot. * * Bring the element H( JOFF, JOFF ) to the position ISTEP, ISTEP, * bring the element H( IOFF, IOFF ) to the position ISTEP1, ISTEP1, * and notify this in vector PIVOT. * 100 CONTINUE CALL SSWAP( NN, H( ISTEP, ISTEP ), 1, H( ISTEP, JOFF ), 1 ) CALL SSWAP( NN, H( ISTEP, ISTEP1 ), 1, H( ISTEP, IOFF ), 1 ) I = PIVOT( ISTEP ) PIVOT( ISTEP ) = PIVOT( JOFF ) PIVOT( JOFF ) = I CALL SSWAP( N, H( ISTEP, 1 ), LDA, H( JOFF, 1 ), LDA ) CALL SSWAP( N, H( ISTEP1, 1 ), LDA, H( IOFF, 1 ), LDA ) I = PIVOT( ISTEP1 ) PIVOT( ISTEP1 ) = PIVOT( IOFF ) PIVOT( IOFF ) = I * * Calculate the eigenvalue decomposition of the pivot submatrix, * and initialize ISTEP2 and NN2. * A = H( ISTEP, ISTEP ) B = H( ISTEP1, ISTEP1 ) TEMP = H( ISTEP1, ISTEP ) CALL SROTJJ( A, B, TEMP, C, S, T, -ONE ) A = A - T * TEMP B = B + T * TEMP ISTEP2 = ISTEP1 + 1 NN2 = NN - 2 * * Update JJ( ISTEP ), JJ( ISTEP1 ) and NPOS. * IF( A .GT. ZERO ) JJ( ISTEP ) = ONE IF( B .GT. ZERO ) JJ( ISTEP1 ) = ONE NPOS = NPOS + 1 * * Update the elements H( ISTEP:N, ISTEP:ISTEP1) of the columns * ISTEP and ISTEP1 of the current matrix H. * A = SQRT( ABS( A ) ) B = SQRT( ABS( B ) ) H( ISTEP, ISTEP ) = C * A H( ISTEP1, ISTEP) = -S * A H( ISTEP, ISTEP1) = S * B H( ISTEP1, ISTEP1 ) = C * B CALL SROT( NN2, H( ISTEP2,ISTEP ), 1, H( ISTEP2,ISTEP1 ), 1,C,-S) CALL SSCAL( NN2, JJ( ISTEP ) / A, H( ISTEP2, ISTEP ), 1) CALL SSCAL( NN2, JJ( ISTEP1 ) / B, H( ISTEP2, ISTEP1 ), 1) * * Set H( ISTEP : ISTEP1, ISTEP2 : N ) to zero. * CALL SSCAL( NN2, ZERO, H( ISTEP, ISTEP2 ), LDA ) CALL SSCAL( NN2, ZERO, H( ISTEP1, ISTEP2 ), LDA ) * * Update the remaining part of the matrix H, * H( ISTEP2 : N, ISTEP2 : N ). * DO 60 I = ISTEP2, N DO 60 J = ISTEP2, I H( I,J ) = H( I,J ) - H( I, ISTEP ) * H( J, ISTEP )*JJ(ISTEP) $ - H( I, ISTEP1 ) * H( J, ISTEP1 ) * JJ( ISTEP1 ) H( J,I ) = H( I,J ) 60 CONTINUE * * Go to the next step, if necessary. * NN = NN2 ISTEP = ISTEP2 IF( NN .GT. 0) GO TO 20 RANK = N * * Finally, permute the columns of H, in order to sort the matrix J. * 1000 CONTINUE J = NPOS + 1 DO 70 I = 1, NPOS IF( JJ( I ) .GT. ZERO ) GO TO 70 75 IF (JJ( J ) .GT. ZERO ) GO TO 80 J = J + 1 GO TO 75 80 CONTINUE CALL SSWAP( N, H( 1,I ), 1, H( 1,J ), 1 ) J = J + 1 70 CONTINUE RETURN * * End of SGJGT. * END CUT HERE.......... cat > sjgjf.f <<'CUT HERE..........' SUBROUTINE SJGJF( N, M, LDA, G, NPOS, EV, V, JJ, JOB ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, NPOS CHARACTER JOB * .. * .. Array Arguments INTEGER JJ( * ) REAL G( LDA, * ), EV( * ), V( * ) * .. * * Purpose * ======= * * Subroutine SJGJF (implicit J-orthogonal Jacobi method on the * pair G, J) computes the rank, the non-zero eigenvalues, and the * corresponding eigenvectors of the implicitly defined N x N real * symmetric matrix * * H = G * J * G' , * * where ' denotes the transpose, G is a N x M real matrix, * and J is direct sum of I( NPOS ) and ( -I( M - NPOS ) ). Here * I denotes an identity matrix. * * * ! If the matrix H is highly conditioned, but the matrix * ! Scal( G ), where G = Scal( G ) * D, has small condition, then * ! SJGJF computes all eigenvalues with high relative accuracy. * ! Here D is diagonal positive definite, and the columns of * ! Scal( G ) have unit norms. More precisely, the computed * ! eigenvalues have approximately * ! LOG10( 1 / EPS ) - LOG10( Cond( B ) ) accurate decimal digits, * ! where EPS is the machine precision. For more details see [1]. * * * On entry, the parameter M is the number of columns of G. On exit, * M returns the rank of G, which is also the rank of H. If, on exit, * M is smaller than its starting value (deflation), then the initial * G was (almost) singular, and the results may be inaccurate. * * * SJGJF implements the implicit (one-sided) J-orthogonal Jacobi method * for solving Det( G' * G - Lambda * J ) = 0, originally defined in * [2]. The formal algorithm and its error analysis are given in [1]. * The global and the quadratic convergence have been proved in [2] and * [3], respectively. * * The implicit method works only on the matrix G, thus obtaining the * eigenvectors automatically. The diagonal of the matrix G' * G, which * is needed to determine the transformation parameters, is kept and * updated in a separate vector. * * SJGJF performs hyperbolic plane rotation if the pivot indices I < J * satisfy ( I <= NPOS ) .AND. ( J > NPOS ), and trigonometric * rotation otherwise. E.g. for N = M = 5 and NPOS = 3 the rotation * scheme looks like * * ( x T T H H ) * ( x T H H ) T = trigonometric rotation * ( x H H ) * ( x T ) H = hyperbolic rotation * ( x ) * * SJGJF uses either fast or fast self-scaling rotations. * * Treshold and convergence tests are both made after the relative * error analysis of [1]. Let * * ( A CC ) * ( CC B ) * * be the pivot submatrix. We rotate if * * SQRT( ABS (CC) ) * ------------------------ > N * EPS , * SQRT( A ) * SQRT ( B ) * * where EPS is the machine precision. The algorithm has converged * if M * ( M - 1 ) / 2 succesive rotations were not performed. * * * References * ========== * * [1] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * [2] K. Veselic: * An Eigenreduction Algorithm for Definite Matrix Pairs, preprint, * FB Mathematik, Fernuniversitaet Hagen, Hagen 1986, 1989, 1992. * * [3] Z. Drmac, V. Hari, On quadratic convergence bounds for the * J-symmetric Jacobi method, to appear in Numer. Math. * * * Arguments * ========= * * N (input) INTEGER * Number of rows of the matrix G. N >= 0. * * M (input/output) INTEGER * On entry, number of columns of the matrix G. N >= M >= 0. * On exit, M = rank( G ) = rank ( H ) = number of non-zero * eigenvalues of H. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * G (input/output) REAL array, dimension (LDA,M). * On entry a N x M real matrix. * On exit, the first M (M can change!) columns of G contain * the normalized eigenvectors corresponding to the non-zero * eigenvalues of the matrix H = G * J * G'. * * NPOS (input/output) INTEGER, M >= NPOS >= 0 * On entry, number of +1's in the matrix J. * On exit, number of +1's in the matrix J and number of the * positive eigenvalues of the matrix H. * (On entry and on exit, M and NPOS define implicitly the * matrix J.) * * EV (output) REAL array, dimension (M) * On exit, EV contains the non-zero eigenvalues of the matrix * H = G J G'. During the computation EV contains the diagonal * of the matrix G' * G. * * V (workspace) REAL array, dimension (M) * accumulates the diagonal of the fast rotations. * * JJ (workspace) INTEGER array, dimension (M) * JJ( I ) = 0 if the I-th column of G is non-zero, * JJ( I ) = 1 if the I-th column of G is zero. * The rotations which involve zero columns are skipped. * * JOB (input) CHARACTER * If JOB = 'f' or 'F', then SJGJ uses fast rotations. * If JOB = 's' or 'S', then SJGJ uses fast self-scaling * rotations. (These rotations should be used if * underflow/owerflow occurs in keeping the diagonal of the * fast rotations, which is unlikely to happen.) * * * Further Details * =============== * * SJGJF uses BLAS1 subroutines SSCAL and SCOPY, * BLAS1 function SDOT, * BLASJ subroutines SROTJJ, SROTF and SROTFS, and * LAPACK auxiliary function SLAMCH. * * BLASJ contains BLAS-type routines for use by J-orthogonal Jacobi * methods. SLAMCH determines the machine precision. * * ============================================================= * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. LOGICAL NOTSS INTEGER INFO, RANK, NP, I, J, ISWEEP INTEGER NTEMP, NMAX, SINGUL, IDI, IDJ, IDIDJ REAL EPS, ONEEPS, A, B, CC, HH, HYP, C, S, T * .. * .. External Subroutines .. EXTERNAL SROTJJ, SROTF, SROTFS, SCOPY, SSCAL * .. * .. External Functions .. REAL SDOT, SLAMCH EXTERNAL SDOT, SLAMCH * .. * .. Intrinsic Functions .. INTRINSIC MAX, SQRT, ABS, FLOAT, INT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( ( M .LT. 0 ) .OR. ( M .GT. N ) ) THEN INFO = 2 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 3 ELSE IF( ( NPOS. LT. 0 ) .OR. ( NPOS .GT. M ) ) THEN INFO = 5 ELSE IF( ( JOB .NE. 'f') .AND. ( JOB .NE. 'F' ) .AND. $ ( JOB .NE. 's' ) .AND. ( JOB .NE. 'S' ) ) THEN INFO = 9 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to SJGJF parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( ( N .EQ. 0 ) .OR. ( M .EQ. 0 ) ) RETURN IF( M .EQ. 1 ) THEN EV( 1 ) = SDOT( N, G( 1, 1 ), 1, G( 1, 1 ), 1) IF( EV( 1 ) .LE. ZERO ) THEN * * We have a zero matrix. * M = 0 ELSE * * Compute the only eigenvector and set the sign of the eigenvalue. * CALL SSCAL( N, ONE / SQRT( EV( 1 ) ), G( 1,1 ), 1 ) IF( NPOS .EQ. 0 ) EV( 1 ) = - EV( 1 ) ENDIF RETURN ENDIF * * Calculate the treshold N * EPS, and ONEEPS = ONE - SQRT( EPS ). * EPS = SLAMCH('E') ONEEPS = ONE - SQRT( EPS ) EPS = FLOAT( N ) * EPS * * Set NMAX to size of the sweep, and set NTEMP to 0. * NMAX = M * ( M - 1 ) / 2 NTEMP = 0 * * Set parameter NOTSS to .FALSE. if self-scaling rotations are to * be used. * NOTSS = .TRUE. IF( ( JOB .EQ. 's' ) .OR. ( JOB .EQ. 'S' ) ) NOTSS = .FALSE. * * Set vector V where we accumulate the diagonal of fast rotations * to 1, and set vector JJ to zero (at the beginning all columns of * G are treated as non-zero columns). * DO 9 I = 1, M V(I) = ONE JJ( I ) = 0 9 CONTINUE * * Set RANK to M, NP to NPOS, and SINGUL to 0. * SINGUL = 0 means that no singularity is detected in the current * rotation. * SINGUL = 1 means that the updated diagonal may be inaccurate. * SINGUL = 2 means that the singularity is detected. * RANK = M NP = NPOS SINGUL = 0 * * Start the iterations. * DO 40 ISWEEP = 1, 30 * * Refresh the diagonal of the matrix G' * G. If EV( I ) = 0, then * set JJ( I ) = 1 and decrease RANK by 1. If I-th column * corresponds to a positive eigenvalue, decrease NP by 1. * DO 30 I = 1, M IF( JJ( I ) .EQ. 1 ) GO TO 30 EV(I) = SDOT( N, G( 1,I ), 1, G( 1,I ), 1) * V( I ) ** 2 IF( EV( I ) . LE. ZERO ) THEN JJ( I ) = 1 RANK = RANK - 1 IF( I .LE. NPOS ) NP = NP - 1 END IF 30 CONTINUE * * Start the main loop. * DO 50 I = 1, M - 1 DO 50 J = I + 1, M * * If I-th or J-th columns are zero, skip this rotation and test * for convergence. * IF( ( JJ( I ) + JJ( J ) ) .GT. 0 ) GO TO 60 * * Choose trigonometric or hyperbolic rotation. * HYP = - ONE IF( ( I .LE. NPOS ) .AND. ( J .GT. NPOS ) ) HYP = ONE * * Calculate the pivot submatrix. * A = EV( I ) B = EV( J ) CC = SDOT( N, G( 1,I ), 1, G( 1,J ), 1 ) * V( I ) * V( J ) * * Test whether the matrix Scal( G )' * Scal( G ) is positive * definite. This means that HH must be smaller than 1. If * HH >= 1, we compute the diagonal elements A and B from the * columns of G once more and repeat the test. If the test fails * again, G is singular and the computed eigensolution might be * inaccurate. * HH = ABS( CC ) / SQRT( A ) / SQRT( B ) IF( HH .GE. ONE ) THEN A = SDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) * V(I) ** 2 B = SDOT( N, G( 1,J ), 1, G( 1,J ), 1 ) * V(J) ** 2 HH = ABS( CC ) / SQRT( A ) / SQRT( B ) IF( HH .GE. ONE ) THEN * * G' * G is singular. Set SINGUL = 2. * SINGUL = 2 GO TO 55 END IF END IF * * If the following test is true, the update of the diagonal may be * inaccurate. In this case set SINGUL = 1, which means that the * diagonal elements are later calculted from the columns of G. * IF( HH .GE. ONEEPS ) SINGUL = 1 * * Test whether to perform this rotation. If the rotation is * not performed, set NTEMP = NTEMP + 1. If NTEMP = NMAX, * the algorithm has converged. * IF( HH .GT. EPS ) GO TO 55 60 NTEMP = NTEMP + 1 IF( NTEMP .GT. NMAX ) GO TO 100 GO TO 50 * * Calculate the rotation parameters and set NTEMP = 0. * 55 CALL SROTJJ( A, B, CC, C, S, T, HYP ) NTEMP = 0 * * If in the hyperbolic case, SROTJJ could not calculate the * transformation parameters, it returns C = - 1. In this case * I-th and J-th columns of G both correspond to zero * eigenvalues. Set JJ( I ) and JJ( J ) to 1, decrease RANK * by 2, decrease NP by 1, and go to the next rotation. On exit * NPOS = NP is the number of the positive eigenvalues of H. * IF( C .LT. ZERO ) THEN JJ( I ) = 1 JJ( J ) = 1 RANK = RANK - 2 NP = NP - 1 SINGUL = 0 GO TO 50 END IF * * Update the columns I and J of G, and the elements of the * accumulated diagonal of fast rotations, V( I ) and V( J ). * SROTF performs fast J-orthogonal Jacobi rotation. * SROTFS performs fast self-scaling J-orthogonal Jacobi rotation. * * Self--scaling rotation, if used, pushes the element of V, * which is further away from 1, towards 1. For details see [1]. * IF ( NOTSS ) GO TO 44 IDI = 1 IDJ = 1 IDIDJ = 1 IF( V( I ) .LT. ONE ) IDI = - 1 IF( V( J ) .LT. ONE ) IDJ = - 1 IF( V( I ) .LE. V( J ) ) IDIDJ = - 1 IF( ( IDI + IDJ + 2 * INT( HYP ) ) * IDIDJ ) 43, 44, 42 * 42 CALL SROTFS ( N, G( 1,I ), 1, G( 1,J ), 1, T * V(I) / V(J), $ HYP * C * S * V( J ) / V( I ) ) V( I ) = V( I ) / C V( J ) = V( J ) * C GO TO 45 * 43 CALL SROTFS( N, G( 1,J ), 1, G( 1,I ), 1, $ HYP * T * V( J ) / V( I ), C * S * V( I ) / V( J ) ) V( I ) = V( I ) * C V( J ) = V( J ) / C GO TO 45 * 44 CALL SROTF( N, G( 1,I ), 1, G( 1,J ), 1, T * V( I ) / V( J ), $ HYP * T * V( J ) / V( I ) ) V( I ) = V( I ) * C V( J ) = V( J ) * C * * Both, hyperbolic and trigonometric rotations decrease the smaller * diagonal element. If the matrix is singular (SINGUL = 2), then * the column corresponding to the smaller diagonal element becomes * after rotation a zero-column. Set the corresponding element of JJ * to 1 and decrease RANK by 1. If this column corresponds to a * positive eigenvalue, decrease NP by 1. The non-zero diagonal * element is computed from the corresponding column of G. * 45 CONTINUE IF( SINGUL .EQ. 2 ) THEN IF( A .LE. B ) THEN JJ( I ) = 1 IF( I .LE. NPOS ) NP = NP - 1 EV( J ) = SDOT( N, G( 1,J ), 1, G( 1,J ), 1 ) * V(J) ** 2 ELSE JJ( J ) = 1 IF( J .LE. NPOS ) NP = NP - 1 EV( I ) = SDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) * V(I) ** 2 END IF RANK = RANK - 1 SINGUL = 0 GO TO 50 END IF * * Update the diagonal elements of G' * G, EV( I ) and EV( J ). * IF( SINGUL .EQ. 0 ) THEN EV( I ) = A + HYP * CC * T EV( J ) = B + CC * T ELSE * * This point is reached if SINGUL = 1. This means that the * updates of the diagonal elements might be inaccurate. The * diagonal elements are therefore computed from the I-th and * J-th column of G. Although no singularity was detected, one * of the diagonal elements may have become zero. In this case, * set the corresponding element of JJ to 1 and decrease RANK * by 1. If this column corresponds to a positive eigenvalue, * decrease NP by 1. * EV( I ) = SDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) * V(I) ** 2 IF( EV( I ) .LE. ZERO ) THEN JJ( I ) = 1 RANK = RANK - 1 IF( I .LE. NPOS ) NP = NP - 1 END IF EV( J ) = SDOT( N, G( 1,J ), 1, G( 1,J ), 1 ) * V(J) ** 2 IF( EV( J ) .LE. ZERO ) THEN JJ( J ) = 1 RANK = RANK - 1 IF( J .LE. NPOS ) NP = NP - 1 END IF SINGUL = 0 ENDIF * * End of the main loop. * 50 CONTINUE * * End of the ISWEEP loop. * 40 CONTINUE * * If this point is reached, the algorithm did not converge. * WRITE(*,*) ' SJGJF executed ', ISWEEP, ' sweeps without', $ ' convergence. ' 100 CONTINUE * * If no deflation occured, skip the next part. * IF( M .EQ. RANK ) GO TO 200 * * The deflation has occured. The results might be inaccurate! * Move NP non-zero columns of G to the first NP columns, * and do the same for the corresponding elements of the * vectors EV and V. * DO 110 I = 1, NP IF( JJ( I ) .EQ. 1 ) THEN J = I + 1 120 CONTINUE IF( JJ( J ) .EQ. 0 ) THEN CALL SCOPY( N, G( 1, J ), 1, G( 1, I ), 1 ) EV( I ) = EV( J ) V( I ) = V( J ) JJ( J ) = 1 GO TO 110 END IF J = J + 1 GO TO 120 END IF 110 CONTINUE * * Move the rest of the RANK - NP non-zero columns of G to * the columns NP + 1, ... , RANK, and do the same for the * corresponding elements of the vectors EV and V. * DO 140 I = NP + 1, RANK IF( JJ( I ) .EQ. 1 ) THEN J = I + 1 150 CONTINUE IF( JJ( J ) .EQ. 0 ) THEN CALL SCOPY( N, G( 1, J ), 1, G( 1, I ), 1 ) EV( I ) = EV( J ) V( I ) = V( J ) JJ( J ) = 1 GO TO 140 END IF J = J + 1 GO TO 150 END IF 140 CONTINUE * * Set M to rank of G, and NPOS to the number of positive * eigenvalues of H. * M = RANK NPOS = NP * * Calculate the eigenvalues and the corresponding eigenvectors. * Absolute values of the eigenvalues are the squares of the norms * of the columns of G. The corresponding eigenvectors are the * columns of G divided by their norms. * 200 CONTINUE DO 80 I = 1, M C = SDOT( N, G( 1,I ), 1, G( 1,I ), 1 ) CALL SSCAL( N, ONE / SQRT( C ), G( 1,I ), 1 ) EV( I ) = C * V( I ) ** 2 80 CONTINUE * * Change the signs for the negative eigenvalues. * DO 70 I = NPOS + 1, M EV( I ) = - EV( I ) 70 CONTINUE RETURN * * End of SJGJF * END CUT HERE.......... cat > djac.f <<'CUT HERE..........' SUBROUTINE DJAC( N, LDA, H, EVEC, EV, W ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 30, 1992 * * .. Scalar Arguments INTEGER N, LDA * .. * .. Array Arguments DOUBLE PRECISION H( LDA, * ), EVEC( LDA, * ), EV( * ), W( * ) * .. * * Purpose * ======= * * DJAC computes all eigenvalues and and eigenvectors of a real * symmetric N x N matrix H using the standard Jacobi method * with delayed updates of the diagonal. For more information see * * H. Rutishauser: * The Jacobi Method for Real Symmetric Matrices, * Numer. Math. 9, 1-10 (1966) * * Only upper triangle of H is referenced. * * Treshold is relative. Let * * ( A CC ) * ( CC B ) * * be the pivot submatrix. We perform the rotation if * * SQRT( ABS( CC ) ) > EPS * SQRT( ( ABS( A ) ) * SQRT( ABS( B ) ) * * where EPS is the machine precision. The algorithm has converged * if N * ( N - 1 ) / 2 succesive rotations were not performed. * * * Arguments * ========= * * N (input) INTEGER * Dimension of the symmetric matrix H. N >= 0. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * H (input) DOUBLE PRECISION array, dimension (LDA,N). * N x N real symmetric matrix. * * EVEC (output) DOUBLE PRECISION array, dimension (LDA,N). * Contains the eigenvectors. * * EV (output) DOUBLE PRECISION array, dimansion (N). * contains eigenvalues of H. * * W (workspace) DOUBLE PRECISION array, dimansion (N) * * * Further Details * =============== * * DJAC uses BLAS1 subroutines DROT, DCOPY and DAXPY, * BLASJ subroutine DROTJJ, and * LAPACK auxiliary function DLAMCH. * * BLASJ contains BLAS-type routines for use by J-orthogonal Jacobi * methods. DLAMCH determines the machine precision. * * ============================================================= * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, J, ISWEEP, NTEMP, NMAX DOUBLE PRECISION EPS, A, B, CC, C, S, T * .. * .. External Subroutines .. EXTERNAL DROTJJ, DROT, DCOPY, DAXPY * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC DSQRT, DABS * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 2 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to DJAC parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( N .EQ. 0 ) RETURN IF( N .EQ. 1 ) THEN EVEC( 1, 1 ) = ONE EV( 1 ) = H( 1, 1 ) RETURN END IF * * Calculate the machine precision EPS, set NMAX to size of the * sweep, and set NTEMP to 0. * EPS = DLAMCH('E') NMAX = N * ( N - 1 ) / 2 NTEMP = 0 * * Set the array EVEC to an identity matrix. * DO 30 I = 1, N EVEC( I, I ) = ONE DO 30 J = I + 1, N EVEC( I, J ) = ZERO EVEC( J, I ) = ZERO 30 CONTINUE * * Set EV to diagonal of H. * CALL DCOPY( N, H( 1,1 ), LDA + 1, EV, 1 ) * * Start the iterations. * DO 40 ISWEEP = 1, 30 * * Set W to zero. * DO 45 I = 1, N W( I ) = ZERO 45 CONTINUE * * Begin the main loop. * DO 50 I=1,N-1 DO 50 J=I+1,N * * Calculate the pivot submatrix. * A = H( I, I ) B = H( J, J ) CC = H( I, J ) * * Test whether to perform this rotation. If the rotation is * performed, set NTEMP to 0, otherwise NTEMP = NTEMP + 1. * IF NTEMP = NMAX, the algorithm has converged. * IF( DABS( CC ) .LE. ( EPS * DSQRT( DABS( A ) ) * $ DSQRT( DABS( B ) ) ) ) THEN NTEMP = NTEMP + 1 IF( NTEMP .GT. NMAX ) GO TO 100 GO TO 50 END IF NTEMP = 0 * * Calculate the Jacobi rotation parameters. * CALL DROTJJ( A, B, CC, C, S, T, - ONE ) * * Update the upper triangle of H except the elements of the * pivot submatrix. * CALL DROT( I-1, H( 1,I ), 1, H( 1,J ), 1, C, - S ) CALL DROT( J-I-1, H( I, I + 1 ), LDA, H( I + 1, J ), 1, C, - S ) CALL DROT( N-J, H( I, J + 1 ), LDA, H( J, J + 1 ), LDA, C, - S ) * * Update the elements of the pivot submatrix, and add the updates * of the diagonal to elements W( I ) and W( J ). * CC = CC * T H( I, I ) = A - CC H( J, J ) = B + CC H( I, J ) = ZERO W( I ) = W( I ) - CC W( J ) = W( J ) + CC * * Update the eigenvalue matrix. * CALL DROT( N, EVEC( 1,I ), 1, EVEC( 1,J ), 1, C, - S ) * * End of the main loop. * 50 CONTINUE * * Add the updates of the diagonal stored in W to the diagonal * which was stored at the end of the previous sweep EV. * CALL DAXPY( N, ONE, W, 1, EV, 1 ) * * Set diagonal of H to EV, and set W to zero. * CALL DCOPY( N, EV, 1, H( 1,1 ), LDA + 1) DO 60 I = 1, N W( I ) = ZERO 60 CONTINUE * * End of the ISWEEP loop. * 40 CONTINUE * * If this point is reached, the algorithm did not converge. * WRITE(*,*) ' DJAC executed ', ISWEEP, ' sweeps without', $ ' convergence. ' RETURN * * Add the updates stored in W to EV. * 100 CONTINUE CALL DAXPY( N, ONE, W, 1, EV, 1 ) RETURN * * End of DJAC. * END CUT HERE.......... cat > sjac.f <<'CUT HERE..........' SUBROUTINE SJAC( N, LDA, H, EVEC, EV, W ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 30, 1992 * * .. Scalar Arguments INTEGER N, LDA * .. * .. Array Arguments REAL H( LDA, * ), EVEC( LDA, * ), EV( * ), W( * ) * .. * * Purpose * ======= * * SJAC computes all eigenvalues and and eigenvectors of a real * symmetric N x N matrix H using the standard Jacobi method * with delayed updates of the diagonal. For more information see * * H. Rutishauser: * The Jacobi Method for Real Symmetric Matrices, * Numer. Math. 9, 1-10 (1966) * * Only the upper triangle of H is referenced. * * Treshold is relative. Let * * ( A CC ) * ( CC B ) * * be the pivot submatrix. We perform the rotation if * * SQRT( ABS( CC ) ) > EPS * SQRT( ABS( A ) ) * SQRT( ABS( B ) ) * * where EPS is the machine precision. The algorithm has converged * if N * ( N - 1 ) / 2 succesive rotations were not performed. * * * Arguments * ========= * * N (input) INTEGER * Dimension of the symmetric matrix H. N >= 0. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * H (input) REAL array, dimension (LDA,N). * N x N real symmetric matrix. * * EVEC (output) REAL array, dimension (LDA,N). * Contains the eigenvectors. * * EV (output) REAL array, dimansion (N). * contains eigenvalues of H. * * W (workspace) REAL array, dimansion (N) * * * Further Details * =============== * * SJAC uses BLAS1 subroutines SROT, SCOPY and SAXPY, * BLASJ subroutine SROTJJ, and * LAPACK auxiliary function SLAMCH. * * BLASJ contains BLAS-type routines for use by J-orthogonal Jacobi * methods. SLAMCH determines the machine precision. * * ============================================================= * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, J, ISWEEP, NTEMP, NMAX REAL EPS, A, B, CC, C, S, T * .. * .. External Subroutines .. EXTERNAL SROTJJ, SROT, SCOPY, SAXPY * .. * .. External Functions .. REAL SLAMCH EXTERNAL SLAMCH * .. * .. Intrinsic Functions .. INTRINSIC SQRT, ABS * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 2 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to SJAC parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( N .EQ. 0 ) RETURN IF( N .EQ. 1 ) THEN EVEC( 1, 1 ) = ONE EV( 1 ) = H( 1, 1 ) RETURN END IF * * Calculate the machine precision EPS, set NMAX to size of the * sweep, and set NTEMP to 0. * EPS = SLAMCH('E') NMAX = N * ( N - 1 ) / 2 NTEMP = 0 * * Set the array EVEC to an identity matrix. * DO 30 I = 1, N EVEC( I, I ) = ONE DO 30 J = I + 1, N EVEC( I, J ) = ZERO EVEC( J, I ) = ZERO 30 CONTINUE * * Set EV to diagonal of H. * CALL SCOPY( N, H( 1,1 ), LDA + 1, EV, 1 ) * * Start the iterations. * DO 40 ISWEEP = 1, 30 * * Set W to zero. * DO 45 I = 1, N W( I ) = ZERO 45 CONTINUE * * Begin the main loop. * DO 50 I=1,N-1 DO 50 J=I+1,N * * Calculate the pivot submatrix. * A = H( I, I ) B = H( J, J ) CC = H( I, J ) * * Test whether to perform this rotation. If the rotation is * performed, set NTEMP to 0, otherwise NTEMP = NTEMP + 1. * IF NTEMP = NMAX, the algorithm has converged. * IF( ABS( CC ) .LE. ( EPS * SQRT( ABS( A ) ) * $ SQRT( ABS( B ) ) ) ) THEN NTEMP = NTEMP + 1 IF( NTEMP .GT. NMAX ) GO TO 100 GO TO 50 END IF NTEMP = 0 * * Calculate the Jacobi rotation parameters. * CALL SROTJJ( A, B, CC, C, S, T, - ONE ) * * Update the upper triangle of H except the elements of the * pivot submatrix. * CALL SROT( I-1, H( 1,I ), 1, H( 1,J ), 1, C, - S ) CALL SROT( J-I-1, H( I, I + 1 ), LDA, H( I + 1, J ), 1, C, - S ) CALL SROT( N-J, H( I, J + 1 ), LDA, H( J, J + 1 ), LDA, C, - S ) * * Update the elements of the pivot submatrix, and add the updates * of the diagonal to elements W( I ) and W( J ). * CC = CC * T H( I, I ) = A - CC H( J, J ) = B + CC H( I, J ) = ZERO W( I ) = W( I ) - CC W( J ) = W( J ) + CC * * Update the eigenvalue matrix. * CALL SROT( N, EVEC( 1,I ), 1, EVEC( 1,J ), 1, C, - S ) * * End of the main loop. * 50 CONTINUE * * Add the updates of the diagonal stored in W to the diagonal * which was stored at the end of the previous sweep EV. * CALL SAXPY( N, ONE, W, 1, EV, 1 ) * * Set diagonal of H to EV, and set W to zero. * CALL SCOPY( N, EV, 1, H( 1,1 ), LDA + 1) DO 60 I = 1, N W( I ) = ZERO 60 CONTINUE * * End of the ISWEEP loop. * 40 CONTINUE * * If this point is reached, the algorithm did not converge. * WRITE(*,*) ' SJAC executed ', ISWEEP, ' sweeps without', $ ' convergence. ' RETURN * * Add the updates stored in W to EV. * 100 CONTINUE CALL SAXPY( N, ONE, W, 1, EV, 1 ) RETURN * * End of SJAC. * END CUT HERE.......... cat > blasj.f <<'CUT HERE..........' * * This file contains following BLAS-type subroutines * for use with the J-orthogonal Jacobi methods: * * DROTJJ, SROTJJ - construct a J-orthogonal Jacobi rotation * DROTF, SROTF - apply a fast J-orthogonal Jacobi rotation * DROTFS, SROTFS - apply a fast self-scaling J-orthogonal * Jacobi rotation * DIPERR, SIPERR - permute rows of a matrix * SUBROUTINE DROTJJ( A, B, CC, C, S, T, HYP) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments DOUBLE PRECISION A, B, CC, C, S, T, HYP * * Purpose * ======= * * DROTJJ constructs trigonometric or hyperbolic Jacobi plane rotation * of the form * * ( C S ) * R = ( ) * ( HYP * S C ) * * such that the matrix * * ( A CC ) * R' * ( ) * R * ( CC B ) * * is diagonal. If HYP = - 1, the rotation is trigonometric, that is, * * C ** 2 + S ** 2 = 1 . * * If HYP = 1, the rotation is hyperbolic, that is * * C ** 2 - S ** 2 = 1 . * * Arguments * ========= * * A, B, CC (input) DOUBLE PRECISION * define the 2 x 2 input matrix. * * C, S (output) DOUBLE PRECISION * C and S define the plane rotation. * If C = - 1, DROTJJ failed to construct the rotation. * This can only happen in the hyperbolic case. Maximal * C that can be calculated is of order of the fourth root * of the machine precision. * * T (output) DOUBLE PRECISION * T = S / C * * HYP (input) DOUBLE PRECISION * If HYP = 1.0D0, DROTJJ computes a hyperbolic rotation. * If HYP = -1.0D0, DROTJJ computes a trigonometric rotation. * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. DOUBLE PRECISION ABSCC, TMP, TMPBA * .. * .. Intrinsic Functions .. INTRINSIC DSQRT, DABS * .. * .. Executable Statements .. * ABSCC = DABS( CC ) TMP = 100.0D0 * ABSCC TMPBA = B + HYP * A * * If TMP is small compared to TMPBA, use the approximate formulae. * IF( ( TMPBA + TMP ) .EQ. TMPBA ) THEN T = - HYP * CC / TMPBA C = ONE S = T RETURN END IF * * Calculate the tangent. If the rotation is hyperbolic, we have * to check whether it can be calculated. * TMP = - HYP * 0.5D0 * TMPBA / CC TMPBA = TMP * TMP - HYP IF( TMPBA .GT. ZERO ) THEN * * Standard formula for tangent. * T = ONE / ( DABS( TMP ) + DSQRT( TMPBA ) ) ELSE * * Hyperbolic rotation cannot be calculated. Set C to - 1 * and return. * C = - ONE RETURN END IF * * Set the signum of the tangent. * IF( TMP .LT. ZERO ) T = - T * * Calculate cosine and sine. * TMP = DSQRT( ONE - HYP * T * T ) C = ONE / TMP S = T / TMP RETURN * * End of DROTJJ. * END * * * SUBROUTINE SROTJJ( A, B, CC, C, S, T, HYP) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments REAL A, B, CC, C, S, T, HYP * * Purpose * ======= * * SROTJJ constructs trigonometric or hyperbolic Jacobi plane rotation * of the form * * ( C S ) * R = ( ) * ( HYP * S C ) * * such that the matrix * * ( A CC ) * R' * ( ) * R * ( CC B ) * * is diagonal. If HYP = - 1, the rotation is trigonometric, that is, * * C ** 2 + S ** 2 = 1 . * * If HYP = 1, the rotation is hyperbolic, that is * * C ** 2 - S ** 2 = 1 . * * * Arguments * ========= * * A, B, CC (input) REAL * define the 2 x 2 input matrix. * * C, S (output) REAL * C and S define the plane rotation. * If C = - 1, SROTJJ failed to construct the rotation. * This can only happen in the hyperbolic case. Maximal * C that can be calculated is of order of the fourth root * of the machine precision. * * T (output) REAL * T = S / C * * HYP (input) REAL * If HYP = 1.0E0, SROTJJ computes a hyperbolic rotation. * If HYP = -1.0E0, SROTJJ computes a trigonometric rotation. * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. REAL ABSCC, TMP, TMPBA * .. * .. Intrinsic Functions .. INTRINSIC SQRT, ABS * .. * .. Executable Statements .. * ABSCC = ABS( CC ) TMP = 100.0E0 * ABSCC TMPBA = B + HYP * A * * If TMP is small compared to TMPBA, use the approximate formulae. * IF( ( TMPBA + TMP ) .EQ. TMPBA ) THEN T = - HYP * CC / TMPBA C = ONE S = T RETURN END IF * * Calculate the tangent. If the rotation is hyperbolic, we have * to check whether it can be calculated. * TMP = - HYP * 0.5E0 * TMPBA / CC TMPBA = TMP * TMP - HYP IF( TMPBA .GT. ZERO ) THEN * * Standard formula for tangent. * T = ONE / ( ABS( TMP ) + SQRT( TMPBA ) ) ELSE * * Hyperbolic rotation cannot be calculated. Set C to - 1 * and return. * C = - ONE RETURN END IF * * Set the signum of the tangent. * IF( TMP .LT. ZERO ) T = - T * * Calculate cosine and sine. * TMP = SQRT( ONE - HYP * T * T ) C = ONE / TMP S = T / TMP RETURN * * End of SROTJJ. * END * * * SUBROUTINE DROTF( N, X, INX, Y, INY, ALPHA, BETA) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, INX, INY DOUBLE PRECISION ALPHA, BETA * .. * .. Array Arguments DOUBLE PRECISION X( * ), Y( * ) * * Purpose * ======= * * DROTF applies a fast plane rotation of the form * * ( 1 ALPHA ) * ( BETA 1 ) * * Arguments * ========= * * N (input) INTEGER * lenght of the vectors X and Y which are rotated * * X, Y (input) DOUBLE PRECISION array * * INX, INY (input) INTEGER * Increments of the arrays X and Y, respectively. * * ALPHA (input) DOUBLE PRECISION * * BETA (input) DOUBLE PRECISION * * .. * .. Local Scalars .. INTEGER NN, IX, IY, INCX, INCY, I DOUBLE PRECISION TEMP * .. * .. Executable Statements .. * IF( N .LE. 0 ) RETURN NN = N INCX = INX INCY = INY IF( INCX .EQ. 1 .AND. INCY .EQ. 1 ) GO TO 20 * * Code for unequal increments or equal increments not equal to 1. * IX = 1 IY = 1 IF( INCX .LT. 0 ) IX = ( - NN + 1 ) * INCX + 1 IF( INCY .LT. 0 ) IY = ( - NN + 1 ) * INCY + 1 DO 10 I = 1, NN TEMP = X( IX ) + BETA * Y( IY ) Y( IY ) = Y( IY ) + ALPHA * X( IX ) X( IX ) = TEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN * * Code for both increments equal to 1. * 20 CONTINUE DO 30 I = 1, NN TEMP = X( I ) + BETA * Y( I ) Y( I ) = Y( I ) + ALPHA * X( I ) X( I ) = TEMP 30 CONTINUE RETURN * * End of DROTF. * END * * * SUBROUTINE SROTF( N, X, INX, Y, INY, ALPHA, BETA) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, INX, INY REAL ALPHA, BETA * .. * .. Array Arguments REAL X( * ), Y( * ) * * Purpose * ======= * * SROTF applies a fast plane rotation of the form * * ( 1 ALPHA ) * ( BETA 1 ) * * Arguments * ========= * * N (input) INTEGER * lenght of the vectors X and Y which are rotated * * X, Y (input) REAL array * * INX, INY (input) INTEGER * Increments of the arrays X and Y, respectively. * * ALPHA (input) REAL * * BETA (input) REAL * * .. * .. Local Scalars .. INTEGER NN, IX, IY, INCX, INCY, I REAL TEMP * .. * .. Executable Statements .. * IF( N .LE. 0 ) RETURN NN = N INCX = INX INCY = INY IF( INCX .EQ. 1 .AND. INCY .EQ. 1 ) GO TO 20 * * Code for unequal increments or equal increments not equal to 1. * IX = 1 IY = 1 IF( INCX .LT. 0 ) IX = ( - NN + 1 ) * INCX + 1 IF( INCY .LT. 0 ) IY = ( - NN + 1 ) * INCY + 1 DO 10 I = 1, NN TEMP = X( IX ) + BETA * Y( IY ) Y( IY ) = Y( IY ) + ALPHA * X( IX ) X( IX ) = TEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN * * Code for both increments equal to 1. * 20 CONTINUE DO 30 I = 1, NN TEMP = X( I ) + BETA * Y( I ) Y( I ) = Y( I ) + ALPHA * X( I ) X( I ) = TEMP 30 CONTINUE RETURN * * End of SROTF. * END * * * SUBROUTINE DROTFS( N, X, INX, Y, INY, ALPHA, BETA) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, INX, INY DOUBLE PRECISION ALPHA, BETA * .. * .. Array Arguments DOUBLE PRECISION X( * ), Y( * ) * * Purpose * ======= * * DROTF applies a fast self-scaling plane rotation of the form * * ( 1 ALPHA ) * ( 1 0 ) * ( 0 1 ) ( BETA 1 ) * * Interchenging the vectors X and Y , "changes" the * transformation matrix to * * ( 1 0 ) * ( 1 BETA ) * ( ALPHA 1 ) ( 0 1 ) * * * Arguments * ========= * * N (input) INTEGER * lenght of the vectors X and Y which are rotated * * X, Y (input) DOUBLE PRECISION array * * INX, INY (input) INTEGER * Increments of the arrays X and Y, respectively. * * ALPHA (input) DOUBLE PRECISION * * BETA (input) DOUBLE PRECISION * * .. * .. Local Scalars .. INTEGER NN, IX, IY, INCX, INCY, I * .. * .. Executable Statements .. * IF( N .LE. 0 ) RETURN NN = N INCX = INX INCY = INY IF( INCX .EQ. 1 .AND. INCY .EQ. 1 ) GO TO 20 * * Code for unequal increments or equal increments not equal to 1. * IX = 1 IY = 1 IF( INCX .LT. 0 ) IX = ( - NN + 1 ) * INCX + 1 IF( INCY .LT. 0 ) IY = ( - NN + 1 ) * INCY + 1 DO 10 I = 1, NN Y( IY ) = Y( IY ) + ALPHA * X( IX ) X( IX ) = X( IX ) + BETA * Y( IY ) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN * * Code for both increments equal to 1. * 20 CONTINUE DO 30 I = 1, NN Y( I ) = Y( I ) + ALPHA * X( I ) X( I ) = X( I ) + BETA * Y( I ) 30 CONTINUE RETURN * * End of DROTFS. * END * * * SUBROUTINE SROTFS( N, X, INX, Y, INY, ALPHA, BETA) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, INX, INY REAL ALPHA, BETA * .. * .. Array Arguments REAL X( * ), Y( * ) * * Purpose * ======= * * SROTF applies a fast self-scaling plane rotation of the form * * ( 1 ALPHA ) * ( 1 0 ) * ( 0 1 ) ( BETA 1 ) * * Interchenging the vectors X and Y , "changes" the * transformation matrix to * * ( 1 0 ) * ( 1 BETA ) * ( ALPHA 1 ) ( 0 1 ) * * * Arguments * ========= * * N (input) INTEGER * lenght of the vectors X and Y which are rotated * * X, Y (input) REAL array * * INX, INY (input) INTEGER * Increments of the arrays X and Y, respectively. * * ALPHA (input) REAL * * BETA (input) REAL * * .. * .. Local Scalars .. INTEGER NN, IX, IY, INCX, INCY, I * .. * .. Executable Statements .. * IF( N .LE. 0 ) RETURN NN = N INCX = INX INCY = INY IF( INCX .EQ. 1 .AND. INCY .EQ. 1 ) GO TO 20 * * Code for unequal increments or equal increments not equal to 1. * IX = 1 IY = 1 IF( INCX .LT. 0 ) IX = ( - NN + 1 ) * INCX + 1 IF( INCY .LT. 0 ) IY = ( - NN + 1 ) * INCY + 1 DO 10 I = 1, NN Y( IY ) = Y( IY ) + ALPHA * X( IX ) X( IX ) = X( IX ) + BETA * Y( IY ) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN * * Code for both increments equal to 1. * 20 CONTINUE DO 30 I = 1, NN Y( I ) = Y( I ) + ALPHA * X( I ) X( I ) = X( I ) + BETA * Y( I ) 30 CONTINUE RETURN * * End of SROTFS. * END * * * SUBROUTINE DIPERR( N, M, LDA, A, PIVOT ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 29, 1992 * * .. Scalar Arguments INTEGER N, M, LDA * .. * .. Array Arguments DOUBLE PRECISION A( LDA, * ) INTEGER PIVOT( * ) * .. * * Purpose * ======= * * DIPERR applies inverse of the permutation given by the vector * PIVOT to rows of the N x M matrix A. * * * Arguments * ========= * * N (input) INTEGER * Number of rows of the matrix A. N >= 0. * * M (input) INTEGER * Number of columns of the matrix A. M >= 0. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N) * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N). * On entry, real N x M matrix. * On exit, rows of A are permuted. * * PIVOT (input) INTEGER array, dimension (N) * defines the permutation whise inverse is applied to rows * of A. * * Further Details * =============== * * DIPERR uses BLAS1 subroutine DSWAP. * * * .. Local Scalars .. INTEGER NN, MM, I, J, TEMP * .. * .. External Subroutines .. EXTERNAL DSWAP * .. * .. Executable Statements .. * * Quick return, if possible. * IF( ( N .LE. 0 ) .OR. ( M .LE. 0 ) ) RETURN * NN = N MM = M DO 10 I = 1, NN DO 10 J = I + 1, NN IF( PIVOT( J ) .EQ. I ) THEN CALL DSWAP( MM, A( I,1 ), LDA, A( J,1 ), LDA ) TEMP = PIVOT( I ) PIVOT( I ) = PIVOT( J ) PIVOT( J ) = TEMP END IF 10 CONTINUE * * End of DIPERR. * RETURN END * * * SUBROUTINE SIPERR( N, M, LDA, A, PIVOT ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 29, 1992 * * .. Scalar Arguments INTEGER N, M, LDA * .. * .. Array Arguments REAL A( LDA, * ) INTEGER PIVOT( * ) * .. * * Purpose * ======= * * SIPERR applies inverse of the permutation given by the vector * PIVOT to rows of the N x M matrix A. * * * Arguments * ========= * * N (input) INTEGER * Number of rows of the matrix A. N >= 0. * * M (input) INTEGER * Number of columns of the matrix A. M >= 0. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N) * * A (input/output) REAL array, dimension (LDA,N). * On entry, real N x M matrix. * On exit, rows of A are permuted. * * PIVOT (input) INTEGER array, dimension (N) * defines the permutation whise inverse is applied to rows * of A. * * Further Details * =============== * * SIPERR uses BLAS1 subroutine SSWAP. * * * .. Local Scalars .. INTEGER NN, MM, I, J, TEMP * .. * .. External Subroutines .. EXTERNAL SSWAP * .. * .. Executable Statements .. * * Quick return, if possible. * IF( ( N .LE. 0 ) .OR. ( M .LE. 0 ) ) RETURN * NN = N MM = M DO 10 I = 1, NN DO 10 J = I + 1, NN IF( PIVOT( J ) .EQ. I ) THEN CALL SSWAP( MM, A( I,1 ), LDA, A( J,1 ), LDA ) TEMP = PIVOT( I ) PIVOT( I ) = PIVOT( J ) PIVOT( J ) = TEMP END IF 10 CONTINUE * * End of SIPERR. * RETURN END CUT HERE.......... cat > gensym.f <<'CUT HERE..........' SUBROUTINE GENSYM( N, LDA, SEED, ASCAL, HSCAL, H, A, EV, D, P ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * Juli 3, 1992 * * .. Scalar Arguments INTEGER N, LDA, SEED DOUBLE PRECISION ASCAL, HSCAL * .. * .. Array Arguments DOUBLE PRECISION H( LDA, * ), A( LDA, * ), EV( * ), D( * ) INTEGER P( * ) * .. * * Purpose * ======= * * GENSYM generates a non-singular symmetric matrix H of order N, * with spectral condition approximately equal to 10 ** HSCAL and * the measure C( A, A ) from [1] approximately equal to * 10 ** ASCAL. * * GENSYM is used with program COMPAR which compares several double * and single precision routines. * * * We first generate a random diagonal matrix A, whose entries' * logarithms are uniformly distributed between * [ - ASCAL / 2, ASCAL / 2 ]. We then apply 5 sweeps of random * trigonometric rotations on A (the spectral condition remains * unchanged). We then apply 5 sweeps of "anti-Jacobi" method on * A (the spectral condition remain unchanged again). This step * ensures that an arbitrary diagonal scaling, D * A * D, cannot * improve the spectral condition of A considerably. The * "anti_Jacobi" method is due to Veselic. For details see [1]. * We then scale A so that all diagonal elements become 1. This * step changes the spectral condition, but, as already mentioned, * not by much. Thus, the condition of A is now approximately the * condition of the starting diagonal matrix. We further generate * diagonal matrix D whose entries' logarithms are uniformly * distributed between [ - HSCAL / 4, HSCAL / 4 ], and form the * matrix A = D * A * D. The spectral condition of A is now * approximately 10 ** HSCAL. We then use DSYEVJ to compute the * eigendecomposition A = U' * EV * U. We then change signs of * some randomly chosen eigenvalues, thus obtaining EV'. * The final random symmetric indefinite matrix is then * H = U' * EV' * U. * * More details about this generator are found in [1]. * * References * ========== * * [1] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix H. N >= 0. * * LDA (input) INTEGER * The leading dimension of the array H. LDA >= max(1,N) * * SEED (input) INTEGER, SEED < 0. * Seed for the random number generator RAN1. * * ASCAL (input) DOUBLE PRECISION * LOG10 of the measure C( A, A ) (approximately). * * HSCAL (input) DOUBLE PRECISION * LOG10 of the spectral condition of H (approximately). * * H (output) DOUBLE PRECISION array, dimension (LDA,N). * Randomly genenerated real symmetric matrix of order N. * * A (workspace) DOUBLE PRECISION array, dimension (LDA,N). * * NPOS (output) INTEGER * number of positive eigenvalues of H. * (R and NPOS define implicitly the matrix J.) * * EV (workspace) DOUBLE PRECISION array, dimension (N) * * D (workspace) DOUBLE PRECISION array, dimension (N) * * P (workspace) INTEGER array, dimension (N) * * Further Details * =============== * * GENSYM uses BLAS1 subroutines DCOPY and DROT, * BLAS1 function DDOT, * subroutine DSYEVJ, and * subroutine DSCALM and function RAN1 from * the file auxil.f. * * ============================================================= * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, HALF PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER INFO, I, J, RANK DOUBLE PRECISION D10, TEMP, C, S, T, DM5 * .. * .. External Subroutines .. EXTERNAL DROT, DCOPY, DSCALM, DSYEVJ * .. * .. External Functions .. DOUBLE PRECISION RAN1, DDOT EXTERNAL RAN1, DDOT * .. * .. Intrinsic Functions .. INTRINSIC DSQRT, DABS * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N .LT. 0 ) THEN INFO = 1 ELSE IF( LDA .LT. MAX( 1,N ) ) THEN INFO = 2 ELSE IF( SEED .GE. 0 ) THEN INFO = 3 END IF IF( INFO .NE. 0 ) THEN WRITE(*,*) ' On entry to GENSYM parameter number ', INFO WRITE(*,*) ' had an illegal value. ' STOP END IF * * Quick return, if possible. * IF( N .EQ. 0 ) RETURN * * Set the array A to zero. * DO 10 I = 1, N DO 10 J = 1, N A( I, J ) = ZERO 10 CONTINUE * * Initialize the random number generator and define the parameter * D10. * T = RAN1( SEED ) D10 = DLOG( 10.0D0 ) * * Generate the diagonal of A. * DO 20 I = 1, N A( I,I ) = DEXP( D10 * ASCAL * ( RAN1( 0 ) - HALF ) ) 20 CONTINUE * * Apply 5 random sweeps of trigonometric rotations on A. * DO 30 ISWEEP = 1, 5 DO 30 I = 1, N - 1 DO 30 J = I + 1, N * * Compute a random tangent. * T = TWO * RAN1( 0 ) - ONE * * Compute sine and cosine. * TEMP = DSQRT( ONE + T * T ) C = ONE / TEMP S = T / TEMP * * Perform the rotation on A. * CALL DROT( N, A( 1,I ), 1, A( 1,J ), 1, C, - S ) A( I,I ) = C * A(I,I)- S * A( J,I ) TEMP = S * A( I,J ) + C * A( J,J ) A( I,J ) = C * A( I,J ) - S * A( J,J ) A( J,J ) = TEMP DO 40 K = 1, N A( I,K ) = A( K,I ) A( J,K ) = A( K,J ) 40 CONTINUE 30 CONTINUE * * Apply 5 sweeps of "anti-Jacobi" to A. The rotation is performed * only if ABS( A( I,I ) - A( J,J ) ) > 1.0D-5 * ABS( A( I,J ) ). * DM5 = 1.0D-5 DO 50 ISWEEP = 1, 5 DO 60 I = 1, N - 1 DO 60 J = I + 1, N * * Compute tangent. * TEMP = A( J,J ) - A( I,I ) IF( DABS( TEMP ) .LE. ( 1E-5 * DABS( A( I,J ) ) ) ) GOTO 60 TEMP = 2 * A( I,J ) / TEMP T = ONE /( DABS( TEMP ) + DSQRT( TEMP * TEMP + ONE ) ) IF ( TEMP .LT. ZERO) T = -T * * Compute sine and cosine. * TEMP = DSQRT( 1 + T * T ) C = ONE / TEMP S = - T/ TEMP * * Perform the rotation on A. * CALL DROT( N, A( 1,I ) , 1, A( 1,J ), 1, C, - S ) A( I,I ) = C * A( I,I ) - S * A( J,I ) TEMP = S * A( I,J ) + C * A( J,J ) A( I,J ) = C * A( I,J ) - S * A( J,J ) A( J,J ) = TEMP DO 65 K = 1, N A( I,K ) = A( K,I ) A( J,K ) = A( K,J ) 65 CONTINUE 60 CONTINUE 50 CONTINUE * * Scale the matrix A so that the diagonal elements become 1. * CALL DCOPY( N, A( 1,1 ), LDA + 1, D, 1 ) DO 70 I = 1, N D(I) = ONE / DSQRT( D( I ) ) 70 CONTINUE CALL DSCALM( N, LDA, A, D ) * * Generate a random diagonal matrix D. * DO 80 I = 1, N D( I ) = DEXP( HALF * D10 * HSCAL * ( RAN1( 0 ) - HALF ) ) 80 CONTINUE * * Scale A from both sides by with D. * CALL DSCALM( N, LDA, A, D ) * * Compute the eigendecomposition of A. Eigenvalues are stored in * EV, corresponding eigenvectors are stored in A. * CALL DSYEVJ( N, RANK, LDA, A, EV, D, P, 'F' ) * * Change the signs of some randomly selected eigenvalues. * DO 90 I = 1, RANK IF( ( RAN1( 0 ) - HALF ) .LT. ZERO ) EV( I ) = - EV( I ) 90 CONTINUE * * Compute the final matrix H = A' * EV * A. * DO 100 I = 1, N DO 110 K = 1, N D( K ) = A( I,K ) * EV( K ) 110 CONTINUE DO 120 J = I, N H( I,J ) = DDOT( N, D, 1, A( J,1 ), LDA ) H( J,I ) = H( I,J ) 120 CONTINUE 100 CONTINUE RETURN * * End of GENSYM. * END CUT HERE.......... cat > auxil.f <<'CUT HERE..........' * * This file contains following auxiliary routines: * * DDISPM, SDISPM - display a N x M matrix. * DSTORM, SSTORM - store a N x M matrix in a disk file. * DREADM, SREADM - read a N x M matrix from a disk file. * DSCALM - multiplies a symmetric matrix from both sides * with a diagonal matrix (scaling). * DOUT, SOUT - write sorted eigenvalues and corresponding * eigenvectors to a disc file. * COMP - compares eigensolutions computed by several * routines. * RAN1 - double precision random number generator. * * SUBROUTINE DDISPM( N, M, LDA, A, DEV, FNAME) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, DEV CHARACTER*(*) FNAME * .. * .. Array Arguments DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DDISPM displays a N x M DOUBLE PRECISION matrix A in a file * called FNAME on the device number DEV. It must be M <= 6 . * * LDA is the leading dimension of the array A, LDA >= max( 1, N ) * If DEV is the standard output (here it is set to 6), FNAME can * be blank. * * .. * .. Local Scalars .. INTEGER I, J, STOUT * .. * .. Executable Statements .. * IF( M .GT. 6 ) GO TO 100 STOUT = 6 IF ( DEV .NE. STOUT ) THEN OPEN( UNIT = DEV, FILE = FNAME ) CLOSE( IDEV, STATUS = 'DELETE' ) OPEN( UNIT = IDEV, FILE = FNAME ) END IF DO 10 I = 1, N WRITE( DEV, 20 ) ( A( I, J ), J = 1, M ) 20 FORMAT( 6 ( ' ', D11.5 ) ) 10 CONTINUE WRITE( DEV, * ) ' ' 100 CONTINUE IF ( DEV .NE. STOUT ) CLOSE( DEV ) RETURN * * End of DDISPM. * END * * * SUBROUTINE SDISPM( N, M, LDA, A, DEV, FNAME) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, DEV CHARACTER*(*) FNAME * .. * .. Array Arguments REAL A( LDA, * ) * .. * * Purpose * ======= * * SDISPM displays a N x M REAL matrix A in a file * called FNAME on the device number DEV. It must be M <= 6 . * * LDA is the leading dimension of the array A, LDA >= max( 1, N ) * If DEV is the standard output (here it is set to 6), FNAME can * be blank. * * .. * .. Local Scalars .. INTEGER I, J, STOUT * .. * .. Executable Statements .. * IF( M .GT. 6 ) GO TO 100 STOUT = 6 IF ( DEV .NE. STOUT ) THEN OPEN( UNIT = DEV, FILE = FNAME ) CLOSE( IDEV, STATUS = 'DELETE' ) OPEN( UNIT = IDEV, FILE = FNAME ) END IF DO 10 I = 1, N WRITE( DEV, 20 ) ( A( I, J ), J = 1, M ) 20 FORMAT( 6 ( ' ', E11.5 ) ) 10 CONTINUE WRITE( DEV, * ) ' ' 100 CONTINUE IF ( DEV .NE. STOUT ) CLOSE( DEV ) RETURN * * End of SDISPM. * END * * * SUBROUTINE DSTORM( N, M, LDA, A, DEV, FNAME) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, DEV CHARACTER*(*) FNAME * .. * .. Array Arguments DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DSTORM stores a N x M DOUBLE PRECISION matrix A in a file * called FNAME on the device number DEV after the following scheme: * N * M * A( 1, 1 ) * A( 2, 1 ) * ... * A( N, 1 ) * A( 1, 2 ) * ... * A( N, M ) * * LDA is the leading dimension of the array A, LDA >= max( 1, N ) * If DEV is the standard output (here it is set to 6), FNAME can * be blank. * * .. Local Scalars .. INTEGER I, J, STOUT * .. * .. Executable Statements .. * STOUT = 6 IF( DEV .NE. STOUT ) THEN OPEN( UNIT = DEV , FILE = FNAME ) CLOSE( DEV, STATUS = 'DELETE' ) OPEN( UNIT = DEV, FILE = FNAME ) END IF WRITE( DEV, * ) N WRITE( DEV, * ) M DO 10 I = 1, M DO 10 J = 1, N 10 WRITE( DEV, * ) A( J, I ) IF ( DEV .NE. STOUT ) CLOSE( DEV ) RETURN * * End of DSTORM. * END * * * SUBROUTINE SSTORM( N, M, LDA, A, DEV, FNAME) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, DEV CHARACTER*(*) FNAME * .. * .. Array Arguments REAL A( LDA, * ) * .. * * Purpose * ======= * * SSTORM stores a N x M DOUBLE PRECISION matrix A in a file * called FNAME on the device number DEV after the following scheme: * N * M * A( 1, 1 ) * A( 2, 1 ) * ... * A( N, 1 ) * A( 1, 2 ) * ... * A( N, M ) * * LDA is the leading dimension of the array A, LDA >= max( 1, N ) * If DEV is the standard output (here it is set to 6), FNAME can * be blank. * * * .. Local Scalars .. INTEGER I, J, STOUT * .. * .. Executable Statements .. * STOUT = 6 IF( DEV .NE. STOUT ) THEN OPEN( UNIT = DEV , FILE = FNAME ) CLOSE( DEV, STATUS = 'DELETE' ) OPEN( UNIT = DEV, FILE = FNAME ) END IF WRITE( DEV, * ) N WRITE( DEV, * ) M DO 10 I = 1, M DO 10 J = 1, N 10 WRITE( DEV, * ) A( J, I ) IF ( DEV .NE. STOUT ) CLOSE( DEV ) RETURN * * End of SSTORM. * END * * * SUBROUTINE DREADM( N, M, LDA, A, DEV, FNAME) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, DEV CHARACTER*(*) FNAME * .. * .. Array Arguments DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DREADM reads a N x M DOUBLE PRECISION matrix A from a file * called FNAME on the device number DEV. The matrix A is stored * after the following scheme: * N * M * A( 1, 1 ) * A( 2, 1 ) * ... * A( N, 1 ) * A( 1, 2 ) * ... * A( N, M ) * * LDA is the leading dimension of the array A, LDA >= max( 1, N ) * If DEV is the standard input (here it is set to 5), FNAME can be * blank. * * * .. Local Scalars .. INTEGER I, J, STIN * .. * .. Executable Statements .. * STIN = 5 IF( DEV .NE. STIN ) OPEN( UNIT = DEV , FILE = FNAME ) READ( DEV, * ) N READ( DEV, * ) M DO 10 I = 1, M DO 10 J = 1, N 10 READ( DEV, * ) A( J, I ) IF ( DEV .NE. STIN ) CLOSE( DEV ) RETURN * * End of DREADM. * END * * * SUBROUTINE SREADM( N, M, LDA, A, DEV, FNAME) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, DEV CHARACTER*(*) FNAME * .. * .. Array Arguments REAL A( LDA, * ) * .. * * Purpose * ======= * * SREADM reads a N x M REAL matrix A from a file * called FNAME on the device number DEV. The matrix A is * stored after the following scheme: * N * M * A( 1, 1 ) * A( 2, 1 ) * ... * A( N, 1 ) * A( 1, 2 ) * ... * A( N, M ) * * LDA is the leading dimension of the array A, LDA >= max( 1, N ) * If DEV is the standard input (here it is set to 5), FNAME can be * blank. * * * .. Local Scalars .. INTEGER I, J, STIN * .. * .. Executable Statements .. * STIN = 5 IF( DEV .NE. STIN ) OPEN( UNIT = DEV , FILE = FNAME ) READ( DEV, * ) N READ( DEV, * ) M DO 10 I = 1, M DO 10 J = 1, N 10 READ( DEV, * ) A( J, I ) IF ( DEV .NE. STIN ) CLOSE( DEV ) RETURN * * End of SREADM. * END * * * SUBROUTINE DSCALM( N, LDA, A, D) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Arguments INTEGER N, LDA * .. * .. Array Arguments DOUBLE PRECISION A( LDA, * ), D( * ) * .. * * Purpose * ======= * * DSCALM multiplies a symmetric N x N matrix A from both sides * with a diagonal matrix stored in vector D. * * LDA is the leading dimension of the array A, LDA >= max( 1, N ) * * .. Local Scalars .. INTEGER I, J * .. * .. Executable Statements .. * DO 10 I = 1, N DO 10 J = I, N A( I,J ) = A( I,J ) * D( I ) * D( J ) A( J,I ) = A( I,J) 10 CONTINUE RETURN * * End of DSCALM. * END * * * SUBROUTINE DOUT( N, M, LDA, EIG, VEC, SORT, NAME ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 29, 1992 * * .. Scalar Arguments CHARACTER*(*) NAME INTEGER N, M, LDA * .. * .. Array Arguments DOUBLE PRECISION EIG( * ), VEC( LDA, * ) INTEGER SORT( * ) * .. * * Purpose * ======= * * DOUT sorts M eigenvalues of a N x N matrix in descending order. * Sorted eigenvalues and the corresponding eigenvectors are stored * in the file NAME on unit number DEV. * * Arguments * ========= * * N (input) INTEGER * Length of the eigenvectors. N >= 0. * * M (input) INTEGER * Number of eigenvalues and eigenvectors. N >= M >= 0. * * LDA (input) INTEGER * The leading dimension of the array VEC. LDA >= max(1,N) * * EIG (input) DOUBLE PRECISION array, dimension (M). * Contains the eigenvalues. * * VEC (input) DOUBLE PRECISION array, dimension (LDA,M). * Contains the corresponding eigenvectors. * * SORT (workspace) INTEGER array, dimension (M) * Used for sorting the eigenvalues. * * NAME (input) CHARACTER*(*) * Name of the file where eigenvalues and eigenvectors * are to be stored. * * ============================================================= * * .. Local Scalars .. INTEGER I, J, AUX, DEV * .. * .. Executable Statements .. * * Open the file NAME on the unit DEV, and make sure the file * is empty. If the following unit number does not work, change it. * DEV = 7 OPEN( UNIT = DEV, FILE = NAME ) CLOSE( DEV, STATUS = 'DELETE' ) OPEN( UNIT = DEV, FILE = NAME ) * * Sort the eigenvalues. * DO 10 I = 1, M SORT( I ) = I 10 CONTINUE DO 20 I = 1, M - 1 DO 20 J = I + 1, M IF( EIG( SORT( I ) ) .GE. EIG( SORT( J ) ) ) GO TO 20 AUX = SORT( I ) SORT( I ) = SORT( J ) SORT( J ) = AUX 20 CONTINUE * * Store the eigenvalues. * WRITE( DEV, * ) ' EIGENVALUES ' DO 30 I = 1, M WRITE( DEV, * ) EIG( SORT( I ) ) 30 CONTINUE * * Store the eigenvectors. * DO 40 I = 1, M AUX=SORT( I ) WRITE( DEV, 45 ) I 45 FORMAT( I5, '. EIGENVECTOR ' ) DO 40 J=1,N WRITE( DEV, * ) VEC( J, AUX ) 40 CONTINUE CLOSE( DEV ) RETURN * * End of DOUT. * END * * * SUBROUTINE SOUT( N, M, LDA, EIG, VEC, SORT, NAME ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 29, 1992 * * .. Scalar Arguments CHARACTER*(*) NAME INTEGER N, M, LDA * .. * .. Array Arguments REAL EIG( * ), VEC( LDA, * ) INTEGER SORT( * ) * .. * * Purpose * ======= * * SOUT sorts M eigenvalues of a N x N matrix in descending order. * Sorted eigenvalues and the corresponding eigenvectors are stored * in the file NAME on unit number DEV. * * Arguments * ========= * * N (input) INTEGER * Length of the eigenvectors. N >= 0. * * M (input) INTEGER * Number of eigenvalues and eigenvectors. N >= M >= 0. * * LDA (input) INTEGER * The leading dimension of the array VEC. LDA >= max(1,N) * * EIG (input) REAL array, dimension (M). * Contains the eigenvalues. * * VEC (input) REAL array, dimension (LDA,M). * Contains the corresponding eigenvectors. * * SORT (workspace) INTEGER array, dimension (M) * Used for sorting the eigenvalues. * * NAME (input) CHARACTER*(*) * Name of the file where eigenvalues and eigenvectors * are to be stored. * * ============================================================= * * .. Local Scalars .. INTEGER I, J, AUX, DEV * .. * .. Executable Statements .. * * Open the file NAME on the unit DEV, and make sure the file * is empty. If the following unit number does not work, change it. * DEV = 7 OPEN( UNIT = DEV, FILE = NAME ) CLOSE( DEV, STATUS = 'DELETE' ) OPEN( UNIT = DEV, FILE = NAME ) * * Sort the eigenvalues. * DO 10 I = 1, M SORT( I ) = I 10 CONTINUE DO 20 I = 1, M - 1 DO 20 J = I + 1, M IF( EIG( SORT( I ) ) .GE. EIG( SORT( J ) ) ) GO TO 20 AUX = SORT( I ) SORT( I ) = SORT( J ) SORT( J ) = AUX 20 CONTINUE * * Store the eigenvalues. * WRITE( DEV, * ) ' EIGENVALUES ' DO 30 I = 1, M WRITE( DEV, * ) EIG( SORT( I ) ) 30 CONTINUE * * Store the eigenvectors. * DO 40 I = 1, M AUX=SORT( I ) WRITE( DEV, 45 ) I 45 FORMAT( I5, '. EIGENVECTOR ' ) DO 40 J=1,N WRITE( DEV, * ) VEC( J, AUX ) 40 CONTINUE CLOSE( DEV ) RETURN * * End of SOUT. * END * * * SUBROUTINE COMP( N, M, LDA, PNUM, A, NAME, DEVNUM ) * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 30, 1992 * * .. Scalar Arguments INTEGER N, M, LDA, PNUM * .. * .. Array Arguments CHARACTER*(*) NAME( 7 ) DOUBLE PRECISION A( LDA, * ) INTEGER DEVNUM( * ) * .. * * Purpose * ======= * * COMP compares eigensolutions computed by PNUM different routines. * The eigensolutions are read from files NAME( I ) on units * DEVNUM( I ), where 2 <= I <= PNUM. The first file contains the * eigensolution of reference. * * Eigenvalues are compared relatively: * * | Lambda - Lambda' | * ---------------------- . * | Lambda | * * Eigenvectors are compared in 2-norm || v - v' ||. Here * Here ( Lambda, v ) denotes the eigenpairs of reference, and * ( Lambda', v' ) denotes the corresponding eigenpair computed by * other routines. Therefore, the eigenvalues of reference must be * non-zero. * * Results of the comparison are written in file 'COMP' on unit * DEV. * * Arguments * ========= * * N (input) INTEGER * Length of the eigenvectors. N >= 0. * * M (input) INTEGER * Number of eigenvalues and eigenvectors. N >= M >= 0. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N) * * PNUM (input) INTEGER. PNUM < = 7. * number of compared eigensolutions. * * A (workspace) DOUBLE PRECISION array, dimension (LDA,4). * * NAME (input) CHARACTER*(*) array, dimension (7) * names of the files where the eigensolutions are stored * (by DOUT or SOUT). * * DEVNUM (workspace) INTEGER array, dimension (7) * the numbers of input units. * * * Further Details * =============== * * COMP uses BLAS1 subroutine DAXPY and function DNRM2. * * Although it works if M < N, COMP should only be used if M = N, * that is if the matrix whose eigensolutions are compared is * non-singular. The 2-norm errors in eigenvectors which correspond * to multiple eigenvalues may naturally be large. * * ============================================================= * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER DEV, I, J, K DOUBLE PRECISION TEMP, TEMP1 * .. * .. External Subroutines .. EXTERNAL DAXPY, DCOPY * .. * .. External Functions .. DOUBLE PRECISION DNRM2 EXTERNAL DNRM2 * .. * .. Intrinsic Functions .. INTRINSIC DABS, DMIN1 * .. * .. Executable Statements .. * * Give a unit number to the output file, open it, and write the * heading. If the unit numbers are inconvenient, they con be * changed. * DEV = 7 OPEN( UNIT = DEV, FILE = 'COMP' ) CLOSE( DEV, STATUS = 'DELETE' ) OPEN( UNIT = DEV, FILE = 'COMP' ) WRITE( DEV, 2 ) NAME(1) 2 FORMAT(' COMPARING EIGENSOLUTIONS '// $ ' EIGENSOLUTION OF REFERENCE IS FROM ',A8// $ ' EIGENVALUES ' / ) * * Give unit numbers to input files and open them. * DO 10 I = 1, PNUM DEVNUM( I ) = DEV + I OPEN( UNIT = DEVNUM( I ), FILE = NAME( I ) ) READ( DEVNUM( I ), * ) 10 CONTINUE * * Read the eigenvalues of reference into the first column of * array A, and write them to the output file. * DO 20 J = 1, M READ( DEVNUM( 1 ), * ) A( J, 1 ) WRITE( DEV, 3 ) J, A( J, 1 ) 20 CONTINUE 3 FORMAT( 1X, I4, '. ', E27.18 ) * * Write names of other files. * WRITE( DEV, 4 ) ( NAME( I ), I = 2, PNUM ) 4 FORMAT(/' RELATIVE ERRORS IN EIGENVALUES ' // 8X, 6( A8, 3X ) ) WRITE( DEV, * ) * * Start the main loop. * DO 40 I = 1, M TEMP = A( I, 1 ) DO 50 J = 2, PNUM READ( DEVNUM( J ), * ) TEMP1 A( J,2 ) = DABS( ( TEMP - TEMP1 ) / TEMP ) 50 CONTINUE WRITE( DEV, 5 ) I, ( A( J,2 ), J = 2, PNUM ) 5 FORMAT( 1X, I4 , '. ' , 6( E7.1, 4X ) ) 40 CONTINUE * * End of the main loop. * * WRITE(DEV,7) ( NAME(I),I = 2, PNUM ) 7 FORMAT(/ ' NORM ERRORS IN EIGENVECTORS '// 8X, 6( A8, 3X ) ) WRITE( DEV, * ) * * Start the main loop. * DO 60 I = 1, M * * Read the I-th eigenvector of reference into the first column * of array A. * READ( DEVNUM( 1 ), * ) DO 70 J = 1, N READ( DEVNUM( 1 ), * ) A( J, 1 ) 70 CONTINUE DO 80 K = 2, PNUM * * Read the I-th eigenvector from the K-th file into the second * column of A, and copy it into the third column. * READ( DEVNUM( K ), * ) DO 90 J = 1, N READ( DEVNUM( K ), * ) A( J, 2 ) 90 CONTINUE CALL DCOPY( N, A( 1,2 ), 1, A( 1,3 ), 1 ) * * Form the sum and the difference of the eigenvector of reference * and the corresponding eigenvector from the K-th file. * CALL DAXPY( N, ONE, A( 1,1 ), 1, A( 1,2 ), 1 ) CALL DAXPY( N, - ONE, A( 1,1 ), 1, A( 1,3 ), 1 ) * * We take the smaller of the 2-norms. * A( K,4 ) = DMIN1( DNRM2( N, A( 1,2 ), 1 ), $ DNRM2( N, A( 1,3 ), 1 ) ) 80 CONTINUE * Write the 2-nrom errors. * WRITE( DEV, 5 ) I, ( A( K,4 ), K = 2, PNUM ) * * End of the main loop. * 60 CONTINUE * * Close units. * CLOSE( DEV ) DO 100 K = 1, PNUM CLOSE( DEVNUM( K ) ) 100 CONTINUE RETURN * * End of COMP. * END * * * DOUBLE PRECISION FUNCTION RAN1(IDUM) *---------------------------------------------------------------- * THIS FUNCTION IS A STANDART RANDOM NUMBER GENERATOR * IT WAS SUGGESTED IN 'NUMERICAL RECIPIES'. *---------------------------------------------------------------- DIMENSION R(97) PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6) PARAMETER (M3=243000,IA3=4561,IC3=51349) DATA IFF /0/ IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF=1 IX1=MOD(IC1-IDUM,M1) IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IX1,M2) IX1=MOD(IA1*IX1+IC1,M1) IX3=MOD(IX1,M3) DO 11 J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 11 CONTINUE IDUM=1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 IF(J.GT.97.OR.J.LT.1)PAUSE 'IN RAN1' RAN1=R(J) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 RETURN * * End of RAN1. * END CUT HERE.......... cat > dtestj.f <<'CUT HERE..........' PROGRAM DTESTJ * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Variables INTEGER N, RANK, LDA, I, J CHARACTER JOB * .. * .. Array Variables DOUBLE PRECISION H( 50, 50 ), EV( 50 ), V( 50 ) INTEGER PIVOT( 50 ) * .. * .. External Subroutines .. EXTERNAL DSYEVJ, DDISPM * .. * * Purpose * ======= * * DTESTJ is a small example of how to use the subroutine DSYEVJ * which computes the non-zero eigenvalues and the corresponding * eigenvectors of a real symmetric matrix. * * * Further Details * =============== * * The object file dtestj.o must be linked with the level one BLAS, * and with the object files * dsyevj.o * dgjgt.o * djgjf.o * blasj.o (we need the subroutines DROTJJ, DROTF, DROTFS * and DIPERR) * lamch.o (we need the LAPACK auxiliary function DLAMCH, * which determines the machine precision) * auxil.o (we need the subroutine DDISPM) * * .. * .. Executable Statements .. * * Set LDA to the leading dimension of the arrays, and choose * fast rotations. * LDA = 50 JOB = 'F' * * Generate the test matrix, and display it on the screen. * * ! It is interesting to compare the computed eigensolution * ! with the eigensolutions obtained by the single precision * ! routine SSYEVJ (the program STESTJ.F is supplied with the * ! same test matrix), and the LAPACK routines DSYEV and SSYEV. * N = 4 H( 1, 1 ) = 1600.0D0 H( 1, 2 ) = -300.0D0 H( 2, 1 ) = H( 1, 2 ) H( 1, 3 ) = 14.0D0 H( 3, 1 ) = H( 1, 3 ) H( 1, 4 ) = 300000.0D0 H( 4, 1 ) = H( 1, 4 ) H( 2, 2 ) = 43.5D0 H( 2, 3 ) = -4.75D0 H( 3, 2 ) = H( 2, 3 ) H( 2, 4 ) = -423212.0D0 H( 4, 2 ) = H( 2, 4 ) H( 3, 3 ) = 0.1875D0 H( 3, 4 ) = 19800.0D0 H( 4, 3 ) = H( 3, 4 ) H( 4, 4 ) = 3207938D3 WRITE(*,*) ' Test matrix: ' CALL DDISPM( N, N, LDA, H, 6, ' ' ) * * Calculate the non-zero eigenvalues and the corresponding * eigenvectors of H. * CALL DSYEVJ( N, RANK, LDA, H, EV, V, PIVOT, JOB ) * * Display on the screen. * WRITE(*,*) ' Non-zero eigenvalues computed by DSYEVJ: ' CALL DDISPM( RANK, 1, LDA, EV, 6, ' ' ) WRITE(*,*) ' Corresponding eigenvectors: ' CALL DDISPM( N, RANK, LDA, H, 6, ' ' ) * STOP * * End of DTESTJ. * END CUT HERE.......... cat > stestj.f <<'CUT HERE..........' PROGRAM STESTJ * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * June 21, 1992 * * .. Scalar Variables INTEGER N, RANK, LDA, I, J CHARACTER JOB * .. * .. Array Variables REAL H( 50, 50 ), EV( 50 ), V( 50 ) INTEGER PIVOT( 50 ) * .. * .. External Subroutines .. EXTERNAL SSYEVJ, SDISPM * .. * * Purpose * ======= * * STESTJ is a small example of how to use the subroutine SSYEVJ * which computes the non-zero eigenvalues and the corresponding * eigenvectors of a real symmetric matrix. * * * Further Details * =============== * * The object file stestj.o must be linked with the level one BLAS, * and with the object files * ssyevj.o * sgjgt.o * sjgjf.o * blasj.o (we need the subroutines SROTJJ, SROTF, SROTFS * and SIPERR) * lamch.o (we need the LAPACK auxiliary function SLAMCH, * which determines the machine precision) * auxil.o (we need the subroutine SDISPM) * * .. * .. Executable Statements .. * * Set LDA to the leading dimension of the arrays, and choose * fast rotations. * LDA = 50 JOB = 'F' * * Generate the test matrix, and display it on the screen. * * ! It is interesting to compare the computed eigensolution * ! with the eigensolutions obtained by the double precision * ! routine DSYEVJ (the program DTESTJ.F is supplied with the * ! same test matrix), and the LAPACK routines DSYEV and SSYEV. * N = 4 H( 1, 1 ) = 1600.0E0 H( 1, 2 ) = -300.0E0 H( 2, 1 ) = H( 1, 2 ) H( 1, 3 ) = 14.0E0 H( 3, 1 ) = H( 1, 3 ) H( 1, 4 ) = 300000.0E0 H( 4, 1 ) = H( 1, 4 ) H( 2, 2 ) = 43.5E0 H( 2, 3 ) = -4.75E0 H( 3, 2 ) = H( 2, 3 ) H( 2, 4 ) = -423212.0E0 H( 4, 2 ) = H( 2, 4 ) H( 3, 3 ) = 0.1875E0 H( 3, 4 ) = 19800.0E0 H( 4, 3 ) = H( 3, 4 ) H( 4, 4 ) = 3207938E3 WRITE(*,*) ' Test matrix: ' CALL SDISPM( N, N, LDA, H, 6, ' ' ) * * Calculate the non-zero eigenvalues and the corresponding * eigenvectors of H. * CALL SSYEVJ( N, RANK, LDA, H, EV, V, PIVOT, JOB ) * * Display on the screen. * WRITE(*,*) ' Non-zero eigenvalues computed by SSYEVJ: ' CALL SDISPM( RANK, 1, LDA, EV, 6, ' ' ) WRITE(*,*) ' Corresponding eigenvectors: ' CALL SDISPM( N, RANK, LDA, H, 6, ' ' ) * STOP * * End of STESTJ. * END CUT HERE.......... cat > compar.f <<'CUT HERE..........' PROGRAM COMPAR * * Ivan Slapnicar, Faculty of Electrical Engineering, Mechanical * Engineering and Naval Architecture, R. Boskovica b.b, * 58000 Split, Croatia, e-mail islapnicar at uni-zg.ac.mail.yu * This version: * July 6, 1992 * * .. Scalar Variables INTEGER SEED, N, LDA, I, J, DEV, RANKD, RANKS INTEGER INFOD, INFOS, PNUM DOUBLE PRECISION ASCAL, HSCAL CHARACTER JOB * .. * .. Array Variables CHARACTER*(8) NAME( 7 ) DOUBLE PRECISION H( 50, 50 ), EV( 50 ), V( 50 ) DOUBLE PRECISION EVEC( 50, 50 ) REAL SH( 50, 50 ), SEV( 50 ), SV( 50 ) REAL SEVEC( 50, 50 ) INTEGER PIVOT( 50 ) * .. * .. External Subroutines .. EXTERNAL GENSYM, DSYEVJ, SSYEVJ, DSYEV, SSYEV EXTERNAL DJAC, SJAC, DOUT, SOUT, COMP * .. * * Purpose * ======= * * Program COMPAR compares eigensolutions obtained by several different * routines. User can add one more routine or substitute some routines * with his own ones. * * COMPAR first calls the subroutine GENSYM to generate a non-singular * symmetric matrix. User can, at its will, change this generator or read * a symmetric matrix. * * GENSYM generates a non-singular symmetric matrix H of order N, * with condition approximately equal to 10 ** HSCAL and the measure * C( A, A ) from [1] approximately equal to 10 ** ASCAL. * * For meaningful experiments, ASCAL should be smaller than * LOG10( 1 / SEPS ), where SEPS is machine precision in single precision, * e.g. 0 <= ASCAL <= 4. HSCAL should be smaller than logarithm of the * overflow in single precision, e.g. 0 <= HSCAL <= 30. * * * COMPAR then solves the eigenproblem by double precision routines * * DSYEVJ - implements symmetric indefinite decomposition followed * by implicit J-orthogonal Jacobi, * * DSYEV - implements tridiagonalization followed by QR iteration * (this is a LAPACK driver routine which can be obtained * from NETLIB. If DSYEV is not readily available, * delete the part where it is being called, and change * the character array NAME, and the number of programs * PNUM at the end of COMPAR. This remark holds for * LAPACK driver routine SSYEV, as well.), * * DJAC - implements the standard Jacobi algorithm, * * and by the corresponding single precision routines SSYEVJ, * SSYEV (LAPACK), and SJAC. * * * Finally, COMPAR calls the subroutine COMP to compare the eigensolution * computed by DSYEVJ with the eigensolutions computed by other routines. * COMP computes relativle errors in eigenvalues and 2-nrom errors in * eigenvectors. Output is in file 'COMP'. For more details see comments * of the subroutine COMP in file auxil.f . * * * ! For matrices generated by GENSYM, our single precision algorithm * ! SSYEVJ is likely to be more accurate than SSYEV and SJAC, and * ! sometimes even more accurate than DSYEV. More precisely, * ! eigenvalues computed by SSYEVJ have approximately * ! LOG10( 1 / SEPS ) - ASCAL accurate decimal digits. For other * ! types of matrices our algorithms are not worse than other * ! algorithms. * * * For more theoretical and practical details see [1]. * * References * ========== * * [1] I. Slapnicar: * Accurate Symmetric Eigenreduction by a Jacobi Method, Ph.D. Thesis, * FB Mathematik, Fernuniversitaet Hagen, Hagen, 1992. * * * Further Details * =============== * * The object file compar.o must be linked with the object files * * dsyevj.o , ssyevj.o , * dgjgt.o , sgjgt.o , * djgjf.o , sjgjf.o , * djac.o , sjac.o , * blasj.o , * gensym.o , auxil.o . * * We also need to link LAPACK driver routines DSYEV and SSYEV and * all routines which they use. DSYEV and SSYEV can be obtained from * NETLIB. Since they already contain the auxiliary functions DLAMCH * and SLAMCH, respectively, there is no need to link the file * lamch.o. * * ================================================================ * .. * .. Executable Statements .. * * Set LDA to the leading dimension of the arrays, choose fast * rotations, and set the unit number for hard disc. * LDA = 50 JOB = 'F' DEV = 7 * * Define order of the test matrix N < = 50. (For larger N's, just * increase dimensions of arrays and LDA.) * WRITE( *, * ) ' Order of the matrix N <= 50 ? ' READ( *, * ) N * * Define parameters ASCAL and HSCAL for the random matrix generator. * WRITE( *, * ) ' Parameter ASCAL (0 <= ASCAL <= 4) ? ' READ( *, * ) ASCAL WRITE( *, * ) ' Parameter HSCAL (0 <= HSCAL <= 30) ? ' READ( *, * ) HSCAL * * Define the seed for the random number generator. * WRITE( *, * ) ' Seed for the random number generator ', $ ' (integer > 0) ' READ( *, * ) SEED * * Generate random symmetric matrix H of order N using * SEED, ASCAL and HSCAL. * CALL GENSYM( N, LDA, - SEED, ASCAL, HSCAL, H, EVEC, EV, V, PIVOT ) * * If N <= 6, display the test matrix. * IF( N .LE. 6 ) THEN WRITE(*,*) ' Test matrix: ' CALL DDISPM( N, N, LDA, H, 6, ' ' ) END IF * * Copy double precision matrix into real matrix. * DO 10 I = 1, N DO 10 J = 1, N SH( I,J ) = H( I,J ) 10 CONTINUE * * Store both matrices. * CALL DSTORM( N, N, LDA, H, DEV, 'DSYMM' ) CALL SSTORM( N, N, LDA, SH, DEV, 'SSYMM' ) WRITE( *, * ) ' Test matrix is stored in files DSYMM ', $ ' and SSYMM. ' * * Solve the eigenproblem using DSYEVJ and store the computed * eigensolution in the file 'DSYEVJ'. * CALL DSYEVJ( N, RANKD, LDA, H, EV, V, PIVOT, JOB ) CALL DOUT( N, RANKD, LDA, EV, H, PIVOT, 'DSYEVJ' ) * * Solve the eigenproblem using DSYEV and store the computed * eigensolution in the file 'DSYEV'. * CALL DREADM( N, N, LDA, H, DEV, 'DSYMM' ) CALL DSYEV( 'V', 'L', N, H, LDA, EV, V, 3 * LDA, INFOD ) CALL DOUT( N, N, LDA, EV, H, PIVOT, 'DSYEV' ) * * Solve the eigenproblem using DJAC and store the computed * eigensolution in the file 'DJAC'. * CALL DREADM( N, N, LDA, H, DEV, 'DSYMM' ) CALL DJAC( N, LDA, H, EVEC, EV, V ) CALL DOUT( N, N, LDA, EV, EVEC, PIVOT, 'DJAC' ) * * Solve the eigenproblem using SSYEVJ and store the computed * eigensolution in the file 'SSYEVJ'. * CALL SSYEVJ( N, RANKS, LDA, SH, SEV, SV, PIVOT, JOB ) CALL SOUT( N, RANKS, LDA, SEV, SH, PIVOT, 'SSYEVJ' ) * * Solve the eigenproblem using SSYEV and store the computed * eigensolution in the file 'SSYEV'. * CALL SREADM( N, N, LDA, SH, DEV, 'SSYMM' ) CALL SSYEV( 'V', 'L', N, SH, LDA, SEV, SV, 3 * LDA, INFOS ) CALL SOUT( N, N, LDA, SEV, SH, PIVOT, 'SSYEV' ) * * Solve the eigenproblem using SJAC and store the computed * eigensolution in the file 'SJAC'. * CALL SREADM( N, N, LDA, SH, DEV, 'SSYMM' ) CALL SJAC( N, LDA, SH, SEVEC, SEV, SV ) CALL SOUT( N, N, LDA, SEV, SEVEC, PIVOT, 'SJAC' ) * * Compare eigensolution computed by DSYEVJ with eigensolutions * computed by DSYEV, DJAC, SSYEVJ, SSYEV and SJAC. The comparison * is made only if all subroutines returned all N eigenvalues. * IF( ( RANKD .EQ. N ) .AND. ( INFOD .EQ. 0 ) .AND. $ ( RANKS .EQ. N ) .AND. ( INFOS .EQ. 0 ) ) THEN NAME( 1 ) = 'DSYEVJ' NAME( 2 ) = 'DSYEV' NAME( 3 ) = 'DJAC' NAME( 4 ) = 'SSYEVJ' NAME( 5 ) = 'SSYEV' NAME( 6 ) = 'SJAC' PNUM = 6 CALL COMP( N, N, LDA, PNUM, H, NAME, PIVOT ) WRITE( *, * ) ' Output of the comparison is in file COMP. ' ELSE WRITE( *, * ) ' Some subroutines did not return all N ', $ 'eigenvalues.' WRITE( *, * ) ' The comparison is not performed. ' END IF * STOP * * End of COMPAR. * END CUT HERE.......... cat > lamch.f <<'CUT HERE..........' REAL FUNCTION SLAMCH( CMACH ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER CMACH * .. * * Purpose * ======= * * SLAMCH determines single precision machine parameters. * * Arguments * ========= * * CMACH (input) CHARACTER*1 * Specifies the value to be returned by SLAMCH: * = 'E' or 'e', SLAMCH := eps * = 'S' or 's , SLAMCH := sfmin * = 'B' or 'b', SLAMCH := base * = 'P' or 'p', SLAMCH := eps*base * = 'N' or 'n', SLAMCH := t * = 'R' or 'r', SLAMCH := rnd * = 'M' or 'm', SLAMCH := emin * = 'U' or 'u', SLAMCH := rmin * = 'L' or 'l', SLAMCH := emax * = 'O' or 'o', SLAMCH := rmax * * where * * eps = relative machine precision * sfmin = safe minimum, such that 1/sfmin does not overflow * base = base of the machine * prec = eps*base * t = number of (base) digits in the mantissa * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise * emin = minimum exponent before (gradual) underflow * rmin = underflow threshold - base**(emin-1) * emax = largest exponent before overflow * rmax = overflow threshold - (base**emax)*(1-eps) * * * .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. LOGICAL FIRST, LRND INTEGER BETA, IMAX, IMIN, IT REAL BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, $ RND, SFMIN, SMALL, T * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL SLAMC2 * .. * .. Save statement .. SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, $ EMAX, RMAX, PREC * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) BASE = BETA T = IT IF( LRND ) THEN RND = ONE EPS = ( BASE**( 1-IT ) ) / 2 ELSE RND = ZERO EPS = BASE**( 1-IT ) END IF PREC = EPS*BASE EMIN = IMIN EMAX = IMAX SFMIN = RMIN SMALL = ONE / RMAX IF( SMALL.GE.SFMIN ) THEN * * Use SMALL plus a bit, to avoid the possibility of rounding * causing overflow when computing 1/sfmin. * SFMIN = SMALL*( ONE+EPS ) END IF END IF * IF( LSAME( CMACH, 'E' ) ) THEN RMACH = EPS ELSE IF( LSAME( CMACH, 'S' ) ) THEN RMACH = SFMIN ELSE IF( LSAME( CMACH, 'B' ) ) THEN RMACH = BASE ELSE IF( LSAME( CMACH, 'P' ) ) THEN RMACH = PREC ELSE IF( LSAME( CMACH, 'N' ) ) THEN RMACH = T ELSE IF( LSAME( CMACH, 'R' ) ) THEN RMACH = RND ELSE IF( LSAME( CMACH, 'M' ) ) THEN RMACH = EMIN ELSE IF( LSAME( CMACH, 'U' ) ) THEN RMACH = RMIN ELSE IF( LSAME( CMACH, 'L' ) ) THEN RMACH = EMAX ELSE IF( LSAME( CMACH, 'O' ) ) THEN RMACH = RMAX END IF * SLAMCH = RMACH RETURN * * End of SLAMCH * END * ************************************************************************ * SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T * .. * * Purpose * ======= * * SLAMC1 determines the machine parameters given by BETA, T, RND, and * IEEE1. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * IEEE1 (output) LOGICAL * Specifies whether rounding appears to be done in the IEEE * 'round to nearest' style. * * Further Details * =============== * * The routine is based on the routine ENVRON by Malcolm and * incorporates suggestions by Gentleman and Marovich. See * * Malcolm M. A. (1972) Algorithms to reveal properties of * floating-point arithmetic. Comms. of the ACM, 15, 949-951. * * Gentleman W. M. and Marovich S. B. (1974) More on algorithms * that reveal properties of floating point arithmetic units. * Comms. of the ACM, 17, 276-277. * * * .. Local Scalars .. LOGICAL FIRST, LIEEE1, LRND INTEGER LBETA, LT REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2 * .. * .. External Functions .. REAL SLAMC3 EXTERNAL SLAMC3 * .. * .. Save statement .. SAVE FIRST, LIEEE1, LBETA, LRND, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ONE = 1 * * LBETA, LIEEE1, LT and LRND are the local values of BETA, * IEEE1, T and RND. * * Throughout this routine we use the function SLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * Compute a = 2.0**m with the smallest positive integer m such * that * * fl( a + 1.0 ) = a. * A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 10 CONTINUE IF( C.EQ.ONE ) THEN A = 2*A C = SLAMC3( A, ONE ) C = SLAMC3( C, -A ) GO TO 10 END IF *+ END WHILE * * Now compute b = 2.0**m with the smallest positive integer m * such that * * fl( a + b ) .gt. a. * B = 1 C = SLAMC3( A, B ) * *+ WHILE( C.EQ.A )LOOP 20 CONTINUE IF( C.EQ.A ) THEN B = 2*B C = SLAMC3( A, B ) GO TO 20 END IF *+ END WHILE * * Now compute the base. a and c are neighbouring floating point * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so * their difference is beta. Adding 0.25 to c is to ensure that it * is truncated to beta and not ( beta - 1 ). * QTR = ONE / 4 SAVEC = C C = SLAMC3( C, -A ) LBETA = C + QTR * * Now determine whether rounding or chopping occurs, by adding a * bit less than beta/2 and a bit more than beta/2 to a. * B = LBETA F = SLAMC3( B / 2, -B / 100 ) C = SLAMC3( F, A ) IF( C.EQ.A ) THEN LRND = .TRUE. ELSE LRND = .FALSE. END IF F = SLAMC3( B / 2, B / 100 ) C = SLAMC3( F, A ) IF( ( LRND ) .AND. ( C.EQ.A ) ) $ LRND = .FALSE. * * Try and decide whether rounding is done in the IEEE 'round to * nearest' style. B/2 is half a unit in the last place of the two * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit * zero, and SAVEC is odd. Thus adding B/2 to A should not change * A, but adding B/2 to SAVEC should change SAVEC. * T1 = SLAMC3( B / 2, A ) T2 = SLAMC3( B / 2, SAVEC ) LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND * * Now find the mantissa, t. It should be the integer part of * log to the base beta of a, however it is safer to determine t * by powering. So we find t as the smallest positive integer for * which * * fl( beta**t + 1.0 ) = 1.0. * LT = 0 A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 30 CONTINUE IF( C.EQ.ONE ) THEN LT = LT + 1 A = A*LBETA C = SLAMC3( A, ONE ) C = SLAMC3( C, -A ) GO TO 30 END IF *+ END WHILE * END IF * BETA = LBETA T = LT RND = LRND IEEE1 = LIEEE1 RETURN * * End of SLAMC1 * END * ************************************************************************ * SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL RND INTEGER BETA, EMAX, EMIN, T REAL EPS, RMAX, RMIN * .. * * Purpose * ======= * * SLAMC2 determines the machine parameters specified in its argument * list. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * EPS (output) REAL * The smallest positive number such that * * fl( 1.0 - EPS ) .LT. 1.0, * * where fl denotes the computed value. * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow occurs. * * RMIN (output) REAL * The smallest normalized number for the machine, given by * BASE**( EMIN - 1 ), where BASE is the floating point value * of BETA. * * EMAX (output) INTEGER * The maximum exponent before overflow occurs. * * RMAX (output) REAL * The largest positive number for the machine, given by * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point * value of BETA. * * Further Details * =============== * * The computation of EPS is based on a routine PARANOIA by * W. Kahan of the University of California at Berkeley. * * * .. Local Scalars .. LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, $ NGNMIN, NGPMIN REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, $ SIXTH, SMALL, THIRD, TWO, ZERO * .. * .. External Functions .. REAL SLAMC3 EXTERNAL SLAMC3 * .. * .. External Subroutines .. EXTERNAL SLAMC1, SLAMC4, SLAMC5 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Save statement .. SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, $ LRMIN, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / , IWARN / .FALSE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ZERO = 0 ONE = 1 TWO = 2 * * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of * BETA, T, RND, EPS, EMIN and RMIN. * * Throughout this routine we use the function SLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. * CALL SLAMC1( LBETA, LT, LRND, LIEEE1 ) * * Start to find EPS. * B = LBETA A = B**( -LT ) LEPS = A * * Try some tricks to see whether or not this is the correct EPS. * B = TWO / 3 HALF = ONE / 2 SIXTH = SLAMC3( B, -HALF ) THIRD = SLAMC3( SIXTH, SIXTH ) B = SLAMC3( THIRD, -HALF ) B = SLAMC3( B, SIXTH ) B = ABS( B ) IF( B.LT.LEPS ) $ B = LEPS * LEPS = 1 * *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 10 CONTINUE IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN LEPS = B C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) C = SLAMC3( HALF, -C ) B = SLAMC3( HALF, C ) C = SLAMC3( HALF, -B ) B = SLAMC3( HALF, C ) GO TO 10 END IF *+ END WHILE * IF( A.LT.LEPS ) $ LEPS = A * * Computation of EPS complete. * * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). * Keep dividing A by BETA until (gradual) underflow occurs. This * is detected when we cannot recover the previous A. * RBASE = ONE / LBETA SMALL = ONE DO 20 I = 1, 3 SMALL = SLAMC3( SMALL*RBASE, ZERO ) 20 CONTINUE A = SLAMC3( ONE, SMALL ) CALL SLAMC4( NGPMIN, ONE, LBETA ) CALL SLAMC4( NGNMIN, -ONE, LBETA ) CALL SLAMC4( GPMIN, A, LBETA ) CALL SLAMC4( GNMIN, -A, LBETA ) IEEE = .FALSE. * IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN IF( NGPMIN.EQ.GPMIN ) THEN LEMIN = NGPMIN * ( Non twos-complement machines, no gradual underflow; * e.g., VAX ) ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN LEMIN = NGPMIN - 1 + LT IEEE = .TRUE. * ( Non twos-complement machines, with gradual underflow; * e.g., IEEE standard followers ) ELSE LEMIN = MIN( NGPMIN, GPMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) * ( Twos-complement machines, no gradual underflow; * e.g., CYBER 205 ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. $ ( GPMIN.EQ.GNMIN ) ) THEN IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT * ( Twos-complement machines with gradual underflow; * no known machine ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF *** * Comment out this if block if EMIN is ok IF( IWARN ) THEN FIRST = .TRUE. WRITE( 6, FMT = 9999 )LEMIN END IF *** * * Assume IEEE arithmetic if we found denormalised numbers above, * or if arithmetic seems to round in the IEEE style, determined * in routine SLAMC1. A true IEEE machine should have both things * true; however, faulty machines may have one or the other. * IEEE = IEEE .OR. LIEEE1 * * Compute RMIN by successive division by BETA. We could compute * RMIN as BASE**( EMIN - 1 ), but some machines underflow during * this computation. * LRMIN = 1 DO 30 I = 1, 1 - LEMIN LRMIN = SLAMC3( LRMIN*RBASE, ZERO ) 30 CONTINUE * * Finally, call SLAMC5 to compute EMAX and RMAX. * CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) END IF * BETA = LBETA T = LT RND = LRND EPS = LEPS EMIN = LEMIN RMIN = LRMIN EMAX = LEMAX RMAX = LRMAX * RETURN * 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', $ ' EMIN = ', I8, / $ ' If, after inspection, the value EMIN looks', $ ' acceptable please comment out ', $ / ' the IF block as marked within the code of routine', $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / ) * * End of SLAMC2 * END * ************************************************************************ * REAL FUNCTION SLAMC3( A, B ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. REAL A, B * .. * * Purpose * ======= * * SLAMC3 is intended to force A and B to be stored prior to doing * the addition of A and B , for use in situations where optimizers * might hold one of these in a register. * * Arguments * ========= * * A, B (input) REAL * The values A and B. * * * .. Executable Statements .. * SLAMC3 = A + B * RETURN * * End of SLAMC3 * END * ************************************************************************ * SUBROUTINE SLAMC4( EMIN, START, BASE ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. INTEGER BASE, EMIN REAL START * .. * * Purpose * ======= * * SLAMC4 is a service routine for SLAMC2. * * Arguments * ========= * * EMIN (output) EMIN * The minimum exponent before (gradual) underflow, computed by * setting A = START and dividing by BASE until the previous A * can not be recovered. * * START (input) REAL * The starting point for determining EMIN. * * BASE (input) INTEGER * The base of the machine. * * * .. Local Scalars .. INTEGER I REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO * .. * .. External Functions .. REAL SLAMC3 EXTERNAL SLAMC3 * .. * .. Executable Statements .. * A = START ONE = 1 RBASE = ONE / BASE ZERO = 0 EMIN = 1 B1 = SLAMC3( A*RBASE, ZERO ) C1 = A C2 = A D1 = A D2 = A *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 10 CONTINUE IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. $ ( D2.EQ.A ) ) THEN EMIN = EMIN - 1 A = B1 B1 = SLAMC3( A / BASE, ZERO ) C1 = SLAMC3( B1*BASE, ZERO ) D1 = ZERO DO 20 I = 1, BASE D1 = D1 + B1 20 CONTINUE B2 = SLAMC3( A*RBASE, ZERO ) C2 = SLAMC3( B2 / RBASE, ZERO ) D2 = ZERO DO 30 I = 1, BASE D2 = D2 + B2 30 CONTINUE GO TO 10 END IF *+ END WHILE * RETURN * * End of SLAMC4 * END * ************************************************************************ * SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL IEEE INTEGER BETA, EMAX, EMIN, P REAL RMAX * .. * * Purpose * ======= * * SLAMC5 attempts to compute RMAX, the largest machine floating-point * number, without overflow. It assumes that EMAX + abs(EMIN) sum * approximately to a power of 2. It will fail on machines where this * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, * EMAX = 28718). It will also fail if the value supplied for EMIN is * too large (i.e. too close to zero), probably with overflow. * * Arguments * ========= * * BETA (input) INTEGER * The base of floating-point arithmetic. * * P (input) INTEGER * The number of base BETA digits in the mantissa of a * floating-point value. * * EMIN (input) INTEGER * The minimum exponent before (gradual) underflow. * * IEEE (input) LOGICAL * A logical flag specifying whether or not the arithmetic * system is thought to comply with the IEEE standard. * * EMAX (output) INTEGER * The largest exponent before overflow * * RMAX (output) REAL * The largest machine floating-point number. * * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) * .. * .. Local Scalars .. INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP REAL OLDY, RECBAS, Y, Z * .. * .. External Functions .. REAL SLAMC3 EXTERNAL SLAMC3 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * * First compute LEXP and UEXP, two powers of 2 that bound * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum * approximately to the bound that is closest to abs(EMIN). * (EMAX is the exponent of the required number RMAX). * LEXP = 1 EXBITS = 1 10 CONTINUE TRY = LEXP*2 IF( TRY.LE.( -EMIN ) ) THEN LEXP = TRY EXBITS = EXBITS + 1 GO TO 10 END IF IF( LEXP.EQ.-EMIN ) THEN UEXP = LEXP ELSE UEXP = TRY EXBITS = EXBITS + 1 END IF * * Now -LEXP is less than or equal to EMIN, and -UEXP is greater * than or equal to EMIN. EXBITS is the number of bits needed to * store the exponent. * IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN EXPSUM = 2*LEXP ELSE EXPSUM = 2*UEXP END IF * * EXPSUM is the exponent range, approximately equal to * EMAX - EMIN + 1 . * EMAX = EXPSUM + EMIN - 1 NBITS = 1 + EXBITS + P * * NBITS is the total number of bits needed to store a * floating-point number. * IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN * * Either there are an odd number of bits used to store a * floating-point number, which is unlikely, or some bits are * not used in the representation of numbers, which is possible, * (e.g. Cray machines) or the mantissa has an implicit bit, * (e.g. IEEE machines, Dec Vax machines), which is perhaps the * most likely. We have to assume the last alternative. * If this is true, then we need to reduce EMAX by one because * there must be some way of representing zero in an implicit-bit * system. On machines like Cray, we are reducing EMAX by one * unnecessarily. * EMAX = EMAX - 1 END IF * IF( IEEE ) THEN * * Assume we are on an IEEE machine which reserves one exponent * for infinity and NaN. * EMAX = EMAX - 1 END IF * * Now create RMAX, the largest machine number, which should * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . * * First compute 1.0 - BETA**(-P), being careful that the * result is less than 1.0 . * RECBAS = ONE / BETA Z = BETA - ONE Y = ZERO DO 20 I = 1, P Z = Z*RECBAS IF( Y.LT.ONE ) $ OLDY = Y Y = SLAMC3( Y, Z ) 20 CONTINUE IF( Y.GE.ONE ) $ Y = OLDY * * Now multiply by BETA**EMAX to get RMAX. * DO 30 I = 1, EMAX Y = SLAMC3( Y*BETA, ZERO ) 30 CONTINUE * RMAX = Y RETURN * * End of SLAMC5 * END * * * DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER CMACH * .. * * Purpose * ======= * * DLAMCH determines double precision machine parameters. * * Arguments * ========= * * CMACH (input) CHARACTER*1 * Specifies the value to be returned by DLAMCH: * = 'E' or 'e', DLAMCH := eps * = 'S' or 's , DLAMCH := sfmin * = 'B' or 'b', DLAMCH := base * = 'P' or 'p', DLAMCH := eps*base * = 'N' or 'n', DLAMCH := t * = 'R' or 'r', DLAMCH := rnd * = 'M' or 'm', DLAMCH := emin * = 'U' or 'u', DLAMCH := rmin * = 'L' or 'l', DLAMCH := emax * = 'O' or 'o', DLAMCH := rmax * * where * * eps = relative machine precision * sfmin = safe minimum, such that 1/sfmin does not overflow * base = base of the machine * prec = eps*base * t = number of (base) digits in the mantissa * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise * emin = minimum exponent before (gradual) underflow * rmin = underflow threshold - base**(emin-1) * emax = largest exponent before overflow * rmax = overflow threshold - (base**emax)*(1-eps) * * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL FIRST, LRND INTEGER BETA, IMAX, IMIN, IT DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, $ RND, SFMIN, SMALL, T * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLAMC2 * .. * .. Save statement .. SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, $ EMAX, RMAX, PREC * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) BASE = BETA T = IT IF( LRND ) THEN RND = ONE EPS = ( BASE**( 1-IT ) ) / 2 ELSE RND = ZERO EPS = BASE**( 1-IT ) END IF PREC = EPS*BASE EMIN = IMIN EMAX = IMAX SFMIN = RMIN SMALL = ONE / RMAX IF( SMALL.GE.SFMIN ) THEN * * Use SMALL plus a bit, to avoid the possibility of rounding * causing overflow when computing 1/sfmin. * SFMIN = SMALL*( ONE+EPS ) END IF END IF * IF( LSAME( CMACH, 'E' ) ) THEN RMACH = EPS ELSE IF( LSAME( CMACH, 'S' ) ) THEN RMACH = SFMIN ELSE IF( LSAME( CMACH, 'B' ) ) THEN RMACH = BASE ELSE IF( LSAME( CMACH, 'P' ) ) THEN RMACH = PREC ELSE IF( LSAME( CMACH, 'N' ) ) THEN RMACH = T ELSE IF( LSAME( CMACH, 'R' ) ) THEN RMACH = RND ELSE IF( LSAME( CMACH, 'M' ) ) THEN RMACH = EMIN ELSE IF( LSAME( CMACH, 'U' ) ) THEN RMACH = RMIN ELSE IF( LSAME( CMACH, 'L' ) ) THEN RMACH = EMAX ELSE IF( LSAME( CMACH, 'O' ) ) THEN RMACH = RMAX END IF * DLAMCH = RMACH RETURN * * End of DLAMCH * END * ************************************************************************ * SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T * .. * * Purpose * ======= * * DLAMC1 determines the machine parameters given by BETA, T, RND, and * IEEE1. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * IEEE1 (output) LOGICAL * Specifies whether rounding appears to be done in the IEEE * 'round to nearest' style. * * Further Details * =============== * * The routine is based on the routine ENVRON by Malcolm and * incorporates suggestions by Gentleman and Marovich. See * * Malcolm M. A. (1972) Algorithms to reveal properties of * floating-point arithmetic. Comms. of the ACM, 15, 949-951. * * Gentleman W. M. and Marovich S. B. (1974) More on algorithms * that reveal properties of floating point arithmetic units. * Comms. of the ACM, 17, 276-277. * * * .. Local Scalars .. LOGICAL FIRST, LIEEE1, LRND INTEGER LBETA, LT DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Save statement .. SAVE FIRST, LIEEE1, LBETA, LRND, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ONE = 1 * * LBETA, LIEEE1, LT and LRND are the local values of BETA, * IEEE1, T and RND. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * Compute a = 2.0**m with the smallest positive integer m such * that * * fl( a + 1.0 ) = a. * A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 10 CONTINUE IF( C.EQ.ONE ) THEN A = 2*A C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 10 END IF *+ END WHILE * * Now compute b = 2.0**m with the smallest positive integer m * such that * * fl( a + b ) .gt. a. * B = 1 C = DLAMC3( A, B ) * *+ WHILE( C.EQ.A )LOOP 20 CONTINUE IF( C.EQ.A ) THEN B = 2*B C = DLAMC3( A, B ) GO TO 20 END IF *+ END WHILE * * Now compute the base. a and c are neighbouring floating point * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so * their difference is beta. Adding 0.25 to c is to ensure that it * is truncated to beta and not ( beta - 1 ). * QTR = ONE / 4 SAVEC = C C = DLAMC3( C, -A ) LBETA = C + QTR * * Now determine whether rounding or chopping occurs, by adding a * bit less than beta/2 and a bit more than beta/2 to a. * B = LBETA F = DLAMC3( B / 2, -B / 100 ) C = DLAMC3( F, A ) IF( C.EQ.A ) THEN LRND = .TRUE. ELSE LRND = .FALSE. END IF F = DLAMC3( B / 2, B / 100 ) C = DLAMC3( F, A ) IF( ( LRND ) .AND. ( C.EQ.A ) ) $ LRND = .FALSE. * * Try and decide whether rounding is done in the IEEE 'round to * nearest' style. B/2 is half a unit in the last place of the two * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit * zero, and SAVEC is odd. Thus adding B/2 to A should not change * A, but adding B/2 to SAVEC should change SAVEC. * T1 = DLAMC3( B / 2, A ) T2 = DLAMC3( B / 2, SAVEC ) LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND * * Now find the mantissa, t. It should be the integer part of * log to the base beta of a, however it is safer to determine t * by powering. So we find t as the smallest positive integer for * which * * fl( beta**t + 1.0 ) = 1.0. * LT = 0 A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 30 CONTINUE IF( C.EQ.ONE ) THEN LT = LT + 1 A = A*LBETA C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 30 END IF *+ END WHILE * END IF * BETA = LBETA T = LT RND = LRND IEEE1 = LIEEE1 RETURN * * End of DLAMC1 * END * ************************************************************************ * SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL RND INTEGER BETA, EMAX, EMIN, T DOUBLE PRECISION EPS, RMAX, RMIN * .. * * Purpose * ======= * * DLAMC2 determines the machine parameters specified in its argument * list. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * EPS (output) DOUBLE PRECISION * The smallest positive number such that * * fl( 1.0 - EPS ) .LT. 1.0, * * where fl denotes the computed value. * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow occurs. * * RMIN (output) DOUBLE PRECISION * The smallest normalized number for the machine, given by * BASE**( EMIN - 1 ), where BASE is the floating point value * of BETA. * * EMAX (output) INTEGER * The maximum exponent before overflow occurs. * * RMAX (output) DOUBLE PRECISION * The largest positive number for the machine, given by * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point * value of BETA. * * Further Details * =============== * * The computation of EPS is based on a routine PARANOIA by * W. Kahan of the University of California at Berkeley. * * * .. Local Scalars .. LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, $ NGNMIN, NGPMIN DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, $ SIXTH, SMALL, THIRD, TWO, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. External Subroutines .. EXTERNAL DLAMC1, DLAMC4, DLAMC5 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Save statement .. SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, $ LRMIN, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / , IWARN / .FALSE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ZERO = 0 ONE = 1 TWO = 2 * * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of * BETA, T, RND, EPS, EMIN and RMIN. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. * CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) * * Start to find EPS. * B = LBETA A = B**( -LT ) LEPS = A * * Try some tricks to see whether or not this is the correct EPS. * B = TWO / 3 HALF = ONE / 2 SIXTH = DLAMC3( B, -HALF ) THIRD = DLAMC3( SIXTH, SIXTH ) B = DLAMC3( THIRD, -HALF ) B = DLAMC3( B, SIXTH ) B = ABS( B ) IF( B.LT.LEPS ) $ B = LEPS * LEPS = 1 * *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 10 CONTINUE IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN LEPS = B C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) C = DLAMC3( HALF, -C ) B = DLAMC3( HALF, C ) C = DLAMC3( HALF, -B ) B = DLAMC3( HALF, C ) GO TO 10 END IF *+ END WHILE * IF( A.LT.LEPS ) $ LEPS = A * * Computation of EPS complete. * * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). * Keep dividing A by BETA until (gradual) underflow occurs. This * is detected when we cannot recover the previous A. * RBASE = ONE / LBETA SMALL = ONE DO 20 I = 1, 3 SMALL = DLAMC3( SMALL*RBASE, ZERO ) 20 CONTINUE A = DLAMC3( ONE, SMALL ) CALL DLAMC4( NGPMIN, ONE, LBETA ) CALL DLAMC4( NGNMIN, -ONE, LBETA ) CALL DLAMC4( GPMIN, A, LBETA ) CALL DLAMC4( GNMIN, -A, LBETA ) IEEE = .FALSE. * IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN IF( NGPMIN.EQ.GPMIN ) THEN LEMIN = NGPMIN * ( Non twos-complement machines, no gradual underflow; * e.g., VAX ) ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN LEMIN = NGPMIN - 1 + LT IEEE = .TRUE. * ( Non twos-complement machines, with gradual underflow; * e.g., IEEE standard followers ) ELSE LEMIN = MIN( NGPMIN, GPMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) * ( Twos-complement machines, no gradual underflow; * e.g., CYBER 205 ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. $ ( GPMIN.EQ.GNMIN ) ) THEN IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT * ( Twos-complement machines with gradual underflow; * no known machine ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF *** * Comment out this if block if EMIN is ok IF( IWARN ) THEN FIRST = .TRUE. WRITE( 6, FMT = 9999 )LEMIN END IF *** * * Assume IEEE arithmetic if we found denormalised numbers above, * or if arithmetic seems to round in the IEEE style, determined * in routine DLAMC1. A true IEEE machine should have both things * true; however, faulty machines may have one or the other. * IEEE = IEEE .OR. LIEEE1 * * Compute RMIN by successive division by BETA. We could compute * RMIN as BASE**( EMIN - 1 ), but some machines underflow during * this computation. * LRMIN = 1 DO 30 I = 1, 1 - LEMIN LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) 30 CONTINUE * * Finally, call DLAMC5 to compute EMAX and RMAX. * CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) END IF * BETA = LBETA T = LT RND = LRND EPS = LEPS EMIN = LEMIN RMIN = LRMIN EMAX = LEMAX RMAX = LRMAX * RETURN * 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', $ ' EMIN = ', I8, / $ ' If, after inspection, the value EMIN looks', $ ' acceptable please comment out ', $ / ' the IF block as marked within the code of routine', $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) * * End of DLAMC2 * END * ************************************************************************ * DOUBLE PRECISION FUNCTION DLAMC3( A, B ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. DOUBLE PRECISION A, B * .. * * Purpose * ======= * * DLAMC3 is intended to force A and B to be stored prior to doing * the addition of A and B , for use in situations where optimizers * might hold one of these in a register. * * Arguments * ========= * * A, B (input) DOUBLE PRECISION * The values A and B. * * * .. Executable Statements .. * DLAMC3 = A + B * RETURN * * End of DLAMC3 * END * ************************************************************************ * SUBROUTINE DLAMC4( EMIN, START, BASE ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. INTEGER BASE, EMIN DOUBLE PRECISION START * .. * * Purpose * ======= * * DLAMC4 is a service routine for DLAMC2. * * Arguments * ========= * * EMIN (output) EMIN * The minimum exponent before (gradual) underflow, computed by * setting A = START and dividing by BASE until the previous A * can not be recovered. * * START (input) DOUBLE PRECISION * The starting point for determining EMIN. * * BASE (input) INTEGER * The base of the machine. * * * .. Local Scalars .. INTEGER I DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Executable Statements .. * A = START ONE = 1 RBASE = ONE / BASE ZERO = 0 EMIN = 1 B1 = DLAMC3( A*RBASE, ZERO ) C1 = A C2 = A D1 = A D2 = A *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 10 CONTINUE IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. $ ( D2.EQ.A ) ) THEN EMIN = EMIN - 1 A = B1 B1 = DLAMC3( A / BASE, ZERO ) C1 = DLAMC3( B1*BASE, ZERO ) D1 = ZERO DO 20 I = 1, BASE D1 = D1 + B1 20 CONTINUE B2 = DLAMC3( A*RBASE, ZERO ) C2 = DLAMC3( B2 / RBASE, ZERO ) D2 = ZERO DO 30 I = 1, BASE D2 = D2 + B2 30 CONTINUE GO TO 10 END IF *+ END WHILE * RETURN * * End of DLAMC4 * END * ************************************************************************ * SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. LOGICAL IEEE INTEGER BETA, EMAX, EMIN, P DOUBLE PRECISION RMAX * .. * * Purpose * ======= * * DLAMC5 attempts to compute RMAX, the largest machine floating-point * number, without overflow. It assumes that EMAX + abs(EMIN) sum * approximately to a power of 2. It will fail on machines where this * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, * EMAX = 28718). It will also fail if the value supplied for EMIN is * too large (i.e. too close to zero), probably with overflow. * * Arguments * ========= * * BETA (input) INTEGER * The base of floating-point arithmetic. * * P (input) INTEGER * The number of base BETA digits in the mantissa of a * floating-point value. * * EMIN (input) INTEGER * The minimum exponent before (gradual) underflow. * * IEEE (input) LOGICAL * A logical flag specifying whether or not the arithmetic * system is thought to comply with the IEEE standard. * * EMAX (output) INTEGER * The largest exponent before overflow * * RMAX (output) DOUBLE PRECISION * The largest machine floating-point number. * * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP DOUBLE PRECISION OLDY, RECBAS, Y, Z * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * * First compute LEXP and UEXP, two powers of 2 that bound * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum * approximately to the bound that is closest to abs(EMIN). * (EMAX is the exponent of the required number RMAX). * LEXP = 1 EXBITS = 1 10 CONTINUE TRY = LEXP*2 IF( TRY.LE.( -EMIN ) ) THEN LEXP = TRY EXBITS = EXBITS + 1 GO TO 10 END IF IF( LEXP.EQ.-EMIN ) THEN UEXP = LEXP ELSE UEXP = TRY EXBITS = EXBITS + 1 END IF * * Now -LEXP is less than or equal to EMIN, and -UEXP is greater * than or equal to EMIN. EXBITS is the number of bits needed to * store the exponent. * IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN EXPSUM = 2*LEXP ELSE EXPSUM = 2*UEXP END IF * * EXPSUM is the exponent range, approximately equal to * EMAX - EMIN + 1 . * EMAX = EXPSUM + EMIN - 1 NBITS = 1 + EXBITS + P * * NBITS is the total number of bits needed to store a * floating-point number. * IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN * * Either there are an odd number of bits used to store a * floating-point number, which is unlikely, or some bits are * not used in the representation of numbers, which is possible, * (e.g. Cray machines) or the mantissa has an implicit bit, * (e.g. IEEE machines, Dec Vax machines), which is perhaps the * most likely. We have to assume the last alternative. * If this is true, then we need to reduce EMAX by one because * there must be some way of representing zero in an implicit-bit * system. On machines like Cray, we are reducing EMAX by one * unnecessarily. * EMAX = EMAX - 1 END IF * IF( IEEE ) THEN * * Assume we are on an IEEE machine which reserves one exponent * for infinity and NaN. * EMAX = EMAX - 1 END IF * * Now create RMAX, the largest machine number, which should * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . * * First compute 1.0 - BETA**(-P), being careful that the * result is less than 1.0 . * RECBAS = ONE / BETA Z = BETA - ONE Y = ZERO DO 20 I = 1, P Z = Z*RECBAS IF( Y.LT.ONE ) $ OLDY = Y Y = DLAMC3( Y, Z ) 20 CONTINUE IF( Y.GE.ONE ) $ Y = OLDY * * Now multiply by BETA**EMAX to get RMAX. * DO 30 I = 1, EMAX Y = DLAMC3( Y*BETA, ZERO ) 30 CONTINUE * RMAX = Y RETURN * * End of DLAMC5 * END * * * LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 1.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END CUT HERE..........