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