C ALGORITHM 776, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 23,NO. 4, December, 1997, P. 494--513. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc # Drivers # Src # This archive created: Fri Jul 31 08:49:19 1998 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'README' then echo shar: will not over-write existing file "'README'" else cat << \SHAR_EOF > 'README' Document for SRRIT - A Fortran subroutine to calculate the dominant invariant subspace of a nonsymmetric matrix Z. Bai and G. W. Stewart Addresses of authors: Zhaojun Bai G. W. Stewart Department of Mathematics Department of Computer Science University of Kentucky University of Maryland Lexington, KY 40506 College Park, MD 20742 phone: 606-257-3167 phone: email: bai@ms.uky.edu email: stewart@cs.umd.edu ---------------------------------------------------------------------------- In the file test.f, it includes the following F77 test routines: ............................................................................ TEST test driver routine SCHECK test routine for SLAQR3 ATQxxx matrix-vector multiplication subroutines. They are problem-dependent routines, and in general, supplied by user, but the calling sequence has to be as described in the subroutine. ............................................................................ In the source code file srrit.f, it includes the following F77 subroutines: ............................................................................ SRRIT main routine, computes a nested sequence of orthonormal bases for the dominant invariant subspaces of a real matrix A of order N. SRRSTP performs an Schur-Rayleigh-Ritz iteration step. ORTH orthonormalizes columns of a matrix. RESID computes the each column norm of residual vectors GROUP finds a cluster of a set of complex numbers. COND estimates the infty-norm condition number with respect to inversion of an upper Hessenberg matrix. SLAQR3 computes the Schur factorization of a real upper Hessenberg matrix, the eigenvalues of Schur form appear in descending order of magnitude along its diagonal. This subroutine is a variant of LAPACK subroutine SLAHQR for computing the Schur decomposition. SLARAN generates a random real number from a uniform (0,1) distribution. It is in LAPACK test matrix generator suit. ............................................................................ In addition, aux.f contains the standard BLAS and LAPACK routines used in SRRIT. If BLAS and LAPACK have been installed on your local machines, you may replace aux.f by linking to your local BLAS and LAPACK libraries. On a Unix system, by typing ``make'', it will automatically generate executive files ``subspiter'' and test the single precision routines. The output file is ``soutput''. A real double precision version of SRRIT is also available. The main subroutines names are called DSRRIT, DSRRSP, DORTH, DRESID, DGROUP, DLAQR3, DCOND, and DLARAN. The rest of routines are from real double precision BLAS and LAPACK routines. Using the provided makefile, the double precision version will be automatically tested. The output file is ``doutput''. SHAR_EOF fi # end of overwriting check if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << \SHAR_EOF > 'Makefile' # ################################################################### # # This is the makefile for SRRIT - A FORTRAN subroutine to calculate # the Dominant Invariant Subspace of a Nonsymmetric Matrix developed # by Z. Bai at University of Kentucky, and G. W. Stewart at University # of Maryland. # # If you have installed BLAS and LAPACK on your local machine, you may # replace lapacksub.o and blas.o by your BLAS and LAPACK libraries. # # For any comments, please send to # bai@ms.uky.edu or stewart@cs.umd.edu # # ################################################################### # .SUFFIXES: .f .o FFLAGS = -sloppy -C -d1 -g -temp=/tmp -u F77 = epcf90 .f.o: $(F77) -c $(FFLAGS) $*.f BLAS = LAPACK = # SRRIT = srrit.o sblas.o slapack.o DSRRIT = dsrrit.o dblas.o dlapack.o MAKE = make single double # #all: # ( $(MAKE); \ # subspiter < data; \ # dsubspiter < data ) # single: driver.o $(SRRIT) ;\ $(F77) -g -o subspiter driver.o $(SRRIT) $(LAPACK) $(BLAS) double: driver.o $(DSRRIT) ;\ $(F77) -g -o dsubspiter driver.o $(DSRRIT) $(LAPACK) $(BLAS) clean: rm -f *.o SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'DRES' then echo shar: will not over-write existing file "'DRES'" else cat << \SHAR_EOF > 'DRES' ...... Type of Test problem: RWK ...... The test matrix size = 496 NV and M = 2 6 IT = 0 WR = 9.547D-01 8.114D-02 -4.440D-02 3.793D-02 -1.019D-02 -1.019D-02 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 6.322D-03 -6.322D-03 RSD = 2.534D-01 5.941D-01 5.856D-01 6.041D-01 6.015D-01 6.015D-01 NGRP = 1 CTR = 9.547D-01 AE = 9.547D-01 ARSD = 2.534D-01 NXTSRR = 5 IDORT = 1 IT = 5 WR = 9.969D-01 3.230D-01 3.230D-01 -2.954D-01 -1.914D-01 2.835D-02 WI = 0.000D+00 5.205D-02 -5.205D-02 0.000D+00 0.000D+00 0.000D+00 RSD = 8.032D-02 8.593D-01 8.593D-01 8.675D-01 8.878D-01 9.221D-01 NGRP = 1 CTR = 9.969D-01 AE = 9.969D-01 ARSD = 8.032D-02 NXTSRR = 7 IDORT = 1 IT = 7 WR = 9.985D-01 -4.240D-01 3.563D-01 -2.810D-01 2.145D-01 6.637D-02 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 6.545D-02 8.244D-01 9.003D-01 8.937D-01 8.756D-01 8.881D-01 NGRP = 1 CTR = 9.985D-01 AE = 9.985D-01 ARSD = 6.545D-02 NXTSRR = 10 IDORT = 1 IT = 10 WR = 9.983D-01 4.873D-01 -4.041D-01 -3.772D-01 3.235D-01 1.557D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 6.003D-02 8.040D-01 8.982D-01 8.660D-01 9.193D-01 9.263D-01 NGRP = 1 CTR = 9.983D-01 AE = 9.983D-01 ARSD = 6.003D-02 NXTSRR = 15 IDORT = 2 IT = 15 WR = 9.973D-01 -6.102D-01 5.842D-01 4.607D-01 2.890D-01 -2.379D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 4.633D-02 7.755D-01 7.253D-01 8.774D-01 8.990D-01 9.637D-01 NGRP = 1 CTR = 9.973D-01 AE = 9.973D-01 ARSD = 4.633D-02 NXTSRR = 22 IDORT = 2 IT = 22 WR = 9.981D-01 -8.038D-01 7.285D-01 6.324D-01 3.781D-01 -3.540D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 4.356D-02 5.885D-01 6.990D-01 7.443D-01 9.302D-01 8.976D-01 NGRP = 1 CTR = 9.981D-01 AE = 9.981D-01 ARSD = 4.356D-02 NXTSRR = 33 IDORT = 2 IT = 33 WR = 1.000D+00 -9.398D-01 9.217D-01 8.784D-01 -5.927D-01 1.048D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 2.765D-02 3.097D-01 3.701D-01 4.638D-01 7.542D-01 9.283D-01 NGRP = 1 CTR = 1.000D+00 AE = 1.000D+00 ARSD = 2.765D-02 NXTSRR = 49 IDORT = 1 IT = 49 WR = 1.000D+00 -9.870D-01 9.717D-01 9.000D-01 -7.130D-01 -4.267D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 6.739D-03 1.101D-01 1.273D-01 4.398D-01 5.964D-01 8.724D-01 NGRP = 1 CTR = 1.000D+00 AE = 1.000D+00 ARSD = 6.739D-03 NXTSRR = 73 IDORT = 2 IT = 73 WR = 1.000D+00 -9.917D-01 9.867D-01 -9.631D-01 9.392D-01 -8.748D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 2.695D-03 2.716D-02 3.938D-02 2.733D-01 2.978D-01 4.538D-01 NGRP = 1 CTR = 1.000D+00 AE = 1.000D+00 ARSD = 2.695D-03 NXTSRR = 109 IDORT = 2 IT = 109 WR = 1.000D+00 -9.974D-01 -9.974D-01 9.924D-01 9.712D-01 -9.691D-01 WI = 0.000D+00 8.067D-04 -8.067D-04 0.000D+00 0.000D+00 0.000D+00 RSD = 1.207D-03 3.388D-02 3.388D-02 1.288D-02 1.060D-01 1.298D-01 NGRP = 1 CTR = 1.000D+00 AE = 1.000D+00 ARSD = 1.207D-03 NXTSRR = 163 IDORT = 2 IT = 163 WR = 1.000D+00 -1.000D+00 -9.938D-01 9.934D-01 -9.753D-01 9.752D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 9.891D-05 3.319D-03 2.206D-03 1.268D-03 2.643D-02 2.649D-02 NGRP = 1 CTR = 1.000D+00 AE = 1.000D+00 ARSD = 9.891D-05 NXTSRR = 244 IDORT = 2 IT = 244 WR = 1.000D+00 -1.000D+00 -9.935D-01 9.935D-01 -9.755D-01 9.755D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 1.674D-06 5.403D-05 6.562D-05 3.625D-05 3.097D-03 3.296D-03 NGRP = 1 CTR = 1.000D+00 AE = 1.000D+00 ARSD = 1.674D-06 NXTSRR = 366 IDORT = 2 IT = 366 WR = 1.000D+00 -1.000D+00 -9.935D-01 9.935D-01 -9.755D-01 9.755D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 3.497D-09 1.129D-07 3.047D-07 1.685D-07 1.322D-04 1.421D-04 NGRP = 2 CTR = 1.000D+00 AE = 1.107D-12 ARSD = 7.989D-08 NXTSRR = 549 IDORT = 2 IT = 549 WR = 1.000D+00 -1.000D+00 -9.935D-01 9.935D-01 -9.755D-01 9.755D-01 WI = 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 RSD = 3.345D-13 1.083D-11 9.747D-11 5.276D-11 1.173D-06 1.275D-06 NGRP = 2 CTR = 1.000D+00 AE = 3.886D-16 ARSD = 7.664D-12 NGRP = 2 CTR = 9.935D-01 AE = -8.062D-12 ARSD = 7.837D-11 NGRP = 2 CTR = 9.755D-01 AE = -1.410D-12 ARSD = 1.225D-06 The INFO from SRRIT = 0 Total number of SRRIT iterations = 549 .. Converged Eigenvalues .. WR = 1.000000000000D+00 WI = 0.000000000000D+00 WR = -1.000000000000D+00 WI = 0.000000000000D+00 WR = -9.934621902459D-01 WI = 0.000000000000D+00 WR = 9.934621902298D-01 WI = 0.000000000000D+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 3.034D-14 1.001D-12 9.216D-12 4.924D-12 norm(Q*Q - I) = 1.998401444325D-15 ...... Type of Test problem: LRW ...... The test matrix size = 66 NV and M = 2 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 163 .. Converged Eigenvalues .. WR = -1.000000000000D+00 WI = 0.000000000000D+00 WR = 1.000000000000D+00 WI = 0.000000000000D+00 WR = 9.475580155401D-01 WI = 0.000000000000D+00 WR = -9.475580155401D-01 WI = 0.000000000000D+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 3.053D-16 3.331D-16 1.896D-12 7.013D-13 norm(Q*Q - I) = 1.332267629550D-15 ...... Type of Test problem: DRW ...... The test matrix size = 66 NV and M = 2 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 109 .. Converged Eigenvalues .. WR = -1.000000000001D+00 WI = 0.000000000000D+00 WR = 1.000000000000D+00 WI = 0.000000000000D+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 7.255D-12 1.023D-12 norm(Q*Q - I) = 9.992007221626D-16 ...... Type of Test problem: ODE ...... The test matrix size = 301 NV and M = 4 6 IT = 0 WR = -1.528D-01 3.013D-04 3.013D-04 1.455D-04 1.004D-04 2.828D-05 WI = 0.000D+00 1.064D-04 -1.064D-04 0.000D+00 0.000D+00 0.000D+00 RSD = 1.434D-01 1.094D-02 1.094D-02 2.703D-03 6.455D-03 6.423D-03 NGRP = 1 CTR = 1.528D-01 AE = -1.528D-01 ARSD = 1.434D-01 NXTSRR = 5 IDORT = 1 IT = 5 WR = -1.264D-02 -1.264D-02 4.429D-03 4.429D-03 2.542D-03 2.542D-03 WI = 2.313D-02 -2.313D-02 7.288D-03 -7.288D-03 2.575D-03 -2.575D-03 RSD = 1.093D-07 1.093D-07 2.913D-05 2.913D-05 6.126D-04 6.126D-04 NGRP = 2 CTR = 2.636D-02 AE = -1.264D-02 ARSD = 1.093D-07 NXTSRR = 7 IDORT = 1 IT = 7 WR = -1.264D-02 -1.264D-02 4.446D-03 4.446D-03 2.904D-03 2.904D-03 WI = 2.313D-02 -2.313D-02 7.309D-03 -7.309D-03 2.316D-03 -2.316D-03 RSD = 1.451D-09 1.451D-09 1.988D-06 1.988D-06 2.358D-04 2.358D-04 NGRP = 2 CTR = 2.636D-02 AE = -1.264D-02 ARSD = 1.451D-09 NXTSRR = 9 IDORT = 1 IT = 9 WR = -1.264D-02 -1.264D-02 4.447D-03 4.447D-03 2.886D-03 2.886D-03 WI = 2.313D-02 -2.313D-02 7.308D-03 -7.308D-03 2.242D-03 -2.242D-03 RSD = 1.178D-11 1.178D-11 7.568D-08 7.568D-08 6.098D-05 6.098D-05 NGRP = 2 CTR = 2.636D-02 AE = -1.264D-02 ARSD = 1.178D-11 NXTSRR = 10 IDORT = 1 IT = 10 WR = -1.264D-02 -1.264D-02 4.447D-03 4.447D-03 2.897D-03 2.897D-03 WI = 2.313D-02 -2.313D-02 7.308D-03 -7.308D-03 2.233D-03 -2.233D-03 RSD = 5.713D-13 5.713D-13 2.530D-08 2.530D-08 3.334D-05 3.334D-05 NGRP = 2 CTR = 2.636D-02 AE = -1.264D-02 ARSD = 5.713D-13 NGRP = 2 CTR = 8.555D-03 AE = 4.447D-03 ARSD = 2.530D-08 NXTSRR = 15 IDORT = 1 IT = 15 WR = -1.264D-02 -1.264D-02 4.447D-03 4.447D-03 2.897D-03 2.897D-03 WI = 2.313D-02 -2.313D-02 7.308D-03 -7.308D-03 2.204D-03 -2.204D-03 RSD = 5.713D-13 5.713D-13 1.042D-11 1.042D-11 1.467D-06 1.467D-06 NGRP = 2 CTR = 8.555D-03 AE = 4.447D-03 ARSD = 1.042D-11 NXTSRR = 16 IDORT = 1 IT = 16 WR = -1.264D-02 -1.264D-02 4.447D-03 4.447D-03 2.897D-03 2.897D-03 WI = 2.313D-02 -2.313D-02 7.308D-03 -7.308D-03 2.204D-03 -2.204D-03 RSD = 5.713D-13 5.713D-13 3.569D-12 3.569D-12 9.619D-07 9.619D-07 NGRP = 2 CTR = 8.555D-03 AE = 4.447D-03 ARSD = 3.569D-12 NXTSRR = 17 IDORT = 1 IT = 17 WR = -1.264D-02 -1.264D-02 4.447D-03 4.447D-03 2.897D-03 2.897D-03 WI = 2.313D-02 -2.313D-02 7.308D-03 -7.308D-03 2.204D-03 -2.204D-03 RSD = 5.713D-13 5.713D-13 7.349D-13 7.349D-13 5.041D-07 5.041D-07 NGRP = 2 CTR = 8.555D-03 AE = 4.447D-03 ARSD = 7.349D-13 NGRP = 2 CTR = 3.640D-03 AE = 2.897D-03 ARSD = 5.041D-07 The INFO from SRRIT = 0 Total number of SRRIT iterations = 17 .. Converged Eigenvalues .. WR = -1.264346493143D-02 WI = 2.312526126815D-02 WR = -1.264346493143D-02 WI = -2.312526126815D-02 WR = 4.446825337358D-03 WI = 7.308366846422D-03 WR = 4.446825337358D-03 WI = -7.308366846422D-03 Final residual norms of (A*Q - Q*T) of converged RR pairs = 7.651D-14 2.445D-14 4.922D-14 9.934D-14 norm(Q*Q - I) = 7.771561172376D-16 ...... Type of Test problem: PDE ...... The test matrix size = 961 NV and M = 3 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 2085 .. Converged Eigenvalues .. WR = 7.977818149246D+00 WI = 0.000000000000D+00 WR = 7.949033322109D+00 WI = 0.000000000000D+00 WR = 7.949033322098D+00 WI = 0.000000000000D+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 4.295D-12 5.855D-12 4.953D-12 norm(Q*Q - I) = 7.771561172376D-16 ...... Type of Test problem: PDE ...... The test matrix size = 225 NV and M = 1 1 The INFO from SRRIT = 0 Total number of SRRIT iterations = 1513 .. Converged Eigenvalues .. WR = 7.911564989140D+00 WI = 0.000000000000D+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 8.987D-12 norm(Q*Q - I) = 4.440892098501D-16 ...... Type of Test problem: BWM ...... The test matrix size = 200 NV and M = 2 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 2000 .. Converged Eigenvalues .. WR = -1.235506919564D+03 WI = 0.000000000000D+00 WR = -1.234607256326D+03 WI = 0.000000000000D+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 3.596D-09 1.335D-08 norm(Q*Q - I) = 8.881784197001D-16 ...... Type of Test problem: N44 ...... The test matrix size = 4 NV and M = 1 2 DORTH subroutine fails, INFO = 1 ...... Type of Test problem: F44 ...... The test matrix size = 4 NV and M = 1 2 Reduced matrix T is singular, INFO = 3 ...... Type of Test problem: M44 ...... The test matrix size = 4 NV and M = 1 2 DSRRIT fails to converge after MAXIT, INF0 = 4 2-norms of Residual Vectors 6.594D-03 6.594D-03 Iteration numbers where the res. are computed 40 40 WR = 2.063278042415D+00 WI = 6.130769118509D-02 WR = 2.063278042415D+00 WI = -6.130769118509D-02 ...... Type of Test problem: SUB ...... Passed all DLAQR3 tests ..... END OF DSRRIT TEST ..... SHAR_EOF fi # end of overwriting check if test -f 'data' then echo shar: will not over-write existing file "'data'" else cat << \SHAR_EOF > 'data' RWK # standard test (istart = -1) 30 2 6 LRW # initial vectors as input (istart = 0) 10 2 6 DRW # initial vectors as input (istart = 1) 10 2 6 ODE # standard test 300 4 6 PDE # standard test 31 3 6 PDE # special case NV = M = 1 15 1 1 BWM # standard test 200 2 6 N44 # initial zero vectors (info = 1) 4 1 2 F44 # will make T to be singular, and SRRIT fails. 4 1 2 M44 # method fails convergence after MAXIT iteration 4 1 2 SUB # test SGEHD2, SORGN2, SLAQR3 10 # number of dimensions 1 2 4 9 10 20 25 30 35 50 # MAX different dimension of test matrices END # end of all tests SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << \SHAR_EOF > 'driver.f' PROGRAM DTEST * * ====================================================================== * * Test routine of SRRIT - a subroutine to calculate the dominant * invariant subspace of a nonsymmetric matrix. The SRRIT user guide is * documented in * * Z. Bai and G. W. Stewart, SRRIT - A FORTRAN Subroutine to * calculate the dominant invariant subspace of a nonsymmetric * matrix, submitted to ACM TOMS. * * With the file ``input'' as input, this test driver will exercise * major features of SRRIT. The output is going to be in the file * named ``output''. * * The current array sizes for testing SRRIT are set to handle the order * of a test matrix less than 1000, and the dimension of subspace less * than 100. For running a larger test or a larger subspace, user needs * to change the dimension parameters LDA and MAXSPC. * In addition, see the comments of DCHECK for the limit of the test * problem sizes for the testing subroutine DLAQR3. * * ---------------- * * Copyright (C) 1994, Zhaojun Bai and G. W. Stewart * All rights reserved. Use at your own risk. * Please send comments to bai@ms.uky.edu or stewart@cs.umd.edu * * ====================================================================== * * .. Parameters .. INTEGER LDA, LDQ, MAXSPC, LDT PARAMETER ( LDA = 1000, LDQ = LDA, MAXSPC = 100, $ LDT = MAXSPC ) INTEGER LWORK PARAMETER ( LWORK = MAXSPC*MAXSPC+5*MAXSPC ) DOUBLE PRECISION ZERO, HALF, ONE PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) * .. * .. Local Scalars .. CHARACTER*3 PTYPE INTEGER I, INFO, ISTART, J, K, M, MAXIT, N, NV DOUBLE PRECISION QNORM, TOL * .. * .. Local Arrays .. INTEGER ITRSD( MAXSPC ), IWORK( LWORK ), ISEED( 4 ), $ NMAX( MAXSPC ) DOUBLE PRECISION AQ( LDA, MAXSPC ), DUMMY( 1 ), $ Q( LDQ, MAXSPC ), RSD( MAXSPC ), T( LDT, LDT ), $ WI( LDT ), WORK( LWORK ), WR( LDT ) * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLANGE, DLARAN EXTERNAL LSAME, DLANGE, DLARAN * .. * .. External Subroutines .. EXTERNAL DGEMM, DLASET, DSRRIT EXTERNAL ATQODE, ATQPDE, ATQRWK, ATQBWM, ATQ4X4, ATQF44 * * .. Intrinsic functions .. INTRINSIC SIGN * .. * .. Executable Statements .. * OPEN( NOUT, FILE = 'doutput' ) * 10 CONTINUE * * Input the Type of Test Problem * READ( NIN, FMT = * )PTYPE * IF( .NOT.LSAME( PTYPE, 'END' ) ) $ WRITE( NOUT, FMT = 9999 )PTYPE * IF( LSAME( PTYPE, 'END' ) ) THEN WRITE( NOUT, FMT = 9998 ) GO TO 80 END IF * * Input Parameter K, which decides the size of test problem * READ( NIN, FMT = * )K * IF( LSAME( PTYPE, 'RWK' ) .OR. LSAME( PTYPE, 'LRW' ) .OR. $ LSAME( PTYPE, 'DWK' ) ) THEN N = ( K+1 )*( K+2 ) / 2 ELSE IF( LSAME( PTYPE, 'ODE' ) ) THEN N = K + 1 ELSE IF( LSAME( PTYPE, 'PDE' ) ) THEN N = K*K ELSE IF( LSAME( PTYPE, 'BWM' ) ) THEN N = K ELSE IF( LSAME( PTYPE, 'M44' ) .OR. LSAME( PTYPE, 'N44' ) ) THEN N = K ELSE IF( LSAME( PTYPE, 'F44' ) ) THEN N = K ELSE IF( LSAME( PTYPE, 'SUB' ) ) THEN GO TO 90 END IF * WRITE( NOUT, FMT = 9996 )N * IF( N.GT.LDA ) THEN WRITE( NOUT, FMT = 9995 ) GO TO 80 END IF * * Input NV, the number of the desired dominant eigenvalues. * M, the subspace size * Note that it is required that M >= NV+2. * READ( NIN, FMT = * )NV, M * WRITE( NOUT, FMT = 1111 )NV, M 1111 FORMAT( ' NV and M = ', 2I5) * * Set stopping criterion TOL and maximal number of iterations. * They may be altered if desired. * TOL = 1D-10 MAXIT = 10*N * * Call main SRRIT routine * IF( LSAME( PTYPE, 'RWK' ) ) THEN INFO = 1 ISTART = -1 CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQRWK ) ELSE IF( LSAME( PTYPE, 'LRW' ) ) THEN INFO = 0 ISTART = 0 * * generate initial vectors (istart = 0) * ISEED( 1 ) = 3192 ISEED( 2 ) = 1221 ISEED( 3 ) = 1322 ISEED( 4 ) = 8097 DO 30 J = 1, M DO 20 I = 1, N Q( I,J ) = SIGN( ONE, DLARAN( ISEED ) - HALF ) 20 CONTINUE 30 CONTINUE * CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQRWK ) ELSE IF( LSAME( PTYPE, 'DRW' ) ) THEN INFO = 0 ISTART = 1 * * generate and orthogonalize initial vectors (istart = 1) * ISEED( 1 ) = 3192 ISEED( 2 ) = 1221 ISEED( 3 ) = 1322 ISEED( 4 ) = 8097 DO 50 J = 1, M DO 40 I = 1, N Q( I,J ) = DLARAN( ISEED ) 40 CONTINUE 50 CONTINUE CALL DORTH( N, 1, M, Q, LDQ, INFO ) IF( INFO.NE.0 )THEN WRITE( NOUT, 1112 )INFO 1112 FORMAT( 'Given starting vectors are linearly dependent' ) GO TO 80 END IF * CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQRWK ) ELSE IF( LSAME( PTYPE, 'ODE' ) ) THEN INFO = 1 ISTART = -1 CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQODE ) ELSE IF( LSAME( PTYPE, 'PDE' ) ) THEN INFO = 0 ISTART = -1 CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQPDE ) ELSE IF( LSAME( PTYPE, 'BWM' ) ) THEN INFO = 0 ISTART = -1 CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQBWM ) ELSE IF( LSAME( PTYPE, 'M44' ) ) THEN INFO = 0 ISTART = -1 CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQ4X4 ) ELSE IF( LSAME( PTYPE, 'N44' ) ) THEN INFO = 0 ISTART = 0 DO 55 J = 1, M DO 45 I = 1, N Q( I,J ) = ONE 45 CONTINUE 55 CONTINUE CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQ4X4 ) ELSE IF( LSAME( PTYPE, 'F44' ) ) THEN INFO = 0 ISTART = 1 DO 75 J = 1, M DO 65 I = 1, N Q( I,J ) = ZERO 65 CONTINUE Q( J, J ) = ONE 75 CONTINUE CALL DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQF44 ) END IF * IF( INFO.EQ.1 )THEN WRITE( NOUT, FMT = 1113 )INFO 1113 FORMAT( ' DORTH subroutine fails, INFO = ', I2 ) ELSE IF( INFO.EQ.2 )THEN WRITE( NOUT, FMT = 1114 )INFO 1114 FORMAT( ' DSRRSP subroutine fails, INFO = ', I2 ) ELSE IF( INFO.EQ.3 )THEN WRITE( NOUT, FMT = 1115 )INFO 1115 FORMAT( ' Reduced matrix T is singular, INFO = ', I2 ) ELSE IF( INFO.EQ.4 )THEN WRITE( NOUT, FMT = 1116 )INFO 1116 FORMAT( ' DSRRIT fails to converge after MAXIT, INF0 = ', I2) WRITE( NOUT, FMT = 1118 ) 1118 FORMAT( ' 2-norms of Residual Vectors ' ) WRITE( NOUT, FMT = 9989 )( RSD( I ), I = 1, M ) WRITE( NOUT, FMT = 1119 ) 1119 FORMAT( ' Iteration numbers where the res. are computed ' ) WRITE( NOUT, FMT = 9985 )( ITRSD( I ), I = 1, M ) DO 61 K = 1, M WRITE( NOUT, FMT = 9991 )WR( K ), WI( K ) 61 CONTINUE ELSE WRITE( NOUT, FMT = 9993 )INFO WRITE( NOUT, FMT = 1117 )MAXIT 1117 FORMAT( ' Total number of SRRIT iterations = ', I5 ) END IF * IF( INFO.EQ.0 .AND. NV.GT.0 )THEN WRITE( NOUT, FMT = 9992 ) DO 60 K = 1, NV WRITE( NOUT, FMT = 9991 )WR( K ), WI( K ) 60 CONTINUE * * Test residual: R = A*Q - Q*T * CALL DGEMM( 'No transpose', 'No transpose', N, NV, NV, -ONE, $ Q, LDQ, T, LDT, ONE, AQ, LDA ) * DO 70 I = 1, NV WORK( I ) = DLANGE( 'MAX', N, 1, AQ( 1, I ), 1, DUMMY ) 70 CONTINUE * WRITE( NOUT, FMT = 9990 ) WRITE( NOUT, FMT = 9989 )( WORK( I ), I = 1, NV ) * CALL DLASET( 'Full', NV, NV, ZERO, ONE, WORK, NV ) CALL DGEMM( 'Transpose', 'No transpose', NV, NV, N, -ONE, Q, $ LDQ, Q, LDQ, ONE, WORK, NV ) QNORM = DLANGE( 'Max', NV, NV, WORK, NV, DUMMY ) * WRITE( NOUT, FMT = 9988 )QNORM * END IF * GO TO 10 * 9999 FORMAT( / 10X, ' ...... Type of Test problem: ', A3 , ' ......' ) 9998 FORMAT( / ' ..... END OF DSRRIT TEST .....' ) 9996 FORMAT( / ' The test matrix size = ', I5 ) 9995 FORMAT( ' N is larger than array leading dimension LDA' ) 9993 FORMAT( / ' The INFO from SRRIT = ', I2 ) 9992 FORMAT( ' .. Converged Eigenvalues ..' ) 9991 FORMAT( ' WR = ', 1P,D20.12, ' WI = ', 1P,D20.12 ) 9990 FORMAT( ' Final residual norms ', $ 'of (A*Q - Q*T) of converged RR pairs = ' ) 9989 FORMAT( 1P,6D10.3 ) 9985 FORMAT( 6I10 ) 9988 FORMAT( ' norm(Q*Q - I) = ', 1P,D20.12 ) * * Test subroutine DLAQR3 * 90 CONTINUE * READ( NIN, FMT = * )( NMAX( I ), I = 1, K ) * CALL DCHECK( K, NMAX, INFO ) IF( INFO .EQ. 0 )THEN WRITE( NOUT, FMT = 9987 ) 9987 FORMAT( ' Passed all DLAQR3 tests ' ) ELSE WRITE( NOUT, FMT = 9986 )INFO 9986 FORMAT( ' Test routine DCHECK INFO = ', I5 ) END IF * GO TO 10 * 80 CONTINUE * CLOSE ( NOUT ) * END * ****** begin of dcheck.f ******************************************** * SUBROUTINE DCHECK( MAX, NMAX, INFO ) * * .. Scalar Arguments .. INTEGER MAX, INFO * * .. Arrary Arguments .. INTEGER NMAX( * ) * * Purpose * ======= * * DCHECK tests the subroutine DLAQR3. * On normal return, INFO is to be zero, otherwise, if * * INFO = 1, DLAQR3 has illegal parameter in the calling sequence * or it fails to compute the complete Schur decomposition. * = 2, the backward errors of the Schur decomposition fails * test. * = 3, the order of the computed eigenvalues has not been * sorted correctly with option WANTT = .TRUE. and * WANTZ = .TRUE. in DLAQR3 * = 4, the order of the computed eigenvalues has not been * sorted correctly with option WANTT = .FALSE. and * WANTZ = .FALSE. in DLAQR3 * * Note: the internal parameter LDA is set so that the order of any * test matrix is smaller than 50. * * =================================================================== * * .. Parameters .. INTEGER LDA, LDB, LDH, LDU, LWORK PARAMETER ( LDA = 50, LDB = LDA, LDH = LDA, LDU = LDA, $ LWORK = 2*LDA*LDA ) * DOUBLE PRECISION ZERO, TOL PARAMETER ( ZERO = 0.0D+0, TOL = 2.0D+01 ) * * .. Local Scalars .. INTEGER I, J, K, N, IERR * * .. Local Arrays .. INTEGER ISEED( 4 ) DOUBLE PRECISION A( LDA,LDA ), B( LDB,LDB ), U( LDU,LDU ), $ H( LDH, LDH ), WR( LDA ), WI(LDA), $ WORK( LWORK ), RESULT( 2 ) * * .. External functions .. DOUBLE PRECISION DLAPY2, DLARAN EXTERNAL DLAPY2, DLARAN * * .. External Subroutines .. EXTERNAL DLACPY, DGEHRD, DORGHR, DLAQR3, DGET21 * INFO = 0 ISEED( 1 ) = 8321 ISEED( 2 ) = 9735 ISEED( 3 ) = 8673 ISEED( 4 ) = 2437 * DO 10 K = 1, MAX N = NMAX( K ) * * Generate the random test matrices * DO 30 J = 1,N DO 20 I = 1, N A( I,J ) = DLARAN( ISEED ) 20 CONTINUE 30 CONTINUE * * Save a copy of the matrix * CALL DLACPY( 'FULL', N, N, A, LDA, B, LDB ) * * Hessenberg reduction * CALL DGEHRD( N, 1, N, A, LDA, WR, WORK, LWORK, IERR ) DO 50 J = 1, N - 1 DO 40 I = J + 2, N U( I, J ) = A( I, J ) A( I, J ) = ZERO 40 CONTINUE 50 CONTINUE CALL DORGHR( N, 1, N, U, LDU, WR, WORK, LWORK, IERR ) * CALL DLACPY( 'FULL', N, N, A, LDA, H, LDH ) * * Compute the Schur decomposition of upper Hessenberg matrix * CALL DLAQR3( .TRUE., .TRUE., N, 1, N, A, LDA, WR, WI, 1, N, $ U, LDU, WORK, IERR ) IF( IERR .NE. 0 )THEN INFO = 1 RETURN END IF * CALL DGET21( N, B, LDB, A, LDA, U, LDU, WORK, RESULT ) IF( RESULT( 1 ) .GT. TOL .OR. RESULT( 2 ) .GT. TOL )THEN INFO = 2 RETURN END IF * * Check the ordering * DO 60 I = 1,N WORK( I ) = DLAPY2( WR( I ), WI( I ) ) 60 CONTINUE DO 80 I = N-1,-1,1 DO 70 J = I+1,-1,1 IF( WORK(I) .LT. WORK(J) )THEN INFO = 3 RETURN END IF 70 CONTINUE 80 CONTINUE * * Compute the eigenvalues only. * CALL DLAQR3( .FALSE., .FALSE., N, 1, N, H, LDH, WR, WI, $ 1, N, U, LDU, WORK, IERR ) * * Check the ordering * DO 90 I = 1,N WORK( I ) = DLAPY2( WR( I ), WI( I ) ) 90 CONTINUE DO 110 I = N-1,-1,1 DO 100 J = I+1,-1,1 IF( WORK(I) .LT. WORK(J) )THEN INFO = 4 RETURN END IF 100 CONTINUE 110 CONTINUE 10 CONTINUE * RETURN END * ********* begin of dget21.f ***************************************** * SUBROUTINE DGET21( N, A, LDA, B, LDB, U, LDU, WORK, RESULT ) * .. * .. Scalar Arguments .. INTEGER LDA, LDB, LDU, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), RESULT( * ), $ U( LDU, * ), WORK( * ) * * Purpose * ======= * * SGET21 generally checks a decomposition of the form * * A = U B U' * * where ' means transpose and U is orthogonal. Specifically, * * RESULT(1) = | A - U B U' | / ( |A| n ulp ) * and * RESULT(2) = | I - UU' | / ( n ulp ) * * Arguments * ========== * * N - INTEGER * The size of the matrix. If it is zero, SGET21 does nothing. * It must be at least zero. * Not modified. * * A - DOUBLE PRECISION array of dimension ( LDA , N ) * The original (unfactored) matrix. * Not modified. * * LDA - INTEGER * The leading dimension of A. It must be at least 1 * and at least N. * Not modified. * * B - DOUBLE PRECISION array of dimension ( LDB , N ) * The factored matrix. * Not modified. * * LDB - INTEGER * The leading dimension of B. It must be at least 1 * and at least N. * Not modified. * * U - DOUBLE PRECISION array of dimension ( LDU, N ). * The orthogonal matrix in the decomposition. * Not modified. * * LDU - INTEGER * The leading dimension of U. LDU must be at least N and * at least 1. * Not modified. * * WORK - DOUBLE PRECISION array of dimension ( 2*N**2 ) * Workspace. * Modified. * * RESULT - DOUBLE PRECISION array of dimension ( 2 ) * The values computed by the two tests described above. The * values are currently limited to 1/ulp, to avoid overflow. * Errors are flagged by RESULT(1)=10/ulp. * *----------------------------------------------------------------------- * * .. Parameters .. * DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. * INTEGER JDIAG DOUBLE PRECISION ANORM, ULP, UNFL, WNORM * .. * .. External Functions .. * DOUBLE PRECISION DLAMCH, DLANGE EXTERNAL DLAMCH, DLANGE * .. * .. External Subroutines .. * EXTERNAL DGEMM, DLACPY * .. * .. Intrinsic Functions .. * INTRINSIC MAX, MIN, DBLE * .. * .. Executable Statements .. * RESULT( 1 ) = ZERO RESULT( 2 ) = ZERO IF( N.LE.0 ) $ RETURN * * Constants * UNFL = DLAMCH( 'Safe minimum' ) ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) * * Do Test 1 * * Norm of A: * ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), UNFL ) * * Norm of A - UBU' * CALL DLACPY( ' ', N, N, A, LDA, WORK, N ) CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO, $ WORK( N**2+1 ), N ) * CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, U, LDU, $ ONE, WORK, N ) * WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ) * IF( ANORM.GT.WNORM ) THEN RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) ELSE IF( ANORM.LT.ONE ) THEN RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) ELSE RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) END IF END IF * * Do Test 2 * * Compute UU' - I * CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, $ N ) * DO 40 JDIAG = 1, N WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+ $ 1 ) - ONE 40 CONTINUE * RESULT( 2 ) = MIN( DLANGE( '1', N, N, WORK, N, $ WORK( N**2+1 ) ), DBLE( N ) ) / ( N*ULP ) * RETURN * * End of DGET21 * END * * ****** begin of atqrwk.f ********************************************* * SUBROUTINE ATQRWK( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ) * .. * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is for the random walk example in * the manuscript ``SRRIT - A Fortran Subroutine to Calculate the * Dominant Invariant subspace of a Nonsymmetric matrix'' by Z. Bai * and G. W. Stewart. * * The order of the random walk test matrix is * * N = (K+1)*(K+2)/2 * * where K is the size of triangular random walk grid. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) DOUBLE PRECISION array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) DOUBLE PRECISION array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, HALF PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER I, J, K, KK, LK, NN DOUBLE PRECISION UPPER * .. * .. Intrinsic Functions .. INTRINSIC DBLE, SQRT * .. * .. Statement Functions .. INTEGER IINDEX * .. * .. Statement Function definitions .. IINDEX( K, I, J ) = I*K - ( ( I-1 )*( I-2 ) ) / 2 + 1 + J + 1 * .. * .. Executable Statements .. * NN = ( -3+SQRT( DBLE( 8*N+1 ) ) ) / 2 DO 30 KK = L, M * * Compute the vector AQ(1:N,KK) = A*Q(1:N,KK) * DO 20 I = 0, NN DO 10 J = 0, NN - I * K = IINDEX( NN, I, J ) UPPER = DBLE( I+J ) / DBLE( 2*NN ) AQ( K, KK ) = ZERO * IF( I.EQ.0 .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * ELSE IF( I.EQ.0 .AND. J.LE.NN-1 ) THEN LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * ELSE IF( I.EQ.0 .AND. J.EQ.NN ) THEN LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) * ELSE IF( I.LE.NN-1 .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * ELSE IF( I.EQ.NN .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) * ELSE IF( ( I+J ).EQ.NN ) THEN LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) * ELSE LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * END IF 10 CONTINUE 20 CONTINUE * 30 CONTINUE * RETURN * * End of ATQRWK END * ****** begin of atqpde.f ********************************************** SUBROUTINE ATQPDE( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ) * .. * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is for the convection-diffusion * PDE example in the manuscript ``SRRIT - A Fortran Subroutine to * Calculate the Dominant Invariant subspace of a Nonsymmetric * matrix'' by Z. Bai and G. W. Stewart. * * The convection-diffusion operator is discretized on a K x K square * grid, and resulted a PDE test matrix of order * * N = K^2. * * The constant parameters p1, p2, and p3 in the convection-diffusion * operator can changed in this subroutine. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) DOUBLE PRECISION array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) DOUBLE PRECISION array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, FOUR PARAMETER ( ONE = 1.0D+0, FOUR = 4.0D+0 ) DOUBLE PRECISION P1, P2, P3 PARAMETER ( P1 = 1.0D+0, P2 = 1.0D+0, P3 = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, II, K, NS DOUBLE PRECISION BETA, GAMMA, HSTEP, SIGMA * .. * .. Intrinsic Functions .. INTRINSIC DBLE, SQRT * .. * .. Executable Statements .. * NS = SQRT( DBLE( N ) ) HSTEP = ONE / ( DBLE( NS )+ONE ) SIGMA = P3*HSTEP*HSTEP GAMMA = P2*HSTEP BETA = P1*HSTEP * DO 70 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * DO 60 I = 1, NS * DO 10 II = ( I-1 )*NS + 1, I*NS AQ( II, K ) = ( FOUR-SIGMA )*Q( II, K ) 10 CONTINUE * DO 20 II = ( I-1 )*NS + 1, I*NS - 1 AQ( II, K ) = AQ( II, K ) + ( GAMMA-ONE )*Q( II+1, K ) 20 CONTINUE * DO 30 II = ( I-1 )*NS + 1 + 1, I*NS AQ( II, K ) = AQ( II, K ) + ( -GAMMA-ONE )*Q( II-1, K ) 30 CONTINUE * IF( I.GT.1 ) THEN DO 40 II = ( I-1 )*NS + 1, I*NS AQ( II, K ) = AQ( II, K ) - ( BETA+ONE )*Q( II-NS, K ) 40 CONTINUE END IF IF( I.LT.NS ) THEN DO 50 II = ( I-1 )*NS + 1, I*NS AQ( II, K ) = AQ( II, K ) + ( BETA-ONE )*Q( II+NS, K ) 50 CONTINUE END IF * 60 CONTINUE * 70 CONTINUE * RETURN * * End of ATQPDE END * ********* begin of atqode.f ******************************************* SUBROUTINE ATQODE( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ) * .. * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is for the ODE boundary value * problem example in the manuscript ``SRRIT - A Fortran Subroutine * to Calculate the Dominant Invariant subspace of a Nonsymmetric * matrix'' by Z. Bai and G. W. Stewart. * * The order of the ODE test matrix is * * N = K + 1, * * where K is the number of subintervals of [0,1]. * * NOTE: N should be smaller than 1001, if one wants to test a larger * size of N, the size of working arrays A and U in this subroutine * should be increased. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) DOUBLE PRECISION array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) DOUBLE PRECISION array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) DOUBLE PRECISION THREE, FOUR PARAMETER ( THREE = 3.0D+0, FOUR = 4.0D+0 ) DOUBLE PRECISION GAMMA PARAMETER ( GAMMA = 1.0D-2 ) INTEGER LDU PARAMETER ( LDU = 1001 ) * .. * .. Local Scalars .. INTEGER I, INFO, J, K DOUBLE PRECISION DD, DELTA, HSTEP, HSTEP2 * .. * .. Local Arrays .. DOUBLE PRECISION A( 2, LDU ), U( LDU, 1 ) * .. * .. Intrinsic Functions .. INTRINSIC DBLE * .. * .. External Subroutines .. EXTERNAL DPBTF2, DPBTRS * .. * .. Executable Statements .. * HSTEP = ONE / DBLE( N ) HSTEP2 = HSTEP*HSTEP * A( 1, 1 ) = ZERO DO 10 J = 2, N - 1 A( 1, J ) = -ONE 10 CONTINUE * DO 20 J = 1, N - 1 A( 2, J ) = TWO 20 CONTINUE * * Cholesky decomposition of Symmetric tridiagonal matrix * CALL DPBTF2( 'Upper', N-1, 1, A, 2, INFO ) * DO 30 J = 1, N - 2 U( J, 1 ) = ZERO 30 CONTINUE U( N-1, 1 ) = ONE * CALL DPBTRS( 'Upper', N-1, 1, 1, A, 2, U, 301, INFO ) * DELTA = THREE*GAMMA + ( FOUR*U( 1, 1 )-U( 2, 1 )+GAMMA* $ U( N-2, 1 )-FOUR*GAMMA*U( N-1, 1 ) ) * DO 60 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * DO 40 I = 1, N - 1 AQ( I, K ) = HSTEP2*Q( I, K ) 40 CONTINUE AQ( N, K ) = ZERO * CALL DPBTRS( 'Upper', N-1, 1, 1, A, 2, AQ( 1, K ), LDA, INFO ) * DD = FOUR*AQ( 1, K ) - AQ( 2, K ) + GAMMA*AQ( N-2, K ) - $ FOUR*GAMMA*AQ( N-1, K ) AQ( N, K ) = DD / DELTA * DO 50 I = 1, N - 1 AQ( I, K ) = -AQ( I, K ) + U( I, 1 )*AQ( N, K ) 50 CONTINUE * 60 CONTINUE * DO 80 K = L, M DO 70 I = 1, N AQ( I, K ) = -AQ( I, K ) 70 CONTINUE 80 CONTINUE * RETURN * * End of ATQODE END ********* begin of atqbwm.f ****************************************** SUBROUTINE ATQBWM( N, L, M, X, LDX, Y, LDY ) * .. * .. Scalar Arguments .. INTEGER LDY, LDX, L, M, N * .. * .. Array Arguments .. DOUBLE PRECISION Y( LDY, * ), X( LDX, * ) * .. * * Purpose * ======= * * Compute * * Y(:,L:M) = A*X(:,L:M) * * The matrix A is a 2 by 2 block matrix resulted from the finite * difference discretization of the Brusselator wave model. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N has to be an even number. * * L,M (input) INTEGERS * The number of the first and the last column of Q to * multiply by A. * * X (input) DOUBLE PRECISION array, dimension ( LDX, M ) * contains the vectors X. * * LDX (input) INTEGER * The leading dimension of the array X, LDX >= max( 1, N ) * * Y (output) DOUBLE PRECISION array, dimension (LDY, M ) * contains the product of the matrix op(A) with X. * * LDY (input) INTEGER * The leading dimension of the array Y, LDY >= max( 1, N ) * * =================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, TWO PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0 ) DOUBLE PRECISION ALPHA, BETA, DELT1, DELT2, LL PARAMETER ( ALPHA = 2.0D+0, BETA = 5.45D+0, $ DELT1 = 8.0D-3, DELT2 = 4.0D-3, $ LL = 0.51302D+0 ) * * .. Local Variables .. INTEGER I, K, NHALF DOUBLE PRECISION HSTEP, TAU1, TAU2, DIAG1, DIAG2, ALPHA2 * * set up middle parameters * NHALF = N / 2 HSTEP = ONE / ( DBLE( NHALF ) + 1 ) TAU1 = DELT1 / (HSTEP*LL)**2 TAU2 = DELT2 / (HSTEP*LL)**2 DIAG1 = -TWO*TAU1 + BETA - ONE ALPHA2 = ALPHA**2 DIAG2 = -TWO*TAU2 - ALPHA2 * * Compute Y(:,L:M) = A*X(:,L:M) * DO 30 K = L, M * * the 1st row * Y( 1, K ) = DIAG1*X( 1, K ) + TAU1*X( 2, K ) + $ ALPHA2*X( NHALF+1, K ) * * the 2nd row to (NHALF-1)th row * DO 10 I = 2, NHALF - 1 Y( I, K ) = TAU1*X( I-1, K ) + DIAG1*X( I, K ) + $ TAU1*X( I+1, K ) + ALPHA2*X( NHALF+I, K ) 10 CONTINUE * * the NHALFth row * Y( NHALF, K ) = TAU1*X( NHALF-1, K ) + DIAG1*X( NHALF, K ) $ + ALPHA2*X( N,K ) * * the (NHALF+1)th row * Y( NHALF+1, K ) = -BETA*X( 1, K ) + DIAG2*X( NHALF+1, K ) $ + TAU2*X( NHALF+2, K ) * * the (NHALF+2)th row to (N-1)th row * DO 20 I = NHALF+2, N - 1 Y( I, K ) = -BETA*X( I-NHALF, K ) + TAU2*X( I-1, K ) + $ DIAG2*X( I, K ) + TAU2*X( I+1, K ) 20 CONTINUE * * the last (Nth) row * Y( N, K ) = -BETA*X( N-NHALF, K ) + TAU2*X( N-1,K ) + $ DIAG2*X( N, K ) * 30 CONTINUE * RETURN * * End of ATQBWM * END ********* begin of atq4x4.f ****************************************** SUBROUTINE ATQ4X4( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ) * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is a 3 by 3 matrix: * * A = [ 2 1 1 1 ] * [ 0 2 1 1 ] * [ 0 0 2 1 ] * [ 0 0 1 ] * * which has exact eigenvalues 2, 2, 2, 1. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) DOUBLE PRECISION array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) DOUBLE PRECISION array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. DOUBLE PRECISION TWO PARAMETER ( TWO = 2.0D+0 ) * * .. * .. Local Scalars .. INTEGER K * N = 4 DO 10 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * AQ( 1, K ) = TWO*Q( 1, K ) + Q( 2, K ) + Q( 3, K ) + Q( 4, K ) AQ( 2, K ) = TWO*Q( 2, K ) + Q( 3, K ) + Q( 4, K ) AQ( 3, K ) = TWO*Q( 3, K ) + Q( 4, K ) AQ( 4, K ) = Q( 4, K ) * 10 CONTINUE * RETURN END ********* begin of atqf44.f ****************************************** SUBROUTINE ATQF44( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ) * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is a 3 by 3 matrix: * * A = [ 2 0 1 2 ] * [ 0 0 1 1 ] * [ 0 0 2 1 ] * [ 0 0 1 1 ] * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) DOUBLE PRECISION array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) DOUBLE PRECISION array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. DOUBLE PRECISION TWO PARAMETER ( TWO = 2.0D+0 ) * * .. * .. Local Scalars .. INTEGER K * N = 4 DO 10 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * AQ( 1, K ) = TWO*Q( 1, K ) + Q( 3, K ) + TWO*Q( 4, K ) AQ( 2, K ) = Q( 3, K ) + Q( 4, K ) AQ( 3, K ) = TWO*Q( 3, K ) + Q( 4, K ) AQ( 4, K ) = Q( 3, K ) + Q( 4, K ) * 10 CONTINUE * RETURN END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'SRES' then echo shar: will not over-write existing file "'SRES'" else cat << \SHAR_EOF > 'SRES' ...... Type of Test problem: RWK ...... The test matrix size = 496 NV and M = 2 6 IT = 0 WR = 9.547E-01 8.114E-02 -4.440E-02 3.793E-02 -1.019E-02 -1.019E-02 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 6.322E-03 -6.322E-03 RSD = 2.534E-01 5.941E-01 5.856E-01 6.041E-01 6.015E-01 6.015E-01 NGRP = 1 CTR = 9.547E-01 AE = 9.547E-01 ARSD = 2.534E-01 NXTSRR = 5 IDORT = 1 IT = 5 WR = 9.969E-01 3.230E-01 3.230E-01 -2.954E-01 -1.914E-01 2.835E-02 WI = 0.000E+00 5.205E-02 -5.205E-02 0.000E+00 0.000E+00 0.000E+00 RSD = 8.032E-02 8.593E-01 8.593E-01 8.675E-01 8.878E-01 9.221E-01 NGRP = 1 CTR = 9.969E-01 AE = 9.969E-01 ARSD = 8.032E-02 NXTSRR = 7 IDORT = 1 IT = 7 WR = 9.985E-01 -4.240E-01 3.563E-01 -2.810E-01 2.145E-01 6.637E-02 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 6.545E-02 8.244E-01 9.003E-01 8.937E-01 8.756E-01 8.881E-01 NGRP = 1 CTR = 9.985E-01 AE = 9.985E-01 ARSD = 6.545E-02 NXTSRR = 10 IDORT = 1 IT = 10 WR = 9.983E-01 4.873E-01 -4.041E-01 -3.772E-01 3.235E-01 1.557E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 6.004E-02 8.040E-01 8.982E-01 8.660E-01 9.193E-01 9.263E-01 NGRP = 1 CTR = 9.983E-01 AE = 9.983E-01 ARSD = 6.004E-02 NXTSRR = 15 IDORT = 2 IT = 15 WR = 9.973E-01 -6.102E-01 5.842E-01 4.607E-01 2.890E-01 -2.379E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 4.633E-02 7.755E-01 7.253E-01 8.774E-01 8.990E-01 9.637E-01 NGRP = 1 CTR = 9.973E-01 AE = 9.973E-01 ARSD = 4.633E-02 NXTSRR = 22 IDORT = 2 IT = 22 WR = 9.981E-01 -8.038E-01 7.285E-01 6.324E-01 3.781E-01 -3.540E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 4.356E-02 5.885E-01 6.990E-01 7.443E-01 9.302E-01 8.976E-01 NGRP = 1 CTR = 9.981E-01 AE = 9.981E-01 ARSD = 4.356E-02 NXTSRR = 33 IDORT = 2 IT = 33 WR = 1.000E+00 -9.398E-01 9.217E-01 8.784E-01 -5.927E-01 1.048E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 2.765E-02 3.097E-01 3.701E-01 4.638E-01 7.542E-01 9.283E-01 NGRP = 1 CTR = 1.000E+00 AE = 1.000E+00 ARSD = 2.765E-02 NXTSRR = 49 IDORT = 1 IT = 49 WR = 1.000E+00 -9.870E-01 9.717E-01 9.000E-01 -7.130E-01 -4.267E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 6.739E-03 1.101E-01 1.273E-01 4.398E-01 5.964E-01 8.724E-01 NGRP = 1 CTR = 1.000E+00 AE = 1.000E+00 ARSD = 6.739E-03 NXTSRR = 73 IDORT = 2 IT = 73 WR = 1.000E+00 -9.917E-01 9.867E-01 -9.631E-01 9.392E-01 -8.748E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 2.695E-03 2.716E-02 3.938E-02 2.733E-01 2.978E-01 4.538E-01 NGRP = 1 CTR = 1.000E+00 AE = 1.000E+00 ARSD = 2.695E-03 NXTSRR = 109 IDORT = 2 IT = 109 WR = 1.000E+00 -9.974E-01 -9.974E-01 9.924E-01 9.712E-01 -9.691E-01 WI = 0.000E+00 8.061E-04 -8.061E-04 0.000E+00 0.000E+00 0.000E+00 RSD = 1.207E-03 3.388E-02 3.388E-02 1.288E-02 1.060E-01 1.298E-01 NGRP = 1 CTR = 1.000E+00 AE = 1.000E+00 ARSD = 1.207E-03 NXTSRR = 163 IDORT = 2 IT = 163 WR = 1.000E+00 -1.000E+00 -9.938E-01 9.934E-01 -9.753E-01 9.752E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 9.870E-05 3.319E-03 2.206E-03 1.268E-03 2.643E-02 2.649E-02 NGRP = 2 CTR = 1.000E+00 AE = 8.017E-06 ARSD = 2.348E-03 NXTSRR = 244 IDORT = 2 IT = 244 WR = 1.000E+00 -1.000E+00 -9.935E-01 9.935E-01 -9.755E-01 9.755E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 1.738E-06 5.403E-05 6.572E-05 3.626E-05 3.097E-03 3.296E-03 NGRP = 2 CTR = 1.000E+00 AE = 5.960E-08 ARSD = 3.823E-05 NXTSRR = 274 IDORT = 2 IT = 274 WR = 1.000E+00 -1.000E+00 9.935E-01 -9.935E-01 -9.755E-01 9.755E-01 WI = 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 RSD = 2.440E-06 1.188E-05 9.631E-06 1.754E-05 1.422E-03 1.522E-03 NGRP = 2 CTR = 1.000E+00 AE = 1.371E-06 ARSD = 8.576E-06 NGRP = 2 CTR = 9.935E-01 AE = 3.576E-07 ARSD = 1.415E-05 The INFO from SRRIT = 0 Total number of SRRIT iterations = 274 .. Converged Eigenvalues .. WR = 1.000002E+00 WI = 0.000000E+00 WR = -9.999996E-01 WI = 0.000000E+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 1.565E-07 1.095E-06 norm(Q*Q - I) = 9.536743E-07 ...... Type of Test problem: LRW ...... The test matrix size = 66 NV and M = 2 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 73 .. Converged Eigenvalues .. WR = -1.000000E+00 WI = 0.000000E+00 WR = 9.999996E-01 WI = 0.000000E+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 4.470E-08 1.714E-07 norm(Q*Q - I) = 4.768372E-07 ...... Type of Test problem: DRW ...... The test matrix size = 66 NV and M = 2 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 49 .. Converged Eigenvalues .. WR = -1.000001E+00 WI = 0.000000E+00 WR = 1.000000E+00 WI = 0.000000E+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 4.880E-06 6.557E-07 norm(Q*Q - I) = 5.960464E-07 ...... Type of Test problem: ODE ...... The test matrix size = 301 NV and M = 4 6 IT = 0 WR = -1.528E-01 3.013E-04 3.013E-04 1.455E-04 1.004E-04 2.828E-05 WI = 0.000E+00 1.064E-04 -1.064E-04 0.000E+00 0.000E+00 0.000E+00 RSD = 1.434E-01 1.094E-02 1.094E-02 2.704E-03 6.453E-03 6.423E-03 NGRP = 1 CTR = 1.528E-01 AE = -1.528E-01 ARSD = 1.434E-01 NXTSRR = 5 IDORT = 1 IT = 5 WR = -1.264E-02 -1.264E-02 4.429E-03 4.429E-03 2.542E-03 2.542E-03 WI = 2.312E-02 -2.312E-02 7.288E-03 -7.288E-03 2.575E-03 -2.575E-03 RSD = 1.115E-07 1.115E-07 2.912E-05 2.912E-05 6.126E-04 6.126E-04 NGRP = 2 CTR = 2.636E-02 AE = -1.264E-02 ARSD = 1.115E-07 NXTSRR = 7 IDORT = 1 IT = 7 WR = -1.264E-02 -1.264E-02 4.446E-03 4.446E-03 2.904E-03 2.904E-03 WI = 2.312E-02 -2.312E-02 7.309E-03 -7.309E-03 2.316E-03 -2.316E-03 RSD = 2.000E-08 2.000E-08 1.994E-06 1.994E-06 2.358E-04 2.358E-04 NGRP = 2 CTR = 2.636E-02 AE = -1.264E-02 ARSD = 2.000E-08 NGRP = 2 CTR = 8.555E-03 AE = 4.446E-03 ARSD = 1.994E-06 NXTSRR = 8 IDORT = 1 IT = 8 WR = -1.264E-02 -1.264E-02 4.447E-03 4.447E-03 2.888E-03 2.888E-03 WI = 2.312E-02 -2.312E-02 7.308E-03 -7.308E-03 2.254E-03 -2.254E-03 RSD = 2.000E-08 2.000E-08 4.426E-07 4.426E-07 1.292E-04 1.292E-04 NGRP = 2 CTR = 8.555E-03 AE = 4.447E-03 ARSD = 4.426E-07 NXTSRR = 9 IDORT = 1 IT = 9 WR = -1.264E-02 -1.264E-02 4.447E-03 4.447E-03 2.886E-03 2.886E-03 WI = 2.312E-02 -2.312E-02 7.309E-03 -7.309E-03 2.242E-03 -2.242E-03 RSD = 2.000E-08 2.000E-08 8.069E-08 8.069E-08 6.098E-05 6.098E-05 NGRP = 2 CTR = 8.555E-03 AE = 4.447E-03 ARSD = 8.069E-08 NGRP = 2 CTR = 3.654E-03 AE = 2.886E-03 ARSD = 6.098E-05 The INFO from SRRIT = 0 Total number of SRRIT iterations = 9 .. Converged Eigenvalues .. WR = -1.264390E-02 WI = 2.312457E-02 WR = -1.264390E-02 WI = -2.312457E-02 WR = 4.446721E-03 WI = 7.308519E-03 WR = 4.446721E-03 WI = -7.308519E-03 Final residual norms of (A*Q - Q*T) of converged RR pairs = 2.212E-09 6.752E-09 2.066E-08 1.240E-08 norm(Q*Q - I) = 2.384186E-07 ...... Type of Test problem: PDE ...... The test matrix size = 961 NV and M = 3 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 759 .. Converged Eigenvalues .. WR = 7.977825E+00 WI = 0.000000E+00 WR = 7.949039E+00 WI = 0.000000E+00 WR = 7.949021E+00 WI = 0.000000E+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 8.345E-07 2.235E-06 2.205E-06 norm(Q*Q - I) = 7.152557E-07 ...... Type of Test problem: PDE ...... The test matrix size = 225 NV and M = 1 1 The INFO from SRRIT = 0 Total number of SRRIT iterations = 549 .. Converged Eigenvalues .. WR = 7.911541E+00 WI = 0.000000E+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 1.013E-05 norm(Q*Q - I) = 1.788139E-07 ...... Type of Test problem: BWM ...... The test matrix size = 200 NV and M = 2 6 The INFO from SRRIT = 0 Total number of SRRIT iterations = 1234 .. Converged Eigenvalues .. WR = -1.235506E+03 WI = 0.000000E+00 WR = -1.234607E+03 WI = 0.000000E+00 WR = -1.233109E+03 WI = 0.000000E+00 WR = -1.231013E+03 WI = 0.000000E+00 WR = -1.228322E+03 WI = 0.000000E+00 Final residual norms of (A*Q - Q*T) of converged RR pairs = 5.646E-04 9.918E-05 4.120E-04 4.044E-04 3.052E-04 norm(Q*Q - I) = 5.960464E-07 ...... Type of Test problem: N44 ...... The test matrix size = 4 NV and M = 1 2 ORTH subroutine fails, INFO = 1 ...... Type of Test problem: F44 ...... The test matrix size = 4 NV and M = 1 2 Reduced matrix T is singular, INFO = 3 ...... Type of Test problem: M44 ...... The test matrix size = 4 NV and M = 1 2 SRRIT fails to converge after MAXIT, INF0 = 4 2-norms of Residual Vectors 6.594E-03 6.594E-03 Iteration numbers where the res. are computed 40 40 WR = 2.063278E+00 WI = 6.130946E-02 WR = 2.063278E+00 WI = -6.130946E-02 ...... Type of Test problem: SUB ...... Passed all SLAQR3 tests ..... END OF SRRIT TEST ..... SHAR_EOF fi # end of overwriting check if test -f 'data' then echo shar: will not over-write existing file "'data'" else cat << \SHAR_EOF > 'data' RWK # standard test (istart = -1) 30 2 6 LRW # initial vectors as input (istart = 0) 10 2 6 DRW # initial vectors as input (istart = 1) 10 2 6 ODE # standard test 300 4 6 PDE # standard test 31 3 6 PDE # special case NV = M = 1 15 1 1 BWM # standard test 200 2 6 N44 # initial zero vectors (info = 1) 4 1 2 F44 # will make T to be singular, and SRRIT fails. 4 1 2 M44 # method fails convergence after MAXIT iteration 4 1 2 SUB # test SGEHD2, SORGN2, SLAQR3 10 # number of dimensions 1 2 4 9 10 20 25 30 35 50 # MAX different dimension of test matrices END # end of all tests SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << \SHAR_EOF > 'driver.f' PROGRAM TEST * * ====================================================================== * * Test routine of SRRIT - a subroutine to calculate the dominant * invariant subspace of a nonsymmetric matrix. The SRRIT user guide is * documented in * * Z. Bai and G. W. Stewart, SRRIT - A FORTRAN Subroutine to * calculate the dominant invariant subspace of a nonsymmetric * matrix, submitted to ACM TOMS. * * With the file ``input'' as input, this test driver will exercise * major features of SRRIT. The output is going to be in the file * named ``output''. * * The current array sizes for testing SRRIT are set to handle the order * of a test matrix less than 1000, and the dimension of subspace less * than 100. For running a larger test or a larger subspace, user needs * to change the dimension parameters LDA and MAXSPC. * In addition, see the comments of SCHECK for the limit of the test * problem sizes for the testing subroutine SLAQR3. * * ---------------- * * Copyright (C) 1994, Zhaojun Bai and G. W. Stewart * All rights reserved. Use at your own risk. * Please send comments to bai@ms.uky.edu or stewart@cs.umd.edu * * ====================================================================== * * .. Parameters .. INTEGER LDA, LDQ, MAXSPC, LDT PARAMETER ( LDA = 1000, LDQ = LDA, MAXSPC = 100, $ LDT = MAXSPC ) INTEGER LWORK PARAMETER ( LWORK = MAXSPC*MAXSPC+5*MAXSPC ) REAL ZERO, HALF, ONE PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 ) INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) * .. * .. Local Scalars .. CHARACTER*3 PTYPE INTEGER I, INFO, ISTART, J, K, M, MAXIT, N, NV REAL QNORM, TOL * .. * .. Local Arrays .. INTEGER ITRSD( MAXSPC ), IWORK( LWORK ), ISEED( 4 ), $ NMAX( MAXSPC ) REAL AQ( LDA, MAXSPC ), DUMMY( 1 ), $ Q( LDQ, MAXSPC ), RSD( MAXSPC ), T( LDT, LDT ), $ WI( LDT ), WORK( LWORK ), WR( LDT ) * .. * .. External Functions .. LOGICAL LSAME REAL SLANGE, SLARAN EXTERNAL LSAME, SLANGE, SLARAN * .. * .. External Subroutines .. EXTERNAL SGEMM, SLASET, SRRIT EXTERNAL ATQODE, ATQPDE, ATQRWK, ATQBWM, ATQ4X4, ATQF44 * * .. Intrinsic functions .. INTRINSIC SIGN * .. * .. Executable Statements .. * OPEN( NOUT, FILE = 'soutput' ) * 10 CONTINUE * * Input the Type of Test Problem * READ( NIN, FMT = * )PTYPE * IF( .NOT.LSAME( PTYPE, 'END' ) ) $ WRITE( NOUT, FMT = 9999 )PTYPE * IF( LSAME( PTYPE, 'END' ) ) THEN WRITE( NOUT, FMT = 9998 ) GO TO 80 END IF * * Input Parameter K, which decides the size of test problem * READ( NIN, FMT = * )K * IF( LSAME( PTYPE, 'RWK' ) .OR. LSAME( PTYPE, 'LRW' ) .OR. $ LSAME( PTYPE, 'DWK' ) ) THEN N = ( K+1 )*( K+2 ) / 2 ELSE IF( LSAME( PTYPE, 'ODE' ) ) THEN N = K + 1 ELSE IF( LSAME( PTYPE, 'PDE' ) ) THEN N = K*K ELSE IF( LSAME( PTYPE, 'BWM' ) ) THEN N = K ELSE IF( LSAME( PTYPE, 'M44' ) .OR. LSAME( PTYPE, 'N44' ) ) THEN N = K ELSE IF( LSAME( PTYPE, 'F44' ) ) THEN N = K ELSE IF( LSAME( PTYPE, 'SUB' ) ) THEN GO TO 90 END IF * WRITE( NOUT, FMT = 9996 )N * IF( N.GT.LDA ) THEN WRITE( NOUT, FMT = 9995 ) GO TO 80 END IF * * Input NV, the number of the desired dominant eigenvalues. * M, the subspace size * Note that it is required that M >= NV+2. * READ( NIN, FMT = * )NV, M * WRITE( NOUT, FMT = 1111 )NV, M 1111 FORMAT( ' NV and M = ', 2I5) * * Set stopping criterion TOL and maximal number of iterations. * They may be altered if desired. * TOL = 1E-5 MAXIT = 10*N * * Call main SRRIT routine * IF( LSAME( PTYPE, 'RWK' ) ) THEN INFO = 1 ISTART = -1 CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQRWK ) ELSE IF( LSAME( PTYPE, 'LRW' ) ) THEN INFO = 0 ISTART = 0 * * generate initial vectors (istart = 0) * ISEED( 1 ) = 3192 ISEED( 2 ) = 1221 ISEED( 3 ) = 1322 ISEED( 4 ) = 8097 DO 30 J = 1, M DO 20 I = 1, N Q( I,J ) = SIGN( ONE, SLARAN( ISEED ) - HALF ) 20 CONTINUE 30 CONTINUE * CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQRWK ) ELSE IF( LSAME( PTYPE, 'DRW' ) ) THEN INFO = 0 ISTART = 1 * * generate and orthogonalize initial vectors (istart = 1) * ISEED( 1 ) = 3192 ISEED( 2 ) = 1221 ISEED( 3 ) = 1322 ISEED( 4 ) = 8097 DO 50 J = 1, M DO 40 I = 1, N Q( I,J ) = SLARAN( ISEED ) 40 CONTINUE 50 CONTINUE CALL ORTH( N, 1, M, Q, LDQ, INFO ) IF( INFO.NE.0 )THEN WRITE( NOUT, 1112 )INFO 1112 FORMAT( 'Given starting vectors are linearly dependent' ) GO TO 80 END IF * CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQRWK ) ELSE IF( LSAME( PTYPE, 'ODE' ) ) THEN INFO = 1 ISTART = -1 CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQODE ) ELSE IF( LSAME( PTYPE, 'PDE' ) ) THEN INFO = 0 ISTART = -1 CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQPDE ) ELSE IF( LSAME( PTYPE, 'BWM' ) ) THEN INFO = 0 ISTART = -1 CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQBWM ) ELSE IF( LSAME( PTYPE, 'M44' ) ) THEN INFO = 0 ISTART = -1 CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQ4X4 ) ELSE IF( LSAME( PTYPE, 'N44' ) ) THEN INFO = 0 ISTART = 0 DO 55 J = 1, M DO 45 I = 1, N Q( I,J ) = ONE 45 CONTINUE 55 CONTINUE CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQ4X4 ) ELSE IF( LSAME( PTYPE, 'F44' ) ) THEN INFO = 0 ISTART = 1 DO 75 J = 1, M DO 65 I = 1, N Q( I,J ) = ZERO 65 CONTINUE Q( J, J ) = ONE 75 CONTINUE CALL SRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, LDT, $ WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, INFO, TOL, $ ATQF44 ) END IF * IF( INFO.EQ.1 )THEN WRITE( NOUT, FMT = 1113 )INFO 1113 FORMAT( ' ORTH subroutine fails, INFO = ', I2 ) ELSE IF( INFO.EQ.2 )THEN WRITE( NOUT, FMT = 1114 )INFO 1114 FORMAT( ' SRRSTP subroutine fails, INFO = ', I2 ) ELSE IF( INFO.EQ.3 )THEN WRITE( NOUT, FMT = 1115 )INFO 1115 FORMAT( ' Reduced matrix T is singular, INFO = ', I2 ) ELSE IF( INFO.EQ.4 )THEN WRITE( NOUT, FMT = 1116 )INFO 1116 FORMAT( ' SRRIT fails to converge after MAXIT, INF0 = ', I2) WRITE( NOUT, FMT = 1118 ) 1118 FORMAT( ' 2-norms of Residual Vectors ' ) WRITE( NOUT, FMT = 9989 )( RSD( I ), I = 1, M ) WRITE( NOUT, FMT = 1119 ) 1119 FORMAT( ' Iteration numbers where the res. are computed ' ) WRITE( NOUT, FMT = 9985 )( ITRSD( I ), I = 1, M ) DO 61 K = 1, M WRITE( NOUT, FMT = 9991 )WR( K ), WI( K ) 61 CONTINUE ELSE WRITE( NOUT, FMT = 9993 )INFO WRITE( NOUT, FMT = 1117 )MAXIT 1117 FORMAT( ' Total number of SRRIT iterations = ', I5 ) END IF * IF( INFO.EQ.0 .AND. NV.GT.0 )THEN WRITE( NOUT, FMT = 9992 ) DO 60 K = 1, NV WRITE( NOUT, FMT = 9991 )WR( K ), WI( K ) 60 CONTINUE * * Test residual: R = A*Q - Q*T * CALL SGEMM( 'No transpose', 'No transpose', N, NV, NV, -ONE, $ Q, LDQ, T, LDT, ONE, AQ, LDA ) * DO 70 I = 1, NV WORK( I ) = SLANGE( 'MAX', N, 1, AQ( 1, I ), 1, DUMMY ) 70 CONTINUE * WRITE( NOUT, FMT = 9990 ) WRITE( NOUT, FMT = 9989 )( WORK( I ), I = 1, NV ) * CALL SLASET( 'Full', NV, NV, ZERO, ONE, WORK, NV ) CALL SGEMM( 'Transpose', 'No transpose', NV, NV, N, -ONE, Q, $ LDQ, Q, LDQ, ONE, WORK, NV ) QNORM = SLANGE( 'Max', NV, NV, WORK, NV, DUMMY ) * WRITE( NOUT, FMT = 9988 )QNORM * END IF * GO TO 10 * 9999 FORMAT( / 10X, ' ...... Type of Test problem: ', A3 , ' ......' ) 9998 FORMAT( / ' ..... END OF SRRIT TEST .....' ) 9996 FORMAT( / ' The test matrix size = ', I5 ) 9995 FORMAT( ' N is larger than array leading dimension LDA' ) 9993 FORMAT( / ' The INFO from SRRIT = ', I2 ) 9992 FORMAT( ' .. Converged Eigenvalues ..' ) 9991 FORMAT( ' WR = ', 1P,E15.6, ' WI = ', 1P,E15.6 ) 9990 FORMAT( ' Final residual norms ', $ 'of (A*Q - Q*T) of converged RR pairs = ' ) 9989 FORMAT( 1P,6E10.3 ) 9985 FORMAT( 6I10 ) 9988 FORMAT( ' norm(Q*Q - I) = ', 1P,E15.6 ) * * Test subroutines SLAQR3 * 90 CONTINUE * READ( NIN, FMT = * )( NMAX( I ), I = 1, K ) * CALL SCHECK( K, NMAX, INFO ) IF( INFO .EQ. 0 )THEN WRITE( NOUT, FMT = 9987 ) 9987 FORMAT( ' Passed all SLAQR3 tests ' ) ELSE WRITE( NOUT, FMT = 9986 )INFO 9986 FORMAT( ' Test routine SCHECK INFO = ', I5 ) END IF * GO TO 10 * 80 CONTINUE * CLOSE ( NOUT ) * END * ****** begin of scheck.f ******************************************** * SUBROUTINE SCHECK( MAX, NMAX, INFO ) * * .. Scalar Arguments .. INTEGER MAX, INFO * * .. Arrary Arguments .. INTEGER NMAX( * ) * * Purpose * ======= * * SCHECK tests the subroutines SLAQR3. * On normal return, INFO is to be zero, otherwise, if * * INFO = 1, SLAQR3 has illegal parameter in the calling sequence * or it fails to compute the complete Schur decomposition. * = 2, the backward errors of the Schur decomposition fails * test. * = 3, the order of the computed eigenvalues has not been * sorted correctly with option WANTT = .TRUE. and * WANTZ = .TRUE. in SLAQR3 * = 4, the order of the computed eigenvalues has not been * sorted correctly with option WANTT = .FALSE. and * WANTZ = .FALSE. in SLAQR3 * * Note: the internal parameter LDA is set so that the order of any * test matrix is smaller than 50. * * =================================================================== * * .. Parameters .. INTEGER LDA, LDB, LDH, LDU, LWORK PARAMETER ( LDA = 50, LDB = LDA, LDH = LDA, LDU = LDA, $ LWORK = 2*LDA*LDA ) * REAL ZERO, TOL PARAMETER ( ZERO = 0.0E+0, TOL = 2.0E+01 ) * * .. Local Scalars .. INTEGER I, J, K, N, IERR * * .. Local Arrays .. INTEGER ISEED( 4 ) REAL A( LDA,LDA ), B( LDB,LDB ), U( LDU,LDU ), $ H( LDH, LDH ), WR( LDA ), WI(LDA), $ WORK( LWORK ), RESULT( 2 ) * * .. External functions .. REAL SLAPY2, SLARAN EXTERNAL SLAPY2, SLARAN * * .. External Subroutines .. EXTERNAL SLACPY, SGEHRD, SORGHR, SLAQR3, SGET21 * INFO = 0 ISEED( 1 ) = 8321 ISEED( 2 ) = 9735 ISEED( 3 ) = 8673 ISEED( 4 ) = 2437 * DO 10 K = 1, MAX N = NMAX( K ) * * Generate the random test matrices * DO 30 J = 1,N DO 20 I = 1, N A( I,J ) = SLARAN( ISEED ) 20 CONTINUE 30 CONTINUE * * Save a copy of the matrix * CALL SLACPY( 'FULL', N, N, A, LDA, B, LDB ) * * Hessenberg reduction * CALL SGEHRD( N, 1, N, A, LDA, WR, WORK, LWORK, IERR ) DO 50 J = 1, N - 1 DO 40 I = J + 2, N U( I, J ) = A( I, J ) A( I, J ) = ZERO 40 CONTINUE 50 CONTINUE CALL SORGHR( N, 1, N, U, LDU, WR, WORK, LWORK, IERR ) * CALL SLACPY( 'FULL', N, N, A, LDA, H, LDH ) * * Compute the Schur decomposition of upper Hessenberg matrix * CALL SLAQR3( .TRUE., .TRUE., N, 1, N, A, LDA, WR, WI, 1, N, $ U, LDU, WORK, IERR ) IF( IERR .NE. 0 )THEN INFO = 1 RETURN END IF * CALL SGET21( N, B, LDB, A, LDA, U, LDU, WORK, RESULT ) IF( RESULT( 1 ) .GT. TOL .OR. RESULT( 2 ) .GT. TOL )THEN INFO = 2 RETURN END IF * * Check the ordering * DO 60 I = 1,N WORK( I ) = SLAPY2( WR( I ), WI( I ) ) 60 CONTINUE DO 80 I = N-1,-1,1 DO 70 J = I+1,-1,1 IF( WORK(I) .LT. WORK(J) )THEN INFO = 3 RETURN END IF 70 CONTINUE 80 CONTINUE * * Compute the eigenvalues only. * CALL SLAQR3( .FALSE., .FALSE., N, 1, N, H, LDH, WR, WI, $ 1, N, U, LDU, WORK, IERR ) * * Check the ordering * DO 90 I = 1,N WORK( I ) = SLAPY2( WR( I ), WI( I ) ) 90 CONTINUE DO 110 I = N-1,-1,1 DO 100 J = I+1,-1,1 IF( WORK(I) .LT. WORK(J) )THEN INFO = 4 RETURN END IF 100 CONTINUE 110 CONTINUE 10 CONTINUE * RETURN END * ********* begin of sget21.f ***************************************** * SUBROUTINE SGET21( N, A, LDA, B, LDB, U, LDU, WORK, RESULT ) * .. * .. Scalar Arguments .. INTEGER LDA, LDB, LDU, N * .. * .. Array Arguments .. REAL A( LDA, * ), B( LDB, * ), RESULT( * ), $ U( LDU, * ), WORK( * ) * * Purpose * ======= * * SGET21 generally checks a decomposition of the form * * A = U B U' * * where ' means transpose and U is orthogonal. Specifically, * * RESULT(1) = | A - U B U' | / ( |A| n ulp ) * and * RESULT(2) = | I - UU' | / ( n ulp ) * * Arguments * ========== * * N - INTEGER * The size of the matrix. If it is zero, SGET21 does nothing. * It must be at least zero. * Not modified. * * A - REAL array of dimension ( LDA , N ) * The original (unfactored) matrix. * Not modified. * * LDA - INTEGER * The leading dimension of A. It must be at least 1 * and at least N. * Not modified. * * B - REAL array of dimension ( LDB , N ) * The factored matrix. * Not modified. * * LDB - INTEGER * The leading dimension of B. It must be at least 1 * and at least N. * Not modified. * * U - REAL array of dimension ( LDU, N ). * The orthogonal matrix in the decomposition. * Not modified. * * LDU - INTEGER * The leading dimension of U. LDU must be at least N and * at least 1. * Not modified. * * WORK - REAL array of dimension ( 2*N**2 ) * Workspace. * Modified. * * RESULT - REAL array of dimension ( 2 ) * The values computed by the two tests described above. The * values are currently limited to 1/ulp, to avoid overflow. * Errors are flagged by RESULT(1)=10/ulp. * *----------------------------------------------------------------------- * * .. Parameters .. * REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. * INTEGER JDIAG REAL ANORM, ULP, UNFL, WNORM * .. * .. External Functions .. * REAL SLAMCH, SLANGE EXTERNAL SLAMCH, SLANGE * .. * .. External Subroutines .. * EXTERNAL SGEMM, SLACPY * .. * .. Intrinsic Functions .. * INTRINSIC MAX, MIN, REAL * .. * .. Executable Statements .. * RESULT( 1 ) = ZERO RESULT( 2 ) = ZERO IF( N.LE.0 ) $ RETURN * * Constants * UNFL = SLAMCH( 'Safe minimum' ) ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) * * Do Test 1 * * Norm of A: * ANORM = MAX( SLANGE( '1', N, N, A, LDA, WORK ), UNFL ) * * Norm of A - UBU' * CALL SLACPY( ' ', N, N, A, LDA, WORK, N ) CALL SGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO, $ WORK( N**2+1 ), N ) * CALL SGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, U, LDU, $ ONE, WORK, N ) * WNORM = SLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ) * IF( ANORM.GT.WNORM ) THEN RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) ELSE IF( ANORM.LT.ONE ) THEN RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) ELSE RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP ) END IF END IF * * Do Test 2 * * Compute UU' - I * CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, $ N ) * DO 40 JDIAG = 1, N WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+ $ 1 ) - ONE 40 CONTINUE * RESULT( 2 ) = MIN( SLANGE( '1', N, N, WORK, N, $ WORK( N**2+1 ) ), REAL( N ) ) / ( N*ULP ) * RETURN * * End of SGET21 * END * * ****** begin of atqrwk.f ********************************************* * SUBROUTINE ATQRWK( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. REAL AQ( LDA, * ), Q( LDQ, * ) * .. * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is for the random walk example in * the manuscript ``SRRIT - A Fortran Subroutine to Calculate the * Dominant Invariant subspace of a Nonsymmetric matrix'' by Z. Bai * and G. W. Stewart. * * The order of the random walk test matrix is * * N = (K+1)*(K+2)/2 * * where K is the size of triangular random walk grid. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) REAL array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) REAL array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * ===================================================================== * * .. Parameters .. REAL ZERO, HALF PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER I, J, K, KK, LK, NN REAL UPPER * .. * .. Intrinsic Functions .. INTRINSIC REAL, SQRT * .. * .. Statement Functions .. INTEGER IINDEX * .. * .. Statement Function definitions .. IINDEX( K, I, J ) = I*K - ( ( I-1 )*( I-2 ) ) / 2 + 1 + J + 1 * .. * .. Executable Statements .. * NN = ( -3+SQRT( REAL( 8*N+1 ) ) ) / 2 DO 30 KK = L, M * * Compute the vector AQ(1:N,KK) = A*Q(1:N,KK) * DO 20 I = 0, NN DO 10 J = 0, NN - I * K = IINDEX( NN, I, J ) UPPER = REAL( I+J ) / REAL( 2*NN ) AQ( K, KK ) = ZERO * IF( I.EQ.0 .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * ELSE IF( I.EQ.0 .AND. J.LE.NN-1 ) THEN LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * ELSE IF( I.EQ.0 .AND. J.EQ.NN ) THEN LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) * ELSE IF( I.LE.NN-1 .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * ELSE IF( I.EQ.NN .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + 2*UPPER*Q( LK, KK ) * ELSE IF( ( I+J ).EQ.NN ) THEN LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) * ELSE LK = IINDEX( NN, I-1, J ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) LK = IINDEX( NN, I, J-1 ) AQ( K, KK ) = AQ( K, KK ) + UPPER*Q( LK, KK ) LK = IINDEX( NN, I, J+1 ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) LK = IINDEX( NN, I+1, J ) AQ( K, KK ) = AQ( K, KK ) + ( HALF-UPPER )*Q( LK, KK ) * END IF 10 CONTINUE 20 CONTINUE * 30 CONTINUE * RETURN * * End of ATQRWK END * ****** begin of atqpde.f ********************************************** SUBROUTINE ATQPDE( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. REAL AQ( LDA, * ), Q( LDQ, * ) * .. * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is for the convection-diffusion * PDE example in the manuscript ``SRRIT - A Fortran Subroutine to * Calculate the Dominant Invariant subspace of a Nonsymmetric * matrix'' by Z. Bai and G. W. Stewart. * * The convection-diffusion operator is discretized on a K x K square * grid, and resulted a PDE test matrix of order * * N = K^2. * * The constant parameters p1, p2, and p3 in the convection-diffusion * operator can changed in this subroutine. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) REAL array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) REAL array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. REAL ONE, FOUR PARAMETER ( ONE = 1.0E+0, FOUR = 4.0E+0 ) REAL P1, P2, P3 PARAMETER ( P1 = 1.0E+0, P2 = 1.0E+0, P3 = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, II, K, NS REAL BETA, GAMMA, HSTEP, SIGMA * .. * .. Intrinsic Functions .. INTRINSIC REAL, SQRT * .. * .. Executable Statements .. * NS = SQRT( REAL( N ) ) HSTEP = ONE / ( REAL( NS )+ONE ) SIGMA = P3*HSTEP*HSTEP GAMMA = P2*HSTEP BETA = P1*HSTEP * DO 70 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * DO 60 I = 1, NS * DO 10 II = ( I-1 )*NS + 1, I*NS AQ( II, K ) = ( FOUR-SIGMA )*Q( II, K ) 10 CONTINUE * DO 20 II = ( I-1 )*NS + 1, I*NS - 1 AQ( II, K ) = AQ( II, K ) + ( GAMMA-ONE )*Q( II+1, K ) 20 CONTINUE * DO 30 II = ( I-1 )*NS + 1 + 1, I*NS AQ( II, K ) = AQ( II, K ) + ( -GAMMA-ONE )*Q( II-1, K ) 30 CONTINUE * IF( I.GT.1 ) THEN DO 40 II = ( I-1 )*NS + 1, I*NS AQ( II, K ) = AQ( II, K ) - ( BETA+ONE )*Q( II-NS, K ) 40 CONTINUE END IF IF( I.LT.NS ) THEN DO 50 II = ( I-1 )*NS + 1, I*NS AQ( II, K ) = AQ( II, K ) + ( BETA-ONE )*Q( II+NS, K ) 50 CONTINUE END IF * 60 CONTINUE * 70 CONTINUE * RETURN * * End of ATQPDE END * ********* begin of atqode.f ******************************************* SUBROUTINE ATQODE( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. REAL AQ( LDA, * ), Q( LDQ, * ) * .. * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is for the ODE boundary value * problem example in the manuscript ``SRRIT - A Fortran Subroutine * to Calculate the Dominant Invariant subspace of a Nonsymmetric * matrix'' by Z. Bai and G. W. Stewart. * * The order of the ODE test matrix is * * N = K + 1, * * where K is the number of subintervals of [0,1]. * * NOTE: N should be smaller than 1001, if one wants to test a larger * size of N, the size of working arrays A and U in this subroutine * should be increased. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) REAL array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) REAL array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. REAL ZERO, ONE, TWO PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) REAL THREE, FOUR PARAMETER ( THREE = 3.0E+0, FOUR = 4.0E+0 ) REAL GAMMA PARAMETER ( GAMMA = 1.0E-2 ) INTEGER LDU PARAMETER ( LDU = 1001 ) * .. * .. Local Scalars .. INTEGER I, INFO, J, K REAL DD, DELTA, HSTEP, HSTEP2 * .. * .. Local Arrays .. REAL A( 2, LDU ), U( LDU, 1 ) * .. * .. Intrinsic Functions .. INTRINSIC REAL * .. * .. External Subroutines .. EXTERNAL SPBTF2, SPBTRS * .. * .. Executable Statements .. * HSTEP = ONE / REAL( N ) HSTEP2 = HSTEP*HSTEP * A( 1, 1 ) = ZERO DO 10 J = 2, N - 1 A( 1, J ) = -ONE 10 CONTINUE * DO 20 J = 1, N - 1 A( 2, J ) = TWO 20 CONTINUE * * Cholesky decomposition of Symmetric tridiagonal matrix * CALL SPBTF2( 'Upper', N-1, 1, A, 2, INFO ) * DO 30 J = 1, N - 2 U( J, 1 ) = ZERO 30 CONTINUE U( N-1, 1 ) = ONE * CALL SPBTRS( 'Upper', N-1, 1, 1, A, 2, U, 301, INFO ) * DELTA = THREE*GAMMA + ( FOUR*U( 1, 1 )-U( 2, 1 )+GAMMA* $ U( N-2, 1 )-FOUR*GAMMA*U( N-1, 1 ) ) * DO 60 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * DO 40 I = 1, N - 1 AQ( I, K ) = HSTEP2*Q( I, K ) 40 CONTINUE AQ( N, K ) = ZERO * CALL SPBTRS( 'Upper', N-1, 1, 1, A, 2, AQ( 1, K ), LDA, INFO ) * DD = FOUR*AQ( 1, K ) - AQ( 2, K ) + GAMMA*AQ( N-2, K ) - $ FOUR*GAMMA*AQ( N-1, K ) AQ( N, K ) = DD / DELTA * DO 50 I = 1, N - 1 AQ( I, K ) = -AQ( I, K ) + U( I, 1 )*AQ( N, K ) 50 CONTINUE * 60 CONTINUE * DO 80 K = L, M DO 70 I = 1, N AQ( I, K ) = -AQ( I, K ) 70 CONTINUE 80 CONTINUE * RETURN * * End of ATQODE END ********* begin of atqbwm.f ****************************************** SUBROUTINE ATQBWM( N, L, M, X, LDX, Y, LDY ) * .. * .. Scalar Arguments .. INTEGER LDY, LDX, L, M, N * .. * .. Array Arguments .. REAL Y( LDY, * ), X( LDX, * ) * .. * * Purpose * ======= * * Compute * * Y(:,L:M) = A*X(:,L:M) * * The matrix A is a 2 by 2 block matrix resulted from the finite * difference discretization of the Brusselator wave model. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N has to be an even number. * * L,M (input) INTEGERS * The number of the first and the last column of Q to * multiply by A. * * X (input) REAL array, dimension ( LDX, M ) * contains the vectors X. * * LDX (input) INTEGER * The leading dimension of the array X, LDX >= max( 1, N ) * * Y (output) REAL array, dimension (LDY, M ) * contains the product of the matrix op(A) with X. * * LDY (input) INTEGER * The leading dimension of the array Y, LDY >= max( 1, N ) * * =================================================================== * * .. Parameters .. REAL ONE, TWO PARAMETER ( ONE = 1.0E+0, TWO = 2.0E+0 ) REAL ALPHA, BETA, DELT1, DELT2, LL PARAMETER ( ALPHA = 2.0E+0, BETA = 5.45E+0, $ DELT1 = 8.0E-3, DELT2 = 4.0E-3, $ LL = 0.51302E+0 ) * * .. Local Variables .. INTEGER I, K, NHALF REAL HSTEP, TAU1, TAU2, DIAG1, DIAG2, ALPHA2 * * set up middle parameters * NHALF = N / 2 HSTEP = ONE / ( REAL( NHALF ) + 1 ) TAU1 = DELT1 / (HSTEP*LL)**2 TAU2 = DELT2 / (HSTEP*LL)**2 DIAG1 = -TWO*TAU1 + BETA - ONE ALPHA2 = ALPHA**2 DIAG2 = -TWO*TAU2 - ALPHA2 * * Compute Y(:,L:M) = A*X(:,L:M) * DO 30 K = L, M * * the 1st row * Y( 1, K ) = DIAG1*X( 1, K ) + TAU1*X( 2, K ) + $ ALPHA2*X( NHALF+1, K ) * * the 2nd row to (NHALF-1)th row * DO 10 I = 2, NHALF - 1 Y( I, K ) = TAU1*X( I-1, K ) + DIAG1*X( I, K ) + $ TAU1*X( I+1, K ) + ALPHA2*X( NHALF+I, K ) 10 CONTINUE * * the NHALFth row * Y( NHALF, K ) = TAU1*X( NHALF-1, K ) + DIAG1*X( NHALF, K ) $ + ALPHA2*X( N,K ) * * the (NHALF+1)th row * Y( NHALF+1, K ) = -BETA*X( 1, K ) + DIAG2*X( NHALF+1, K ) $ + TAU2*X( NHALF+2, K ) * * the (NHALF+2)th row to (N-1)th row * DO 20 I = NHALF+2, N - 1 Y( I, K ) = -BETA*X( I-NHALF, K ) + TAU2*X( I-1, K ) + $ DIAG2*X( I, K ) + TAU2*X( I+1, K ) 20 CONTINUE * * the last (Nth) row * Y( N, K ) = -BETA*X( N-NHALF, K ) + TAU2*X( N-1,K ) + $ DIAG2*X( N, K ) * 30 CONTINUE * RETURN * * End of ATQBWM * END ********* begin of atq4x4.f ****************************************** SUBROUTINE ATQ4X4( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. REAL AQ( LDA, * ), Q( LDQ, * ) * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is a 3 by 3 matrix: * * A = [ 2 1 1 1 ] * [ 0 2 1 1 ] * [ 0 0 2 1 ] * [ 0 0 1 ] * * which has exact eigenvalues 2, 2, 2, 1. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) REAL array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) REAL array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. REAL TWO PARAMETER ( TWO = 2.0E+0 ) * * .. * .. Local Scalars .. INTEGER K * N = 4 DO 10 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * AQ( 1, K ) = TWO*Q( 1, K ) + Q( 2, K ) + Q( 3, K ) + Q( 4, K ) AQ( 2, K ) = TWO*Q( 2, K ) + Q( 3, K ) + Q( 4, K ) AQ( 3, K ) = TWO*Q( 3, K ) + Q( 4, K ) AQ( 4, K ) = Q( 4, K ) * 10 CONTINUE * RETURN END SUBROUTINE ATQF44( N, L, M, Q, LDQ, AQ, LDA ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, M, N * .. * .. Array Arguments .. REAL AQ( LDA, * ), Q( LDQ, * ) * * Purpose * ======= * * Compute * * AQ(:,L:M) = A*Q(:,L:M) * * Note that we only need to access the matrix A in question via * matrix-vector product form. Therefore, any special structure * and/or sparsity of matrix A can be exploited in forming such * matrix-vector products. * * The matrix A in this subroutine is a 3 by 3 matrix: * * A = [ 2 0 1 2 ] * [ 0 0 1 1 ] * [ 0 0 2 1 ] * [ 0 0 1 1 ] * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * L, M (input) INTEGERS * The first and last columns of Q to multiply by the matrix Q. * * Q (input) REAL array, dimension ( LDQ, M ) * Q contains the matrix Q. * * AQ (output) REAL array, dimension (LDQ, M ) * On return, columns L through M of AQ should contain the product * of the matrix A with columns L through M of Q. * * =================================================================== * * .. Parameters .. REAL TWO PARAMETER ( TWO = 2.0E+0 ) * * .. * .. Local Scalars .. INTEGER K * N = 4 DO 10 K = L, M * * Compute the vector AQ(1:N,K) = A*Q(1:N,K) * AQ( 1, K ) = TWO*Q( 1, K ) + Q( 3, K ) + TWO*Q( 4, K ) AQ( 2, K ) = Q( 3, K ) + Q( 4, K ) AQ( 3, K ) = TWO*Q( 3, K ) + Q( 4, K ) AQ( 4, K ) = Q( 3, K ) + Q( 4, K ) * 10 CONTINUE * RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'dsrrit.f' then echo shar: will not over-write existing file "'dsrrit.f'" else cat << \SHAR_EOF > 'dsrrit.f' SUBROUTINE DSRRIT( N, NV, M, MAXIT, ISTART, Q, LDQ, AQ, LDA, T, $ LDT, WR, WI, RSD, ITRSD, IWORK, WORK, LWORK, $ INFO, EPS, ATQ ) * * ==================================================================== * * Copyright (C) 1994, Zhaojun Bai and G. W. Stewart * All rights reserved. Use at your own risk. * Please send comments to bai@ms.uky.edu or stewart@cs.umd.edu * * ==================================================================== * .. * .. Scalar Arguments .. INTEGER INFO, ISTART, LDA, LDQ, LDT, LWORK, M, MAXIT, $ N, NV DOUBLE PRECISION EPS * .. * .. Array Arguments .. INTEGER ITRSD( * ), IWORK( * ) DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ), RSD( * ), $ T( LDT, * ), WI( * ), WORK( * ), WR( * ) * .. * .. Subroutine Arguments .. EXTERNAL ATQ * .. * * Purpose * ======= * * SRRIT computes a nested sequence of orthonormal bases for the * dominant invariant subspaces of a real matrix A of order N. * Specifically, the program returns an N x NV matrix Q with * orthonormal columns and an NV x NV quasi-triangular matrix T * satisfying * * A*Q = Q*T + O(eps) * * The eigenvalues of T in 1 x 1 diagonal blocks are real, in the 2 x 2 * diagonal blocks are complex conjugate pairs. The eigenvalues of T: * W(J) = WR(J) + i*WI(J), J = 1, 2, ... , NV, are ordered so that * * |W(1)| >= |W(2)| >= ... >= |W(NV)|. * * These eigenvalues approximate the dominant eigenvalues of A, i.e., * the eigenvalues of largest absolute magnitudes. These facts have the * following consequences * * 1. If W(L).ne.W(L+1) and W(L).ne.conj(W(L+1)), then columns * 1,2,...,L of Q form an approximate basis for the invariant * subspace corresponding to the L dominant eigenvalues of A. * The L by L leading principal submatrix of T is a representation * of A in that subspace with respect to the basis Q. * * 2. If z is an eigenvector of T corresponding to W(J), then Q*z is * an approximate eigenvector of A corresponding to W(J). z may * be computed by calling LAPACK subroutine STREVC. * * The program actually iterates with an N x M matrix Q and an M x M * matrix T. The rate of convergence of the L-th column of Q is * essentially linear with ratio |W(M+1)/W(L)|. The program may fail to * converge after the user specified the maximal number of iterations, * MAXIT, see the arguments NV, MAXIT and INFO. In general, it pays that * the user sets M larger than the desired number, NV, of eigenvalues, * say M = NV + 5 or larger. * * The user must furnish a subroutine to compute the product A*Q. * The calling sequence is * * ATQ( N, L, M, Q, LDQ, AQ, LDA ) * * For J = L, L+1, ..., M, the program must place the product A*Q(:,J) * in AQ(:,J). * * The method is based on the paper by G. W. Stewart 'Simultaneous * iteration for computing invariant subspaces of non-hermitian * matrices', Numer, Math. 25(1976), pp.123-136. * * For the further details of algorithm implementation and usage, see * * [1] Z. Bai and G. W. Stewart, SRRIT - A Fortran subroutine to * calculate the dominant invariant subspace of a nonsymmetric matrix, * Technical Report TR-2908, Department of Computer Science, * University of Maryland, 1992 (Submitted to ACM TOMS for * publication). * * Arguments * ========= * * N (input) INTEGER * N is the order of the matrix A. * * NV (input/output) INTEGER * One entry, NV is the size of the dominant invariant subspace * of A that the user desired. On exit, NV is the actual size of * converged invariant subspace. * Note that if input NV is larger than output NV, say, on exit, * NV = 0, then it means that the method fails to converge to * the number of desired eigenvalues after the MAXIT number of * iterations. * * M (input) INTEGER * M is the size of iteration space. NV <= M <= N. * * MAXIT (input/output) INTEGER * On entry, MAXIT is an upper bound on the number of iterations * the program is to execute. On exit, MAXIT is the actual number * of iteration the program executed (also see NV). * * ISTART (input) INTEGER * ISTART specifies whether user supply initial Q. * < 0, Q is initialized by the program, random vectors with * uniform (0,1) distribution. * = 0, starting Q has been set in the input, but * it is not orthonormal. * > 0, starting Q has been set in the input, it is also * orthonormal. * * Q (input/output) DOUBLE PRECISION array, dimension( LDQ,M ) * On entry, if ISTART > 0, Q contains the starting Q which will * be used in the simultaneous iteration. * On exit, Q contains the orthonormal vectors described above. * * LDQ (input) INTEGER * The leading dimension of Q, LDQ >= max(1,N). * * AQ (output) DOUBLE PRECISION array, dimension( LDA, M) * On exit, AQ contains the product A*Q. * * LDA (input) INTEGER * The leading dimension of AQ, LDA >= max(1,N). * * T (output) DOUBLE PRECISION array, dimension( LDT,M ) * On exit, T contains of representation of A described above. * * LDT (input) INTEGER * The leading dimension of T, LDT >= max(1,M). * * WR,WI (output) DOUBLE PRECISION arrays, dimension (M) * On exit, WR and WI contain the real and imaginary parts, resp. * of the eigenvalues of T, which are the approximations of the * dominant eigenvalues of matrix A. The eigenvalues are ordered * in decreasing. * * RSD (output) DOUBLE PRECISION arrays, dimension( M ) * On exit, RSD contains the 2-norm of the residual vectors * R(:,J) = A*Q(:,J) - Q*T(:,J) for J = 1,...,M. * * ITRSD (output) INTEGER array, dimension( M ) * On exit, ITRSD contains the iteration numbers at which * the residuals were computed. * * IWORK (workspace) INTEGER array, dimension( 2*M ) * * WORK (workspace) DOUBLE PRECISION array, WORK( LWORK ) * * LWORK (input) INTEGER * The length of workspace WORK, LWORK >= M*M + 5*M * Note: the first M*M workspace is used for the matrix U * in DSRRSP, next 2*M is the workspace for DSRRSP, additional * 2*M for storing old WR and WI, and finally, M for storing * old RSD in this program. * * INFO (input/output) INTEGER * On entry, if INFO is set to be a nonzero integer, then * the information about the course of the iteration will be * printed, otherwise, it is omitted. * On exit, if INFO is set to * 0, normal return. * -k, if input argument number k is illegal. * 1, Orthonormalization subroutine DORTH fails. * 2, error from subroutine DSRRSP. * 3, matrix T is singular, see DCOND. * 4, method fails to converge the number of desired * eigenvalues (NV) after MAXIT iterations. RSD * contains the 2-norm of residual vectors. * * EPS (Input) DOUBLE PRECISION * A convergence criterion supplied by user. * * ATQ External procedure name, which passes the matrix-vectors * product operation. * * * Further details: * ================ * * The internal control parameters INIT, STPFAC, ALPHA, BETA, GRPTOL, * CNVTOL and ORTTOL are described in reference [1]. Currently, they * are set as default initial values. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) DOUBLE PRECISION STPFAC, ALPHA, BETA PARAMETER ( STPFAC = 1.5, ALPHA = ONE, BETA = 1.1D+0 ) DOUBLE PRECISION ORTTOL PARAMETER ( ORTTOL = TWO ) DOUBLE PRECISION GRPTOL, CNVTOL PARAMETER ( GRPTOL = 1.0D-8, CNVTOL = 1.0D-6 ) * INTEGER INIT, NOUT PARAMETER ( INIT = 5, NOUT = 6 ) * .. * .. Local Scalars .. LOGICAL PINFO INTEGER I, IDORT, IDSRR, IERR, IT, J, L, M2, NGRP, $ NOGRP, NV0, NXTORT, NXTSRR DOUBLE PRECISION AE, AQMAX, ARSD, CTR, OAE, OARSD, OCTR, REC, $ TCOND * .. * .. Local Arrays .. INTEGER ISEED( 4 ) DOUBLE PRECISION DUMMY( 1 ) * .. * .. External Functions .. DOUBLE PRECISION DCOND, DLANGE, DLARAN EXTERNAL DCOND, DLANGE, DLARAN * .. * .. External Subroutines .. EXTERNAL DGROUP, DORTH, DRESID, DSRRSP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, LOG10, MAX, MIN, DBLE * .. * .. Executable Statements .. * * Decide whether print the middle computed results. * PINFO = .TRUE. IF( INFO.EQ.0 ) $ PINFO = .FALSE. * * Test the input arguments * INFO = 0 IF( N.LT.0 ) THEN INFO = -1 ELSE IF( NV.GT.N .OR. NV.GT.M ) THEN INFO = -2 ELSE IF( M.GT.N ) THEN INFO = -3 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -7 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -9 ELSE IF( LDT.LT.MAX( 1, M ) ) THEN INFO = -11 ELSE IF( LWORK.LT.( M*M+5*M ) ) THEN INFO = -18 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSRRIT', -INFO ) RETURN END IF * * Initialize * IT = 0 L = 1 M2 = M*M + 2*M DO 10 J = 1, M WR( J ) = ZERO WI( J ) = ZERO RSD( J ) = ZERO ITRSD( J ) = -1 10 CONTINUE NV0 = NV * IF( ISTART.GE.0 ) $ GO TO 40 * * Set initial random vectors with uniform (0,1) distribution. * ISEED( 1 ) = 1978 ISEED( 2 ) = 1989 ISEED( 3 ) = 1992 ISEED( 4 ) = 1993 DO 30 J = 1, M DO 20 I = 1, N Q( I, J ) = DLARAN( ISEED ) 20 CONTINUE 30 CONTINUE * * Use the input vectors as initial vectors. * 40 CONTINUE IF( ISTART.GT.0 ) $ GO TO 50 * * Orthogonalize Q if necessary * CALL DORTH( N, 1, M, Q, LDQ, IERR ) IF( IERR.NE.0 )THEN INFO = 1 RETURN END IF * 50 CONTINUE * * The main Schur-Rayleigh-Ritz (SRR) loop * 60 CONTINUE * DO 70 I = L, M WORK( M2+I ) = WR( I ) WORK( M2+M+I ) = WI( I ) 70 CONTINUE * IF( PINFO ) THEN WRITE( NOUT, FMT = 9999 )IT 9999 FORMAT( 30X, 'IT = ', I4 ) END IF * CALL DSRRSP( N, L, M, Q, LDQ, AQ, LDA, T, LDT, WR, WI, WORK, M, $ WORK( M*M+1 ), IERR, ATQ ) IF( IERR.NE.0 )THEN INFO = 2 RETURN END IF * DO 80 I = L, M WORK( M2+2*M+I ) = RSD( I ) IWORK( M+I ) = ITRSD( I ) ITRSD( I ) = IT 80 CONTINUE * * Compute the 2-norms of residual vectors: * || R(1:M,L:M) || = || A*Q(1:M,L:M) - Q*T(1:M,L:M) || * CALL DRESID( N, L, M, Q, LDQ, AQ, LDA, T, LDT, RSD ) * IF( PINFO ) THEN WRITE( NOUT, FMT = 9998 )( WR( I ), I = 1, M ) 9998 FORMAT( ' WR =', 1P,6D12.3 ) WRITE( NOUT, FMT = 9997 )( WI( I ), I = 1, M ) 9997 FORMAT( ' WI =', 1P,6D12.3 ) WRITE( NOUT, FMT = 9996 )( RSD( I ), I = 1, M ) 9996 FORMAT( ' RSD =', 1P,6D12.3 ) END IF * 90 CONTINUE * * Find clusters of computed eigenvalues * CALL DGROUP( L, M, WR, WI, RSD, NGRP, CTR, AE, ARSD, GRPTOL ) CALL DGROUP( L, M, WORK( M2+1 ), WORK( M2+M+1 ), WORK( M2+2*M+1 ), $ NOGRP, OCTR, OAE, OARSD, GRPTOL ) * IF( PINFO ) THEN WRITE( NOUT, FMT = 9995 )NGRP 9995 FORMAT( ' NGRP =', I5 ) WRITE( NOUT, FMT = 9994 )CTR, AE, ARSD 9994 FORMAT( ' CTR =', 1P,D12.3, ' AE =', 1P,D12.3, $' ARSD = ', 1P,D12.3 ) END IF * IF( NGRP.NE.NOGRP ) $ GO TO 100 IF( NGRP.EQ.0 ) $ GO TO 100 * IF( ABS( AE-OAE ).GT.CTR*CNVTOL*DBLE( ITRSD( L )-IWORK( M+L ) ) ) $ GO TO 100 * * Test for the convergence. * IF( ARSD.GT.CTR*EPS ) $ GO TO 100 L = L + NGRP IF( L.GT.M ) $ GO TO 100 * GO TO 90 100 CONTINUE * * Exit if the required number of eigenvalues have converged. * IF( L.GT.NV ) $ GO TO 210 * * Exit if iteration count exceeds the maximum number of iterations * IF( IT.GE.MAXIT ) $ GO TO 210 * * Determine when the next SRR step (NXTSRR) is to be taken. * NXTSRR = MIN( MAXIT, MAX( INT( STPFAC*DBLE( IT ) ), INIT ) ) IDSRR = NXTSRR - IT IF( NGRP.NE.NOGRP ) $ GO TO 110 IF( NGRP.EQ.0 ) $ GO TO 110 IF( ARSD.GE.OARSD ) $ GO TO 110 REC = ALPHA + BETA*DBLE( IWORK( M+L )-ITRSD( L ) )* $ LOG( ARSD / EPS ) / LOG( ARSD / OARSD ) IDSRR = MAX( 1, INT( REC ) ) * 110 CONTINUE NXTSRR = MIN( NXTSRR, IT+IDSRR ) * * Determine the interval between orthogonalizations. * The next orthogonalization step at NXTORT * DO 130 J = 1, M DO 120 I = 1, M WORK( I+( J-1 )*M ) = T( I, J ) 120 CONTINUE 130 CONTINUE TCOND = DCOND( M, WORK, M, IWORK, WORK( M*M+1 ), IERR ) IF( IERR.NE.0 )THEN INFO = 3 RETURN END IF IDORT = MAX( 1, INT( ORTTOL / MAX( ONE, LOG10( TCOND ) ) ) ) NXTORT = MIN( IT+IDORT, NXTSRR ) * IF( PINFO ) THEN WRITE( NOUT, FMT = 9993 )NXTSRR, IDORT 9993 FORMAT( ' NXTSRR =', I4, ' IDORT =', I4 ) END IF * DO 150 J = L, M DO 140 I = 1, N Q( I, J ) = AQ( I, J ) 140 CONTINUE 150 CONTINUE IT = IT + 1 * 160 CONTINUE * * Orthogonalization loop * * Power loop. After each matrix-vector step, we scale the * product to prevent the possible overflow. * 170 CONTINUE IF( IT.EQ.NXTORT ) $ GO TO 200 * CALL ATQ( N, L, M, Q, LDQ, AQ, LDA ) AQMAX = DLANGE( 'M', N, M-L+1, AQ( 1, L ), LDA, DUMMY ) REC = ONE / AQMAX DO 190 J = L, M DO 180 I = 1, N Q( I, J ) = AQ( I, J )*REC 180 CONTINUE 190 CONTINUE IT = IT + 1 GO TO 170 * 200 CONTINUE CALL DORTH( N, L, M, Q, LDQ, IERR ) IF( IERR.NE.0 )THEN INFO = 1 RETURN END IF NXTORT = MIN( IT+IDORT, NXTSRR ) IF( IT.LT.NXTSRR ) $ GO TO 160 * GO TO 60 * 210 CONTINUE * * Final number of converged approx. eigenvalues, and number of * total iterations (= the number of calls to ATQ). * NV = L - 1 MAXIT = IT * IF( NV.LT.NV0 ) $ INFO = 4 * RETURN * * End of DSRRIT END * ******** begin of dsrrsp.f ******************************************** SUBROUTINE DSRRSP( N, L, M, Q, LDQ, AQ, LDA, T, LDT, WR, WI, U, $ LDU, WORK, INFO, ATQ ) * .. * .. Scalar Arguments .. INTEGER INFO, L, LDA, LDQ, LDT, LDU, M, N * .. * .. Array Arguments .. DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ), T( LDT, * ), $ U( LDU, * ), WI( * ), WORK( * ), WR( * ) * .. * .. Subroutine Arguments .. EXTERNAL ATQ * .. * * Purpose * ======= * * DRRSP performs a Schur-Rayleigh-Ritz refinement on the set of * M-L+1 orthogonal N-vectors contained in the array Q(:,L:M) in * the following three steps: * * 1). Compute T = Q(:,L:M)'*A*Q(:,L:M); * 2). Compute Schur decomposition of T: T := U'*T*U; * 3). Update Q(:,L:M) = Q(:,L:M)*U; * AQ(:,L:M) = AQ(:,L:M)*U. * * Note that in the Schur decomposition of T, the eigenvalues in the * Schur form are ordered in decreasing in terms of absolute magnitude. * * Arguments * ========= * * N (input) INTEGER * N is the order of the matrix A. * * L, M (input) INTEGER * The first and last column indices of Q for SRR iteration. * * Q (input/output) DOUBLE PRECISION array, dimension( LDQ,M ) * On entry, Q contains approximate invariant subpace matrix. * On exit, Q is updated by Schur vectors of matrix T. * * LDQ (input) INTEGER * The leading dimension of Q, LDQ >= max(1,N). * * AQ (output) DOUBLE PRECISION array, dimension( LDA, M) * On exit, AQ contains the product A*Q. * * LDA (input) INTEGER * The leading dimension of AQ, LDA >= max(1,N). * * T (output) DOUBLE PRECISION array, dimension( LDT,M ) * On exit, T contains of Schur form computed by SRR iteration. * * LDT (input) INTEGER * The leading dimension of T, LDT >= max(1,M). * * WR,WI (output) DOUBLE PRECISION arrays, dimension (M) * On exit, WR and WI contain the real and imaginary parts of * the eigenvalues of T, respectively. The eigenvalues are * ordered in decreasing. * * U (output) DOUBLE PRECISION array, dimension(M,M) * On exit, U contains the orthogonal matrix of the Schur * decomposition. * * LDU (input) INTEGER * The leading dimension of U, LDU >= max(1,M). * * WORK (workspace) DOUBLE PRECISION array, dimension( 2*M ) * * INFO (output) INTEGER * On exit, if INFO is set to * 0, normal return * 1, error from reduction T to Hessenberg form (DGEHD2). * 2, error from forming orthogonal matrix U (DORGN2). * 3, error from Schur decomposition routine (DLAQR3). * * ATQ External procedure name, which passes the matrix-vectors * product operation. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, IERR, J, K, NC * .. * .. External Subroutines .. EXTERNAL DGEHRD, DGEMM, DLAQR3, DLASET, DORGHR * .. * .. Executable Statements .. * INFO = 0 * * Calculate T(1:M,L:M) = Q(:,1:M)'*A*Q(:,L:M) * NC = M - L + 1 CALL ATQ( N, L, M, Q, LDQ, AQ, LDA ) CALL DGEMM( 'Transpose', 'No transpose', M, NC, N, ONE, Q, LDQ, $ AQ( 1, L ), LDA, ZERO, T( 1, L ), LDT ) * * Hessenberg reduction * CALL DGEHRD( M, L, M, T, LDT, WORK, WORK( M+1 ), M, IERR ) IF( IERR.NE.0 )THEN INFO = 1 RETURN END IF DO 20 J = 1, M - 1 DO 10 I = J + 2, M U( I, J ) = T( I, J ) T( I, J ) = ZERO 10 CONTINUE 20 CONTINUE CALL DORGHR( M, L, M, U, LDU, WORK, WORK( M+1 ), M, IERR ) * * Compute Schur decomposition * CALL DLAQR3( .TRUE., .TRUE., M, L, M, T, LDT, WR, WI, 1, M, $ U, LDU, WORK, IERR ) IF( IERR.NE.0 )THEN INFO = 3 RETURN END IF * * Update: Q(:,L:M) := Q(:,L:M)*U * AQ(:,L:M) := AQ(:,L:M)*U * DO 60 I = 1, N DO 40 J = L, M WORK( J ) = ZERO WORK( J+M ) = ZERO DO 30 K = 1, M WORK( J ) = WORK( J ) + Q( I, K )*U( K, J ) WORK( J+M ) = WORK( J+M ) + AQ( I, K )*U( K, J ) 30 CONTINUE 40 CONTINUE DO 50 J = L, M Q( I, J ) = WORK( J ) AQ( I, J ) = WORK( J+M ) 50 CONTINUE 60 CONTINUE * RETURN * * End of DSRRSP END * ******** begin of dorth.f ********************************************** SUBROUTINE DORTH( N, L, M, Q, LDQ, INFO ) * .. * .. Scalar Arguments .. INTEGER INFO, L, LDQ, M, N * .. * .. Array Arguments .. DOUBLE PRECISION Q( LDQ, * ) * .. * * Purpose * ======= * * DORTH orthonormalizes columns L through M of the array Q with respect * to columns 1 to M. Column 1 through L-1 are assumed to be * orthonormal. The method is the modified Gram-Schmidt method with * reorthogonalization. A column is accepted when an orthogonalization * does not reduce its 2-norm by a factor of more than TOL (0.5). If * this is not done in MAXTRY (5) attempts the program stops. DORTH also * stops if it encounters a zero vector. * * Arguments * ========= * * N (input) INTEGER * The number of rows of Q. * * L,M (input) INTEGERs * The first and last column indices of Q which is * going to be orthonormalized. * * Q (input/output) DOUBLE PRECISION array, dimension(LDQ,M) * On entry, Q contains the vectors which is going to be * orthonormalized. * On exit, Q is overwritten by the orthonormalized vectors. * * LDQ (input) INTEGER * The leading dimension of array Q, LDQ >=max(1,N) * * INFO (output) INTEGER * On exit, if INFO is set to * 0, successful return. * 1, zero vectors in Q. * 2, reorthogonalization iteration fails after MAXTRY * attempts. * * ===================================================================== * * .. Parameters .. INTEGER MAXTRY PARAMETER ( MAXTRY = 5 ) DOUBLE PRECISION ZERO, ONE, TOL PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TOL = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER ITRY, J, JM1, K DOUBLE PRECISION NORM, QQ, REC * .. * .. External Functions .. DOUBLE PRECISION DDOT, DNRM2 EXTERNAL DDOT, DNRM2 * .. * .. External Subroutines .. EXTERNAL DAXPY, DSCAL * .. * .. Executable Statements .. * INFO = 0 * DO 30 J = L, M JM1 = J - 1 ITRY = 0 10 CONTINUE * * Compute the 2-norm * NORM = DNRM2( N, Q( 1, J ), 1 ) IF( NORM.EQ.ZERO )THEN INFO = 1 RETURN END IF * * Scale the vector * REC = ONE / NORM CALL DSCAL( N, REC, Q( 1, J ), 1 ) * * Test to see if the J-th vector is orthogonal * IF( J.EQ.1 ) $ GO TO 30 IF( ITRY.EQ.0 ) $ NORM = ZERO IF( NORM.GT.TOL ) $ GO TO 30 ITRY = ITRY + 1 IF( ITRY.GT.MAXTRY )THEN INFO = 2 RETURN END IF * * Perform one step of the modified Gram-Schmidt process: * Q(:,J) = ( I - Q(:,1:J-1)*Q(:,1:J-1)')*Q(:,J) * DO 20 K = 1, JM1 QQ = DDOT( N, Q( 1, K ), 1, Q( 1, J ), 1 ) CALL DAXPY( N, -QQ, Q( 1, K ), 1, Q( 1, J ), 1 ) 20 CONTINUE GO TO 10 * 30 CONTINUE * RETURN * * End of DORTH END * ******** begin of dresid.f ******************************************* SUBROUTINE DRESID( N, L, M, Q, LDQ, AQ, LDA, T, LDT, RSD ) * .. * .. Scalar Arguments .. INTEGER L, LDA, LDQ, LDT, M, N * .. * .. Array Arguments .. DOUBLE PRECISION AQ( LDA, * ), Q( LDQ, * ), RSD( * ), $ T( LDT, * ) * .. * * Purpose * ======= * * Computes the column norms of residual vectors * * A*Q(1:N,L:M) - Q*T(1:M,L:M) * * where on entry, A*Q has been computed and stored in the array AQ. * * Arguments * ========= * * N (input) INTEGER * N is the order of the matrix A. * * L, M (input) INTEGER * The first and last column indices of Q. * * Q (input) DOUBLE PRECISION array, dimension( LDQ,M ) * On entry, Q contains approximate subspace matrix. * * LDQ (input) INTEGER * The leading dimension of Q, LDQ >= max(1,N). * * AQ (input) DOUBLE PRECISION array, dimension(LDA, M) * On entry, AQ contains the product A*Q. * * LDA (input) INTEGER * The leading dimension of AQ, LDA >= max(1,N). * * T (input) DOUBLE PRECISION array, dimension( LDT,M ) * On entry, T contains of Schur form. * * LDT (input) INTEGER * The leading dimension of T, LDT >= max(1,M). * * RSD (output) DOUBLE PRECISION array, dimension(M) * On exit, RSD(L) to RSD(M) contains the 2-norm of the residual * vectors. * * ====================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, TWO PARAMETER ( ZERO = 0.0D+0, TWO = 2.0D+0 ) * .. * .. Local Scalars .. INTEGER I, J, JNEXT, K, KU DOUBLE PRECISION S * .. * .. Intrinsic Functions .. INTRINSIC MIN, SQRT * .. * .. Executable Statements .. * * Quick return if necessary * IF( L.GT.M ) $ RETURN * * Compute the residual R = A*Q - Q*T, where A*Q has been * computed and stored in the array AQ. * DO 30 J = L, M KU = MIN( J+1, M ) IF( J.LT.M )THEN IF( T( J+1, J ).EQ.ZERO ) $ KU = J END IF * RSD( J ) = ZERO DO 20 I = 1, N S = ZERO DO 10 K = 1, KU S = S