C ALGORITHM 841, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 31, NO. 1, March, 2005, P. 166--185. #! /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/ # Doc/Makefile # Doc/README # Fortran/ # Fortran/Dp/ # Fortran/Dp/Drivers/ # Fortran/Dp/Drivers/Res1 # Fortran/Dp/Drivers/Res2 # Fortran/Dp/Drivers/balanc.f # Fortran/Dp/Drivers/data # Fortran/Dp/Drivers/driver1.f # Fortran/Dp/Drivers/driver2.f # Fortran/Dp/Drivers/interactive.f # Fortran/Dp/Src/ # Fortran/Dp/Src/dblas1.f # Fortran/Dp/Src/src.f # Matlab/ # Matlab/abh3.m # Matlab/compare.m # Matlab/mat.m # Matlab/mat2.m # This archive created: Thu May 12 16:46:57 2005 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' INCSRC = ../Src src.o: $(INCSRC)/src.f $(F77) $(F77OPTS) -c $(INCSRC)/src.f driver2.o: driver2.f $(F77) $(F77OPTS) -c driver2.f driver1.o: driver1.f $(F77) $(F77OPTS) -c driver1.f DRIVERS= driver1 driver2 Libs1= $(BLAS) Objs1= driver1.o src.o Objs2= driver2.o src.o all: res1 res2 driver1: $(Objs1) $(F77) $(F77OPTS) -o driver1 $(Objs1) $(Libs1) res1: driver1 data driver1 < data > res1 driver2: $(Objs2) $(F77) $(F77OPTS) -o driver2 $(Objs2) $(Libs1) res2: driver2 data driver2 < data > res2 diffres:Res1 res1 Res2 res2 echo "Differences in results from driver1" $(DIFF) Res1 res1 echo "Differences in results from driver2" $(DIFF) Res2 res2 clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'README' then echo shar: will not over-write existing file "'README'" else cat << "SHAR_EOF" > 'README' Please report bugs to howell@zach.fit.edu or gary.howell@hp.com The codes given here have no implicit or explicit warranty. They are free for personal or academic use. To extract the tar file, type gunzip -9c bhess.tar.gz | tar xvf - CONTENTS OF README 1. CONTENTS OF TAR FILE 2. CREATING AN EXECUTABLE TESTING PROGRAM TESTBH 3. MATLAB SCRIPTS 4. CHECKING THE EXECUTABLE AGAINST SAMPLE FILES 5. CREATING MATLAB OUTPUT THAT CAN BE CHECKED AGAINST YOUR EXECUTABLE 6. CHECKING OUTPUT FROM THE FORTRAN EXECUTABLE WITH MATLAB 7. TESTS WITHOUT MATLAB 8. OTHER MATERIALS 1. CONTENTS OF TAR FILE the tar file includes the subdirectories SRC and docs The SRC subdirectory contains BHESS.IN ! a sample input matrix README ! current file SAMPOUT ! output to compare to input matrix abh3.m ! a matlab file whose results can be compared ! to the Fortran results. Also the user can ! compare eigenvalues of the banded Hessenberg ! matrix produced by BHESS to those produced ! by LAPACK routines underlying the Matlab ! or octave eig routines. balanc.f ! balancing routine from EISPACK -- recommended ! preprocessing step for finding eigenvalues. bhess.f ! warning BHESS acts on a general square matrix of ! n columns but requires n+1 columns to be dimensioned ! as it uses a sparse column for a work vector. bhtest.f ! fortran files for testing program (use described below) compare.m ! a matlab script file (use described below) dblas1.f ! blas 1 routines. tuned libraries will give faster ! performance. makbhes ! a unix makefile to get an executable to test BHESS The docs subdirectory contains esub2acm.cls ! TOMS latex auxiliary file esub2acm.bst ! TOMS latex auxiliary file howell.ps -- A preprint of a TOMS paper on bhess (and figures) bandwidt.eps ! figures bdiff.eps blaserr.eps blastime.eps howell01.eps howell02.eps howell03.eps lanczos.eps randerr.eps randerr2.eps The tex file howell.tex may also be included. 2. CREATING AN EXECUTABLE TESTING PROGRAM BHTEST Users with UNIX can try to build a version of bhess by moving to the SRC subdirectory and typing make -f makbhes Other users will need to link the included FORTRAN files. These files are in FORTRAN 77. One minor extension used is the use of an ! mark to start comments in any column. This extension is allowable on most FORTRAN 77 compilers and on all FORTRAN 90 compilers. Performance will be greatly improved if tuned BLAS are available instead of the dblas1.f file. Also, typically a -O2 or -O3 compiler option will speed execution. For using BHESS with other programs, only the dblas1.f file and the bhess.f files are required. All of the bhess and auxiliary routines are in the file bhess.f. All start with the letters BH The code as given does not call balanc.f (the diagonal rebalancing routine from EISPACK) but it is a good idea to use such a routine as a preconditioner. The corresponding LAPACK routine is likely to be more efficient on modern machines but the time to call balanc.f is not very significant. 3. MATLAB SCRIPTS. abh3.m and compare.m are included. abh3.m runs from octave or matlab. abh3.m requires tol and n and an n times n matrix a. abh3 performs the bhess reduction to banded Hessenberg form b (albeit in a very inefficient way). Running abh3 from inside matlab or octave produces a banded Hessenberg matrix b by bhess. It also produces z such that b = float( inv(z) * a * z ) (A note on running the script files. When the script stalls in matlab, hit enter, when the script stalls in octave, hit ctrl c once, then f for forward). 4. CHECKING THE EXECUTABLE AGAINST SAMPLE FILES. Having created an executable, it is desirable to check the correctness of the output. BHESS.IN is a sample output file from octave. The file SAMPOUT has the input matrix a, the matrix b returned by BHESS, z the accumulated similarity transformation, and the results of multiplying z and inv(z) by the operations of the BHAP1, BHAP2, and BHAP3 respectively. A user can test (to 8 decimal digit accuracy) SAMPOUT against the screen output from the driver routine with input data from BHESS.IN 5. CREATING MATLAB OUTPUT THAT CAN BE CHECKED AGAINST YOUR EXECUTABLE First save your BHESS.IN file under a different name. Though very inefficient abh3 allows a user to check the proper operation of bhess. There is also documentation in abh3.m which can be accessed by typing help abh3 from inside octave or matlab. From octave for input to the Fortran exeutable save BHESS.IN tol n a (then edit the resulting file to get rid of extraneous nonnumeric lines) From matlab for input to the Fortran executable save BHESS.IN tol n a -ASCII Running the compiled file allows input from the BHESS.IN file. The defaults for matlab or octave put the data in single precision, so you have to declare format long (octave) to actually check accuracies. In matlab save with save BHESS.IN tol n a -ASCII -DOUBLE The included make file does not use optimized compiler options and does not link to an optimized blas library. For fast performance of this code or of LAPACK you need to link to a vendor supplied blas library. Optimizing compiler options are also a good idea. 6. CHECKING OUTPUT FROM THE FORTRAN EXECUTABLE WITH MATLAB There is an interactive option to produce a matlab file mat.m (though it is saved to only 14 digit accuracy) Before executing the driver, you should delete your old mat.m. The file mat.m will contain a randomly generated matrix a the size n, the tolerance n, and the banded Hessenberg matrix b produced. If you select the options for BHAP1, BHAP2, and BHAP3, mat.m will also contain onezinv = ones(1,n)*inv(z) z1 = z*ones(n,1) onesz = ones(1,n)*z From the MATLAB prompt, the sequence of commands mat abh3 ( need to hit space many times) compare will give some error numbers BackErr = norm(a - z*b*inv(z)) % an error internal to MATLAB ForErr = norm(B-b) % forward error between MATLAB b and FORTRAN B NormZ1Err= norm(z1-z*ones(n,1)) % vector produced by BHAP1 versus MATLAB Norm1ZErr= norm(onez-ones(1,n)*z) % vector produced by BHAP2 versus MATLAB Norm1InvZErr = norm(oneinvz-ones(1,n)*inv(z)) % vector procuced by BHAP3 versus MATLAB CondZ = co % Actual condition number of similarity transformation which you can compare to the estimated CondZ from BHAPC. (In defense of the poor formatting of the matlab-octave output, note that most screen output commands in matlab and octave are not compatible). 7. TESTS WITHOUT MATLAB On the next prompt after the decision whether or not to create a matlab readable file, the input option 2 allows backward error tests on a given matrix while varying the parameter tol. The driver program asks for TOL It then applies BHESS to get a banded Hessenberg matrix B. It then computes z*B*inv(z) and reports norm(a - z*B*inv(z)) (Frobenius norm) norm(a-z*B*inv(z))/norm(a) (Frobenius norm of relative backward error) estimated cond(z) bandwidth of B and prompts for another TOL This option allows the user to experiment with the tradeoff between TOL and bandwidth and rounding error. 8. OTHER MATERIALS For further information, you can access the anonymous ftp site cs.fit.edu cd pub/howell There is a CONTENTS file. A user interested in economical determination of eigenvalues by BHESS can often successfully apply the BR eigenvalue routine (found on the above anonymous ftp site and documented in a SIMAX 1999 paper the BR Eigenvalue Algorithm, by Geist, Howell, and Watkins). SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran' then mkdir 'Fortran' fi cd 'Fortran' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Res1' then echo shar: will not over-write existing file "'Res1'" else cat << "SHAR_EOF" > 'Res1' Maximal upper bandwidth of returned matrix= 1 Testing BHAP1 XRE(I) I BHPIVT(I) KROW(I) 0.1000000E+01 1 1 1 0.2063444E+00 2 2 2 0.5453415E+00 3 6 3 -0.1580819E+00 4 20 4 -0.3105722E-01 5 6 5 -0.1088136E+01 6 11 6 0.1570807E+01 7 20 7 -0.3868240E+00 8 16 8 -0.1327789E+01 9 16 9 0.4032291E+01 10 11 10 -0.3984943E+00 11 18 11 0.9841883E+00 12 13 12 0.1247343E+01 13 20 13 0.1461975E+01 14 16 14 0.2214046E+01 15 17 15 0.1028519E+01 16 20 16 0.2079648E+01 17 19 17 0.1603216E+01 18 20 18 0.1858040E+01 19 19 0 -0.1290508E+01 20 20 0 The Max Norm is = 4.0322913765407806 XRE(2)= 0.2063444389949818 XRE(SIZE/2)= 4.0322913765407806 Testing BHAP2 XRE(I) I BHPIVT(I) KROW(I) 0.1000000E+01 1 1 1 0.2178715E+01 2 2 2 0.1012945E+01 3 6 3 0.3751216E+00 4 20 4 0.3013311E+00 5 6 5 -0.5633626E-01 6 11 6 -0.1822702E-01 7 20 7 -0.4307611E-02 8 16 8 -0.9910564E-03 9 16 9 -0.2142212E-03 10 11 10 -0.4926650E-04 11 18 11 -0.6956075E-05 12 13 12 0.1726307E-05 13 20 13 -0.1325083E-06 14 16 14 0.2077042E-07 15 17 15 0.1279633E-07 16 20 16 -0.2835617E-07 17 19 17 0.8275181E-08 18 20 18 -0.7747061E-09 19 19 0 0.5439936E-08 20 20 0 The Max Norm is 2.1787150446352221 Testing BHAP3 XRE(I) I BHPIVT(I) KROW(I) 0.1000000E+01 1 1 1 0.5155619E-01 2 2 2 0.9042593E+00 3 6 3 0.3565937E+00 4 20 4 0.8889923E-01 5 6 5 -0.1072735E+01 6 11 6 0.1512145E+01 7 20 7 -0.2293525E+01 8 16 8 0.1604982E+00 9 16 9 0.9720993E+00 10 11 10 -0.5842520E+00 11 18 11 0.8639833E+00 12 13 12 0.1259311E+00 13 20 13 -0.1790023E+00 14 16 14 0.1048614E+01 15 17 15 0.2634614E+00 16 20 16 0.1222311E+01 17 19 17 0.1561692E+01 18 20 18 0.7452806E-01 19 19 0 -0.1289074E+01 20 20 0 The Max Norm is 2.2935251704858834 XRE(2)= 5.1556185575171165E-02 XRE(SIZE/2)= 0.9720992684485887 Call to BHAPC Y vector Y( 1)= 0.1000000E+01 Y( 2)=-0.5170456E+00 Y( 3)=-0.5760890E+00 Y( 4)=-0.3613120E+00 Y( 5)=-0.8089988E+00 Y( 6)=-0.1898542E+01 Y( 7)=-0.6501946E+00 Y( 8)=-0.2497253E+01 Y( 9)=-0.1827317E+01 Y( 10)= 0.2249096E+00 Y( 11)=-0.5831571E+01 Y( 12)=-0.3600122E+01 Y( 13)= 0.3735908E+01 Y( 14)= 0.2464788E+01 Y( 15)= 0.5650530E+01 Y( 16)= 0.2661804E+01 Y( 17)= 0.2048346E+01 Y( 18)= 0.2601913E+01 Y( 19)=-0.6495497E+00 Y( 20)=-0.1942904E+01 ESTIMATED CONDITION OF ACCUMULATED TRANSFORMATIONS = 34.0072177918201319 Testing BHUNPK Relevant part of BH array I J BH(I,J) 1 1 0.0000000E+00 1 2 0.1000000E+01 1 3 0.5000000E+00 2 1 0.1192326E+01 2 2 0.1188108E+01 2 3 0.1911253E+00 3 1 0.5314654E+00 3 2 0.2636027E+00 3 3 0.4121962E-01 4 1 0.3535205E-01 4 2 0.2555916E-01 4 3 0.9683867E-02 5 1 0.1319819E-02 5 2 0.2220766E-02 5 3 -0.1971623E-03 6 1 -0.4510177E-03 6 2 0.1704563E-03 6 3 0.2635024E-04 7 1 0.1756061E-04 7 2 0.1132605E-04 7 3 0.1286576E-05 8 1 0.1383043E-05 8 2 0.6482070E-06 8 3 0.7202418E-07 9 1 0.7003363E-07 9 2 0.3188191E-07 9 3 0.3340549E-08 10 1 0.3145448E-08 10 2 0.1343861E-08 10 3 0.1502516E-09 11 1 0.1063815E-09 11 2 0.4833281E-10 11 3 0.3325461E-11 12 1 0.5284099E-11 12 2 0.1473605E-11 12 3 -0.1958233E-12 13 1 -0.7018309E-13 13 2 0.3617281E-13 13 3 -0.1573459E-14 14 1 -0.3895955E-13 14 2 0.2233795E-14 14 3 -0.1364447E-15 15 1 -0.2053139E-14 15 2 0.1175780E-15 15 3 0.8037530E-16 16 1 0.1619328E-15 16 2 -0.2799151E-16 16 3 -0.8774360E-16 17 1 -0.8684909E-16 17 2 0.8395384E-16 17 3 -0.8338280E-16 18 1 0.3849247E-16 18 2 -0.1610948E-15 18 3 -0.1925781E-16 19 1 0.4759077E-16 19 2 0.1241711E-15 19 3 -0.5662216E-16 20 1 0.3423712E-17 20 2 -0.1222361E-15 Testing BHRPCK Relevant part of MAT0 array I J MAT0(I,J) 1 1 0.1000000E+01 1 2 0.5000000E+00 2 1 0.1192326E+01 2 2 0.1188108E+01 2 3 0.1911253E+00 3 2 0.5314654E+00 3 3 0.2636027E+00 3 4 0.4121962E-01 4 3 0.3535205E-01 4 4 0.2555916E-01 4 5 0.9683867E-02 5 4 0.1319819E-02 5 5 0.2220766E-02 5 6 -0.1971623E-03 6 5 -0.4510177E-03 6 6 0.1704563E-03 6 7 0.2635024E-04 7 6 0.1756061E-04 7 7 0.1132605E-04 7 8 0.1286576E-05 8 7 0.1383043E-05 8 8 0.6482070E-06 8 9 0.7202418E-07 9 8 0.7003363E-07 9 9 0.3188191E-07 9 10 0.3340549E-08 10 9 0.3145448E-08 10 10 0.1343861E-08 10 11 0.1502516E-09 11 10 0.1063815E-09 11 11 0.4833281E-10 11 12 0.3325461E-11 12 11 0.5284099E-11 12 12 0.1473605E-11 12 13 -0.1958233E-12 13 12 -0.7018309E-13 13 13 0.3617281E-13 13 14 -0.1573459E-14 14 13 -0.3895955E-13 14 14 0.2233795E-14 14 15 -0.1364447E-15 15 14 -0.2053139E-14 15 15 0.1175780E-15 15 16 0.8037530E-16 16 15 0.1619328E-15 16 16 -0.2799151E-16 16 17 -0.8774360E-16 17 16 -0.8684909E-16 17 17 0.8395384E-16 17 18 -0.8338280E-16 18 17 0.3849247E-16 18 18 -0.1610948E-15 18 19 -0.1925781E-16 19 18 0.4759077E-16 19 19 0.1241711E-15 19 20 -0.5662216E-16 20 19 0.3423712E-17 20 20 -0.1222361E-15 SHAR_EOF fi # end of overwriting check if test -f 'Res2' then echo shar: will not over-write existing file "'Res2'" else cat << "SHAR_EOF" > 'Res2' n = 20 tol = 0.100000E+01 Frobenius Backward Error Norm = 0.257087E-13 Relative Backward Error Norm = 0.583674E-13 Estimated Condition Number = 0.340072E+02 Bandwidth = 1 SHAR_EOF fi # end of overwriting check if test -f 'balanc.f' then echo shar: will not over-write existing file "'balanc.f'" else cat << "SHAR_EOF" > 'balanc.f' subroutine balanc(nm,n,a,low,igh,scale) c c This was modified! (Howell and Watkins June 1995 so that low, igh c not changed) c Statements low = k and igh = l (at end) have been made into comments. c integer i,j,k,l,m,n,jj,nm,igh,low,iexc double precision a(nm,n),scale(n) double precision c,f,g,r,s,b2,radix logical noconv c c this subroutine is a translation of the algol procedure balance, c num. math. 13, 293-304(1969) by parlett and reinsch. c handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). c c this subroutine balances a real matrix and isolates c eigenvalues whenever possible. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c a contains the input matrix to be balanced. c c on output c c a contains the balanced matrix. c c low and igh are two integers such that a(i,j) c is equal to zero if c (1) i is greater than j and c (2) j=1,...,low-1 or i=igh+1,...,n. c c scale contains information determining the c permutations and scaling factors used. c c suppose that the principal submatrix in rows low through igh c has been balanced, that p(j) denotes the index interchanged c with j during the permutation step, and that the elements c of the diagonal matrix used are denoted by d(i,j). then c scale(j) = p(j), for j = 1,...,low-1 c = d(j,j), j = low,...,igh c = p(j) j = igh+1,...,n. c the order in which the interchanges are made is n to igh+1, c then 1 to low-1. c c note that 1 is returned for igh if igh is zero formally. c c the algol procedure exc contained in balance appears in c balanc in line. (note that the algol roles of identifiers c k,l have been reversed.) c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c radix = 16.0d0 c b2 = radix * radix k = 1 l = n go to 100 c .......... in-line procedure for row and c column exchange .......... 20 scale(m) = j if (j .eq. m) go to 50 c do 30 i = 1, l f = a(i,j) a(i,j) = a(i,m) a(i,m) = f 30 continue c do 40 i = k, n f = a(j,i) a(j,i) = a(m,i) a(m,i) = f 40 continue c 50 go to (80,130), iexc c .......... search for rows isolating an eigenvalue c and push them down .......... 80 if (l .eq. 1) go to 280 l = l - 1 c .......... for j=l step -1 until 1 do -- .......... 100 do 120 jj = 1, l j = l + 1 - jj c do 110 i = 1, l if (i .eq. j) go to 110 if (a(j,i) .ne. 0.0d0) go to 120 110 continue c m = l iexc = 1 go to 20 120 continue c go to 140 c .......... search for columns isolating an eigenvalue c and push them left .......... 130 k = k + 1 c 140 do 170 j = k, l c do 150 i = k, l if (i .eq. j) go to 150 if (a(i,j) .ne. 0.0d0) go to 170 150 continue c m = k iexc = 2 go to 20 170 continue c .......... now balance the submatrix in rows k to l .......... do 180 i = k, l 180 scale(i) = 1.0d0 c .......... iterative loop for norm reduction .......... 190 noconv = .false. c do 270 i = k, l c = 0.0d0 r = 0.0d0 c do 200 j = k, l if (j .eq. i) go to 200 c = c + dabs(a(j,i)) r = r + dabs(a(i,j)) 200 continue c .......... guard against zero c or r due to underflow .......... if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 g = r / radix f = 1.0d0 s = c + r 210 if (c .ge. g) go to 220 f = f * radix c = c * b2 go to 210 220 g = r * radix 230 if (c .lt. g) go to 240 f = f / radix c = c / b2 go to 230 c .......... now balance .......... 240 if ((c + r) / f .ge. 0.95d0 * s) go to 270 g = 1.0d0 / f scale(i) = scale(i) * f noconv = .true. c do 250 j = k, n 250 a(i,j) = a(i,j) * g c do 260 j = 1, l 260 a(j,i) = a(j,i) * f c 270 continue c if (noconv) go to 190 c 280 continue c low = k c igh = l return end 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' 1 20 1 0.5 0.333333333333333 0.25 0.2 0.166666666666667 0.142857142857143 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.5 0.333333333333333 0.25 0.2 0.166666666666667 0.142857142857143 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.333333333333333 0.25 0.2 0.166666666666667 0.142857142857143 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.25 0.2 0.166666666666667 0.142857142857143 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.2 0.166666666666667 0.142857142857143 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.166666666666667 0.142857142857143 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.142857142857143 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.125 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.111111111111111 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.1 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0909090909090909 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.0833333333333333 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.0769230769230769 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0714285714285714 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0303030303030303 0.0666666666666667 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0303030303030303 0.0294117647058824 0.0625 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0303030303030303 0.0294117647058824 0.0285714285714286 0.0588235294117647 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0303030303030303 0.0294117647058824 0.0285714285714286 0.0277777777777778 0.0555555555555556 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0303030303030303 0.0294117647058824 0.0285714285714286 0.0277777777777778 0.027027027027027 0.0526315789473684 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0303030303030303 0.0294117647058824 0.0285714285714286 0.0277777777777778 0.027027027027027 0.0263157894736842 0.05 0.0476190476190476 0.0454545454545455 0.0434782608695652 0.0416666666666667 0.04 0.0384615384615385 0.037037037037037 0.0357142857142857 0.0344827586206897 0.0333333333333333 0.032258064516129 0.03125 0.0303030303030303 0.0294117647058824 0.0285714285714286 0.0277777777777778 0.027027027027027 0.0263157894736842 0.0256410256410256 SHAR_EOF fi # end of overwriting check if test -f 'driver1.f' then echo shar: will not over-write existing file "'driver1.f'" else cat << "SHAR_EOF" > 'driver1.f' PROGRAM DRIVER *------------------------------------------------------- * Update History: * July 9, 2004 by G. Howell * Nov. 21, 2002 by N. Diaa * Nov. 20, 2002 by N. Diaa *----------------------------------------------------------------- * * matrix entries are read from the file data * *----------------------------------------------------------------- C .. Parameters .. INTEGER LDX PARAMETER (LDX=21) C .. C .. Scalars in Common .. INTEGER SEED C .. C .. Local Scalars .. DOUBLE PRECISION CONDIT,SIZE2,TOL,XNORM1,XNORM2,XNORM3 INTEGER AVGNU,I,IFLAG,INFO,J,M,MAXNU,NUMAX,SIZE C .. C .. Local Arrays .. DOUBLE PRECISION BH(LDX,LDX),DUMV(LDX),MAT0(LDX,LDX),XIM(LDX), + XRE(LDX) INTEGER BHPIVT(LDX),KROW(LDX),KRVECT(LDX),NU(LDX),NUMON(LDX) C .. C .. External Subroutines .. EXTERNAL BHAP1,BHAP2,BHAP3,BHAPC,BHESS,BHRPCK,BHUNPK C .. C .. Intrinsic Functions .. INTRINSIC ABS,INT,MAX,MIN C .. C .. Common blocks .. COMMON SEED C .. OPEN (33,FILE='data',STATUS='UNKNOWN',ERR=30) READ (33,FMT=*) TOL READ (33,FMT=*) SIZE2 SIZE = INT(SIZE2) IF ((SIZE.GT.LDX) .OR. (SIZE.LE.2) .OR. (TOL.LE.0)) THEN WRITE(*,*)' ERROR: Invalid data in data ' WRITE(*,*) + ' Real number Tol or integer Size in data ' WRITE(*,*)' Tol must be real number > 0 ' WRITE(*,*) + ' Size must be integer: 2 < size <= ',LDX,' ' CLOSE (33) GO TO 30 END IF DO I = 1,SIZE READ (33,FMT=*) (MAT0(I,J),J=1,SIZE) END DO CLOSE (33) *-------------------------------------------------------------- * * A good general policy is to call a balancing routine * (using diagonal similarity transformations) * before BHESS (or any reduction to Hessenberg form). * A standard routine is balanc.f from EISPACK. * *-------------------------------------------------------------- * Now we have the matrix MAT0, * either read from data or randomly generated *-------------------------------------------------------------- 10 CONTINUE CALL BHESS(LDX,MAT0,SIZE,BHPIVT,NU,KRVECT,KROW,DUMV,TOL,INFO) IF (INFO.NE.0) GO TO 30 20 CONTINUE *-------------------------------------------------------------- * * The average bandwidth is computed as AVGNU * *-------------------------------------------------------------- AVGNU = 0 MAXNU = 0 DO I = 1,SIZE AVGNU = AVGNU + NU(I) MAXNU = MAX(MAXNU,NU(I)) END DO AVGNU = (AVGNU+SIZE/2)/SIZE WRITE(*,*)' Maximal upper bandwidth of returned matrix=', + MAXNU *---------------------------------------------------------------- * * Test BHAP1 * WRITE(*,*)' Testing BHAP1' DO I = 1,SIZE XRE(I) = 1.D0 XIM(I) = 0.D0 END DO CALL BHAP1(SIZE,MAT0,LDX,XRE,XIM,BHPIVT,KROW,IFLAG) WRITE(*,*)' XRE(I) I BHPIVT(I) KROW(I)' XNORM1 = 0.0D0 DO I = 1,SIZE WRITE(*,'(1x,e14.7, 2i5,7x,i5)') + XRE(I),I,BHPIVT(I),KROW(I) XNORM1 = MAX(XNORM1,ABS(XRE(I))) END DO WRITE(*,*)' The Max Norm is = ',XNORM1 WRITE(*,*)' XRE(2)= ',XRE(2) WRITE(*,*)' XRE(SIZE/2)= ',XRE(SIZE/2) WRITE(*,*)' ' *---------------------------------------------------------------- * * Test BHAP2 * WRITE(*,*)' Testing BHAP2' DO I = 1,SIZE XRE(I) = 1.D0 END DO CALL BHAP2(SIZE,MAT0,LDX,XRE,BHPIVT,KROW,IFLAG) XNORM2 = 0.D0 WRITE(*,*)' XRE(I) I BHPIVT(I) KROW(I)' DO I = 1,SIZE WRITE(*,'(1x,e14.7, 2i5,7x,i5)') + XRE(I),I,BHPIVT(I),KROW(I) XNORM2 = MAX(XNORM2,ABS(XRE(I))) END DO WRITE(*,*)' ' WRITE(*,*)' The Max Norm is ',XNORM2 WRITE(*,*)' ' WRITE(*,*)' ' *---------------------------------------------------------------- * * Test BHAP3 * *---------------------------------------------------------------- WRITE(*,*) ' Testing BHAP3' DO I = 1,SIZE XRE(I) = 1.D0 XIM(I) = 0.D0 END DO CALL BHAP3(SIZE,MAT0,LDX,XRE,XIM,BHPIVT,KROW,IFLAG) XNORM3 = 0.D0 WRITE(*,*)' XRE(I) I BHPIVT(I) KROW(I)' DO I = 1,SIZE WRITE(*,'(1x,e14.7, 2i5,7x,i5)') + XRE(I),I,BHPIVT(I),KROW(I) XNORM3 = MAX(XNORM3,ABS(XRE(I))) END DO WRITE(*,*)'The Max Norm is ',XNORM3 WRITE(*,*)' XRE(2)= ',XRE(2) WRITE(*,*)' XRE(SIZE/2)= ',XRE(SIZE/2) * --------------------------------------------------------------- * * The next statements call BHAPC -- a LINPACK style * estimator of the infinity norm of inv(Z) * * --------------------------------------------------------------- WRITE(*,'(// ''Call to BHAPC'')') CALL BHAPC(SIZE,MAT0,LDX,XRE,XIM,BHPIVT,KROW,CONDIT,IFLAG) WRITE(*,*) WRITE(*,*)' Y vector' DO I = 1,SIZE WRITE(*,'(1x,a,i5,a,e14.7)')'Y(',I,')=',XRE(I) END DO WRITE(*,*)' ' WRITE(*,*)' ESTIMATED CONDITION OF ACCUMULATED' WRITE(*,*)' TRANSFORMATIONS =',CONDIT WRITE(*,*)' ' * --------------------------------------------------------------- * * Test BHUNPK or BHRPCK * WRITE(*,*) ' Testing BHUNPK' CALL BHUNPK(LDX,MAT0,BH,SIZE,NU,NUMON,NUMAX) WRITE(*,*) ' Relevant part of BH array' WRITE(*,*) WRITE(*,*)' I J BH(I,J)' DO I = 1,SIZE WRITE(*,'(1x,2i5,2x,e14.7)') (I,J,BH(I,J),J=1,NU(I)+2) END DO * -------------------------------------------- * zero out multipliers in atobhes matrix MAT0 * -------------------------------------------- WRITE(*,*) ' Testing BHRPCK' CALL BHRPCK(LDX,SIZE,MAT0,NU) WRITE(*,*) ' Relevant part of MAT0 array' WRITE(*,*) WRITE(*,*)' I J MAT0(I,J)' DO I = 1,SIZE M = MIN(SIZE,NU(I)+2) DO J = 1,M IF (I+J-2.NE.0) THEN WRITE(*,'(1x,2i5,2x,e14.7)') + I,I + J - 2,MAT0(I,I+J-2) END IF END DO END DO 30 CONTINUE END * ----------------------------------------------- * *==================================================== SUBROUTINE BHBACK(LDX,N,A,B,W,V,TOL,NU,PIV,KRVECT,KROW) * This subroutine tests backward error of BHESS * by running BHESS to get a Hessenberg matrix * then using the stored multipliers to regenerate * the original matrix. * * Outputs are the difference between the original matrix * and the reconstruction (computed backward error) * * The upper bandwidth of the returned matrix and the estimated * condition number of the similariy transformation are * also returned. * * Arguments *------------------- * * N -- Integer Number of rows and columns of input matrices * * LDX -- Integer Leading Dimension of Variable dimension arrays * * A -- Double Precision Input matrix unchanged * * B and W -- Double Precision Matrices, B must of at * least N+1 columns. On return B is an approximation * of the input matrix A (to backward error) and W * holds the output from BHESS. * * TOL -- Double precision parameter to control multiplier * size in BHESS. Unchanged. * * V -- Double precision work vector of three columns. * * Local Variables *------------------------------- * * Copy A to B * C .. Scalar Arguments .. DOUBLE PRECISION TOL INTEGER LDX,N C .. C .. Array Arguments .. DOUBLE PRECISION A(LDX,*),B(LDX,*),V(LDX,3),W(LDX,*) INTEGER KROW(*),KRVECT(*),NU(*),PIV(*) C .. C .. Local Scalars .. DOUBLE PRECISION CONDIT,FNORM,RNORM,SUM,SUM2 REAL SUMA INTEGER I,IFLAG,INFO,J,MXBAND C .. C .. External Subroutines .. EXTERNAL BHAP1,BHAP3,BHAPC,BHESS,BHRPCK,DCOPY C .. C .. Intrinsic Functions .. INTRINSIC DSQRT,MAX C .. DO J = 1,N DO I = 1,N B(I,J) = A(I,J) END DO END DO *-------------------------------------------------------- * * Call BHESS with B Warning B requires N + 1 Columns * *-------------------------------------------------------- CALL BHESS(LDX,B,N,PIV,NU,KRVECT,KROW,V,TOL,INFO) * * Copy B to W * DO J = 1,N DO I = 1,N W(I,J) = B(I,J) END DO END DO * * Zero out all but upper Hessenberg in B * CALL BHRPCK(LDX,N,B,NU) * * Calculate backward A as B = Z * B * inv(Z) * * First B <-- Z * B DO I = 1,N - 1,2 CALL BHAP3(N,W,LDX,B(1,I),B(1,I+1),PIV,KROW,IFLAG) END DO IF (N.NE.2* (N/2)) THEN CALL BHAP3(N,W,LDX,B(1,N),V(1,2),PIV,KROW,IFLAG) END IF * * Then B <-- B * inv(Z) So B approximates A * DO I = 1,N - 1,2 CALL DCOPY(N,B(I,1),LDX,V(1,1),1) CALL DCOPY(N,B(I+1,1),LDX,V(1,2),1) CALL BHAP1(N,W,LDX,V(1,1),V(1,2),PIV,KROW,IFLAG) CALL DCOPY(N,V(1,1),1,B(I,1),LDX) CALL DCOPY(N,V(1,2),1,B(I+1,1),LDX) END DO IF (N.NE.2* (N/2)) THEN CALL DCOPY(N,B(N,1),LDX,V(1,1),1) CALL BHAP1(N,W,LDX,V(1,1),V(1,2),PIV,KROW,IFLAG) CALL DCOPY(N,V(1,1),1,B(N,1),LDX) END IF * * Calculate || A - B || * SUM = 0.D0 SUMA = 0.D0 DO I = 1,N DO J = 1,N SUM = SUM + (A(I,J)-B(I,J))**2 SUM2 = SUM2 + A(I,J)**2 END DO END DO * Frobenius norm normalized so that norm(eye(n)) = 1 FNORM = DSQRT(SUM/N) * Relative Frobenius norm RNORM = DSQRT(SUM/SUM2) * * Estimate conditioning of Z * CALL BHAPC(N,W,LDX,V(1,1),V(1,2),PIV,KROW,CONDIT,IFLAG) * * Calculate Maximal Bandwidth * MXBAND = 1 DO I = 1,N MXBAND = MAX(NU(I),MXBAND) END DO WRITE(*,*)' ' WRITE(*,*)' ' WRITE(*,*)' n = ',N WRITE(*,*)' tol = ',TOL WRITE(*,*)' Frobenius Backward Error Norm = ',FNORM WRITE(*,*)' Relative Backward Error Norm = ',RNORM WRITE(*,*)' Estimated Condition Number = ',CONDIT WRITE(*,*)' Bandwidth = ',MXBAND RETURN END SHAR_EOF fi # end of overwriting check if test -f 'driver2.f' then echo shar: will not over-write existing file "'driver2.f'" else cat << "SHAR_EOF" > 'driver2.f' PROGRAM DRIVER *------------------------------------------------------- * Update History: * July 9, 2004 by G. Howell * Nov. 21, 2002 by N. Diaa * Nov. 20, 2002 by N. Diaa *------------------------------------------------------- C .. Parameters .. INTEGER LDX PARAMETER (LDX=500) C .. C .. Scalars in Common .. INTEGER SEED C .. C .. Local Scalars .. DOUBLE PRECISION SIZE2,TOL INTEGER I,J,SIZE C .. C .. Local Arrays .. DOUBLE PRECISION BH(LDX,LDX),MAT0(LDX,LDX),V(LDX,3),WMAT(LDX,LDX) INTEGER BHPIVT(LDX),KROW(LDX),KRVECT(LDX),NU(LDX) C .. C .. External Subroutines .. EXTERNAL BHBACK C .. C .. Intrinsic Functions .. INTRINSIC INT C .. C .. Common blocks .. COMMON SEED C .. OPEN (33,FILE='data',STATUS='UNKNOWN',ERR=10) READ (33,FMT=*) TOL READ (33,FMT=*) SIZE2 SIZE = INT(SIZE2) IF ((SIZE.GT.LDX) .OR. (SIZE.LE.2) .OR. (TOL.LE.0)) THEN WRITE(*,*)' ERROR: Invalid data in data ' WRITE(*,*) +' Real number Tol or integer Size in data ' WRITE(*,*)' Tol must be real number > 0 ' WRITE(*,*) +' Size must be integer: 2 < size <= ',LDX,' ' CLOSE (33) GO TO 10 END IF DO I = 1,SIZE READ (33,FMT=*) (MAT0(I,J),J=1,SIZE) END DO CLOSE (33) *-------------------------------------------------------------- * * A good general policy is to call a balancing routine * (using diagonal similarity transformations) * before BHESS (or any reduction to Hessenberg form). * A standard routine is balanc.f from EISPACK. * *-------------------------------------------------------------- * Now we have the matrix MAT0, * either read from data or randomly generated *-------------------------------------------------------------- CALL BHBACK(LDX,SIZE,MAT0,BH,WMAT,V,TOL,NU,BHPIVT,KRVECT,KROW) 10 END *==================================================== SUBROUTINE BHBACK(LDX,N,A,B,W,V,TOL,NU,PIV,KRVECT,KROW) * This subroutine tests backward error of BHESS * by running BHESS to get a Hessenberg matrix * then using the stored multipliers to regenerate * the original matrix. * * Outputs are the difference between the original matrix * and the reconstruction (computed backward error) * * The upper bandwidth of the returned matrix and the estimated * condition number of the similariy transformation are * also returned. * * Arguments *------------------- * * N -- Integer Number of rows and columns of input matrices * * LDX -- Integer Leading Dimension of Variable dimension arrays * * A -- Double Precision Input matrix unchanged * * B and W -- Double Precision Matrices, B must of at * least N+1 columns. On return B is an approximation * of the input matrix A (to backward error) and W * holds the output from BHESS. * * TOL -- Double precision parameter to control multiplier * size in BHESS. Unchanged. * * V -- Double precision work vector of three columns. * * Local Variables *------------------------------- * * Copy A to B * C .. Scalar Arguments .. DOUBLE PRECISION TOL INTEGER LDX,N C .. C .. Array Arguments .. DOUBLE PRECISION A(LDX,*),B(LDX,*),V(LDX,3),W(LDX,*) INTEGER KROW(*),KRVECT(*),NU(*),PIV(*) C .. C .. Local Scalars .. DOUBLE PRECISION CONDIT,FNORM,RNORM,SUM,SUM2 INTEGER I,IFLAG,INFO,J,MXBAND C .. C .. External Subroutines .. EXTERNAL BHAP1,BHAP3,BHAPC,BHESS,BHRPCK,DCOPY C .. C .. Intrinsic Functions .. INTRINSIC DSQRT,MAX C .. DO J = 1,N DO I = 1,N B(I,J) = A(I,J) END DO END DO *-------------------------------------------------------- * * Call BHESS with B Warning B requires N + 1 Columns * *-------------------------------------------------------- CALL BHESS(LDX,B,N,PIV,NU,KRVECT,KROW,V,TOL,INFO) * * Copy B to W * DO J = 1,N DO I = 1,N W(I,J) = B(I,J) END DO END DO * * Zero out all but upper Hessenberg in B * CALL BHRPCK(LDX,N,B,NU) * * Calculate backward A as B = Z * B * inv(Z) * * First B <-- Z * B DO I = 1,N - 1,2 CALL BHAP3(N,W,LDX,B(1,I),B(1,I+1),PIV,KROW,IFLAG) END DO IF (N.NE.2* (N/2)) THEN CALL BHAP3(N,W,LDX,B(1,N),V(1,2),PIV,KROW,IFLAG) END IF * * Then B <-- B * inv(Z) So B approximates A * DO I = 1,N - 1,2 CALL DCOPY(N,B(I,1),LDX,V(1,1),1) CALL DCOPY(N,B(I+1,1),LDX,V(1,2),1) CALL BHAP1(N,W,LDX,V(1,1),V(1,2),PIV,KROW,IFLAG) CALL DCOPY(N,V(1,1),1,B(I,1),LDX) CALL DCOPY(N,V(1,2),1,B(I+1,1),LDX) END DO IF (N.NE.2* (N/2)) THEN CALL DCOPY(N,B(N,1),LDX,V(1,1),1) CALL BHAP1(N,W,LDX,V(1,1),V(1,2),PIV,KROW,IFLAG) CALL DCOPY(N,V(1,1),1,B(N,1),LDX) END IF * * Calculate || A - B || * SUM = 0.D0 SUM2 = 0.D0 DO I = 1,N DO J = 1,N SUM = SUM + (A(I,J)-B(I,J))**2 SUM2 = SUM2 + A(I,J)**2 END DO END DO * Frobenius norm normalized so that norm(eye(n)) = 1 FNORM = DSQRT(SUM/N) * Relative Frobenius norm RNORM = DSQRT(SUM/SUM2) * * Estimate conditioning of Z * CALL BHAPC(N,W,LDX,V(1,1),V(1,2),PIV,KROW,CONDIT,IFLAG) * * Calculate Maximal Bandwidth * MXBAND = 1 DO I = 1,N MXBAND = MAX(NU(I),MXBAND) END DO WRITE(*,*)' ' WRITE(*,*)' ' WRITE(*,*)' n = ',N WRITE(*,*)' tol = ',TOL WRITE(*,*)' Frobenius Backward Error Norm = ',FNORM WRITE(*,*)' Relative Backward Error Norm = ',RNORM WRITE(*,*)' Estimated Condition Number = ',CONDIT WRITE(*,*)' Bandwidth = ',MXBAND RETURN END SHAR_EOF fi # end of overwriting check if test -f 'interactive.f' then echo shar: will not over-write existing file "'interactive.f'" else cat << "SHAR_EOF" > 'interactive.f' PROGRAM DRIVER *------------------------------------------------------- * Update History: * July 9, 2004 by G. Howell * Nov. 21, 2002 by N. Diaa * Nov. 20, 2002 by N. Diaa *------------------------------------------------------- C .. Parameters .. INTEGER LDX PARAMETER (LDX=500) C .. C .. Scalars in Common .. INTEGER SEED C .. C .. Local Scalars .. DOUBLE PRECISION CONDIT,SIZE2,TOL,XNORM1,XNORM2,XNORM3 INTEGER AVGNU,I,ICHOIC,IDUM,IFLAG,IGO,IMAT,INFO,INIXSIX,ITEMP, + ITEST,ITESTS,J,M,MAXNU,NUMAX,PRLEVEL,SIZE C .. C .. Local Arrays .. DOUBLE PRECISION BH(LDX,LDX),DUMV(LDX),MAT0(LDX,LDX),V(LDX,3), + WMAT(LDX,LDX),XIM(LDX),XRE(LDX) INTEGER BHPIVT(LDX),KROW(LDX),KRVECT(LDX),NU(LDX),NUMON(LDX) C .. C .. External Functions .. DOUBLE PRECISION RANDOM EXTERNAL RANDOM C .. C .. External Subroutines .. EXTERNAL BHAP1,BHAP2,BHAP3,BHAPC,BHBACK,BHESS,BHRPCK,BHUNPK C .. C .. Intrinsic Functions .. INTRINSIC ABS,INT,MAX,MIN C .. C .. Common blocks .. COMMON SEED C .. PRINT *,' ' PRINT *,' BHESS TESTER. ' PRINT *,' Input an integer on each pause ' PRINT *,' ' PRINT *,' BHESS reduces a general square matrix' PRINT *,' to similar small-band Hessenberg form' PRINT *,' using multipliers bounded in root mean' PRINT *,' by a user-specified bound tol .' PRINT *,' ' PRINT *,' Taking tol = 30 typically results in ' PRINT *,' an almost tridiagonal matrix from' PRINT *,' which eigenvalues can be extracted' PRINT *,' by BR iteration.' PRINT *,' ' PRINT *,' DRIVER calls BHESS and some auxiliary' PRINT *,' routines. It can take input files produced' PRINT *,' by Matlab or Octave and produce files readable' PRINT *,' by those programs, allowing verification of' PRINT *,' correct computation.' PRINT *,' ' PRINT *,' In response to system prompts' PRINT *,' enter integers.' PRINT *,' There is additional info in README.' PRINT *,' ' PRINT *,' ' 10 CONTINUE PRINT *,' Input an integer PRINT LEVEL: ' PRINT *,' 0 will give no screen output.' PRINT *,' 1 will give maximal output.' READ *,PRLEVEL IF ((PRLEVEL.NE.0) .AND. (PRLEVEL.NE.1)) GO TO 10 20 CONTINUE PRINT *,' Choose:' PRINT *,' 1) to randomly generate input matrix ' PRINT *,' 2) to read input matrix from file BHESS.IN ' PRINT *,' 3) to quit ' INIXSIX = 0 IMAT = 0 ITEMP = 0 READ *,ICHOIC IF (ICHOIC.EQ.3) GO TO 110 IF ((ICHOIC.NE.1) .AND. (ICHOIC.NE.2)) GO TO 20 PRINT *,' Output is to screen' *---------------------------------------------------------------- * * For ICHOIC = 1, an n x n matrix is generated with a given seed. * using a linear congruential random number generator. Entries of * the matrix are from a uniform distribution between -1 and 1. * *--------------------------------------------------------------- * 1) to randomly generate input matrix IF (ICHOIC.EQ.1) THEN 30 CONTINUE PRINT *,' ' PRINT *,' Input integer size (2 < size <= ',LDX,' ): ' READ *,SIZE IF ((SIZE.GT.LDX) .OR. (SIZE.LE.2)) GO TO 30 40 CONTINUE PRINT *,' ' PRINT *,'Input integer random generator seed > 0:' READ *,SEED IF (SEED.LE.0) GO TO 40 DO J = 1,SIZE DO I = 1,SIZE MAT0(I,J) = (RANDOM(SEED)-.5)*2.D0 END DO END DO END IF *----------------------------------------------------------------- * * If ICHOIC is 2, then matrix entries are read * from the file BHESS.IN * *----------------------------------------------------------------- * * To create the file BHESS.IN from matlab, * assign a tolerance tol, a size n, and a square matrix A. * Use the matlab command * save BHESS.IN tol n A -ascii * * To create the file BHESS.IN from the matlab clone octave, * assign the same variables and use the comman * save BHESS.IN tol n A * But you then have to edit the resulting file * to remove the nonnumeric lines. * * Typing format long inside octave results in more digits * being saved. * * To produce BHESS.IN from matlab * assign tolerance tol, size n, square matrix A * Use the matlab command * Save BHESS.IN tol n A -ascii * Enter an arbitary integer to continue * Accuracy is * constrained by the 8 digits in the ascii * format used to write the file BHESS.IN * from Matlab * In Octave the produced file has extraneous * characters which will need to be deleted * before BHESS.IN can be used as an input file. * Typing * format long * before saving will give more digits in * the saved file. *-------------------------------------------------------------- * ICHOIC=2 read matrix from BHESS.IN IF (ICHOIC.EQ.2) THEN PRINT *,'You chose to read input matrix from file BHESS.IN' PRINT *,'This option assumes you have a file BHESS.IN' PRINT *,'To see how to produce this file see the ' PRINT *,'README file or the comments in abh3.m ' PRINT *,' Input any integer to continue' READ *,IDUM OPEN (33,FILE='BHESS.IN',STATUS='UNKNOWN',ERR=100) READ (33,FMT=*) TOL READ (33,FMT=*) SIZE2 SIZE = INT(SIZE2) IF ((SIZE.GT.LDX) .OR. (SIZE.LE.2) .OR. (TOL.LE.0)) THEN PRINT *,' ERROR: Invalid data in BHESS.IN ' PRINT *, + ' Real number Tol or integer Size in BHESS.IN ' PRINT *,' Tol must be real number > 0 ' PRINT *,' Size must be integer: 2 < size <= ',LDX, + ' ' CLOSE (33) GO TO 110 END IF DO I = 1,SIZE PRINT *,' Reading row ',I,' of ',SIZE READ (33,FMT=*) (MAT0(I,J),J=1,SIZE) END DO CLOSE (33) END IF *-------------------------------------------------------------- * * A good general policy is to call a balancing routine * (using diagonal similarity transformations) * before BHESS (or any reduction to Hessenberg form). * A standard routine is balanc.f from EISPACK. * *-------------------------------------------------------------- * Now we have the matrix MAT0, * either read from BHESS.IN or randomly generated *-------------------------------------------------------------- 50 CONTINUE PRINT *,' Choose set of tests ' PRINT *,' ' PRINT *,' Input 1 to get output that can be compared ' PRINT *,' to Matlab and illustrates the actions of the' PRINT *,' included subroutines' PRINT *,' ' PRINT *,' Input 2 to test the relation of the input ' PRINT *,' parameter TOL on computed backward error,' PRINT *,' resulting bandwidth, and estimated condition' PRINT *,' number of the similarity transformation' READ *,ITESTS IF ((ITESTS.NE.1) .AND. (ITESTS.NE.2)) GO TO 50 IF (ITESTS.EQ.2) THEN 60 CONTINUE PRINT *,' Input TOL as a parameter for calling BHESS' PRINT *,' TOL is the maximal allowable root mean square' PRINT *,' for vectors of paired row-column multipliers' PRINT *,' TOL should be a real number bigger than zero' READ *,TOL IF (TOL.LE.0) GO TO 60 CALL BHBACK(LDX,SIZE,MAT0,BH,WMAT,V,TOL,NU,BHPIVT,KRVECT,KROW) PRINT *,' ' PRINT *,' ' PRINT *,' Input 1 to do another round of tests with the' PRINT *,' same matrix and a new input TOL, ' PRINT *,' Else other integers quit. ' READ *,ITESTS IF (ITESTS.EQ.1) GO TO 60 GO TO 110 END IF * Choice of randomly generated matrix with some explicit tests of * routines. This choice may not be very useful. IF (ICHOIC.EQ.1) THEN PRINT *,' This choice allows you to get a feel for ' PRINT *,' available subroutines, but does not allow' PRINT *,' easy tests of validity. ' PRINT *,' For a test of BHESS validity, either use matrices' PRINT *,' exported from matlab and compare the results here' PRINT *,' to those from running abh3.m ' PRINT *,' OR ' PRINT *,' start over and choose Test set 2' PRINT *,' on the step just before this message.' PRINT *,' ' 70 CONTINUE PRINT *,' Input Tolerance real number greater than zero ' PRINT *,' TOL bounds the product of the root mean squares' PRINT *,' of vectors of column and row multipliers ' READ *,TOL IF (TOL.LE.0) GO TO 70 80 CONTINUE PRINT *,' ' PRINT *,' Input choice: ' PRINT *,' 1 = To create a MATLAB file called mat.m ' PRINT *,' that will contain a copy of your n by n ' PRINT *,' input matrix and the size and tol you are ' PRINT *,' using. NOTE: If you already have a file, ' PRINT *,' mat.m, this program will not overwrite' PRINT *,' it. ' PRINT *,' If you run BHAP1, BHAP2, and BHAP3 ' PRINT *,' from the ensuing menu, their results ' PRINT *,' will also be stored in mat.m ' PRINT *,' ' PRINT *,' 0 = Do not create a file mat.m' PRINT *,' ' READ *,IMAT IF ((IMAT.NE.1) .AND. (IMAT.NE.0)) GO TO 80 * IMAT=1 means write to mat.m IF (IMAT.EQ.1) THEN OPEN (8,FILE='mat.m',STATUS='NEW',ERR=100) WRITE (8,FMT=*) ' n = ',SIZE,';' WRITE (8,FMT=*) ' tol = ',TOL,';' WRITE (8,FMT=*) ' a = eye(n) ;' DO I = 1,SIZE DO J = 1,SIZE WRITE (8,FMT=*) 'a(',I,',',J,')=',MAT0(I,J),';' END DO END DO END IF END IF *---------------------------------------------------------------- * * WARNING!! BHESS requires MAT0 to have SIZE + 1 Columns * though the matrix to be reduced has only SIZE Columns * *---------------------------------------------------------------- CALL BHESS(LDX,MAT0,SIZE,BHPIVT,NU,KRVECT,KROW,DUMV,TOL,INFO) IF (IMAT.EQ.1) THEN DO I = 1,SIZE DO J = 1,SIZE WMAT(I,J) = MAT0(I,J) END DO END DO CALL BHRPCK(LDX,SIZE,WMAT,NU) PRINT *,' The result of BHESS (with multipliers zeroed ' PRINT *,' by a call to BHRPCK) will be stored mat.m ' PRINT *,' as the matrix B ' PRINT *,' ' WRITE (8,FMT=*) ' B = eye(n) ;' WRITE (8,FMT=*) ' tol = ',TOL,' ;' DO I = 1,SIZE DO J = 1,SIZE WRITE (8,FMT=*) 'B(',I,',',J,')=',WMAT(I,J),';' END DO END DO END IF *-------------------------------------------------------------- * * The following continue statement (300) is returned to from * the bottom of the routine, allowing the user to * retry auxiliary routines * *--------------------------------------------------------------- 90 CONTINUE *-------------------------------------------------------------- * * The average bandwidth is computed as AVGNU * *-------------------------------------------------------------- AVGNU = 0 MAXNU = 0 IF (PRLEVEL.EQ.1) THEN PRINT *,' NU(I) is the number of nonzero entries to the right' + ,' of the diagonal in the ith row ' PRINT *,' ' PRINT *,' Input any integer to contiune' READ *,IDUM END IF DO I = 1,SIZE IF (PRLEVEL.EQ.1) THEN PRINT *,'NU(',I,')=',NU(I) IF (I.EQ. (I/20)*20) THEN PRINT *,' Input any integer to continue' READ *,IDUM END IF END IF AVGNU = AVGNU + NU(I) MAXNU = MAX(MAXNU,NU(I)) END DO AVGNU = (AVGNU+SIZE/2)/SIZE PRINT *,' Maximal upper bandwidth of returned matrix=',MAXNU PRINT *,' ' IF (PRLEVEL.EQ.1) THEN PRINT *,' ' PRINT *,' ' PRINT *,' Input 1 to see the nonzero entries of the' PRINT *,' returned Hessenberg matrix.' PRINT *,' Input any other integer to continue.' READ *,IGO IF (IGO.EQ.1) THEN PRINT *,'Row 1', (MAT0(1,J),J=1,1+NU(1)) DO I = 2,SIZE PRINT *,'Row',I, (MAT0(I,J),J=I-1,I+NU(I)) IF (I.EQ. (I/20)*20) THEN PRINT *,' Input any integer to continue' READ *,IDUM END IF END DO END IF END IF *---------------------------------------------------------------- * * The next parts of the code illustrate auxiliary returns * for working with the returned matrix. * *---------------------------------------------------------------- *---------------------------------------------------------------- * * Test BHAP1 * *---------------------------------------------------------------- PRINT *,' ' PRINT *,' ' PRINT *,' Input 1 to test BHAP1 by multiplying ' PRINT *,' by a vector of ones.' PRINT *,' Any other integer will continue to ' PRINT *,' the next option. ' PRINT *,' BHAP1 applies inv(Z) to a row vector where' PRINT *,' inv(Z)*A*Z= H, H the near tridiagonal matrix' PRINT *,' BHAP1 returns a row vector times inv(Z) ' PRINT *,' This operation converts a left eigenvector of' PRINT *,' H to a left eigenvector of A' PRINT *,' The output real and complex parts overwrite the' PRINT *,' input vector' PRINT *,' ' READ *,ITEST IF ((ITEST.EQ.1) .AND. (INIXSIX.EQ.0)) THEN DO I = 1,SIZE XRE(I) = 1.D0 XIM(I) = 0.D0 END DO CALL BHAP1(SIZE,MAT0,LDX,XRE,XIM,BHPIVT,KROW,IFLAG) XNORM1 = 0.D0 IF (PRLEVEL.EQ.0) THEN PRINT *,' ' PRINT *,' Input 1 to see output vector, ' PRINT *,' else input 0.' PRINT *,' ' READ *,ITEMP END IF IF ((ITEMP.EQ.1) .OR. (PRLEVEL.EQ.1)) THEN PRINT *,' XRE(I) I BHPIVT(I) KROW(I) ' DO I = 1,SIZE PRINT *,XRE(I),I,BHPIVT(I),KROW(I) IF (I.EQ. (I/20)*20) THEN PRINT *,' Enter any integer to continue' READ *,IDUM END IF XNORM1 = MAX(XNORM1,ABS(XRE(I))) END DO PRINT *,' The Max Norm is = ',XNORM1 PRINT *,' XRE(2)= ',XRE(2) PRINT *,' XRE(SIZE/2)= ',XRE(SIZE/2) PRINT *,' ' END IF IF (IMAT.EQ.1) THEN PRINT *,' The resulting vector will be stored as ' PRINT *,' the vector oneinvz ' PRINT *,' in the matlab file mat.m ' PRINT *,' ' DO J = 1,SIZE WRITE (8,FMT=*) 'oneinvz(1,',J,')=',XRE(J),';' END DO END IF ELSE IF ((ITEST.EQ.1) .AND. (INIXSIX.NE.0)) THEN PRINT *, + ' Multipliers have been zeroed so multiplication by' PRINT *,' Z is infeasible' END IF END IF *---------------------------------------------------------------- * * Test BHAP2 * *---------------------------------------------------------------- PRINT *,' ' PRINT *,' ' PRINT *,' Input 2 to test BHAP2 by multiplying ' PRINT *,' Z by a vector of ones.' PRINT *,' ' PRINT *,' BHAP2 multiplies an input row vector X' PRINT *,' by Z where ' PRINT *,' inv(Z)*A*Z= H, H the near tridiagonal matrix' PRINT *,' BHAP2 returns xt = xt*Z ' PRINT *,' ' READ *,ITEST IF ((ITEST.EQ.2) .AND. (INIXSIX.EQ.0)) THEN DO I = 1,SIZE XRE(I) = 1.D0 END DO CALL BHAP2(SIZE,MAT0,LDX,XRE,BHPIVT,KROW,IFLAG) IF (PRLEVEL.EQ.0) THEN PRINT *,' ' PRINT *,' Input 1 to see output vector, ' PRINT *,' else input 0 ' PRINT *,' ' READ *,ITEMP END IF XNORM2 = 0.D0 IF ((ITEMP.EQ.1) .OR. (PRLEVEL.EQ.1)) THEN PRINT *,' X(I) I BHPIVT(I) KROW(I) ' DO I = 1,SIZE PRINT *,XRE(I),I,BHPIVT(I),KROW(I) IF (I.EQ. (I/20)*20) THEN PRINT *,' Enter any integer to continue' READ *,IDUM END IF XNORM2 = MAX(XNORM2,ABS(XRE(I))) END DO PRINT *,' ' PRINT *,' The Max Norm is ',XNORM2 PRINT *,' ' PRINT *,' ' END IF IF (IMAT.EQ.1) THEN PRINT *,' The resulting vector will be stored in ' PRINT *, + ' the file mat.m as onez ' PRINT *,' ' DO J = 1,SIZE WRITE (8,FMT=*) 'onez(1,',J,')=',XRE(J),';' END DO END IF ELSE IF ((ITEST.EQ.2) .AND. (INIXSIX.NE.0)) THEN PRINT *,' ' PRINT *, + ' Multipliers have been zeroed so multiplication by' PRINT *,' Z is infeasible' END IF END IF *---------------------------------------------------------------- * * Test BHAP3 * *---------------------------------------------------------------- PRINT *,' ' PRINT *,' ' PRINT *,' Input 3 to test BHAP3 by multiplying ' PRINT *,' by a vector of ones.' PRINT *,' BHAP3 multiplies by Z where' PRINT *,' inv(Z)*A*Z= H, H the near tridiagonal matrix,' PRINT *,' The output real and complex parts overwrite the' PRINT *,' input vector. If the input vector is an' PRINT *,' eigenvector of H, the ouput vector is an ' PRINT *,' eigenvector A. ' READ *,ITEST IF ((ITEST.EQ.3) .AND. (INIXSIX.EQ.0)) THEN DO I = 1,SIZE XRE(I) = 1.D0 XIM(I) = 0.D0 END DO CALL BHAP3(SIZE,MAT0,LDX,XRE,XIM,BHPIVT,KROW,IFLAG) IF (PRLEVEL.EQ.0) THEN PRINT *,' ' PRINT *,' Input 1 to see output vector, ' PRINT *,' else input 0 ' PRINT *,' ' READ *,ITEMP END IF XNORM3 = 0.D0 IF ((ITEMP.EQ.1) .OR. (PRLEVEL.EQ.1)) THEN PRINT *,' XRE(I) I BHPIVT(I) KROW(I) ' DO I = 1,SIZE PRINT *,XRE(I),I,BHPIVT(I),KROW(I) IF (I.EQ. (I/20)*20) THEN PRINT *,' Enter any integer to continue' READ *,IDUM END IF XNORM3 = MAX(XNORM3,ABS(XRE(I))) END DO PRINT *,'The Max Norm is ',XNORM3 PRINT *,' XRE(2)= ',XRE(2) PRINT *,' XRE(SIZE/2)= ',XRE(SIZE/2) PRINT *,' ' PRINT *,' The result given here can be tested ' PRINT *,' by comparing to the Matlab-Octave version. ' PRINT *,' See comments in the script file abh3.m ' PRINT *,' ' PRINT *,' ' END IF IF (IMAT.EQ.1) THEN PRINT *,' The resulting vector will be stored in ' PRINT *, + ' the file mat.m as z1 ' PRINT *,' ' DO J = 1,SIZE WRITE (8,FMT=*) 'z1(',J,',1)=',XRE(J),';' END DO END IF ELSE IF ((ITEST.EQ.3) .AND. (INIXSIX.NE.0)) THEN PRINT *,' ' PRINT *, + ' Multipliers have been zeroed so multiplication by' PRINT *,' Z is infeasible' END IF END IF * --------------------------------------------------------------- * * The next statements call BHAPC -- a LINPACK style * estimator of the infinity norm of inv(Z) * * --------------------------------------------------------------- PRINT *,' ' PRINT *,' ' PRINT *,' Input 4 to estimate COND(Z) by multiplying ' PRINT *,' by a vector W of ones and minus ones.' PRINT *,' BHAPC multiplies by Z where' PRINT *,' inv(Z)*A*Z= H, H the near tridiagonal matrix' READ *,ITEST IF ((ITEST.EQ.4) .AND. (INIXSIX.EQ.0)) THEN CALL BHAPC(SIZE,MAT0,LDX,XRE,XIM,BHPIVT,KROW,CONDIT,IFLAG) IF (PRLEVEL.EQ.1) THEN PRINT *,' Y' DO I = 1,SIZE PRINT *,'Y(',I,')=',XRE(I) END DO END IF PRINT *,' ' PRINT *,' ESTIMATED CONDITION OF ACCUMULATED' PRINT *,' TRANSFORMATIONS =',CONDIT PRINT *,' ' ELSE IF ((ITEST.EQ.4) .AND. (INIXSIX.NE.0)) THEN PRINT *, + ' Multipliers have been zeroed so multiplication by' PRINT *,' Z is infeasible' END IF END IF * --------------------------------------------------------------- * * Test BHUNPK or BHRPCK * * --------------------------------------------------------------- PRINT *,' ' PRINT *,' Input any integer to continue' READ *,IDUM PRINT *,' ' PRINT *,' ' PRINT *,' Input 5 to test BHRPCK ' PRINT *,' BHRPCK converts the retuned matrix H to' PRINT *,' column storage with the first column' PRINT *,' the subdiagonal. This is convenient' PRINT *,' for BR iteration' PRINT *,' ' PRINT *,' ' PRINT *,' Input 6 TO test BHUNPCK ' PRINT *,' BHRPCK zeros the multipliers so further tests of' PRINT *,' multiplication of Z*x or xtr*inv(Z) and estimation' PRINT *,' of cond(Z) will not be possible' PRINT *,' ' PRINT *,' Input any other integer to go on.' PRINT *,' ' READ *,ITEST IF (ITEST.EQ.6) INIXSIX = 1 IF (ITEST.EQ.5) THEN CALL BHUNPK(LDX,MAT0,BH,SIZE,NU,NUMON,NUMAX) IF (PRLEVEL.EQ.1) THEN DO I = 1,SIZE PRINT *, (I,BH(I,J),J=1,NU(I)+2) END DO END IF END IF IF (ITEST.EQ.6) THEN * -------------------------------------------- * zero out multipliers in atobhes matrix MAT0 * -------------------------------------------- CALL BHRPCK(LDX,SIZE,MAT0,NU) IF (PRLEVEL.EQ.1) THEN DO I = 1,SIZE M = MIN(SIZE,NU(I)+2) DO J = 1,M IF (I+J-2.NE.0) THEN PRINT *,' MAT0(',I,',',I + J - 2,')=', + MAT0(I,I+J-2) END IF END DO END DO END IF END IF PRINT *,' Input 1 to reloop through auxiliary routines' PRINT *,' with the matrix returned from BHESS. ' PRINT *,' ' PRINT *,' Input any other integer to quit. ' IF (IMAT.EQ.1) THEN PRINT *,' ON SOME MACHINES, THE FILE mat.m WILL NOT EXIST ' PRINT *,' UNTIL YOU FINISH EXECUTING THE DRIVER ' END IF READ *,IGO IF (IGO.EQ.1) GO TO 90 IF (IMAT.EQ.1) CLOSE (8) GO TO 110 100 CONTINUE PRINT *,' ERROR: One of the following file open errors occured: ' PRINT *,' - error trying to create new file mat.m.' PRINT *,' this file may already exist. ' PRINT *,' - error trying to open file BHESS.IN. ' 110 CONTINUE PRINT *,' ' PRINT *,' ' PRINT *,' END OF DRIVER PROGRAM. ' STOP END * * ---------------------------------------------- DOUBLE PRECISION FUNCTION RANDOM(SEED) * ---------------------------------------------- * Returns a pseudo-random number between 0-1. * ---------------------------------------------- * COMMON SEED C .. Scalar Arguments .. INTEGER SEED C .. C .. Local Scalars .. DOUBLE PRECISION FMD INTEGER I,M,MD C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Data statements .. DATA M/25173/,I/13849/,MD/65536/,FMD/65536.D0/ C .. SEED = MOD(M*SEED+I,MD) RANDOM = SEED/FMD RETURN END * ----------------------------------------------- * *==================================================== SUBROUTINE BHBACK(LDX,N,A,B,W,V,TOL,NU,PIV,KRVECT,KROW) * This subroutine tests backward error of BHESS * by running BHESS to get a Hessenberg matrix * then using the stored multipliers to regenerate * the original matrix. * * Outputs are the difference between the original matrix * and the reconstruction (computed backward error) * * The upper bandwidth of the returned matrix and the estimated * condition number of the similariy transformation are * also returned. * * Arguments *------------------- * * N -- Integer Number of rows and columns of input matrices * * LDX -- Integer Leading Dimension of Variable dimension arrays * * A -- Double Precision Input matrix unchanged * * B and W -- Double Precision Matrices, B must of at * least N+1 columns. On return B is an approximation * of the input matrix A (to backward error) and W * holds the output from BHESS. * * TOL -- Double precision parameter to control multiplier * size in BHESS. Unchanged. * * V -- Double precision work vector of three columns. * * Local Variables *------------------------------- * * Copy A to B * C .. Scalar Arguments .. DOUBLE PRECISION TOL INTEGER LDX,N C .. C .. Array Arguments .. DOUBLE PRECISION A(LDX,*),B(LDX,*),V(LDX,3),W(LDX,*) INTEGER KROW(*),KRVECT(*),NU(*),PIV(*) C .. C .. Local Scalars .. DOUBLE PRECISION CONDIT,FNORM,RNORM,SUM,SUM2 INTEGER I,IFLAG,INFO,J,MXBAND C .. C .. External Subroutines .. EXTERNAL BHAP1,BHAP3,BHAPC,BHESS,BHRPCK,DCOPY C .. C .. Intrinsic Functions .. INTRINSIC DSQRT,MAX C .. DO J = 1,N DO I = 1,N B(I,J) = A(I,J) END DO END DO *-------------------------------------------------------- * * Call BHESS with B Warning B requires N + 1 Columns * *-------------------------------------------------------- CALL BHESS(LDX,B,N,PIV,NU,KRVECT,KROW,V,TOL,INFO) * * Copy B to W * DO J = 1,N DO I = 1,N W(I,J) = B(I,J) END DO END DO * * Zero out all but upper Hessenberg in B * CALL BHRPCK(LDX,N,B,NU) * * Calculate backward A as B = Z * B * inv(Z) * * First B <-- Z * B DO I = 1,N - 1,2 CALL BHAP3(N,W,LDX,B(1,I),B(1,I+1),PIV,KROW,IFLAG) END DO IF (N.NE.2* (N/2)) THEN CALL BHAP3(N,W,LDX,B(1,N),V(1,2),PIV,KROW,IFLAG) END IF * * Then B <-- B * inv(Z) So B approximates A * DO I = 1,N - 1,2 CALL DCOPY(N,B(I,1),LDX,V(1,1),1) CALL DCOPY(N,B(I+1,1),LDX,V(1,2),1) CALL BHAP1(N,W,LDX,V(1,1),V(1,2),PIV,KROW,IFLAG) CALL DCOPY(N,V(1,1),1,B(I,1),LDX) CALL DCOPY(N,V(1,2),1,B(I+1,1),LDX) END DO IF (N.NE.2* (N/2)) THEN CALL DCOPY(N,B(N,1),LDX,V(1,1),1) CALL BHAP1(N,W,LDX,V(1,1),V(1,2),PIV,KROW,IFLAG) CALL DCOPY(N,V(1,1),1,B(N,1),LDX) END IF * * Calculate || A - B || * SUM = 0.D0 SUM2 = 0.D0 DO I = 1,N DO J = 1,N SUM = SUM + (A(I,J)-B(I,J))**2 SUM2 = SUM2 + A(I,J)**2 END DO END DO * Frobenius norm normalized so that norm(eye(n)) = 1 FNORM = DSQRT(SUM/N) * Relative Frobenius norm RNORM = DSQRT(SUM/SUM2) * * Estimate conditioning of Z * CALL BHAPC(N,W,LDX,V(1,1),V(1,2),PIV,KROW,CONDIT,IFLAG) * * Calculate Maximal Bandwidth * MXBAND = 1 DO I = 1,N MXBAND = MAX(NU(I),MXBAND) END DO PRINT *,' ' PRINT *,' ' PRINT *,' n = ',N PRINT *,' tol = ',TOL PRINT *,' Frobenius Backward Error Norm = ',FNORM PRINT *,' Relative Backward Error Norm = ',RNORM PRINT *,' Estimated Condition Number = ',CONDIT PRINT *,' Bandwidth = ',MXBAND RETURN END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'dblas1.f' then echo shar: will not over-write existing file "'dblas1.f'" else cat << "SHAR_EOF" > 'dblas1.f' cat > dblas1.f <<'CUT HERE..........' c subroutine dblas() c d = dasum(n,dx,incx) c call daxpy(n,da,dx,incx,dy,incy) c call dcopy(n,dx,incx,dy,incy) c d = ddot(n,dx,incx,dy,incy) c d = dmach(job) c d = dnrm2 ( n, dx, incx) c call drot (n,dx,incx,dy,incy,c,s) c call drotg(da,db,c,s) c call dscal(n,da,dx,incx) c call dswap (n,dx,incx,dy,incy) c i = idamax(n,dx,incx) c stop c end CUT HERE.......... cat > idamax.f <<'CUT HERE..........' integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified to correct problem with negative increment, 8/21/90. c double precision dx(1),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 dmax = dabs(dx(ix)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end CUT HERE.......... cat > dswap.f <<'CUT HERE..........' subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end CUT HERE.......... cat > dscal.f <<'CUT HERE..........' subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified to correct problem with negative increment, 8/21/90. c double precision da,dx(1) integer i,incx,ix,m,mp1,n c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n dx(ix) = da*dx(ix) ix = ix + incx 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end CUT HERE.......... cat > drotg.f <<'CUT HERE..........' subroutine drotg(da,db,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c modified 9/27/86. c double precision da,db,c,s,roe,scale,r,z c roe = db if( dabs(da) .gt. dabs(db) ) roe = da scale = dabs(da) + dabs(db) if( scale .ne. 0.0d0 ) go to 10 c = 1.0d0 s = 0.0d0 r = 0.0d0 go to 20 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) r = dsign(1.0d0,roe)*r c = da/r s = db/r 20 z = s if( dabs(c) .gt. 0.0d0 .and. dabs(c) .le. s ) z = 1.0d0/c da = r db = z return end CUT HERE.......... cat > drot.f <<'CUT HERE..........' subroutine drot (n,dx,incx,dy,incy,c,s) c c applies a plane rotation. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp,c,s integer i,incx,incy,ix,iy,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = c*dx(ix) + s*dy(iy) dy(iy) = c*dy(iy) - s*dx(ix) dx(ix) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n dtemp = c*dx(i) + s*dy(i) dy(i) = c*dy(i) - s*dx(i) dx(i) = dtemp 30 continue return end CUT HERE.......... cat > dnrm2.f <<'CUT HERE..........' double precision function dnrm2 ( n, dx, incx) integer i, incx, ix, j, n, next double precision dx(1), cutlo, cuthi, hitest, sum, xmax,zero,one data zero, one /0.0d0, 1.0d0/ c c euclidean norm of the n-vector stored in dx() with storage c increment incx . c if n .le. 0 return with result = 0. c if n .ge. 1 then incx must be .ge. 1 c c c.l.lawson, 1978 jan 08 c modified to correct problem with negative increment, 8/21/90. c modified to correct failure to update ix, 1/25/92. c c four phase method using two built-in constants that are c hopefully applicable to all machines. c cutlo = maximum of dsqrt(u/eps) over all known machines. c cuthi = minimum of dsqrt(v) over all known machines. c where c eps = smallest no. such that eps + 1. .gt. 1. c u = smallest positive no. (underflow limit) c v = largest no. (overflow limit) c c brief outline of algorithm.. c c phase 1 scans zero components. c move to phase 2 when a component is nonzero and .le. cutlo c move to phase 3 when a component is .gt. cutlo c move to phase 4 when a component is .ge. cuthi/m c where m = n for x() real and m = 2*n for complex. c c values for cutlo and cuthi.. c from the environmental parameters listed in the imsl converter c document the limiting values are as follows.. c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are c univac and dec at 2**(-103) c thus cutlo = 2**(-51) = 4.44089e-16 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. c thus cuthi = 2**(63.5) = 1.30438e19 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. c thus cutlo = 2**(-33.5) = 8.23181d-11 c cuthi, d.p. same as s.p. cuthi = 1.30438d19 c data cutlo, cuthi / 8.232d-11, 1.304d19 / c data cutlo, cuthi / 4.441e-16, 1.304e19 / data cutlo, cuthi / 8.232d-11, 1.304d19 / c if(n .gt. 0) go to 10 dnrm2 = zero go to 300 c 10 assign 30 to next sum = zero i = 1 if( incx .lt. 0 )i = (-n+1)*incx + 1 ix = 1 c begin main loop 20 go to next,(30, 50, 70, 110) 30 if( dabs(dx(i)) .gt. cutlo) go to 85 assign 50 to next xmax = zero c c phase 1. sum is zero c 50 if( dx(i) .eq. zero) go to 200 if( dabs(dx(i)) .gt. cutlo) go to 85 c c prepare for phase 2. assign 70 to next go to 105 c c prepare for phase 4. c 100 continue ix = j assign 110 to next sum = (sum / dx(i)) / dx(i) 105 xmax = dabs(dx(i)) go to 115 c c phase 2. sum is small. c scale to avoid destructive underflow. c 70 if( dabs(dx(i)) .gt. cutlo ) go to 75 c c common code for phases 2 and 4. c in phase 4 sum is large. scale to avoid overflow. c 110 if( dabs(dx(i)) .le. xmax ) go to 115 sum = one + sum * (xmax / dx(i))**2 xmax = dabs(dx(i)) go to 200 c 115 sum = sum + (dx(i)/xmax)**2 go to 200 c c c prepare for phase 3. c 75 sum = (sum * xmax) * xmax c c c for real or d.p. set hitest = cuthi/n c for complex set hitest = cuthi/(2*n) c 85 hitest = cuthi/float( n ) c c phase 3. sum is mid-range. no scaling. c do 95 j = ix,n if(dabs(dx(i)) .ge. hitest) go to 100 sum = sum + dx(i)**2 i = i + incx 95 continue dnrm2 = dsqrt( sum ) go to 300 c 200 continue ix = ix + 1 i = i + incx if( ix .le. n ) go to 20 c c end of main loop. c c compute square root and adjust for scaling. c dnrm2 = xmax * dsqrt(sum) 300 continue return end CUT HERE.......... cat > dmach.f <<'CUT HERE..........' double precision function dmach(job) integer job c c smach computes machine parameters of floating point c arithmetic for use in testing only. not required by c linpack proper. c c if trouble with automatic computation of these quantities, c they can be set by direct assignment statements. c assume the computer has c c b = base of arithmetic c t = number of base b digits c l = smallest possible exponent c u = largest possible exponent c c then c c eps = b**(1-t) c tiny = 100.0*b**(-l+t) c huge = 0.01*b**(u-t) c c dmach same as smach except t, l, u apply to c double precision. c c cmach same as smach except if complex division c is done by c c 1/(x+i*y) = (x-i*y)/(x**2+y**2) c c then c c tiny = sqrt(tiny) c huge = sqrt(huge) c c c job is 1, 2 or 3 for epsilon, tiny and huge, respectively. c double precision eps,tiny,huge,s c eps = 1.0d0 10 eps = eps/2.0d0 s = 1.0d0 + eps if (s .gt. 1.0d0) go to 10 eps = 2.0d0*eps c s = 1.0d0 20 tiny = s s = s/16.0d0 if (s*1.0 .ne. 0.0d0) go to 20 tiny = (tiny/eps)*100.0 huge = 1.0d0/tiny c if (job .eq. 1) dmach = eps if (job .eq. 2) dmach = tiny if (job .eq. 3) dmach = huge return end CUT HERE.......... cat > ddot.f <<'CUT HERE..........' double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end CUT HERE.......... cat > dcopy.f <<'CUT HERE..........' subroutine dcopy(n,dx,incx,dy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) 50 continue return end CUT HERE.......... cat > daxpy.f <<'CUT HERE..........' subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end CUT HERE.......... cat > dasum.f <<'CUT HERE..........' double precision function dasum(n,dx,incx) c c takes the sum of the absolute values. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified to correct problem with negative increment, 8/21/90. c double precision dx(1),dtemp integer i,incx,ix,m,mp1,n c dasum = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n dtemp = dtemp + dabs(dx(ix)) ix = ix + incx 10 continue dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' *--------1--------2--------3---------4---------5---------6---------7-* * BHESS and Auxiliary Routines * * * * Fed a square matrix N by N matrix dimensioned A(LDA,N+1) * * and with the N by N matrix having upper left corner A(1,1) * * BHESS overwrites the upper N by N left corner with a * * similar banded Hessenberg matrix and the multipliers * * required to accomplish the elimination. * * * * The code is essentially Fortran 77 but contains minor * * extensions that are accepted by all Fortran 90 compilers * * (and indeed by all Fortran 77 compilers we have found). * *--------------------------------------------------------------------* * * * CONTENTS * * * * BHESS -- Performs H = inv(Z)* A * Z * * A a general square matrix. * * H an upper Hessenberg matrix. * * Z is a product of elementary matrices. * * H and the multipliers for Z overwrite A * * * * BHPIVT - auxiliary routine used by BHESS to determine which * * rows can be eliminated, and to choose permutations * * * * BHAP1 - BHESS gives H = inv(Z)*A*Z. * * BHAP1 performs x^T <-- x^T * inv(Z) * * where X is stored as the two vectors XRE and XIM * * * * BHAP2 - x^T <-- x^T * Z where x is a vector * * * * BHAP3 - (XRE,XIM) <-- Z * (XRE,XIM) * * * * BHAPC - Condition number estimator for Z * * * * BHRPCK - Zeros out the multipliers so that only the upper * * Hessenberg matrix remains * * * * BHUNPK - Puts the upper Hessenberg matrix is diagonal * * storage * * * * For some other useful routines, see the anonymous * * ftp site * * * * ftp cs.fit.edu * * cd pub/howell * * * * In particular, the br eigenvalue algorithm * * (SIMAX, summer 1999) is often successful in extracting * * eigenvalues of the matrix returned by BHESS in time * * proportional to the square of the matrix size. * * * *--------------------------------------------------------------------* * * Code authored by Gary Howell and Nadia Diaa * Authors offer it for public use on condition that * our authorship recognized. We hope the code will be useful * but do not warrant its reliability. * * Note: much of the code is a revision of ATOTRI * by Al Geist (SIMAX 1991). The BHESS version differs * in that multipliers are bounded and the returned * Hessenberg matrix is not necessarily tridiagonal. * * Thanks to Charles Fulton, David Watkins, * Michael Saunders, and Al Geist for suggestions. * * First submitted to TOMS in 1994. * *--------------------------------------------------------------------* SUBROUTINE BHESS(LDA,A,N,PIVOTS,NU,KRVECT,KROW,V,TOL,INFO) * ---------------------- * * .. Scalar Arguments .. * ---------------------- INTEGER LDA, N, INFO DOUBLE PRECISION TOL * -------------------------------------------- * * .. Array Arguments .. * -------------------------------------------- INTEGER PIVOTS(*), NU(*), KRVECT(*), KROW(*) DOUBLE PRECISION A(LDA,*), V(*) * -------------------------------------------- * * .. External Functions .. * -------------------------------------------- DOUBLE PRECISION DDOT * --------------------------------------------------------- * * DDOT is from the DBLAS-1 library. * ---------------------------------------------------------- * * Purpose: * This subroutine reduces an n-by-n real general matrix A to * small-band Hessenberg form using elementary Gaussian * similarity transformations. * * At each step k the permutation that minimizes the maximum entry * in the transformation matrix that reduces column k and then row k * is applied. * * July 9, 2004 * Added quick returns for N less than 2 or N greater than LDA * * September 15, 1997 * Cache performance is improved by performing both the row * and column similarity transformations on a given column * in one access. Execution time relative to LAPACK DGEHRD is * (on a RS-6000) 85% for size 200 and 80% for size 1200. * * June 25, 1995 * This version of BHESS does not use BLAS level 2. Outer loops * corresponding to the BLAS 2 operation DGER and DGMEMV are in * the column variable. This version is twice as fast on various * IBM workstations. * * * Arguments: * LDA -integer * leading dimension of A * * A -double precision array of dimension (LDA,N+1) * On entry, the first N columns of A contain the matrix * being reduced. * On exit the first N columns of A are overwritten by * its banded Hessenberg form and the multipliers used * to produce the small band form. * The extra (N+1st) column of A is used as a work * vector. * * N -integer * N specifies the order of the matrix A. * N must be nonnegative. * not modified. * * PIVOTS -integer vector of dimension (LDA) * On exit pivots contains the pivot sequence used during * the reduction (permutation vector). * * NU -integer vector of dimension N. * NU is initialized as all ones * on exit NU(KR) is the number of consecutive * presumed nonzero entries to right of * diagonal in row KR. * * KRVECT -dummy integer vector of dimension N. * on entry: n/a as it is a dummy var passed only so we * can use variable dimension N. * on exit: n/a * KRVECT is initialized as all ones. * Whenever row KR is zeroed, KRVECT(KR) becomes 0. * KRVECT is used internally to keep track of which rows * are eligible to be zeroed. * * KROW -integer vector of dimension N. * On exit KROW(K) contains the row eliminated after * column K was eliminated. If the row elimination was * never performed, then KROW(K) is zero. * * TOL not modified * - when the L2 norms of rows and columns of multipliers are equal, * then TOL corresponds to the maximal standard deviation * of multipliers from zero. More generally TOL is the maximal * allowed root mean square of standard deviations * from zero of row and column multipliers. * * V -double precision dummy vector * * * INFO -integer * On exit, INFO is set to * 0 normal return. * 1 if an error is detected. * -1 if N less than 3 * -2 if N greater than LDA * * .. Local Scalars .. * -------------------------------------------------------------------- INTEGER FIRST1, I, J, K, KR, L, PIVP, ERR, CNT DOUBLE PRECISION TEMP, MINMLT, ONE, INPROD PARAMETER (ONE = 1.0) * -------------------------------------------------------------------- * * FIRST1 is the first 1 in KRVECT, i.e., the index of the first * nonzero row. * * .. External Subroutines .. * -------------------------------------------------------------------- EXTERNAL BHPIVT, DAXPY, DSCAL * DBLAS DDOT,DAXPY,DSCAL * *=============Executable Statements =============================== * -------------- * * ------------------------------------ * Quick returns for wrong sized matrices * ------------------------------------ IF (N. LE. 2) THEN INFO = -1 PRINT*,' Matrix size N must be greater than 2' RETURN ELSEIF (N. GE. LDA) THEN INFO = -2 PRINT*,' Matrix size N must be smaller than', LDA RETURN ENDIF * -------------- * Initialization * -------------- INFO = 0 CNT = 0 L = 1 PIVOTS(1) = 1 PIVOTS(N) = N K = 1 KR = 1 FIRST1 = 1 DO I = 1,N-1 NU(I) = 1 KRVECT(I) = 1 KROW(I) = 0 END DO NU(N) = 0 KRVECT(N) = 1 KROW(N) = 0 * ------------------------------------ * For each column of the matrix * ------------------------------------ DO 1000 K=1, N-2 * ------------------------------------ * We will come to 10 when we use same column with a different row * i.e., do not increment column counter K. * ------------------------------------ 10 CONTINUE * --------------------------------- * Find a suitable pivot. Integer variable ERR returned by * BHPIVT contains information as to whether a row elimination * should be performed. * ------------------------------------------------------ CALL BHPIVT ( LDA,A,N,K,KR,TOL,PIVP,INPROD,MINMLT,ERR ) PIVOTS(K+1) = PIVP * ------------------------------------------------------ * Check for deflation , i.e. both row & column are zero. * Column elimination is not necessary. * (instead of going directly to 1000, update a few * variables first and then fall into 1000 continue statement) * ------------------------------------------------------ IF( ERR .EQ. 1 ) THEN L = K+1 GO TO 888 ENDIF * ----------------------------------------------------------- * Check for (inner product = zero) OR (inner product < tolerance) * ----------------------------------------------------------- IF ((ERR .EQ. -1) .OR. (ERR .EQ. -2)) THEN * -------------------------------------------------------- * Row KR is not eligible so try the next available row * (if any) * First, update NU(KR) to indicate failure to zero out row KR * -------------------------------------------------------- IF (KR .EQ. K) THEN * ------------------------------------------------------- * We have tried all eligible rows with col K, * all resulted in ERR =-1 or -2. * This means we cannot zero out any rows, so just * go ahead and zero out col K anyway. * update subsequent entries of KR to reflect wider * bandwidth resulting from failure to zero out row KR. * Zero out column K if needed * -------------------------------------------------------- IF( ERR .NE. 3 ) THEN * ------------------------------------------------------ * *********** Zero out column K ********************** * but no row is zeroed in conjunction with column K. * * Interchange row and column * ------------------------------------------------------ IF( PIVP .NE. K+1 ) THEN CALL DSWAP (N,A(PIVP,1),LDA,A(K+1,1),LDA) CALL DSWAP (N,A(1,PIVP),1,A(1,K+1),1) ENDIF TEMP = A(K+1,K) CALL DSCAL((N-(K+2)+1),ONE/TEMP,A(K+2,K),1) * ---------------------------------------------------- * The next block of code combines the rank 1 update * and the matrix vector multiply needed to * accomplish an elimination of column K * and preserve similarity. * Each column K+1 to N is accessed only once for the * conbined rank 1 update and matrix vector multiply. * -------------------- * Rank 1 update * -------------------- * Update is addition of -u*vtr where * u is K+2 to N entries of column K * vtr is K+1 to N entries of row K+1 * The upper left corner of updated block is A(K+2,K). * --------------------- * Matrix vector product * --------------------- * For KR(I) not zero * column K+1 of A starting at row FIRST1 * is incremented by A multiplied by the column * K of A (entries K+2 to N) * where the entries of A used consist of the block * starting at column K+2 and the row K. * * ----------------------------------------------------- CALL DAXPY(N-K-1,-A(K+1,K+1),A(K+2,K),1,A(K+2,K+1),1) DO J = K+2,N CALL DAXPY(N-K-1,-A(K+1,J),A(K+2,K),1,A(K+2,J),1) CALL DAXPY(N-K-1,A(J,K),A(K+2,J),1,A(K+2,K+1),1) DO I = FIRST1,K+1 IF ( KRVECT(I) .NE. 0 ) THEN A(I,K+1) = A(I,K+1) + A(J,K)*A(I,J) ENDIF END DO END DO ENDIF * -------------------------------------------------------- * End of (ERR .NE. 3) * Eliminated column K so drop to bottom of main loop * and start work on column K+1 * -------------------------------------------------------- GO TO 999 ELSE * ---------------------------------------------- * KR is smaller than K, * could not zero out row KR * find next nonzero row to try with col K * ---------------------------------------------- DO I = KR+1,K IF (KRVECT(I) .EQ. 1) GO TO 451 ! Row I has ! not been eliminated. END DO 451 CONTINUE KR = I * ---------------------------------------------------- * Only change KR, not FIRST1, * now go find pivot using same col K with new row KR, * so bypass loop increment of K by going to 10 * ---------------------------------------------------- GO TO 10 ! where we will see if Row I can be ! given a paired elimination with Column K ENDIF ENDIF * ------------------------------------------------------------- * Check if row-column inner product large enough. * If so (ERR .NE. 2), then zero out row KR along with column K * ------------------------------------------------------------ IF( ERR .NE. 2 ) THEN * ------------------------------------------------------ * zero out row KR and set KROW(K) as KR * ***************************************************** * Interchange row and column (using pivot chosen in BHPIVT) * ----------------------------------------------------- IF( PIVP .NE. K+1 ) THEN CALL DSWAP (N,A(PIVP,1),LDA,A(K+1,1),LDA) CALL DSWAP (N,A(1,PIVP),1,A(1,K+1),1) ENDIF KROW(K)=KR TEMP = A(KR,K+1) KRVECT(KR) = 0 * ----------------------------------------------- * The vector V(I) is the row multipliers. These overwrite * the KR row of the matrix and are stored as V(I) to * avoid reaccessing columns of the matrix A. * ----------------------------------------------- DO I = K+2,N V(I) = A(KR,I)/TEMP A(KR,I) = V(I) END DO DO I = FIRST1,K IF( KRVECT(I) .NE. 0 ) A(I,N+1)=A(I,K+1) END DO * ------------------------------------------------ * A backup copy of column K+1 of the matrix is used. * The copy stored in column N+1 . * The copy is updated by multiples of columns as part * of the matrix vector multiply associated with the row * elimination. The other is used for the the rank one * update used to eliminate row KR. * ----------------------------------------------------- CALL DCOPY(N-K,A(K+1,K+1),1,A(K+1,N+1),1) * ------------------------------------------------ * Row K+1 of columns K and K+1 * is updated by the matrix vector * multiply associated with * eliminating column K. * ------------------------------------------------ A(K+1,N+1) = A(K+1,N+1) + DDOT(N-K-1,A(K+2,N+1),1,V(K+2),1) A(K+1,K) = A(K+1,K) + DDOT(N-K-1,A(K+2,K),1,V(K+2),1) * ------------------------------------------------ * The multipliers to eliminate column K are stored * in the entry they are used to eliminate * ------------------------------------------------ TEMP = A(K+1,K) IF ( ABS(TEMP) .GT. 1.D-12 ) THEN CALL DSCAL((N-(K+2)+1),ONE/TEMP,A(K+2,K),1) ENDIF * ------------------------------------------------- * The following loop on columns accesses each column * only once in the course of combined a combined * similarity transformation to eliminate both * column K and row KR * ------------------------------------------------- DO J = N,K+1,-1 * ----------------------------------------------- * The rank one update for eliminating row KR. * ----------------------------------------------- IF ( J .NE. K+1 ) THEN DO I = FIRST1,K IF ( KRVECT(I) .NE. 0 ) THEN A(I,J) = A(I,J) - A(I,K+1)*V(J) ENDIF END DO CALL DAXPY(N-K,-V(J),A(K+1,K+1),1,A(K+1,J),1) ENDIF * ---------------------------------------------- * Increment row K+1 by the inner product of multipliers * and succeeding rows. * ---------------------------------------------- IF ( J .NE. K+1 ) THEN TEMP = DDOT(N-K-1,A(K+2,J),1,V(K+2),1) A(K+1,J) = A(K+1,J) + TEMP ENDIF * ---------------------------------------------- * Increment column K+1 by a multiple of column J * and perform rank one update corresponding to * elimination of column K * ---------------------------------------------- IF ( J .NE. K+1 ) THEN CALL DAXPY(N-K-1,A(J,K),A(K+2,J),1,A(K+2,N+1),1) CALL DAXPY(N-K-1,-A(K+1,J),A(K+2,K),1,A(K+2,J),1) DO I = FIRST1,K+1 IF ( KRVECT(I) .NE. 0 ) THEN A(I,N+1) = A(I,N+1) + A(J,K)*A(I,J) ENDIF END DO ENDIF END DO * ---------------------------------------------------- * Overwrite the old column K+1 with the new version which * has been incrementing in column N+1 and do the * last bit of the rank one update for column elimination. * ----------------------------------------------------- DO I = FIRST1,K+1 IF ( KRVECT(I) .NE. 0 ) THEN A(I,K+1) = A(I,N+1) ENDIF END DO CALL DCOPY(N-K-1,A(K+2,N+1),1,A(K+2,K+1),1) CALL DAXPY(N-K-1,-A(K+1,K+1),A(K+2,K),1,A(K+2,K+1),1) ENDIF 888 CONTINUE * --------------------------------------------------- * We have zeroed row KR so reflect this information * by updating KRVECT(KR) and then FIRST1 * --------------------------------------------------- NU(KR)=K+1-KR KRVECT(KR) = 0 DO I = FIRST1, K + 1 IF ( KRVECT(I) .EQ. 1 ) GO TO 895 END DO 895 CONTINUE FIRST1 = I 999 CONTINUE * ----------------------------------------------------- * Going to get a new col K, so we will want to try * it out with first nonzero row * ----------------------------------------------------- KR = FIRST1 1000 CONTINUE DO I = 1,N-1 IF( KRVECT(I) .EQ. 1 ) NU(I) = N-I END DO RETURN * END ! End of BHESS * ------------------------------------- * * * ============================================================ * * ------------------------------------------------------------ SUBROUTINE BHPIVT (LDA,A,N,K,KR,TOL,PIVP,INPROD,MINMLT,ERR ) * ------------------------------------------------------------ * * .. Scalar Arguments .. * ------------------------------------------------------------ INTEGER LDA, N, K, KR, PIVP, ERR DOUBLE PRECISION MINMLT, TOL, INPROD * ------------------------------------ * * .. Array Arguments .. * ------------------------------------ DOUBLE PRECISION A(LDA,*) * ------------------------------------ * * Purpose: * * This subroutine tests whether row KR is eligible to eliminated * at in conjunction with column K. If so, it finds the pivot * that minimizes the maximum entry * in the matrix N, where N=N_rN_c and N_c reduces column k and * inv(N_r) reduces row k. In this event return err=0. * * If the row is not eligible for elimination with the kth column, then the * pivot is chosen to correspond to the largest entry below the subdiagonal * in column k. Then return err=-2. * * Routine also checks if the problem can be deflated at step k * returning err=1 if so. If inner product of the row and column is * zero then a pivot is selected based on maximum row/col entries * and err is set to 2. * * Arguments: * A -double precision array of dimension (n,n) * On entry A contains the partially reduced matrix. * Not modified * * N -integer * N specifies the order of the matrix A. * N must be at least zero. * Not modified. * * K -integer * K specifies the column under consideration. * Not modified * * KR -integer * KR specifies the row under consideration. * Not modified * * TOL -double precision scalar is between 0 and +1 * on entry contains allowable cosine of angle * between col k and row kr * not modified * * PIVP -integer * On return PIVP contains the pivot index * necessary to minimize the maximum entry in N. * * MINMLT -double precision * On return MINMLT contains the minimized maximum entry. * * ERR -integer * On exit ERR is set to *==> -2 inner product is too small. * -1 inner product = 0, recovery step required.(obsolete) * 0 normal return. * 1 problem deflated, row and column are zero * 2 problem deflated, row is zero * 3 problem deflated, column is zero * * .. Local Scalars .. * ----------------------------------------------- INTEGER I, PIVC, PIVR DOUBLE PRECISION MAXCOL, MAXROW, NMXCOL, NMXROW DOUBLE PRECISION MAXNC, MAXNR, MXDIAG, V1 DOUBLE PRECISION TEMC, TEMR, TEMP,SUMC,SUMR DOUBLE PRECISION BIG,TOLK * ----------------------------------------------- * * .. Parameters .. * ----------------------------------------------- PARAMETER (BIG = 1.0D8) * ----------------------- * * * .. External Functions .. * none * * .. External Subroutines .. * none * * .. Intrinsic Functions .. INTRINSIC ABS c INTRINSIC MAX * *=================== Executable Statements ================================= * * Initialize returned values * ------------------------- PIVP = K+1 MINMLT = BIG ERR = 0 TOLK=TOL*(N-K) TOLK=DMIN1(1.D0/TOLK,1.D0) * ----------------------------------------------------- * Find maximum and next-to-max entries in row and col k * and compute inner product. * ----------------------------------------------------- PIVC = K+1 PIVR = K+1 NMXCOL = 0.D0 NMXROW = 0.D0 MAXCOL = A(K+1,K) MAXROW = A(KR,K+1) SUMC = MAXCOL*MAXCOL SUMR = MAXROW*MAXROW INPROD = MAXCOL * MAXROW MAXCOL = DABS(MAXCOL) MAXROW = DABS(MAXROW) DO I=K+2, N TEMC = A(I,K) TEMR = A(KR,I) SUMC = SUMC + TEMC*TEMC SUMR = SUMR + TEMR*TEMR INPROD = INPROD + TEMC * TEMR TEMC = DABS(TEMC) TEMR = DABS(TEMR) IF( TEMC .GT. NMXCOL ) THEN IF( TEMC .GT. MAXCOL ) THEN NMXCOL = MAXCOL MAXCOL = TEMC PIVC = I ELSE NMXCOL = TEMC ENDIF ENDIF IF( TEMR .GT. NMXROW ) THEN IF( TEMR .GT. MAXROW ) THEN NMXROW = MAXROW MAXROW = TEMR PIVR = I ELSE NMXROW = TEMR ENDIF ENDIF END DO * -------------------------------------------------- * Check for deflation * -------------------------------------------------- IF( MAXCOL .EQ. 0.D0 .AND. MAXROW .EQ. 0.D0 ) THEN ERR = 1 RETURN ENDIF IF( MAXROW .EQ. 0.D0 ) THEN ERR = 2 PIVP = PIVC RETURN ENDIF IF( MAXCOL .EQ. 0.D0 ) THEN ERR = 3 PIVP = PIVR RETURN ENDIF * --------------------------------------------------- * Check for inner product too small. * --------------------------------------------------- SUMC = DSQRT(SUMC) SUMR = DSQRT(SUMR) IF( DABS(INPROD/SUMC/SUMR) .LT. TOLK ) THEN * --------------------------------------------------- * Inner product of row and column is too small for elimination * of row so return pivot to correspond to largest column element * --------------------------------------------------- PIVP = PIVC ERR = -2 RETURN ENDIF * --------------------------------------------------- * Calculate maximum multipliers over all permutations i and * determine the minimum triple (maxnc,maxnr,mxdiag)_i. * In this version, the multipliers are selected for the * row elimination first, then column elimination. * --------------------------------------------------- DO I=K+1, N V1 = A(KR,I) IF( V1 .EQ. 0.D0 ) THEN MAXNR = BIG c ELSE IF( I .EQ. PIVC ) THEN ELSE IF( I .EQ. PIVr ) THEN MAXNR = DABS(NMXROW/V1) ELSE MAXNR = DABS(MAXROW/V1) ENDIF c IF( I .EQ. PIVR ) THEN IF( I .EQ. PIVc ) THEN MAXNC = DABS(NMXCOl*V1/INPROD) ELSE MAXNC = DABS(MAXCOl*V1/INPROD) ENDIF MXDIAG = ABS(V1*A(I,K)/INPROD) TEMP = MAXNR IF( TEMP .LT. MAXNC ) TEMP = MAXNC IF( TEMP .LT. MXDIAG ) TEMP = MXDIAG IF( TEMP .LT. MINMLT ) THEN MINMLT = TEMP PIVP = I ENDIF END DO * RETURN END ! End of BHPIVT * ------------------------------------- * *==================================================== *========================================================================= SUBROUTINE BHAP1(N,Z,LDZ,XRE,XIM,IPVT,KROW,IFLAG) * INTEGER N, LDZ, IPVT(*), KROW(*), IFLAG DOUBLE PRECISION Z(LDZ,*), XRE(*), XIM(*) *--------------------------------------------------------- * * * Purpose: * This routine applies the accumulated elementary * transformations inv(Z) to a row vector. The matrix * Z is the product of the elementary transformations * and permutation matrices that near-tridiagonalize the * original matrix A. The computation is of the form * inv(Z)*A*Z = H, where * inv(Z) = Lc(n-2) * inv(Lr(n-2)) * P(n-2) * * Lc(n-3)*inv(Lr(n-3) * P(n-3) * ... Lc(1)*inv(Lr(1))*P(1) * Each of the elementary transformations is stored as a row or * column of the supplied array in this routine. * Multiplying a row eigenvector using BHAP1 results in * in an eigenvector of A. * * July 9, 2004 -- removed unused dummy variable W. * * Arguments: * * N specifies the order of the matrix Z. * N must be nonnegative. * not modified. * * Z -double precision array, dimension (LDZ,N) * Z contains information about the reduction to * tridiagonal form. * * LDZ -integer * The leading dimension of the array Z. * LDZ >= max(1,N). * * XRE -double precision array, dimension (N) * The real part of the vector used in applying the * transformations. * * XIM -double precision array, dimension (N) * The imaginary part of the vector used in applying * the transformations. * * IPIV -integer array, dimension (N) * Contains the pivot sequence used during the reduction * to tridiagonal form. * * KROW -integer array, dimension (N) * Which row (if any) was eliminated after * column elimination k. * * IFLAG -integer * You can put an error flag if you want. * Not currently used. * * .. Parameters .. * * ------------------------------------------ DOUBLE PRECISION TWO,ZERO PARAMETER (TWO = 2.0, ZERO = 0.0) * ------------------------------------------- * .. * .. Local Scalars .. * ------------------------------------------- INTEGER K, L, KP1 DOUBLE PRECISION MRE, MIM, DDOT EXTERNAL DDOT * ------------------------------------------- * .. * .. Executable Statements .. * * Apply the permutations. * ----------------------------- DO K = N-2,1,-1 KP1 = K+1 * ------------------------------------------------------- * Multiply by Lc(k). * ------------------------------------------------------- XRE(KP1) = XRE(KP1) - DDOT(N-K-1,XRE(K+2),1,Z(K+2,K),1) XIM(KP1) = XIM(KP1) - DDOT(N-K-1,XIM(K+2),1,Z(K+2,K),1) * -------------------------- * Multiply by inv(Lr(k)). * -------------------------- IF ( KROW(K).NE.0 ) THEN * ------------------------------------------------------- * The stride of LDZ in the DAXPY is poor for cache * efficiency, but is a consequence of storing the * multipliers in the space they were used to eliminate. * If multiplications by inv(Z)are to be performed many * times, consider packed triangular storage for * upper triangular multipliers. The packed version would * store the row elimination multipliers sequentially. * -------------------------------------------------------- CALL DAXPY(N-K-1,XRE(KP1),Z(KROW(K),K+2),LDZ,XRE(K+2),1) CALL DAXPY(N-K-1,XIM(KP1),Z(KROW(K),K+2),LDZ,XIM(K+2),1) ENDIF END DO DO 70 K = N-2,1,-1 KP1 = K+1 L = IPVT(KP1) MRE = XRE(L) MIM = XIM(L) IF ( L .EQ. KP1 ) GO TO 70 XRE(L) = XRE(KP1) XIM(L) = XIM(KP1) XRE(KP1) = MRE XIM(KP1) = MIM 70 CONTINUE * RETURN END ! End of BHAP1 *======================================================================= SUBROUTINE BHAP2(N,Z,LDZ,X,IPVT,KROW,IFLAG) * ---------------------------------------------- * * Global variables * --------------------------------------- INTEGER N, LDZ, IPVT(*), KROW(*), IFLAG DOUBLE PRECISION Z(LDZ,*), X(*) * --------------------------------------- * * July 9, 2004 -- removed unused dummy variable W * * Purpose: * This routine is used to multiply the row vector X * by the accumulated elementary transformations stored in * the array Z. See routine BHAP1 for further comments about Z. * * Arguments: * N -integer * The number of rows and columns in the matrix Z. * N >= 0. * * Z -double precision array, dimension (LDZ,N) * Z contains information about the reduction to * tridiagonal form. * * LDZ -integer * The leading dimension of the array Z. * LDZ >= max(1,N). * * X -double precision array, dimension (N) * The vector used in applying the transformations. * * IPIV -integer array, dimension (N) * Contains the pivot sequence used during the * reduction to tridiagonal form. * * KROW -integer array, dimension (N) * KR(K) was the row eliminated after the kth column * eliminated. If no row was eliminated KR(K) = 0. * * IFLAG -integer * Not currently used. Could use to signal exceptions. * * .. Parameters .. * -------------------------------- DOUBLE PRECISION TWO, ZERO PARAMETER (TWO = 2.0, ZERO = 0.0) * --------------------------------- * .. * .. Local Scalars .. * --------------------------------- INTEGER K, L, KP1 DOUBLE PRECISION M, DDOT EXTERNAL DDOT * --------------------------------- * .. * .. Executable Statements .. * * Apply the permutations. * --------------------------------- DO 30 K = 1, N-2 KP1 = K+1 L = IPVT(KP1) M = X(L) IF (L .EQ. KP1) GO TO 30 X(L) = X(KP1) X(KP1) = M 30 CONTINUE DO K=1, N-2 KP1 = K+1 * ------------------------------------------------------- * Postmultiply the row vector x by Lr(k). * ------------------------------------------------------- IF(KROW(K).NE.0)THEN CALL DAXPY(N-K-1,-X(KP1),Z(KROW(K),K+2),LDZ,X(K+2),1) ENDIF * -------------------------------------------------------- * Postmultiply the row vector x by inv(Lc(k)). * -------------------------------------------------------- X(KP1) = X(KP1) + DDOT(N-K-1,Z(K+2,K),1,X(K+2),1) END DO * RETURN END ! End of BHAP2 *=================================================================== SUBROUTINE BHAP3(N,Z,LDZ,XRE,XIM,IPVT,KROW,IFLAG) * ---------------------------------------------------- * * Global Variables * ---------------------------------------------------- INTEGER N, LDZ, IPVT(*), KROW(*), IFLAG DOUBLE PRECISION Z(LDZ,*), XRE(*), XIM(*) * ---------------------------------------------------- * * July 9, 2004 -- removed unused dummy variable W. * * Purpose: * This routine used to multiply the accumulated elementary * transformations stored in the array Z by the column vector * (XRE,XIM). See the routine BHAP1 for comments about the * array Z. * * Arguments: * N -integer * The number of rows and columns in the matrix Z. * N >= 0. * * Z -double precision array, dimension (LDZ,N) * Z contains information about the reduction to * tridiagonal form. * * LDZ -integer * The leading dimension of the array Z. * LDZ >= max(1,N). * * XRE -double precision array, dimension (N) * The real part of the vector used in applying the * transformations. * * XIM -double precision array, dimension (N) * The imaginary part of the vector used in applying the * transformations. * * IPIV -integer array, dimension (N) * Contains the pivot sequence used during the reduction * to tridiagonal form. * * KROW -integer array, dimension (N) * KROW(K) was the row eliminated after the kth column * eliminated. If no row was eliminated KROW(K) = 0. * * IFLAG -integer * Not currently used. Can be used to throw * an exception. * * .. Parameters .. * ------------------------------------------ DOUBLE PRECISION TWO,ZERO PARAMETER (TWO = 2.0, ZERO = 0.0) * --