# 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( D