#!/bin/sh # This is a shell archive (produced by shar 3.49) # To extract the files from this archive, save it to a file, remove # everything above the "!/bin/sh" line above, and type "sh file_name". # # existing files will NOT be overwritten unless -c is specified # # This shar contains: # length mode name # ------ ---------- ------------------------------------------ # 24247 -rw------- dble/README # 5212 -rw------- dble/makefile # 5850 -rw------- dble/v.dat # 5850 -rw------- dble/w.dat # 959 -rw------- dble/local.mak # 2412 -rw------- dble/skeleton.mak # 5400 -rw------- dble/x.dat # 20 -rw------- dble/version.h # 770 -rw------- dble/bru/axb.inc # 1579 -rw------- dble/bru/cpyrit.doc # 9536 -rw------- dble/bru/daxb.f # 2013 -rw------- dble/bru/daxbpr.f # 1077 -rw------- dble/bru/makefile # 1579 -rw------- dble/csr/cpyrit.doc # 999 -rw------- dble/csr/csr.inc # 8640 -rw------- dble/csr/dcilut.f # 7610 -rw------- dble/csr/dclilut.f # 6513 -rw------- dble/csr/dclssor.f # 7610 -rw------- dble/csr/dcrilut.f # 6521 -rw------- dble/csr/dcrssor.f # 10293 -rw------- dble/csr/dcsr.f # 9872 -rw------- dble/csr/ilut.f # 2395 -rw------- dble/csr/makefile # 6422 -rw------- dble/csr/sparskit.f # 9603 -rw------- dble/dat/dns.dat # 480 -rw------- dble/dat/dnsv.dat # 480 -rw------- dble/dat/dnsw.dat # 29298 -rw------- dble/dat/csr.dat # 5400 -rw------- dble/dat/csrx.dat # 5850 -rw------- dble/dat/csrv.dat # 5850 -rw------- dble/dat/csrw.dat # 1579 -rw------- dble/dns/cpyrit.doc # 5671 -rw------- dble/dns/ddlssor.f # 6294 -rw------- dble/dns/ddns.f # 5575 -rw------- dble/dns/ddrssor.f # 750 -rw------- dble/dns/dns.inc # 1600 -rw------- dble/dns/makefile # 224 -rw------- dble/inc/dimblk.inc # 155 -rw------- dble/inc/precon.inc # 1579 -rw------- dble/lal/cpyrit.doc # 2725 -rw------- dble/lal/dcoeff.f # 1597 -rw------- dble/lal/dlal.inc # 13813 -rw------- dble/lal/deig.F # 13177 -rw------- dble/lal/deiglal.F # 23494 -rw------- dble/lal/dlal.F # 1266 -rw------- dble/lal/getomg.f # 8607 -rw------- dble/lal/dsys.F # 28685 -rw------- dble/lal/dsysbcg.F # 26481 -rw------- dble/lal/dsyslal.F # 1502 -rw------- dble/lal/makefile # 1558 -rw------- dble/sup/cpyrit.doc # 42496 -rw------- dble/sup/linpack.f # 3728 -rw------- dble/sup/support.f # 38090 -rw------- dble/sup/eispack.f # 1016 -rw------- dble/sup/makefile # # ============= dble/README ============== if test ! -d 'dble'; then echo 'x - creating directory dble' mkdir 'dble' fi if test -f 'dble/README' -a X"$1" != X"-c"; then echo 'x - skipping dble/README (File already exists)' else echo 'x - extracting dble/README (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/README' && X X X December 31, 1993 X Update Info ----------- X The codes in this package have been upgraded and expanded to a new package, called QMRPACK. We retain this package only for historical reasons, and recommend the use of the codes in QMRPACK. The new package is distributed in the form of the compressed tar file "qmrpack.tar.Z" in the "linalg" section of NETLIB. You can obtain QMRPACK by using xnetlib or by anonymous ftp. In the latter case, ftp to "netlib.att.com" and then get the file "qmrpack.tar.Z" from the directory "netlib/linalg". X The package provides two different implementations of QMR, one based on three-term recurrences, and one based on coupled two-term recurrences. Since several people find it desirable to have iterative methods that require only a few lines of code, we have also included versions of QMR without look-ahead. In particular, there is code for "QMR from BCG", which generates QMR by simply adding one extra SAXPY to each BCG iteration. However, in view of their enhanced stability, we recommend the use of the "true" QMR methods with look-ahead. The package also contains a no-look-ahead version of TFQMR, which is a transpose-free variant of QMR. The package comes with two preconditioners: SSOR and a variant of Youcef Saad's ILUT preconditioner. Finally, the package also includes code for computing eigenvalues of nonsymmetric matrices, using the look-ahead Lanczos algorithm. X ------ X X December 31, 1991 X General Info ------------ X This is the final version of the double precision real Lanczos/QMR codes. It is a major rewrite of the previous version. The strategy used to update the QMR iterate has been changed considerably, and the SVD is no longer used to invert the diagonal blocks, as it resulted in too much roundoff being introduced. On the other hand, the code still does not reuse already computed matrix-vector products when blocks are rebuilt. When a block is rebuilt, the matrix-vector products spent in the discarded portion at the end of the block are wasted. In theory, it is possible to reuse these matrix-vector products, so that N steps of the algorithm require exactly N multiplications by A and A^T each. While we know how this would be done, implementing it robustly is a major headache and we deemed the gains to be insignificant (read: too lazy to do this). X As promised, this release does contain a QMR/BCG code which generates the existing BCG iterates stably from the QMR process; see the makefile options for the corresponding entries. X The distribution consists of FORTRAN files in several subdirectories. The code is actually not standard FORTRAN; one has to run the C preprocessor on the individual files to obtain FORTRAN-77 code. This is done so that we could code easier, in that various pieces of the work arrays are given macro names, resulting in code that says, for example, X X CALL AXB (_V_(N), _V_(N+1)) X rather than the FORTRAN version X X CALL AXB (VW(1,Q(N)), VW(1,Q(N+1))) X We'd rather work with the former! Some compilers, such as the Sun compiler, automatically run the C preprocessor on the input file, if it has the extension .F (rather than the default .f extension). In any case, the makefiles provided are set up to run the C preprocessor `cpp' on the .F files to produce standard .f files. Other than having to correctly specify in the `skeleton.mak' file the path to the C preprocessor, this should be transparent to you. If you run into problems, let us know, and we can send you the plain FORTRAN files. X Once you have unpacked the distribution (by now that should be done), you are ready to solve systems in the formats provided, or you may add your own formats. On a Sun-4 running SunOS v4.0 (and possibly on others as well), the makefile should help you get started. You can type `make help' or simply `make', and it will list its options. To compile the package on other systems, you might have to change the names and/or directories for some of the programs used by the makefile (see the "Various programs..." section of skeleton.mak). The current skeleton.mak file is set up for our (possibly non-standard) system setup (f77 in /usr/lang, cpp in /lib, etc). As an example, on the Cray-2 here at NASA Ames, the following settings apply: X CPP = /lib/cpp FC = /bin/cf77 FFLAGS = X The current version of the codes is listed in the file version.h. Minor bug fixes and modifications will be available as diffs relative to the files in a standard distribution, so you might want to keep an original copy around. Please always reference the version number whenever you have a question or comment about the codes. X In this day and age, some legalese is necessary, so here it goes. Please note that the codes are copyrighted. This was done for two purposes: first, we want to make sure that you (and your lawyers!) understand that we do not warrant these codes to do anything at all. We are distributing them for free, and think they don't do anything whatsoever. For all intent and purpose, any description of what the codes are doing should be construed as being a note of what we thought the codes did on our machine on a particular Tuesday of last year. They might do the same thing again someday, and then again they might not (and, if history is any indication, they probably won't). Having said that, the second purpose is that we have nevertheless invested time and energy in these codes, and while we do not mind you using them for research, we do not want you selling them. So, if you want to make any profit from these codes, you have to have our permission. You are certainly allowed to use the codes for your own research. You are also allowed to distribute the codes, as long as you do not charge for this more than the cost of the media and a reasonable handling fee; an example of what we mean by ``reasonable'' is something not more than three times the current U.S. minimum wage (around $5/hr in 1991). But you are not allowed to sell any part of these codes, either alone or incorporated in some product you might design. You are, of course, welcome to code up your own versions of these codes. X Otherwise, we would be interested in hearing about your experience with the codes, especially bugs (what bugs?), smashing success stories, and stunning failures. We can be contacted at na.freund@na-net.ornl.gov (Roland Freund) and na.nachtigal@na-net.ornl.gov (Noel Nachtigal). Also note that we make no claims about the support codes. For example, the ILUT routine is not particularly efficient; it did the job for us, and we did not spend any time improving it. Finally, if you use these codes in your research, please reference one or the other of the following technical reports, as appropriate: X (For eigenvalue computations) Roland W. Freund, Martin H. Gutknecht, and Noel M. Nachtigal. An Implementation of the Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices. Technical Report 91.09, RIACS, NASA Ames Research Center, April 1991. To appear in the SIAM Journal on Statistical and Scientific Computing. X (For linear systems -- QMR) Roland W. Freund and Noel M. Nachtigal. QMR: a Quasi-Minimal Residual Method for Non-Hermitian Linear System. Technical Report 90.51, RIACS, NASA Ames Research Center, December 1990. To appear in Numerische Mathematik. X In case you wish to find out more about how the codes work, the following references describe it in some detail: X (For the heuristics used in the eigenvalue computations) Jane Cullum and Ralph A. Willoughby. A Practical Procedure for Computing Eigenvalues of Large Sparse Nonsymmetric Matrices. In Large Scale Eigenvalue Problems (J. Cullum and R.A. Willoughby, eds), North-Holland, 1986, pp.193--240. X (For a description of the Harwell-Boeing format) I.S. Duff, R.G Grimes, and J.G. Lewis. Sparse Matrix Test Problems. ACM Transactions of Mathematical Software, vol 15, 1989, pp. 1--14. X (For the basic implementation of the look-ahead Lanczos algorithm) Roland W. Freund, Martin H. Gutknecht, and Noel M. Nachtigal. An Implementation of the Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices, Part I. Technical Report 90.45, RIACS, NASA Ames Research Center, November 1990. X (For the QMR algorithm) Roland W. Freund and Noel M. Nachtigal. An Implementation of the Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices, Part II. Technical Report 90.46, RIACS, NASA Ames Research Center, November 1990. X Roland W. Freund and Noel M. Nachtigal. QMR: a Quasi-Minimal Residual Method for Non-Hermitian Linear System. Technical Report 90.51, RIACS, NASA Ames Research Center, December 1990. X (For a description of the SPARSKIT package) Youcef Saad. SPARSKIT: a Basic Tool Kit for Sparse Matrix Computations. Technical Report 90.20, RIACS, NASA Ames Research Center, May 1990. X Finally, the codes in the sup directory are support routines required from the LINPACK and EISPACK libraries. These should be replaced with the local copies of the libraries. In addition, the DNRM2 routine from the linpack.f file has a DATA statement which follows the variable declaration. This is not standard FORTRAN and some compilers might complain about it. We did not modify it ourselves, since the routine was taken from the LINPACK library. X X An Example ---------- X OK. Suppose now you want to solve a system. We have provided two as examples, a small 3x3 dense system to indicate how dense matrices are stored, and an example in Harwell-Boeing format. To solve the example Harwell-Boeing system, type `make scsr', then, when the compilation completes, type `dsys' and answer the prompts as follows (your answers are denoted with `<=='): X Enter sparse data file name: dat/csr.dat <== NDIM : 250 NZMAX : 1500 JOB : 2 NRHS : 0 GUESOL: NROW : 225 NCOL : 225 NNZ : 1065 TITLE : 7-POINT TEST MATRIX FROM SPARSKIT KEY : rua TYPE : 7-P IERR : 0 Enter estimated matrix norm : 1.0e3 <== Enter convergence tolerance : 1.0e-6 <== Maximum number of steps NLIM : 65 <== Preconditioner: ILUT Precondition (1=Yes, 0=No) ? 0 <== X We ran this example on a Sun-4 and on a pair of Crays (Cray-2 and Cray YMP). The results for the Sun-4 were: X 1 1 0.1000E+01 0.1000E+01 0.1000E+01 0.1000E+04 X 2 1 0.1145E+01 0.8705E+00 0.9874E-01 0.1000E+04 X 3 1 0.1402E+01 0.8691E+00 0.6632E-01 0.1000E+04 X 4 1 0.1590E+01 0.8258E+00 0.2191E-01 0.1000E+04 X 5 1 0.1750E+01 0.7988E+00 0.1657E-01 0.1000E+04 X 6 1 0.1887E+01 0.7921E+00 0.1473E-02 0.1000E+04 X 7 1 0.2032E+01 0.7572E+00 0.1525E-02 0.1000E+04 X 8 1 0.2131E+01 0.7577E+00 0.1317E-01 0.1000E+04 X 9 1 0.2156E+01 0.7355E+00 0.1040E-01 0.1000E+04 X 10 1 0.2233E+01 0.7715E+00 0.8337E-02 0.1000E+04 X 11 1 0.2313E+01 0.7937E+00 0.7020E-02 0.1000E+04 X 12 1 0.2379E+01 0.8358E+00 0.7696E-02 0.1000E+04 X 13 1 0.2475E+01 0.8313E+00 0.1887E-02 0.1000E+04 X 14 1 0.2542E+01 0.8907E+00 0.4649E-02 0.1000E+04 X 15 1 0.2629E+01 0.8819E+00 0.2295E-02 0.1000E+04 X 16 1 0.2697E+01 0.8950E+00 0.4872E-02 0.1000E+04 X 17 1 0.2775E+01 0.9051E+00 0.1411E-02 0.1000E+04 X 18 1 0.2838E+01 0.9258E+00 0.1398E-01 0.1000E+04 X 19 1 0.2901E+01 0.9417E+00 0.1173E-02 0.1000E+04 X 20 1 0.2949E+01 0.9872E+00 0.3449E-02 0.1000E+04 X 21 1 0.2995E+01 0.1035E+01 0.1433E-02 0.1000E+04 X 22 1 0.3041E+01 0.1089E+01 0.2032E-02 0.1000E+04 X 23 1 0.3082E+01 0.1144E+01 0.5925E-02 0.1000E+04 X 24 1 0.3113E+01 0.1205E+01 0.1884E-02 0.1000E+04 X 25 1 0.3010E+01 0.1079E+01 0.2231E-02 0.1000E+04 X 26 1 0.2773E+01 0.9390E+00 0.1533E-01 0.1000E+04 X 27 1 0.2444E+01 0.6838E+00 0.1130E-01 0.1000E+04 X 28 1 0.2409E+01 0.6322E+00 0.2135E-01 0.1000E+04 X 29 1 0.2452E+01 0.6323E+00 0.2038E-01 0.1000E+04 X 30 1 0.2059E+01 0.4352E+00 0.4548E-02 0.1000E+04 X 31 1 0.2006E+01 0.3896E+00 0.3308E-06 0.1000E+04 X 32 1 0.1323E+01 0.2407E+00 0.1692E-06 0.1000E+04 X 33 1 0.1335E+01 0.2222E+00 0.6722E-06 0.1000E+04 X 34 1 0.1336E+01 0.1932E+00 0.7657E-08 0.1000E+04 X 35 1 0.1354E+01 0.1985E+00 0.4373E-08 0.1000E+04 X 36 1 0.1370E+01 0.2062E+00 0.2437E-09 0.1000E+04 X 37 1 0.1389E+01 0.2047E+00 0.2836E-09 0.1000E+04 X 38 1 0.1390E+01 0.1935E+00 0.2421E-08 0.1000E+04 X 39 1 0.1323E+01 0.1399E+00 0.8416E-09 0.1000E+04 X 40 1 0.1340E+01 0.1398E+00 0.3249E-08 0.1000E+04 X 41 1 0.1356E+01 0.1397E+00 0.2340E-09 0.1000E+04 X 42 1 0.1303E+01 0.1370E+00 0.3690E-09 0.1000E+04 X 43 1 0.1141E+01 0.1609E+00 0.1123E-08 0.1000E+04 X 44 1 0.1142E+01 0.1783E+00 0.4717E-09 0.1000E+04 X 45 1 0.1008E+01 0.1986E+00 0.2442E-09 0.1000E+04 X 46 1 0.3137E+00 0.3047E-01 0.5367E-09 0.1000E+04 X 47 1 0.2382E+00 0.3412E-01 0.3704E-08 0.1000E+04 X 48 1 0.9779E-01 0.1188E-01 0.7983E-09 0.1000E+04 X 49 1 0.9109E-01 0.1374E-01 0.3691E-08 0.1000E+04 X 50 1 0.5057E-01 0.1008E-01 0.3673E-09 0.1000E+04 X 51 1 0.1005E-01 0.1499E-02 0.5906E-09 0.1000E+04 X 52 1 0.5223E-02 0.3446E-03 0.4439E-08 0.1000E+04 X 53 1 0.1170E-03 0.1614E-04 0.3954E-08 0.1000E+04 X 54 1 0.7184E-04 0.1322E-04 0.1412E-06 0.1000E+04 X 55 1 0.8707E-05 0.1349E-05 0.2749E-06 0.1000E+04 X 56 1 0.6942E-05 0.1346E-05 0.2776E-05 0.1000E+04 X 57 1 0.5192E-05 0.1097E-05 0.2473E-05 0.1000E+04 X 58 1 0.8281E-06 0.1321E-06 0.4710E-06 0.1000E+04 The residual norm has converged. Play it again (Y/N) ? n <== X The results for both Crays were: X 1 1 0.1000E+01 0.1000E+01 0.1000E+01 0.1000E+04 X 2 1 0.1145E+01 0.8705E+00 0.9874E-01 0.1000E+04 X 3 1 0.1402E+01 0.8691E+00 0.6632E-01 0.1000E+04 X 4 1 0.1590E+01 0.8258E+00 0.2191E-01 0.1000E+04 X 5 1 0.1750E+01 0.7988E+00 0.1657E-01 0.1000E+04 X 6 1 0.1887E+01 0.7921E+00 0.1473E-02 0.1000E+04 X 7 1 0.2032E+01 0.7572E+00 0.1525E-02 0.1000E+04 X 8 1 0.2131E+01 0.7577E+00 0.1317E-01 0.1000E+04 X 9 1 0.2156E+01 0.7355E+00 0.1040E-01 0.1000E+04 X 10 1 0.2233E+01 0.7715E+00 0.8337E-02 0.1000E+04 X 11 1 0.2313E+01 0.7937E+00 0.7020E-02 0.1000E+04 X 12 1 0.2379E+01 0.8358E+00 0.7696E-02 0.1000E+04 X 13 1 0.2475E+01 0.8313E+00 0.1887E-02 0.1000E+04 X 14 1 0.2542E+01 0.8907E+00 0.4649E-02 0.1000E+04 X 15 1 0.2629E+01 0.8819E+00 0.2295E-02 0.1000E+04 X 16 1 0.2697E+01 0.8950E+00 0.4872E-02 0.1000E+04 X 17 1 0.2775E+01 0.9051E+00 0.1411E-02 0.1000E+04 X 18 1 0.2838E+01 0.9258E+00 0.1398E-01 0.1000E+04 X 19 1 0.2901E+01 0.9417E+00 0.1173E-02 0.1000E+04 X 20 1 0.2949E+01 0.9872E+00 0.3449E-02 0.1000E+04 X 21 1 0.2995E+01 0.1035E+01 0.1433E-02 0.1000E+04 X 22 1 0.3041E+01 0.1089E+01 0.2032E-02 0.1000E+04 X 23 1 0.3082E+01 0.1144E+01 0.5925E-02 0.1000E+04 X 24 1 0.3113E+01 0.1205E+01 0.1884E-02 0.1000E+04 X 25 1 0.3010E+01 0.1079E+01 0.2231E-02 0.1000E+04 X 26 1 0.2773E+01 0.9390E+00 0.1533E-01 0.1000E+04 X 27 1 0.2444E+01 0.6838E+00 0.1130E-01 0.1000E+04 X 28 1 0.2409E+01 0.6322E+00 0.2135E-01 0.1000E+04 X 29 1 0.2452E+01 0.6323E+00 0.2038E-01 0.1000E+04 X 30 1 0.2059E+01 0.4352E+00 0.4548E-02 0.1000E+04 X 31 1 0.2006E+01 0.3896E+00 0.3308E-06 0.1000E+04 X 32 1 0.1323E+01 0.2407E+00 0.1692E-06 0.1000E+04 X 33 1 0.1335E+01 0.2222E+00 0.6722E-06 0.1000E+04 X 34 1 0.1336E+01 0.1932E+00 0.7657E-08 0.1000E+04 X 35 1 0.1354E+01 0.1985E+00 0.4373E-08 0.1000E+04 X 36 1 0.1370E+01 0.2062E+00 0.2437E-09 0.1000E+04 X 37 1 0.1389E+01 0.2047E+00 0.2836E-09 0.1000E+04 X 38 1 0.1390E+01 0.1935E+00 0.2421E-08 0.1000E+04 X 39 1 0.1323E+01 0.1399E+00 0.8416E-09 0.1000E+04 X 40 1 0.1340E+01 0.1398E+00 0.3249E-08 0.1000E+04 X 41 1 0.1356E+01 0.1397E+00 0.2340E-09 0.1000E+04 X 42 1 0.1303E+01 0.1370E+00 0.3690E-09 0.1000E+04 X 43 1 0.1141E+01 0.1609E+00 0.1123E-08 0.1000E+04 X 44 1 0.1142E+01 0.1783E+00 0.4717E-09 0.1000E+04 X 45 1 0.1008E+01 0.1986E+00 0.2442E-09 0.1000E+04 X 46 1 0.3136E+00 0.3047E-01 0.5367E-09 0.1000E+04 X 47 1 0.2382E+00 0.3412E-01 0.3704E-08 0.1000E+04 X 48 1 0.9779E-01 0.1188E-01 0.7984E-09 0.1000E+04 X 49 1 0.9109E-01 0.1374E-01 0.3691E-08 0.1000E+04 X 50 1 0.5057E-01 0.1008E-01 0.3673E-09 0.1000E+04 X 51 1 0.1005E-01 0.1499E-02 0.5906E-09 0.1000E+04 X 52 1 0.5222E-02 0.3446E-03 0.4439E-08 0.1000E+04 X 53 1 0.1167E-03 0.1611E-04 0.3954E-08 0.1000E+04 X 54 1 0.7310E-04 0.1350E-04 0.1437E-06 0.1000E+04 X 55 1 0.9228E-05 0.1436E-05 0.2661E-06 0.1000E+04 X 56 1 0.7685E-05 0.1361E-05 0.2075E-05 0.1000E+04 X 57 1 0.3302E-05 0.5971E-06 0.7132E-06 0.1000E+04 The residual norm has converged. Play it again (Y/N) ? n <== X In both cases, the exact solution is the vector of all 1's. You can check the solution returned in the file `x.out' and see how it compares. That's it. To solve bigger systems, you will need to adjust the parameters in `dsys.f', `dimblk.inc', and `csr.inc', but the basic idea remains the same. For the dense system provided in the subdirectory `dns', the exact solution is X X 4.441959994855718e+01 X -3.071517860902447e+02 X -1.857585177589905e+02 X 6.304596860549932e+01 X 7.078107514016136e+02 X 2.004530083313206e+02 X -1.133227360733343e+01 X -3.541347098296361e+02 X 4.374106238114270e+02 X -2.400890005920496e+02 X -3.911009425579379e+02 X -2.798321395307324e+01 X -2.918743864606955e+01 X 1.117419236671314e+01 X -1.145266987648990e+02 X 2.823876794263498e+02 X 3.625516005161455e+01 X -4.460202949059404e+01 X -1.039938190024286e+02 X -1.653064614588225e+02 X It is left as an exercise to the reader to set up and solve the dense linear system. For eigenvalues, use `make ecsr' for Harwell-Boeing matrices, `make edns' for dense matrices. Remember that for eigenvalues, you will eventually be solving an eigenvalue problem whose size equals the number of steps taken, so don't get carried away... X X Changes history --------------- X X Version Comments X 1.0 Initial release X 1.1 Corrected bugs reported by Claude Pommerell (pommy@iis.ethz.ch) X in dsys.f, dlal.f, dsyslal.f X 1.2 Corrected bugs reported by Marlis Hochbruck (marlis@riacs.edu) X in dlal.f. X Changed the SVD tolerance from 0 to eps in deiglal.f, dsyslal.f X Changed the initial norm estimation in dlal.f X 1.3 Corrected bug in dlal.f X 1.4 Updated the README and makefile files X Cleaned up the daxbpr.f file X 1.5 Corrected bugs reported by Eric Lucas @ Westinghouse X 2.0 Final release, new strategy for QMR update, switched from SVD X to QR X Distribution Tree for version 2 ------------------------------- X Below is a listing of the distribution tree, with explanations for each of the subdirectories. X The main distribution directory: total 59 -rw-rw-rw- 1 santa 22816 Dec 31 08:00 README drwxrwxrwx 2 santa 512 Dec 31 08:00 bru/ drwxrwxrwx 2 santa 1024 Dec 31 08:00 csr/ drwxrwxrwx 2 santa 512 Dec 31 08:00 dat/ drwxrwxrwx 2 santa 512 Dec 31 08:00 dns/ drwxrwxrwx 2 santa 512 Dec 31 08:00 inc/ drwxrwxrwx 2 santa 1024 Dec 31 08:00 lal/ -rw-rw-rw- 1 santa 959 Dec 31 08:00 local.mak -rw-rw-rw- 1 santa 5212 Dec 31 08:00 makefile -rw-rw-rw- 1 santa 2412 Dec 31 08:00 skeleton.mak drwxrwxrwx 2 santa 512 Dec 31 08:00 sup/ -rw-rw-rw- 2 santa 5850 Dec 31 08:00 v.dat -rw-rw-rw- 1 santa 20 Dec 31 08:00 version.h -rw-rw-rw- 2 santa 5850 Dec 31 08:00 w.dat -rw-rw-rw- 2 santa 5400 Dec 31 08:00 x.dat X This directory contains example codes for a matrix-free problem. `daxb.f' has the code for the matrix-vector routines, and `daxbpr.f' has the code for the preconditioner routines (which in this example are empty). bru: total 17 -rw-rw-rw- 1 santa 770 Dec 31 08:00 axb.inc -rw-rw-rw- 1 santa 1579 Dec 31 08:00 cpyrit.doc -rw-rw-rw- 1 santa 9536 Dec 31 08:00 daxb.F -rw-rw-rw- 1 santa 2013 Dec 31 08:00 daxbpr.F -rw-rw-rw- 1 santa 1077 Dec 31 08:00 makefile X This directory contains support routines for the Harwell-Boeing format. `dcsr.f' contains the main matrix-vector routines (some taken from SPARSKIT), `sparskit.f' contains the remaining support routines from SPARSKIT. The remaining files contain routines for various preconditioners. csr: total 73 -rw-rw-rw- 1 santa 1579 Dec 31 08:00 cpyrit.doc -rw-rw-rw- 1 santa 999 Dec 31 08:00 csr.inc -rw-rw-rw- 1 santa 8640 Dec 31 08:00 dcilut.F -rw-rw-rw- 1 santa 7610 Dec 31 08:00 dclilut.F -rw-rw-rw- 1 santa 6513 Dec 31 08:00 dclssor.F -rw-rw-rw- 1 santa 7610 Dec 31 08:00 dcrilut.F -rw-rw-rw- 1 santa 6521 Dec 31 08:00 dcrssor.F -rw-rw-rw- 1 santa 10293 Dec 31 08:00 dcsr.F -rw-rw-rw- 1 santa 9872 Dec 31 08:00 ilut.F -rw-rw-rw- 1 santa 2638 Dec 31 08:00 makefile -rw-rw-rw- 1 santa 6422 Dec 31 08:00 sparskit.F X This directory contains example data files for the matrices. `dns.dat' has the data for the 3x3 matrix X A = [ 1 2 3; 5 7 11; 13 17 19 ], while `csr.dat' contains the data for a sparse routine in Harwell-Boeing format. The other files are the right-hand sides (`csrv.dat' and `dnsv.dat'), the second starting vector (`csrw.dat' and `dnsw.dat'), and the starting guess (`csrx.dat' for both). dat: total 59 -rw-rw-rw- 1 santa 29298 Dec 31 08:00 csr.dat -rw-rw-rw- 2 santa 5850 Dec 31 08:00 csrv.dat -rw-rw-rw- 2 santa 5850 Dec 31 08:00 csrw.dat -rw-rw-rw- 2 santa 5400 Dec 31 08:00 csrx.dat -rw-rw-rw- 1 santa 9603 Dec 31 08:00 dns.dat -rw-rw-rw- 2 santa 480 Dec 31 08:00 dnsv.dat -rw-rw-rw- 2 santa 480 Dec 31 08:00 dnsw.dat X This directory contains support routines for the dense format. 'ddns.f' has the main matrix-vector routines, while the remaining files contain routines for various preconditioners. dns: total 24 -rw-rw-rw- 1 santa 1579 Dec 31 08:00 cpyrit.doc -rw-rw-rw- 1 santa 5671 Dec 31 08:00 ddlssor.F -rw-rw-rw- 1 santa 6294 Dec 31 08:00 ddns.F -rw-rw-rw- 1 santa 5575 Dec 31 08:00 ddrssor.F -rw-rw-rw- 1 santa 750 Dec 31 08:00 dns.inc -rw-rw-rw- 1 santa 1600 Dec 31 08:00 makefile X This directory contains the main include files which govern the maximum dimension of the matrix (`dimblk.inc') and the common block declaration for the preconditioner name. All other files `dimblk.inc' and `precon.inc' (there is a pair in every other subdirectory) are linked to these two. inc: total 2 -rw-rw-rw- 1 santa 224 Dec 31 08:00 dimblk.inc -rw-rw-rw- 1 santa 155 Dec 31 08:00 precon.inc X This is the directory which contains the Lanczos routines. `dlal.f' is the main low-level Lanczos routine; it calls `dcoeff.f'. For linear systems, `dsys.f' is the example driver code, which calls the QMR routines `dsyslal.f', which in turn calls `dlal.f' and `getomg.f'. For eigenvalues, `deig.f' is the example driver code, which calls `deiglal.f', which in turn calls `dlal.f'. lal: total 127 -rw-rw-rw- 1 santa 1579 Dec 31 08:00 cpyrit.doc -rw-rw-rw- 1 santa 2725 Dec 31 08:00 dcoeff.F -rw-rw-rw- 1 santa 13813 Dec 31 08:00 deig.F -rw-rw-rw- 1 santa 13177 Dec 31 08:00 deiglal.F -rw-rw-rw- 1 santa 23494 Dec 31 08:00 dlal.F -rw-rw-rw- 1 santa 1597 Dec 31 08:00 dlal.inc -rw-rw-rw- 1 santa 8607 Dec 31 08:00 dsys.F -rw-rw-rw- 1 santa 28685 Dec 31 08:00 dsysbcg.F -rw-rw-rw- 1 santa 26481 Dec 31 08:00 dsyslal.F -rw-rw-rw- 1 santa 1266 Dec 31 08:00 getomg.F -rw-rw-rw- 1 santa 1502 Dec 31 08:00 makefile X This directory contains miscellaneous support routines from EISPACK, LINPACK, and a sorting routine in `support.f'. sup: total 87 -rw-rw-rw- 1 santa 1558 Dec 31 08:00 cpyrit.doc -rw-rw-rw- 1 santa 38090 Dec 31 08:00 eispack.f -rw-rw-rw- 1 santa 42496 Dec 31 08:00 linpack.f -rw-rw-rw- 1 santa 1016 Dec 31 08:00 makefile -rw-rw-rw- 1 santa 3728 Dec 31 08:00 support.f SHAR_EOF chmod 0600 dble/README || echo 'restore of dble/README failed' Wc_c="`wc -c < 'dble/README'`" test 24247 -eq "$Wc_c" || echo 'dble/README: original size 24247, current size' "$Wc_c" fi # ============= dble/makefile ============== if test -f 'dble/makefile' -a X"$1" != X"-c"; then echo 'x - skipping dble/makefile (File already exists)' else echo 'x - extracting dble/makefile (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/makefile' && #********************************************************************** # # Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal # All rights reserved. # # This code is part of a copyrighted package. For details, see the # file `cpyrit.doc' in the `lal' directory. # # ***************************************************************** # ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE # COPYRIGHT NOTICE # ***************************************************************** # #********************************************************************** X # # NOTE: GNU make seems to choke and gag on this makefile. # include skeleton.mak X # # Various compilation directories. # DIRS = bru csr dns lal sup AXB = bru X # # Files included in the distribution package. # TOP = dble REL = $(TOP)/README $(TOP)/bcg $(TOP)/bru $(TOP)/csr $(TOP)/dat $(TOP)/dns $(TOP)/gmr $(TOP)/grc $(TOP)/inc $(TOP)/lal $(TOP)/makefile $(TOP)/sup $(TOP)/*.dat $(TOP)/*.mak X # # This is the local help target. # lochelp: X @$(ECHO) " noout - remove output files" X @$(ECHO) " rel - make the release tar file" X @$(ECHO) " " X @$(ECHO) " baxb - generic format, QMR/BCG" X @$(ECHO) " bcsr - compressed sparse row format, QMR/BCG" X @$(ECHO) " bdns - dense format, QMR/BCG" X @$(ECHO) " eaxb - generic format, Lanczos, ew's" X @$(ECHO) " ecsr - compressed sparse row format, Lanczos, ew's" X @$(ECHO) " edns - dense format, Lanczos, ew's" X @$(ECHO) " saxb - generic format, QMR" X @$(ECHO) " scsr - compressed sparse row format, QMR" X @$(ECHO) " sdns - dense format, QMR" X FEIG = eaxb ecsr edns FLIN = baxb bcsr bdns caxb ccsr cdns gaxb gcsr gdns saxb scsr sdns X # # for - make all FORTRAN files # nofor - remove all FORTRAN files # obj - make all object files # noobj - remove all object files # for nofor obj noobj: X @( \ X for f in $(DIRS); \ X do $(CD) $$f; \ X $(MAKE) $@; \ X $(CD) ..; \ X done ) X # # none - remove everything # none: noout X @( \ X for f in $(DIRS); \ X do $(CD) $$f; \ X $(MAKE) nofor noobj; \ X $(CD) ..; \ X done ) X # # Remove all output files. # noout: X @$(RM) $(RMFLAGS) $(FEIG) $(FLIN) deig dsys *.out X # # Distribution files. # rel: nofor noobj X @( \ X $(CD) ..; \ X $(RM) $(RMFLAGS) $(TOP).T; \ X $(ECHO) Copying release files; \ X $(TAR) cFFhf $(TOP).T $(REL); \ X $(ECHO) Compressing the tar file; \ X $(COMPRESS) $(TOP).T; \ X $(MV) $(MVFLAGS) $(TOP).T.Z $(TOP) ) X # # The generic format, Lanczos-QMR/BCG, linear systems. # BAXBF = $(AXB)/daxb.o $(AXB)/daxbpr.o lal/dcoeff.o lal/dsys.o lal/dsysbcg.o \ X lal/dlal.o lal/getomg.o sup/linpack.o sup/support.o baxb: obj X @$(FC) -o dsys $(FFLAGS) $(LDFLAGS) $(BAXBF) X @$(RM) $(RMFLAGS) $(FLIN) X @$(ECHO) " " > baxb X # # The generic format, Lanczos, eigenvalues. # EAXBF = $(AXB)/daxb.o $(AXB)/daxbpr.o lal/dcoeff.o lal/deig.o lal/deiglal.o \ X lal/dlal.o sup/eispack.o sup/linpack.o sup/support.o eaxb: obj X @$(FC) -o deig $(FFLAGS) $(LDFLAGS) $(EAXBF) X @$(RM) $(RMFLAGS) $(FEIG) X @$(ECHO) " " > eaxb X # # The generic format, Lanczos-QMR, linear systems. # SAXBF = $(AXB)/daxb.o $(AXB)/daxbpr.o lal/dcoeff.o lal/dsys.o lal/dsyslal.o \ X lal/dlal.o lal/getomg.o sup/linpack.o sup/support.o saxb: obj X @$(FC) -o dsys $(FFLAGS) $(LDFLAGS) $(SAXBF) X @$(RM) $(RMFLAGS) $(FLIN) X @$(ECHO) " " > saxb X # # The compressed sparse row (CSR) format, Lanczos-QMR/BCG, linear systems. # BCSRF = csr/dcsr.o csr/dcsrpr.o csr/ilut.o csr/sparskit.o lal/dcoeff.o \ X lal/dsys.o lal/dsysbcg.o lal/dlal.o lal/getomg.o sup/linpack.o \ X sup/support.o bcsr: obj X @$(FC) -o dsys $(FFLAGS) $(LDFLAGS) $(BCSRF) X @$(RM) $(RMFLAGS) $(FLIN) X # # The compressed sparse row (CSR) format, Lanczos, eigenvalues. # ECSRF = csr/dcsr.o csr/dcsrpr.o csr/sparskit.o lal/dcoeff.o lal/deig.o \ X lal/deiglal.o lal/dlal.o sup/eispack.o csr/ilut.o sup/linpack.o \ X sup/support.o ecsr: obj X @$(FC) -o deig $(FFLAGS) $(LDFLAGS) $(ECSRF) X @$(RM) $(RMFLAGS) $(FEIG) X @$(ECHO) " " > ecsr X # # The compressed sparse row (CSR) format, Lanczos-QMR, linear systems. # SCSRF = csr/dcsr.o csr/dcsrpr.o csr/ilut.o csr/sparskit.o lal/dcoeff.o \ X lal/dsys.o lal/dsyslal.o lal/dlal.o lal/getomg.o sup/linpack.o \ X sup/support.o scsr: obj X @$(FC) -o dsys $(FFLAGS) $(LDFLAGS) $(SCSRF) X @$(RM) $(RMFLAGS) $(FLIN) X @$(ECHO) " " > scsr X # # The dense (DNS) format, Lanczos-QMR/BCG, linear systems. # BDNSF = dns/ddns.o dns/ddnspr.o lal/dcoeff.o lal/dsys.o lal/dsysbcg.o \ X lal/dlal.o lal/getomg.o sup/linpack.o sup/support.o bdns: obj X @$(FC) -o dsys $(FFLAGS) $(LDFLAGS) $(BDNSF) X @$(RM) $(RMFLAGS) $(FLIN) X @$(ECHO) " " > bdns X # # The dense (DNS) format, Lanczos, eigenvalues. # EDNSF = dns/ddns.o dns/ddnspr.o lal/dcoeff.o lal/deig.o lal/deiglal.o \ X lal/dlal.o sup/eispack.o sup/linpack.o sup/support.o edns: obj X @$(FC) -o deig $(FFLAGS) $(LDFLAGS) $(EDNSF) X @$(RM) $(RMFLAGS) $(FEIG) X @$(ECHO) " " > edns X # # The dense (DNS) format, Lanczos-QMR, linear systems. # SDNSF = dns/ddns.o dns/ddnspr.o lal/dcoeff.o lal/dsys.o lal/dsyslal.o \ X lal/dlal.o lal/getomg.o sup/linpack.o sup/support.o sdns: obj X @$(FC) -o dsys $(FFLAGS) $(LDFLAGS) $(SDNSF) X @$(RM) $(RMFLAGS) $(FLIN) X @$(ECHO) " " > sdns SHAR_EOF chmod 0600 dble/makefile || echo 'restore of dble/makefile failed' Wc_c="`wc -c < 'dble/makefile'`" test 5212 -eq "$Wc_c" || echo 'dble/makefile: original size 5212, current size' "$Wc_c" fi # ============= dble/v.dat ============== if test -f 'dble/v.dat' -a X"$1" != X"-c"; then echo 'x - skipping dble/v.dat (File already exists)' else echo 'x - extracting dble/v.dat (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/v.dat' && X 0.17265625000000000E+01 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.66796875000000000E+00 X 0.78906250000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.66796875000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.39062500000000000E+00 X -0.26953125000000000E+00 X 0.78906250000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.26953125000000000E+00 X -0.14843750000000000E+00 SHAR_EOF chmod 0600 dble/v.dat || echo 'restore of dble/v.dat failed' Wc_c="`wc -c < 'dble/v.dat'`" test 5850 -eq "$Wc_c" || echo 'dble/v.dat: original size 5850, current size' "$Wc_c" fi # ============= dble/w.dat ============== if test -f 'dble/w.dat' -a X"$1" != X"-c"; then echo 'x - skipping dble/w.dat (File already exists)' else echo 'x - extracting dble/w.dat (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/w.dat' && X -2.52670612285392693e+00 X -3.12981492764229430e-01 X -5.93617618531314117e-01 X 3.32322161717665654e-01 X 5.58850703293451434e-01 X 8.99883573546131688e-01 X -2.00898855586523289e-01 X -2.33734974678243862e-01 X 1.44990660185244824e+00 X 1.83613202602166203e+00 X -3.82918259160919427e-01 X 1.55082744517634041e-01 X -9.64648249056895501e-01 X 3.87564312898231872e-02 X 7.65458387074111402e-01 X -5.94524007839467128e-01 X 1.30245975190897978e-01 X 3.50135051326026445e-02 X -6.24674138678893653e-01 X -5.39775240747623841e-01 X 1.87995711256766906e+00 X -1.00384945408178328e+00 X -4.97445877674741299e-01 X -1.50439715265068541e+00 X -9.54492989363980443e-02 X 3.96727053727293111e-01 X -5.27114907886174100e-01 X 3.44571055586827601e-01 X -7.23290526415817259e-01 X 1.26819336319292564e+00 X -3.12426582958191153e-02 X 7.78211737167248230e-01 X 2.18048355295766028e+00 X 4.37813681537625787e-01 X 1.33332898358725349e+00 X 2.51078139110859744e-01 X -3.10470908178119731e-01 X -9.23003723808715093e-01 X -3.84775736018752812e-01 X 1.15818057089115922e+00 X 8.62500188414089375e-01 X -1.03470562492749885e+00 X -1.92672883298641989e-01 X -1.29972277987506701e+00 X 3.06595916028863658e-01 X 9.68992176157492779e-01 X -7.47317126758189398e-01 X -2.79602442297918508e+00 X 6.96731553587481622e-01 X 3.20690754366444031e+00 X 5.36007044749823192e-01 X 2.98450535720106380e-01 X 2.84043160995206323e-01 X 9.59664371348601719e-01 X 2.08759311209632337e+00 X 1.52468053170777473e+00 X -1.95260790032330211e-01 X 1.72603158008979921e-02 X 2.46340438601763606e-01 X -8.54484721089849630e-01 X 1.15778270176474307e+00 X 1.61907723250136332e-01 X 1.55706375548542164e+00 X -1.93543855158581107e-01 X 1.65130117445043179e+00 X -1.89877818089584816e+00 X 1.82252476317014889e+00 X -1.51841513419531315e+00 X -1.05107060879933423e+00 X 4.99305134332340902e-02 X -1.45474886952808413e+00 X 4.66545849753955855e-01 X 5.45436841452860643e-01 X 1.32031907308504048e+00 X -4.04494327876404325e-01 X 4.18468509073848915e-01 X 2.47348749631457726e-01 X 7.04110315408186027e-01 X 6.31938853341632800e-01 X -9.92362112719314848e-01 X 1.76670836879512128e+00 X -3.82103635072939263e-01 X -9.11425420031350630e-01 X -9.96089984118210370e-01 X 1.19514263014411903e+00 X -1.59447782443430802e-01 X 2.70402604824876036e+00 X -1.98499915965496987e-01 X -1.41404614026417491e-01 X 4.11267926557380703e-01 X -1.17905965667048696e+00 X -2.77775505971886494e-01 X -1.58105341380234377e+00 X 1.04902234978445863e+00 X 3.02689036171394199e-01 X -1.22650234105829847e+00 X 6.96000950977379157e-02 X -3.96516210293235416e-01 X 1.38880676152053462e+00 X 1.36442229049003072e+00 X 6.58152637292665821e-01 X 4.91313668926088520e-01 X 8.00733701087079197e-01 X -7.67268996676585435e-01 X 3.64419504046212817e-01 X -3.97913854767016573e-01 X 8.64279576409737516e-01 X -1.77618078276664149e-01 X 1.87438052046940795e+00 X 1.72400234691113430e-01 X 1.27174349438397227e+00 X -3.53443679957601356e-02 X -1.50132883642183579e+00 X 3.65373411191592334e-01 X -1.98659856001020180e-01 X -1.38972170325750999e+00 X 2.29327812227314215e-01 X 2.71190236967230214e-01 X -3.66360220282281213e-01 X 1.37696039157049333e+00 X -7.97532756562797762e-01 X -9.36740611780525367e-01 X -2.43346548885311044e-03 X 3.96086165525257827e-01 X -5.08693172275514027e-01 X -2.68285778746197690e-01 X -1.08214045362096933e+00 X 2.01413372029120419e+00 X 1.94403112557593682e+00 X -1.52152941634797623e+00 X 1.93931842629591666e+00 X -8.95840360657292223e-01 X -3.04157582743064425e-01 X 5.55253123177883778e-01 X -3.24246850701504052e-01 X 1.33881436714640256e+00 X 1.22229851347237273e+00 X -1.59597816278255955e+00 X -1.06773032044442528e+00 X -7.59919212299574154e-01 X 4.20988804468648503e-01 X -4.33373058325242422e-01 X 7.06251990240337246e-01 X 2.27856907314768625e-01 X -1.01699185125668268e+00 X 1.39860372563678254e-01 X -7.48088838235888787e-01 X -6.28974933137321557e-01 X 1.39483065417114704e+00 X -1.64769114004944095e+00 X -2.01498584386620427e+00 X 4.91716880786256527e-01 X -1.55497527509081546e+00 X -1.40609080683032323e-01 X 2.44943668795265301e-01 X -2.67458499968963315e-01 X -5.70245479900343022e-01 X -1.87266786888367814e-01 X 1.20855664796684303e+00 X -6.38854660397775276e-01 X 6.05540298516074937e-01 X -6.24480544088507727e-01 X 5.72228121730056660e-01 X -7.24410495952223288e-01 X 1.19219550553089348e+00 X 1.86746737068575697e-01 X 1.59493888226368430e+00 X 3.21307055691724686e-01 X 8.66840733181726275e-01 X 1.29184357610291878e+00 X 4.34312653452442632e-01 X -3.86206929335472016e-01 X -1.12563759811723021e-01 X -9.64333079249251268e-01 X -2.05725119297093961e+00 X 1.49996068326108345e-01 X 5.42037570810973812e-01 X 2.54408816480612421e-01 X -3.07240693819984811e-01 X -4.17111829581745308e-01 X 1.13680483289389689e+00 X 3.91313809093234544e-01 X 1.60514781867489997e+00 X 8.25892307356857591e-01 X 1.47039035572010768e+00 X -1.37890689233989572e+00 X -2.60172069009687479e-01 X 9.94768172763982217e-01 X 1.83403368186402838e+00 X -1.71591031873495248e+00 X 8.69317058746622712e-02 X 1.95567435281059465e+00 X 1.61453769615341497e-01 X -6.28688359125363805e-01 X -1.43882446533843478e+00 X -6.65959685875582436e-02 X 3.73380862806065694e-01 X 2.17314078186247150e-01 X -1.79456822070788363e-01 X 2.56729070095519964e-02 X 6.42066361973081534e-01 X 9.23086649379001090e-01 X -1.55510777372327769e+00 X 6.63594032788892285e-01 X -6.09499611051491308e-01 X 5.65239403309624411e-01 X -6.10781446255285299e-01 X 1.23111146649210634e+00 X 9.94299745127406931e-01 X -8.03474713644618865e-01 X -5.91204478397532762e-01 X 1.69154640779536125e+00 X 9.53355517329613988e-01 X -1.93005493739851142e+00 X 5.12844987283965770e-01 X 3.93682448572880705e-01 X -9.05426500446262272e-01 X -1.27447327679614819e+00 X 3.46546103379725356e-01 X -1.19523544023497297e+00 X 6.67201442318699822e-01 X -6.77937745985269097e-02 X -1.73566010510706481e+00 X 8.06348573332824392e-01 X -9.14800737775520956e-01 SHAR_EOF chmod 0600 dble/w.dat || echo 'restore of dble/w.dat failed' Wc_c="`wc -c < 'dble/w.dat'`" test 5850 -eq "$Wc_c" || echo 'dble/w.dat: original size 5850, current size' "$Wc_c" fi # ============= dble/local.mak ============== if test -f 'dble/local.mak' -a X"$1" != X"-c"; then echo 'x - skipping dble/local.mak (File already exists)' else echo 'x - extracting dble/local.mak (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/local.mak' && #********************************************************************** # # Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal # All rights reserved. # # This code is part of a copyrighted package. For details, see the # file `cpyrit.doc' in the `lal' directory. # # ***************************************************************** # ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE # COPYRIGHT NOTICE # ***************************************************************** # #********************************************************************** # # Additional skeleton for the local makefiles. # X # # Make all auxiliary files. # for: $(FOR) $(INC) X # # Clean up auxiliary files. # nofor: X @$(RM) $(RMFLAGS) $(FOR) $(INC) X # # Make all object files. # obj: $(OBJ) X # # Clean up object files. # noobj: X @$(RM) $(RMFLAGS) $(OBJ) X # # Clean up everything. # none: nofor noobj SHAR_EOF chmod 0600 dble/local.mak || echo 'restore of dble/local.mak failed' Wc_c="`wc -c < 'dble/local.mak'`" test 959 -eq "$Wc_c" || echo 'dble/local.mak: original size 959, current size' "$Wc_c" fi # ============= dble/skeleton.mak ============== if test -f 'dble/skeleton.mak' -a X"$1" != X"-c"; then echo 'x - skipping dble/skeleton.mak (File already exists)' else echo 'x - extracting dble/skeleton.mak (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/skeleton.mak' && #********************************************************************** # # Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal # All rights reserved. # # This code is part of a copyrighted package. For details, see the # file `cpyrit.doc' in the `lal' directory. # # ***************************************************************** # ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE # COPYRIGHT NOTICE # ***************************************************************** # #********************************************************************** # # Skeleton for all the makefiles. # X # # Make our own suffixes' list. # .SUFFIXES: .SUFFIXES: .f .o .SUFFIXES: .F .f X # # Default command. # .DEFAULT: X @$(ECHO) "Unknown target $@, try: make help" X # # Command to build .o files from .f files. # .f.o: X @$(ECHO) Making $@ from $< X @$(FC) -c $(FFLAGS) $< X .F.f: X @$(ECHO) Making $@ from $< X @$(CPP) -P $(CPPFLAGS) $< $@ X # # Various compilation programs and flags. # You need to make sure these are correct for your system. # CD = cd X CHMOD = chmod CHFLAGS = -f X COMPRESS = compress X CP = cp X # To find the path for cpp, try `man cpp', and it should list the path # at the top, under `Syntax'. It is usually in /lib. CPP = /lib/cpp CPPFLAGS = X ECHO = echo X # You also need to adjust the path to your FORTRAN compiler. FC = /usr/lang/f77 FFLAGS = X LDFLAGS = X LN = ln LNFLAGS = -s X MAKE = /bin/make X MKDIR = mkdir MDFLAGS = -p X MV = mv MVFLAGS = -f X RM = rm RMFLAGS = -f X SHELL = /bin/sh X TAR = tar X # # Default target is help. # help: genhelp lochelp X # # Dependencies for include files in all directories. # cpyrit.inc: ../inc/cpyrit.inc X @$(LN) $(LNFLAGS) $? $@ X dimblk.inc: ../inc/dimblk.inc X @$(LN) $(LNFLAGS) $? $@ X precon.inc: ../inc/precon.inc X @$(LN) $(LNFLAGS) $? $@ X # # This is the general help target. In the individual makefiles, one must # put a `lochelp' target, which may list additional help. The `lochelp' # target must exist even if it is merely a dummy. # genhelp: X @$(ECHO) "usage: make fmt" X @$(ECHO) " where fmt is one of:" X @$(ECHO) " for - make all auxiliary FORTRAN files" X @$(ECHO) " obj - make all object files" X @$(ECHO) " nofor - remove auxiliary FORTRAN files" X @$(ECHO) " noobj - remove object files" X @$(ECHO) " none - remove everything" SHAR_EOF chmod 0600 dble/skeleton.mak || echo 'restore of dble/skeleton.mak failed' Wc_c="`wc -c < 'dble/skeleton.mak'`" test 2412 -eq "$Wc_c" || echo 'dble/skeleton.mak: original size 2412, current size' "$Wc_c" fi # ============= dble/x.dat ============== if test -f 'dble/x.dat' -a X"$1" != X"-c"; then echo 'x - skipping dble/x.dat (File already exists)' else echo 'x - extracting dble/x.dat (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/x.dat' && X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 X 0.00000000000000e+00 SHAR_EOF chmod 0600 dble/x.dat || echo 'restore of dble/x.dat failed' Wc_c="`wc -c < 'dble/x.dat'`" test 5400 -eq "$Wc_c" || echo 'dble/x.dat: original size 5400, current size' "$Wc_c" fi # ============= dble/version.h ============== if test -f 'dble/version.h' -a X"$1" != X"-c"; then echo 'x - skipping dble/version.h (File already exists)' else echo 'x - extracting dble/version.h (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/version.h' && #define VERSION 2.0 SHAR_EOF chmod 0600 dble/version.h || echo 'restore of dble/version.h failed' Wc_c="`wc -c < 'dble/version.h'`" test 20 -eq "$Wc_c" || echo 'dble/version.h: original size 20, current size' "$Wc_c" fi # ============= dble/bru/axb.inc ============== if test ! -d 'dble/bru'; then echo 'x - creating directory dble/bru' mkdir 'dble/bru' fi if test -f 'dble/bru/axb.inc' -a X"$1" != X"-c"; then echo 'x - skipping dble/bru/axb.inc (File already exists)' else echo 'x - extracting dble/bru/axb.inc (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/bru/axb.inc' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C Common block ABLK. C X INTEGER INIT X DOUBLE PRECISION C1, C2, C3, D1, D2, D3 X COMMON /ABLK/C1, C2, C3, D1, D2, D3, INIT SHAR_EOF chmod 0600 dble/bru/axb.inc || echo 'restore of dble/bru/axb.inc failed' Wc_c="`wc -c < 'dble/bru/axb.inc'`" test 770 -eq "$Wc_c" || echo 'dble/bru/axb.inc: original size 770, current size' "$Wc_c" fi # ============= dble/bru/cpyrit.doc ============== if test -f 'dble/bru/cpyrit.doc' -a X"$1" != X"-c"; then echo 'x - skipping dble/bru/cpyrit.doc (File already exists)' else echo 'x - extracting dble/bru/cpyrit.doc (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/bru/cpyrit.doc' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is provided "as is", without any warranty of any kind, C either expressed or implied, including but not limited to, any C implied warranty of merchantibility or fitness for any purpose. C In no event will any party who distributed the code be liable for C damages or for any claim(s) by any other party, including but not C limited to, any lost profits, lost monies, lost data or data C rendered inaccurate, losses sustained by third parties, or any C other special, incidental or consequential damages arising out of C the use or inability to use the program, even if the possibility C of such damages has been advised against. The entire risk as to C the quality, the performance, and the fitness of the program for C any particular purpose lies with the party using the code. C C No derivative of this code may be used in a commercial package C without the prior explicit written permission of all authors or C their legal proxies. Verbatim copies of this code may be made and C distributed in any medium, provided that this copyright notice C is not removed or altered in any way. No fees may be charged for C distribution of the codes, other than a fee to cover the cost of C the media and a reasonable handling fee. C C********************************************************************** SHAR_EOF chmod 0600 dble/bru/cpyrit.doc || echo 'restore of dble/bru/cpyrit.doc failed' Wc_c="`wc -c < 'dble/bru/cpyrit.doc'`" test 1579 -eq "$Wc_c" || echo 'dble/bru/cpyrit.doc: original size 1579, current size' "$Wc_c" fi # ============= dble/bru/daxb.f ============== if test -f 'dble/bru/daxb.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/bru/daxb.f (File already exists)' else echo 'x - extracting dble/bru/daxb.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/bru/daxb.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C This file contains support routines for the Brusselator problem C used to test the eigenvalue codes. The following routines are in C this file: C C SUBROUTINE AXB (X,B) C Computes B = A * X. C SUBROUTINE ATXB (X,B) C Computes B = A^T * X. C SUBROUTINE EXACT (WR,WI) C Computes the exact eigenvalues of the Brusselator matrix. C SUBROUTINE GETMAT C Initializes the data for the Brusselator matrix. C DOUBLE PRECISION FUNCTION GETNRM() C Returns an estimate for the norm of the matrix based on using C Gershgorin disks. C C********************************************************************** C X SUBROUTINE AXB (X,B) C C Purpose: C This subroutine computes B = A * X for the Brusselator matrix. C For more info on this, see Rascham et al., "Waves in distributed C chemical systems: experiments and computations", published in C `New approaches to nonlinear problems in Dynamics', Ed P. Holmes, C SIAM Pub., 1980, pp. 271-288. See also B.N. Parlett and Y. Saad, C _Complex shift and invert strategies for real matrices_, Linear C Algebra Appl. 88/89 (1987), pp. 575-595. C C Parameters: C X = the vector to be multiplied by A (input). C B = the result of the multiplication (output). C C Noel M. Nachtigal C October 2, 1990 C X DOUBLE PRECISION B(*), X(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'axb.inc' C C Local variables. C X INTEGER I, J, NHALF X DOUBLE PRECISION CTMP, DTMP C C Set up the multiplication. C X CTMP = C2 - 2.0 * C1 X DTMP = D3 - 2.0 * D1 X NHALF = NROW / 2 C C Now do the multiplication. First, the main body. C X DO 10 I = 2, NHALF - 1 X J = I + NHALF X B(I) = CTMP * X(I) + C1 * ( X(I+1) + X(I-1) ) + X $ C3 * X(J) X B(J) = DTMP * X(J) + D1 * ( X(J+1) + X(J-1) ) + X $ D2 * X(I) X 10 CONTINUE C C Then, the exceptional cases. C X I = 1 X J = NHALF + 1 X B(I) = CTMP * X(I) + C1 * X(I+1) + C3 * X(J) X B(J) = DTMP * X(J) + D1 * X(J+1) + D2 * X(I) C X I = NHALF X J = NROW X B(I) = CTMP * X(I) + C1 * X(I-1) + C3 * X(J) X B(J) = DTMP * X(J) + D1 * X(J-1) + D2 * X(I) C X RETURN X END C C********************************************************************** C X SUBROUTINE ATXB (X,B) C C Purpose: C This subroutine computes B = A^T * X for the Brusselator matrix. C C Parameters: C X = the vector to be multiplied by A (input). C B = the result of the multiplication (output). C C Noel M. Nachtigal C October 2, 1990 C X DOUBLE PRECISION B(*), X(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'axb.inc' C C Local variables. C X INTEGER I, J, NHALF X DOUBLE PRECISION CTMP, DTMP C C Set up the multiplication. C X CTMP = C2 - 2.0 * C1 X DTMP = D3 - 2.0 * D1 X NHALF = NROW / 2 C C Now do the multiplication. First, the main body. C X DO 10 I = 2, NHALF - 1 X J = I + NHALF X B(I) = CTMP * X(I) + C1 * ( X(I+1) + X(I-1) ) + X $ D2 * X(J) X B(J) = DTMP * X(J) + D1 * ( X(J+1) + X(J-1) ) + X $ C3 * X(I) X 10 CONTINUE C C Then, the exceptional cases. C X I = 1 X J = NHALF + 1 X B(I) = CTMP * X(I) + C1 * X(I+1) + D2 * X(J) X B(J) = DTMP * X(J) + D1 * X(J+1) + C3 * X(I) C X I = NHALF X J = NROW X B(I) = CTMP * X(I) + C1 * X(I-1) + D2 * X(J) X B(J) = DTMP * X(J) + D1 * X(J-1) + C3 * X(I) C X RETURN X END C C********************************************************************** C X SUBROUTINE EXACT (WR,WI) C C Purpose: C This subroutine computes the exact eigenvalues of the Brusselator C matrix. C C Parameters: C WR = contains on output the real parts (output). C WI = contains on output the imaginary parts (output). C C Noel M. Nachtigal C October 2, 1990 C X INTRINSIC ABS, FLOAT, SIN, SQRT C X DOUBLE PRECISION WR(*), WI(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'axb.inc' C C Local variables. C X DOUBLE PRECISION PI X PARAMETER (PI=3.141592653589793238462643383279D0) C X INTEGER I, J X DOUBLE PRECISION XI, XMU, XN, S, S1, S2, T C C Initialize some data. C X J = 0 X XN = 0.5 / FLOAT( NROW / 2 + 1 ) C C Compute the eigenvalues -- they always come in pairs. C X DO 10 I = 1, NROW / 2 X XI = PI * FLOAT(I) X XMU = -( 2.0 * SIN(XI * XN) )**2.0 X S1 = C1*XMU + C2 X S2 = D1*XMU + D3 X S = 0.5 * (S1 + S2) X T = D2*C3 + 0.25 * (S1 - S2)**2 X IF (T.GE.0.0) THEN C C This is the real case. C X T = SQRT(ABS(T)) X J = J+1 X WR(J) = S+T X WI(J) = 0.0 X J = J+1 X WR(J) = S-T X WI(J) = 0.0 X ELSE C C This is the complex case. C X T = SQRT(ABS(T)) X J = J+1 X WR(J) = S X WI(J) = T X J = J+1 X WR(J) = S X WI(J) = -T X ENDIF X 10 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE GETMAT C C Purpose: C This subroutine initializes the Brusselator matrix by reading its C control parameters from the user. C C External routines used: C subroutine exact(wr,wi) C Computes the exact eigenvalues of the Brusselator matrix. C subroutine hsort(n,x,y) C Sorts arrays x and y based on x. C C Noel M. Nachtigal C October 2, 1990 C X INTRINSIC FLOAT C X INCLUDE 'dimblk.inc' X INCLUDE 'axb.inc' C C Local variables. C X CHARACTER*1 ANS X INTEGER I X DOUBLE PRECISION A, B, DX, DY, H, TMP, WR(NDIM), WI(NDIM), XL C C Initialize the variables to defaults. C X NCOL = 200 X NROW = 200 X DX = 0.00800 X DY = 0.00400 X A = 2.00000 X B = 5.45000 X XL = 0.51302 C C Get the data from the user. C X WRITE (6,'(A31,$)') 'Enter the dimension N (200) : ' X READ (5,*) NROW X NCOL = 2 * ( NROW / 2 ) X IF (NCOL.NE.NROW) THEN X WRITE (6,'(A32,I10)') 'N must be even - truncated to: ', NCOL X ENDIF C WRITE (6,'(A31,$)') 'Enter the spacing DX (0.008) : ' C READ (5,*) DX C WRITE (6,'(A31,$)') 'Enter the spacing DY (0.004) : ' C READ (5,*) DY C WRITE (6,'(A31,$)') 'Enter the constant A (2.0) : ' C READ (5,*) A C WRITE (6,'(A31,$)') 'Enter the constant B (5.45) : ' C READ (5,*) B C WRITE (6,'(A31,$)') 'Enter XL (0.51302) : ' C READ (5,*) XL C C Compute the matrix entries. C X NROW = NCOL / 2 X H = 1.0 / (FLOAT(NROW + 1)) X TMP = ( H * XL )**2 X C1 = DX / TMP X D1 = DY / TMP X C2 = B - 1.0 X D2 = -B X C3 = A**2 X D3 = -A**2 X NROW = NCOL C C Output the matrix parameters. C X WRITE (6,'(A5,E25.18)') 'C1 : ', C1 X WRITE (6,'(A5,E25.18)') 'C2 : ', C2 X WRITE (6,'(A5,E25.18)') 'C3 : ', C3 X WRITE (6,'(A5,E25.18)') 'D1 : ', D1 X WRITE (6,'(A5,E25.18)') 'D2 : ', D2 X WRITE (6,'(A5,E25.18)') 'D3 : ', D3 C C Check that it isn't too big. C X IF (NCOL.GT.NDIM) THEN X WRITE (6,'(A36)') 'Matrix dimension exceeds allocation.' X STOP X ENDIF C C Optionally output the exact eigenvalues. C X WRITE (6,'(A37,$)') 'Output the exact eigenvalues (Y/N) ? ' X READ (5,'(A1)') ANS X IF ((ANS.EQ.'Y').OR.(ANS.EQ.'y')) THEN X CALL EXACT (WR,WI) X CALL HSORT (NROW,WR,WI) X OPEN (10,FILE='exact.out') X DO 10 I = 1, NROW X WRITE (10,'(2E30.18)') WR(I), WI(I) X 10 CONTINUE X CLOSE (10) X ENDIF C C Set the initialized flag. C X INIT = 1 C X RETURN X END C C********************************************************************** C X DOUBLE PRECISION FUNCTION GETNRM() C C Purpose: C This function returns an estimate for the norm of the matrix. The C current estimate is based on Gershgorin disks. C C Noel M. Nachtigal C October 2, 1990 C X INTRINSIC MAX, MIN C X INCLUDE 'dimblk.inc' X INCLUDE 'axb.inc' C C Local variables. C X DOUBLE PRECISION DTMP1, DTMP2, DTMP3 C C Get the maximum column norm. C X DTMP2 = 2.0 * ABS(C1) + ABS(C2 - 2.0*C1) + ABS(D2) X DTMP3 = 2.0 * ABS(D1) + ABS(D3 - 2.0*D1) + ABS(C3) X DTMP1 = MAX(DTMP2,DTMP3) C C Get the maximum row norm. C X DTMP2 = 2.0 * ABS(C1) + ABS(C2 - 2.0*C1) + ABS(C3) X DTMP3 = 2.0 * ABS(D1) + ABS(D3 - 2.0*D1) + ABS(D2) X DTMP2 = MAX(DTMP2,DTMP3) C C Get the smallest of the two. C X GETNRM = MIN(DTMP1,DTMP2) C X RETURN X END C C********************************************************************** SHAR_EOF chmod 0600 dble/bru/daxb.f || echo 'restore of dble/bru/daxb.f failed' Wc_c="`wc -c < 'dble/bru/daxb.f'`" test 9536 -eq "$Wc_c" || echo 'dble/bru/daxb.f: original size 9536, current size' "$Wc_c" fi # ============= dble/bru/daxbpr.f ============== if test -f 'dble/bru/daxbpr.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/bru/daxbpr.f (File already exists)' else echo 'x - extracting dble/bru/daxbpr.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/bru/daxbpr.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C These are dummy preconditioner routines. C C There are four routines: C M1I (X,Y) - empty; C M1T (X,Y) - empty; C M2I (X,Y) - empty; C M2T (X,Y) - empty; C PSETUP - empty; C C Noel M. Nachtigal C January 13, 1991 C C********************************************************************** C X SUBROUTINE M1I(X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE M1T(X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE M2I(X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE M2T(X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE PSETUP C X RETURN X END C C********************************************************************** C X BLOCK DATA C C Purpose: C This sets the preconditioner name to 'NONE'. C C Noel M. Nachtigal C January 13, 1991 C X INCLUDE 'precon.inc' C X DATA PNAME/'NONE'/ C X END C C********************************************************************** SHAR_EOF chmod 0600 dble/bru/daxbpr.f || echo 'restore of dble/bru/daxbpr.f failed' Wc_c="`wc -c < 'dble/bru/daxbpr.f'`" test 2013 -eq "$Wc_c" || echo 'dble/bru/daxbpr.f: original size 2013, current size' "$Wc_c" fi # ============= dble/bru/makefile ============== if test -f 'dble/bru/makefile' -a X"$1" != X"-c"; then echo 'x - skipping dble/bru/makefile (File already exists)' else echo 'x - extracting dble/bru/makefile (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/bru/makefile' && #********************************************************************** # # Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal # All rights reserved. # # This code is part of a copyrighted package. For details, see the # file `cpyrit.doc' in the current directory. # # ***************************************************************** # ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE # COPYRIGHT NOTICE # ***************************************************************** # #********************************************************************** # # Makefile for the Brusselator matrix subdirectory. # # Files in this directory: # INC = dimblk.inc precon.inc FOR = OBJ = daxb.o daxbpr.o SRC = daxb.f daxbpr.f axb.inc X # # Include here the skeleton makefile. # include ../skeleton.mak include ../local.mak X # # This is the local help target. # lochelp: X # # Dependencies for files in this directory. # daxb.o: daxb.f axb.inc dimblk.inc X daxbpr.o: daxbpr.f axb.inc dimblk.inc precon.inc SHAR_EOF chmod 0600 dble/bru/makefile || echo 'restore of dble/bru/makefile failed' Wc_c="`wc -c < 'dble/bru/makefile'`" test 1077 -eq "$Wc_c" || echo 'dble/bru/makefile: original size 1077, current size' "$Wc_c" fi # ============= dble/csr/cpyrit.doc ============== if test ! -d 'dble/csr'; then echo 'x - creating directory dble/csr' mkdir 'dble/csr' fi if test -f 'dble/csr/cpyrit.doc' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/cpyrit.doc (File already exists)' else echo 'x - extracting dble/csr/cpyrit.doc (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/cpyrit.doc' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is provided "as is", without any warranty of any kind, C either expressed or implied, including but not limited to, any C implied warranty of merchantibility or fitness for any purpose. C In no event will any party who distributed the code be liable for C damages or for any claim(s) by any other party, including but not C limited to, any lost profits, lost monies, lost data or data C rendered inaccurate, losses sustained by third parties, or any C other special, incidental or consequential damages arising out of C the use or inability to use the program, even if the possibility C of such damages has been advised against. The entire risk as to C the quality, the performance, and the fitness of the program for C any particular purpose lies with the party using the code. C C No derivative of this code may be used in a commercial package C without the prior explicit written permission of all authors or C their legal proxies. Verbatim copies of this code may be made and C distributed in any medium, provided that this copyright notice C is not removed or altered in any way. No fees may be charged for C distribution of the codes, other than a fee to cover the cost of C the media and a reasonable handling fee. C C********************************************************************** SHAR_EOF chmod 0600 dble/csr/cpyrit.doc || echo 'restore of dble/csr/cpyrit.doc failed' Wc_c="`wc -c < 'dble/csr/cpyrit.doc'`" test 1579 -eq "$Wc_c" || echo 'dble/csr/cpyrit.doc: original size 1579, current size' "$Wc_c" fi # ============= dble/csr/csr.inc ============== if test -f 'dble/csr/csr.inc' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/csr.inc (File already exists)' else echo 'x - extracting dble/csr/csr.inc (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/csr.inc' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C Common block ABLK. C X INTEGER NZMAX X PARAMETER (NZMAX=1500) C X INTEGER IA(NDIM+1), ILU(NDIM+1), JA(NZMAX), JLU(NZMAX) X INTEGER IDA(NDIM), ITMP(NDIM) X DOUBLE PRECISION A(NZMAX), LU(NZMAX), AINV(NDIM) X DOUBLE PRECISION DN(NDIM), DR(NDIM), XTMP(NDIM) X COMMON /ABLK/A, LU, AINV, DN, DR, XTMP, IA, JA, ILU, JLU, IDA, X $ ITMP SHAR_EOF chmod 0600 dble/csr/csr.inc || echo 'restore of dble/csr/csr.inc failed' Wc_c="`wc -c < 'dble/csr/csr.inc'`" test 999 -eq "$Wc_c" || echo 'dble/csr/csr.inc: original size 999, current size' "$Wc_c" fi # ============= dble/csr/dcilut.f ============== if test -f 'dble/csr/dcilut.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/dcilut.f (File already exists)' else echo 'x - extracting dble/csr/dcilut.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/dcilut.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C These are the preconditioner routines for the ILUT preconditioner C applied to sparse matrices. C C There are five routines: C M1I (X,Y) - computes X = M_1^{-1} * X, Y a work vector; C M1T (X,Y) - computes X = M_1^{-T} * X, Y a work vector; C M2I (X,Y) - computes X = M_2^{-1} * X, Y a work vector; C M2T (X,Y) - computes X = M_2^{-T} * X, Y a work vector; C PSETUP - sets up the ILUT preconditioner. C C For the ILUT preconditioner, C M = M_1 * M_2 = L * U. C where C M_1 = L, M_2 = U. C C Noel M. Nachtigal C November 23, 1990 C C********************************************************************** C X SUBROUTINE M1I (X,Y) C C Purpose: C Computes X = M_1^{-1} * X. In the case of ILUT, M_1 is L, a lower C triangular matrix with 1's on the diagonal, stored in row major C order. Hence, the elimination is formulated as dot products of C rows of L with the upper parts of X. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C October 31, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the rows. C X DO 20 I = 1, NROW X DTMP = 0.0 C C Loop over the columns, up to the diagonal. C X DO 10 K = ILU(I), IDA(I)-1 X J = JLU(K) X DTMP = DTMP + LU(K) * X(J) X 10 CONTINUE C C Compute X(I). L(I,I) = 1.0. C X X(I) = X(I) - DTMP C X 20 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M1T (X,Y) C C Purpose: C Computes X = M_1^{-T} * X. In the case of ILUT, M_1 is L, a lower C triangular matrix with 1's on the diagonal, stored in row major C order. Hence, the elimination is formulated as sums of columns of C L^T (i.e., rows of L), subtracted from the corresponding parts of C X. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C October 31, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the columns of L^T. C X DO 20 I = NROW, 1, -1 C C We already have X(I), since L(I,I) = 1. Subtract the next column. C X DTMP = X(I) X DO 10 K = ILU(I), IDA(I)-1 X J = JLU(K) X X(J) = X(J) - DTMP * LU(K) X 10 CONTINUE X 20 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M2I (X,Y) C C Purpose: C Computes X = M_2^{-1} * X. In the case of ILUT, M_2 is U, a upper C triangular matrix, stored in row major order. The elimination is C formulated as dot products of rows of U with the lower parts of C X. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C October 31, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the rows. C X DO 20 I = NROW, 1, -1 X DTMP = 0.0 X DO 10 K = IDA(I)+1, ILU(I+1)-1 X J = JLU(K) X DTMP = DTMP + LU(K) * X(J) X 10 CONTINUE C C Compute X(I). LU(IDA(I)) = 1.0 / U(I,I). C X X(I) = ( X(I) - DTMP ) * LU(IDA(I)) C X 20 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M2T (X,Y) C C Purpose: C Computes X = M_2^{-T} * X. In the case of ILUT, M_2 is U, a upper C triangular matrix, stored in row major order. The elimination is C formulated as sums of solumns of U^T (i.e., rows of U) subtracted C from the corresponding parts of X. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C October 31, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the column of U^T. C X DO 20 I = 1, NROW C C Compute X(I). LU(IDA(I)) = 1.0 / U(I,I). C X X(I) = X(I) * LU(IDA(I)) C C Subtract the next column. C X DTMP = X(I) X DO 10 K = IDA(I)+1, ILU(I+1)-1 X J = JLU(K) X X(J) = X(J) - DTMP * LU(K) X 10 CONTINUE X 20 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE PSETUP C C Purpose: C This subroutine sets up the preconditioner. It assumes that the C global arrays A, IA, JA give the matrix in CSR format, and LU, C ILU, and JLU give the matrix in CSC format (transposed). We also C check for small diagonals, where small is relative to the 1-norm C of the row and the column. If a small diagonal element is found, C PRECON is set to -1. C C External routines used: C double precision dadd(dx,dy) C Computes dx + dy. Used to get around optimizers. C double precision dasum(n,dx,incx) C Computes the 1-norm of dx. C C Noel M. Nachtigal C October 31, 1990 C X INTRINSIC MAX X EXTERNAL DADD, DASUM X DOUBLE PRECISION DADD, DASUM C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER FILL, I, J, K X DOUBLE PRECISION DNRM, DTMP, TOL C C Compute the 1-norms of the rows in DN, of the columns in DTMP. We C save the pointers to the diagonal elements in IDA, and check for C small diagonal elements. C X DO 20 I = 1, NROW X IDA(I) = 0 X K = IA(I) X J = IA(I+1) - K X DN(I) = DASUM(J,A(K),1) X K = ILU(I) X J = ILU(I+1) - K X DTMP = DASUM(J,LU(K),1) X DNRM = MAX(DTMP,DN(I)) X DO 10 K = IA(I), IA(I+1)-1 X IF (JA(K).NE.I) GO TO 10 X IDA(I) = K X DTMP = DADD(DNRM,A(K)) X IF (DTMP.NE.DNRM) GO TO 20 X WRITE (6,'(A30,I10,E25.18)') 'Small diagonal on row:', X $ I, A(K) X PRECON = -1 X GO TO 20 X 10 CONTINUE X 20 CONTINUE X IF (PRECON.NE.1) RETURN C C Get the fill-in and tolerance from the user. C X FILL = (NZMAX - IA(NROW+1) + 1) / NROW X FILL = MIN(FILL,NROW) X FILL = FILL / 2 X WRITE (6,'(A17,I5)') 'Maximal fill-in: ', FILL X WRITE (6,'(A17,$)') 'Enter fill-in : ' X READ (5,'(I10)') I X FILL = MIN(FILL,I) X WRITE (6,'(A17,$)') 'Enter tolerance: ' X READ (5,*) TOL C C Compute the ILUT. C X WRITE (6,'(A16,I6)') 'Elements in A : ', IA(NROW+1)-1 X CALL ILUT(NROW,A,JA,IA,IDA,LU,JLU,ILU,TOL,FILL,ITMP, X $ DN,DR,XTMP) X WRITE (6,'(A16,I6)') 'Elements in LU: ', ILU(NROW+1)-1 C C Make a file with it? C X WRITE (6,'(A32,$)') 'Dump LU to file (1=Yes, 0=No) ? ' X READ (5,'(I10)') I X IF (I.EQ.1) THEN X DO 50 I = 1, NROW X WRITE (15,'(3I5,E30.18)') I, I, I, 1.0/LU(I) X DO 40 K = ILU(I), ILU(I+1)-1 X J = JLU(K) X WRITE (15,'(3I5,E30.18)') I, J, K, LU(K) X 40 CONTINUE X 50 CONTINUE X END IF C X RETURN X END C C********************************************************************** C X BLOCK DATA C C Purpose: C This sets the preconditioner name to 'ILUT'. C C Noel M. Nachtigal C October 23, 1990 C X INCLUDE 'precon.inc' C X DATA PNAME/'ILUT'/ C X END C C********************************************************************** SHAR_EOF chmod 0600 dble/csr/dcilut.f || echo 'restore of dble/csr/dcilut.f failed' Wc_c="`wc -c < 'dble/csr/dcilut.f'`" test 8640 -eq "$Wc_c" || echo 'dble/csr/dcilut.f: original size 8640, current size' "$Wc_c" fi # ============= dble/csr/dclilut.f ============== if test -f 'dble/csr/dclilut.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/dclilut.f (File already exists)' else echo 'x - extracting dble/csr/dclilut.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/dclilut.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C These are the preconditioner routines for the ILUT preconditioner C applied to sparse matrices on the left. C C There are four routines: C M1I (X,Y) - computes X = M_1^{-1} * X, Y a work vector; C M1T (X,Y) - computes X = M_1^{-T} * X, Y a work vector; C M2I (X,Y) - computes X = M_2^{-1} * X, Y a work vector; C M2T (X,Y) - computes X = M_2^{-T} * X, Y a work vector; C PSETUP - sets up the ILUT preconditioner. C C For the left ILUT preconditioner, C M = M_1 * M_2 = L * U. C where C M_1 = L * U, M_2 = I. C C Noel M. Nachtigal C November 23, 1990 C C********************************************************************** C X SUBROUTINE M1I (X,Y) C C Purpose: C Computes X = M_1^{-1} * X. In the case of left ILUT, M_1 is L * U C with L a lower triangular matrix with 1's on the diagonal, and U C an upper triangular matrix, both stored in row major order. Hence C the elimination is formulated as dot products. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C November 23, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the rows of L. C X DO 20 I = 1, NROW X DTMP = 0.0 C C Loop over the columns, up to the diagonal. C X DO 10 K = ILU(I), IDA(I)-1 X J = JLU(K) X DTMP = DTMP + LU(K) * X(J) X 10 CONTINUE C C Compute X(I). L(I,I) = 1.0. C X X(I) = X(I) - DTMP C X 20 CONTINUE C C Loop over the rows of U. C X DO 40 I = NROW, 1, -1 X DTMP = 0.0 X DO 30 K = IDA(I)+1, ILU(I+1)-1 X J = JLU(K) X DTMP = DTMP + LU(K) * X(J) X 30 CONTINUE C C Compute X(I). LU(IDA(I)) = 1.0 / U(I,I). C X X(I) = ( X(I) - DTMP ) * LU(IDA(I)) C X 40 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M1T (X,Y) C C Purpose: C Computes X = M_1^{-T} * X. In the case of left ILUT, M_1 is L * U C with L a lower triangular matrix with 1's on the diagonal, and U C an upper triangular matrix, both stored in row major order. Hence C the elimination is formulated as sums of columns subtracted from C the corresponding parts of X. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C November 23, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the column of U^T. C X DO 20 I = 1, NROW C C Compute X(I). LU(IDA(I)) = 1.0 / U(I,I). C X X(I) = X(I) * LU(IDA(I)) C C Subtract the next column. C X DTMP = X(I) X DO 10 K = IDA(I)+1, ILU(I+1)-1 X J = JLU(K) X X(J) = X(J) - DTMP * LU(K) X 10 CONTINUE X 20 CONTINUE C C Loop over the columns of L^T. C X DO 40 I = NROW, 1, -1 C C We already have X(I), since L(I,I) = 1. Subtract the next column. C X DTMP = X(I) X DO 30 K = ILU(I), IDA(I)-1 X J = JLU(K) X X(J) = X(J) - DTMP * LU(K) X 30 CONTINUE X 40 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M2I (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE M2T (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE PSETUP C C Purpose: C This subroutine sets up the preconditioner. It assumes that the C global arrays A, IA, JA give the matrix in CSR format, and LU, C ILU, and JLU give the matrix in CSC format (transposed). We also C check for small diagonals, where small is relative to the 1-norm C of the row and the column. If a small diagonal element is found, C PRECON is set to -1. C C External routines used: C double precision dadd(dx,dy) C Computes dx + dy. Used to get around optimizers. C double precision dasum(n,dx,incx) C Computes the 1-norm of dx. C C Noel M. Nachtigal C October 31, 1990 C X INTRINSIC MAX X EXTERNAL DADD, DASUM X DOUBLE PRECISION DADD, DASUM C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER FILL, I, J, K X DOUBLE PRECISION DNRM, DTMP, TOL C C Compute the 1-norms of the rows in DN, of the columns in DTMP. We C save the pointers to the diagonal elements in IDA, and check for C small diagonal elements. C X DO 20 I = 1, NROW X IDA(I) = 0 X K = IA(I) X J = IA(I+1) - K X DN(I) = DASUM(J,A(K),1) X K = ILU(I) X J = ILU(I+1) - K X DTMP = DASUM(J,LU(K),1) X DNRM = MAX(DTMP,DN(I)) X DO 10 K = IA(I), IA(I+1)-1 X IF (JA(K).NE.I) GO TO 10 X IDA(I) = K X DTMP = DADD(DNRM,A(K)) X IF (DTMP.NE.DNRM) GO TO 20 X WRITE (6,'(A30,I10,E25.18)') 'Small diagonal on row:', X $ I, A(K) X PRECON = -1 X GO TO 20 X 10 CONTINUE X 20 CONTINUE X IF (PRECON.NE.1) RETURN C C Get the fill-in and tolerance from the user. C X FILL = (NZMAX - IA(NROW+1) + 1) / NROW X FILL = MIN(FILL,NROW) X FILL = FILL / 2 X WRITE (6,'(A17,I5)') 'Maximal fill-in: ', FILL X WRITE (6,'(A17,$)') 'Enter fill-in : ' X READ (5,'(I10)') I X FILL = MIN(FILL,I) X WRITE (6,'(A17,$)') 'Enter tolerance: ' X READ (5,*) TOL C C Compute the ILUT. C X WRITE (6,'(A16,I6)') 'Elements in A : ', IA(NROW+1)-1 X CALL ILUT(NROW,A,JA,IA,IDA,LU,JLU,ILU,TOL,FILL,ITMP, X $ DN,DR,XTMP) X WRITE (6,'(A16,I6)') 'Elements in LU: ', ILU(NROW+1)-1 C C Make a file with it? C X WRITE (6,'(A32,$)') 'Dump LU to file (1=Yes, 0=No) ? ' X READ (5,'(I10)') I X IF (I.EQ.1) THEN X DO 50 I = 1, NROW X WRITE (15,'(3I5,E30.18)') I, I, I, 1.0/LU(I) X DO 40 K = ILU(I), ILU(I+1)-1 X J = JLU(K) X WRITE (15,'(3I5,E30.18)') I, J, K, LU(K) X 40 CONTINUE X 50 CONTINUE X END IF C X RETURN X END C C********************************************************************** C X BLOCK DATA C C Purpose: C This sets the preconditioner name to 'ILUT'. C C Noel M. Nachtigal C November 23, 1990 C X INCLUDE 'precon.inc' C X DATA PNAME/'SPARSE ILUT (LEFT)'/ C X END C C********************************************************************** SHAR_EOF chmod 0600 dble/csr/dclilut.f || echo 'restore of dble/csr/dclilut.f failed' Wc_c="`wc -c < 'dble/csr/dclilut.f'`" test 7610 -eq "$Wc_c" || echo 'dble/csr/dclilut.f: original size 7610, current size' "$Wc_c" fi # ============= dble/csr/dclssor.f ============== if test -f 'dble/csr/dclssor.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/dclssor.f (File already exists)' else echo 'x - extracting dble/csr/dclssor.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/dclssor.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C These are the preconditioner routines for the SSOR preconditioner C applied on the right, sparse matrices. C C There are five routines: C M1I (X,Y) - computes X = M_1^{-1} * X, Y a work vector; C M1T (X,Y) - computes X = M_1^{-T} * X, Y a work vector; C M2I (X,Y) - computes X = M_2^{-1} * X, Y a work vector; C M2T (X,Y) - computes X = M_2^{-T} * X, Y a work vector; C PSETUP - sets up the preconditioner (empty for SSOR). C C For the SSOR preconditioner, C M = M_1 * M_2 = ( D + w L ) D^{-1} ( D + w U ), C where C A = D + L + U. C C For left preconditioning, M_2 = I. C C Note: these routines assume that the indices in the arrays JA are C sorted in increasing order. It is the resposibility of the caller C to ensure this is true. C C Noel M. Nachtigal C October 23, 1990 C C********************************************************************** C X SUBROUTINE M1I (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C SSOR parameter OMEGA. C X DOUBLE PRECISION OMEGA X COMMON /OMG/OMEGA C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C If the SSOR parameter is not initialized, initialize it. C X IF (OMEGA.EQ.-1.0) THEN X WRITE (6,'(A30,$)') 'Enter SSOR parameter OMEGA : ' X READ (5,*) OMEGA X ENDIF C C Multiply by ( D + w L )^{-1}. C X DO 20 I = 1, NROW X DTMP = 0.0 X DO 10 K = IA(I), IDA(I)-1 X J = JA(K) X DTMP = DTMP + A(K) * X(J) X 10 CONTINUE X X(I) = ( X(I) - OMEGA * DTMP ) * AINV(I) X 20 CONTINUE C C Multiply by ( D + w U )^{-1} * D. C X DO 40 I = NROW, 1, -1 X DTMP = 0.0 X DO 30 K = IDA(I)+1, IA(I+1)-1 X J = JA(K) X DTMP = DTMP + A(K) * X(J) X 30 CONTINUE X X(I) = X(I) - OMEGA * DTMP * AINV(I) X 40 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M1T (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C SSOR parameter OMEGA. C X DOUBLE PRECISION OMEGA X COMMON /OMG/OMEGA C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C If the SSOR parameter is not initialized, initialize it. C X IF (OMEGA.EQ.-1.0) THEN X WRITE (6,'(A30,$)') 'Enter SSOR parameter OMEGA : ' X READ (5,*) OMEGA X ENDIF C C Multiply by D * ( D + w U )^{-T}. C X DO 20 I = 1, NROW X DTMP = -X(I) * OMEGA * AINV(I) X DO 10 K = IDA(I)+1, IA(I+1)-1 X J = JA(K) X X(J) = X(J) + DTMP * A(K) X 10 CONTINUE X 20 CONTINUE C C Multiply by ( D + w L )^{-T}. C X DO 40 I = NROW, 1, -1 X X(I) = X(I) * AINV(I) X DTMP = -OMEGA * X(I) X DO 30 K = IA(I), IDA(I)-1 X J = JA(K) X X(J) = X(J) + DTMP * A(K) X 30 CONTINUE X 40 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M2I (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE M2T (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE PSETUP C C Purpose: C This subroutine sets up the preconditioner. It assumes that the C global arrays A, IA, JA give the matrix in CSR format, and LU, C ILU, and JLU give the matrix in CSC format (transposed). We store C the diagonal elements inverted in AINV. We also check for small C diagonals, where small is relative to the 1-norm of the row and C the column. If a small diagonal element is found, PRECON is set C to -1. C C External routines used: C double precision dadd(dx,dy) C Computes dx + dy. Used to get around optimizers. C double precision dasum(n,dx,incx) C Computes the 1-norm of dx. C subroutine dzero(n,dx,incx) C Sets dx to zero. C C Noel M. Nachtigal C October 31, 1990 C X INTRINSIC MAX X EXTERNAL DADD, DASUM X DOUBLE PRECISION DADD, DASUM C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C X INTEGER I, J, K X DOUBLE PRECISION DA, DCOL, DROW, DTMP C C Compute the 1-norms of the rows in DROW, of the columns in DCOL. C We save the pointers to the diagonal elements in IDA, and check C for small diagonal elements. C X DO 20 I = 1, NROW X IDA(I) = 0 X AINV(I) = 0.0 X K = IA(I) X J = IA(I+1) - K X DROW = DASUM(J,A(K),1) X K = ILU(I) X J = ILU(I+1) - K X DCOL = DASUM(J,LU(K),1) X DROW = MAX(DCOL,DROW) X DO 10 K = IA(I), IA(I+1)-1 X IF (JA(K).NE.I) GO TO 10 X IDA(I) = K X DA = A(K) X IF (DA.NE.0.0) AINV(I) = 1.0 / DA X DTMP = DADD(DROW,A(K)) X IF (DTMP.NE.DROW) GO TO 20 X WRITE (6,'(A30,I10,E25.18)') 'Small diagonal on row:', X $ I, A(K) X PRECON = -1 X GO TO 20 X 10 CONTINUE X 20 CONTINUE X IF (PRECON.NE.1) RETURN C X RETURN X END C C********************************************************************** C X BLOCK DATA C C Purpose: C This sets the SSOR parameter OMEGA to -1, the preconditioner name C to 'SPARSE SSOR (RIGHT)'. C C Noel M. Nachtigal C October 23, 1990 C X INCLUDE 'precon.inc' C X DOUBLE PRECISION OMEGA X COMMON /OMG/OMEGA C X DATA OMEGA/-1.0/ X DATA PNAME/'SPARSE SSOR (LEFT)'/ C X END C C********************************************************************** SHAR_EOF chmod 0600 dble/csr/dclssor.f || echo 'restore of dble/csr/dclssor.f failed' Wc_c="`wc -c < 'dble/csr/dclssor.f'`" test 6513 -eq "$Wc_c" || echo 'dble/csr/dclssor.f: original size 6513, current size' "$Wc_c" fi # ============= dble/csr/dcrilut.f ============== if test -f 'dble/csr/dcrilut.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/dcrilut.f (File already exists)' else echo 'x - extracting dble/csr/dcrilut.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/dcrilut.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C These are the preconditioner routines for the ILUT preconditioner C applied to sparse matrices on the right. C C There are four routines: C M1I (X,Y) - computes X = M_1^{-1} * X, Y a work vector; C M1T (X,Y) - computes X = M_1^{-T} * X, Y a work vector; C M2I (X,Y) - computes X = M_2^{-1} * X, Y a work vector; C M2T (X,Y) - computes X = M_2^{-T} * X, Y a work vector; C PSETUP - sets up the ILUT preconditioner. C C For the right ILUT preconditioner, C M = M_1 * M_2 = L * U. C where C M_1 = I, M_2 = L * U. C C Noel M. Nachtigal C November 23, 1990 C C********************************************************************** C X SUBROUTINE M1I (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C********************************************************************** C X SUBROUTINE M1T (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C********************************************************************** C X SUBROUTINE M2I (X,Y) C C Purpose: C Computes X = M_2^{-1} * X. In the case of right ILUT M_2 is L * U C with L a lower triangular matrix with 1's on the diagonal, and U C an upper triangular matrix, both stored in row major order. Hence C the elimination is formulated as dot products. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C November 23, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the rows of L. C X DO 20 I = 1, NROW X DTMP = 0.0 C C Loop over the columns, up to the diagonal. C X DO 10 K = ILU(I), IDA(I)-1 X J = JLU(K) X DTMP = DTMP + LU(K) * X(J) X 10 CONTINUE C C Compute X(I). L(I,I) = 1.0. C X X(I) = X(I) - DTMP C X 20 CONTINUE C C Loop over the rows of U. C X DO 40 I = NROW, 1, -1 X DTMP = 0.0 X DO 30 K = IDA(I)+1, ILU(I+1)-1 X J = JLU(K) X DTMP = DTMP + LU(K) * X(J) X 30 CONTINUE C C Compute X(I). LU(IDA(I)) = 1.0 / U(I,I). C X X(I) = ( X(I) - DTMP ) * LU(IDA(I)) C X 40 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M2T (X,Y) C C Purpose: C Computes X = M_2^{-T} * X. In the case of right ILUT M_2 is L * U C with L a lower triangular matrix with 1's on the diagonal, and U C an upper triangular matrix, both stored in row major order. Hence C the elimination is formulated as sums of columns subtracted from C the corresponding parts of X. C C Parameters: C X = the vector to apply the preconditioner to (input/output). C Y = work vector (output). C C Noel M. Nachtigal C November 23, 1990 C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C Loop over the column of U^T. C X DO 20 I = 1, NROW C C Compute X(I). LU(IDA(I)) = 1.0 / U(I,I). C X X(I) = X(I) * LU(IDA(I)) C C Subtract the next column. C X DTMP = X(I) X DO 10 K = IDA(I)+1, ILU(I+1)-1 X J = JLU(K) X X(J) = X(J) - DTMP * LU(K) X 10 CONTINUE X 20 CONTINUE C C Loop over the columns of L^T. C X DO 40 I = NROW, 1, -1 C C We already have X(I), since L(I,I) = 1. Subtract the next column. C X DTMP = X(I) X DO 30 K = ILU(I), IDA(I)-1 X J = JLU(K) X X(J) = X(J) - DTMP * LU(K) X 30 CONTINUE X 40 CONTINUE C X RETURN X END C C**********************************************************************C C X SUBROUTINE PSETUP C C Purpose: C This subroutine sets up the preconditioner. It assumes that the C global arrays A, IA, JA give the matrix in CSR format, and LU, C ILU, and JLU give the matrix in CSC format (transposed). We also C check for small diagonals, where small is relative to the 1-norm C of the row and the column. If a small diagonal element is found, C PRECON is set to -1. C C External routines used: C double precision dadd(dx,dy) C Computes dx + dy. Used to get around optimizers. C double precision dasum(n,dx,incx) C Computes the 1-norm of dx. C C Noel M. Nachtigal C October 31, 1990 C X INTRINSIC MAX X EXTERNAL DADD, DASUM X DOUBLE PRECISION DADD, DASUM C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER FILL, I, J, K X DOUBLE PRECISION DNRM, DTMP, TOL C C Compute the 1-norms of the rows in DN, of the columns in DTMP. We C save the pointers to the diagonal elements in IDA, and check for C small diagonal elements. C X DO 20 I = 1, NROW X IDA(I) = 0 X K = IA(I) X J = IA(I+1) - K X DN(I) = DASUM(J,A(K),1) X K = ILU(I) X J = ILU(I+1) - K X DTMP = DASUM(J,LU(K),1) X DNRM = MAX(DTMP,DN(I)) X DO 10 K = IA(I), IA(I+1)-1 X IF (JA(K).NE.I) GO TO 10 X IDA(I) = K X DTMP = DADD(DNRM,A(K)) X IF (DTMP.NE.DNRM) GO TO 20 X WRITE (6,'(A30,I10,E25.18)') 'Small diagonal on row:', X $ I, A(K) X PRECON = -1 X GO TO 20 X 10 CONTINUE X 20 CONTINUE X IF (PRECON.NE.1) RETURN C C Get the fill-in and tolerance from the user. C X FILL = (NZMAX - IA(NROW+1) + 1) / NROW X FILL = MIN(FILL,NROW) X FILL = FILL / 2 X WRITE (6,'(A17,I5)') 'Maximal fill-in: ', FILL X WRITE (6,'(A17,$)') 'Enter fill-in : ' X READ (5,'(I10)') I X FILL = MIN(FILL,I) X WRITE (6,'(A17,$)') 'Enter tolerance: ' X READ (5,*) TOL C C Compute the ILUT. C X WRITE (6,'(A16,I6)') 'Elements in A : ', IA(NROW+1)-1 X CALL ILUT(NROW,A,JA,IA,IDA,LU,JLU,ILU,TOL,FILL,ITMP, X $ DN,DR,XTMP) X WRITE (6,'(A16,I6)') 'Elements in LU: ', ILU(NROW+1)-1 C C Make a file with it? C X WRITE (6,'(A32,$)') 'Dump LU to file (1=Yes, 0=No) ? ' X READ (5,'(I10)') I X IF (I.EQ.1) THEN X DO 50 I = 1, NROW X WRITE (15,'(3I5,E30.18)') I, I, I, 1.0/LU(I) X DO 40 K = ILU(I), ILU(I+1)-1 X J = JLU(K) X WRITE (15,'(3I5,E30.18)') I, J, K, LU(K) X 40 CONTINUE X 50 CONTINUE X END IF C X RETURN X END C C********************************************************************** C X BLOCK DATA C C Purpose: C This sets the preconditioner name to 'ILUT'. C C Noel M. Nachtigal C Movember 23, 1990 C X INCLUDE 'precon.inc' C X DATA PNAME/'SPARSE ILUT (RIGHT)'/ C X END C C********************************************************************** SHAR_EOF chmod 0600 dble/csr/dcrilut.f || echo 'restore of dble/csr/dcrilut.f failed' Wc_c="`wc -c < 'dble/csr/dcrilut.f'`" test 7610 -eq "$Wc_c" || echo 'dble/csr/dcrilut.f: original size 7610, current size' "$Wc_c" fi # ============= dble/csr/dcrssor.f ============== if test -f 'dble/csr/dcrssor.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/dcrssor.f (File already exists)' else echo 'x - extracting dble/csr/dcrssor.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/dcrssor.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C These are the preconditioner routines for the SSOR preconditioner C applied on the right, sparse matrices. C C There are five routines: C M1I (X,Y) - computes X = M_1^{-1} * X, Y a work vector; C M1T (X,Y) - computes X = M_1^{-T} * X, Y a work vector; C M2I (X,Y) - computes X = M_2^{-1} * X, Y a work vector; C M2T (X,Y) - computes X = M_2^{-T} * X, Y a work vector; C PSETUP - sets up the preconditioner (empty for SSOR). C C For the SSOR preconditioner, C M = M_1 * M_2 = ( D + w L ) D^{-1} ( D + w U ), C where C A = D + L + U. C C For right preconditioning, M_1 = I. C C Note: these routines assume that the indices in the arrays JA are C sorted in increasing order. It is the resposibility of the caller C to ensure this is true. C C Noel M. Nachtigal C October 23, 1990 C C********************************************************************** C X SUBROUTINE M1I (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE M1T (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X RETURN X END C C********************************************************************** C X SUBROUTINE M2I (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C SSOR parameter OMEGA. C X DOUBLE PRECISION OMEGA X COMMON /OMG/OMEGA C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C If the SSOR parameter is not initialized, initialize it. C X IF (OMEGA.EQ.-1.0) THEN X WRITE (6,'(A30,$)') 'Enter SSOR parameter OMEGA : ' X READ (5,*) OMEGA X ENDIF C C Multiply by ( D + w L )^{-1}. C X DO 20 I = 1, NROW X DTMP = 0.0 X DO 10 K = IA(I), IDA(I)-1 X J = JA(K) X DTMP = DTMP + A(K) * X(J) X 10 CONTINUE X X(I) = ( X(I) - OMEGA * DTMP ) * AINV(I) X 20 CONTINUE C C Multiply by ( D + w U )^{-1} * D. C X DO 40 I = NROW, 1, -1 X DTMP = 0.0 X DO 30 K = IDA(I)+1, IA(I+1)-1 X J = JA(K) X DTMP = DTMP + A(K) * X(J) X 30 CONTINUE X X(I) = X(I) - OMEGA * DTMP * AINV(I) X 40 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE M2T (X,Y) C X DOUBLE PRECISION X(*), Y(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C SSOR parameter OMEGA. C X DOUBLE PRECISION OMEGA X COMMON /OMG/OMEGA C C Local variables. C X INTEGER I, J, K X DOUBLE PRECISION DTMP C C If the SSOR parameter is not initialized, initialize it. C X IF (OMEGA.EQ.-1.0) THEN X WRITE (6,'(A30,$)') 'Enter SSOR parameter OMEGA : ' X READ (5,*) OMEGA X ENDIF C C Multiply by D * ( D + w U )^{-T}. C X DO 20 I = 1, NROW X DTMP = -X(I) * OMEGA * AINV(I) X DO 10 K = IDA(I)+1, IA(I+1)-1 X J = JA(K) X X(J) = X(J) + DTMP * A(K) X 10 CONTINUE X 20 CONTINUE C C Multiply by ( D + w L )^{-T}. C X DO 40 I = NROW, 1, -1 X X(I) = X(I) * AINV(I) X DTMP = -OMEGA * X(I) X DO 30 K = IA(I), IDA(I)-1 X J = JA(K) X X(J) = X(J) + DTMP * A(K) X 30 CONTINUE X 40 CONTINUE C X RETURN X END C C********************************************************************** C X SUBROUTINE PSETUP C C Purpose: C This subroutine sets up the preconditioner. It assumes that the C global arrays A, IA, JA give the matrix in CSR format, and LU, C ILU, and JLU give the matrix in CSC format (transposed). We store C the diagonal elements inverted in AINV. We also check for small C diagonals, where small is relative to the 1-norm of the row and C the column. If a small diagonal element is found, PRECON is set C to -1. C C External routines used: C double precision dadd(dx,dy) C Computes dx + dy. Used to get around optimizers. C double precision dasum(n,dx,incx) C Computes the 1-norm of dx. C subroutine dzero(n,dx,incx) C Sets dx to zero. C C Noel M. Nachtigal C October 31, 1990 C X INTRINSIC MAX X EXTERNAL DADD, DASUM X DOUBLE PRECISION DADD, DASUM C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C X INTEGER I, J, K X DOUBLE PRECISION DA, DCOL, DROW, DTMP C C Compute the 1-norms of the rows in DROW, of the columns in DCOL. C We save the pointers to the diagonal elements in IDA, and check C for small diagonal elements. C X DO 20 I = 1, NROW X IDA(I) = 0 X AINV(I) = 0.0 X K = IA(I) X J = IA(I+1) - K X DROW = DASUM(J,A(K),1) X K = ILU(I) X J = ILU(I+1) - K X DCOL = DASUM(J,LU(K),1) X DROW = MAX(DCOL,DROW) X DO 10 K = IA(I), IA(I+1)-1 X IF (JA(K).NE.I) GO TO 10 X IDA(I) = K X DA = A(K) X IF (DA.NE.0.0) AINV(I) = 1.0 / DA X DTMP = DADD(DROW,A(K)) X IF (DTMP.NE.DROW) GO TO 20 X WRITE (6,'(A30,I10,E25.18)') 'Small diagonal on row:', X $ I, A(K) X PRECON = -1 X GO TO 20 X 10 CONTINUE X 20 CONTINUE X IF (PRECON.NE.1) RETURN C X RETURN X END C C********************************************************************** C X BLOCK DATA C C Purpose: C This sets the SSOR parameter OMEGA to -1, the preconditioner name C to 'SPARSE SSOR (RIGHT)'. C C Noel M. Nachtigal C October 23, 1990 C X INCLUDE 'precon.inc' C X DOUBLE PRECISION OMEGA X COMMON /OMG/OMEGA C X DATA OMEGA/-1.0/ X DATA PNAME/'SPARSE SSOR (RIGHT)'/ C X END C C********************************************************************** SHAR_EOF chmod 0600 dble/csr/dcrssor.f || echo 'restore of dble/csr/dcrssor.f failed' Wc_c="`wc -c < 'dble/csr/dcrssor.f'`" test 6521 -eq "$Wc_c" || echo 'dble/csr/dcrssor.f: original size 6521, current size' "$Wc_c" fi # ============= dble/csr/dcsr.f ============== if test -f 'dble/csr/dcsr.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/dcsr.f (File already exists)' else echo 'x - extracting dble/csr/dcsr.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/dcsr.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C This file contains support routines for sparse matrices. Most of C the routines herein use SPARSKIT routines to handle the matrix. C SPARSKIT is Copyright 1990, Youcef Saad. The following routines C are in this file: C C SUBROUTINE AXB (X,B) C Computes B = A * X. A is preconditioned. C SUBROUTINE ATXB (X,B) C Computes B = A^T * X. A is preconditioned. C SUBROUTINE GETMAT C Reads in a sparse matrix from a Boeing-Harwell format file. C DOUBLE PRECISION FUNCTION GETNRM() C Returns an estimate for the norm of the matrix based on using C Gershgorin disks. C SUBROUTINE SORTJA (N,AROW,JROW) C Sorts the indices in JROW in increasing order. C C********************************************************************** C X SUBROUTINE AXB (X,B) C C Purpose: C This subroutine computes B = A * X, for a sparse matrix A. The C code assumes that A is preconditioned to M_1^{-1} A M_2^{-1}. C C Parameters: C X = the vector to be multiplied by A (input). C B = the result of the multiplication (output). C C External routines used: C subroutine amux(nrow,x,b,a,ja,ia) C SPARSKIT routine to compute sparse matrix-vector product. C subroutine dcopy(n,dx,incx,dy,incy) C Computes dy = dx. C subroutine m1i(x,y) C Computes x = M_1^{-1} * x. C subroutine m2i(x,y) C Computes x = M_2^{-1} * x. C C Noel M. Nachtigal C September 27, 1990 C X DOUBLE PRECISION B(*), X(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Copy X to the local vector. C X CALL DCOPY (NROW,X,1,XTMP,1) C C Multiply by preconditioner matrix M_2^{-1}. C X IF (PRECON.EQ.1) CALL M2I (XTMP,B) C C Multiply by A. C X CALL AMUX (NROW,XTMP,B,A,JA,IA) C C Multiply by the preconditioner matrix M_1^{-1}. C X IF (PRECON.EQ.1) CALL M1I (B,XTMP) C X RETURN X END C C********************************************************************** C X SUBROUTINE ATXB (X,B) C C Purpose: C This subroutine computes B = A^T * X, for sparse matrix A. The C code assumes that A is preconditioned to M_1^{-1} A M_2^{-1}. C C Parameters: C X = the vector to be multiplied by A^T (input). C B = the result of the multiplication (output). C C External routines used: C subroutine atmux(nrow,x,b,a,ja,ia) C SPARSKIT routine to compute sparse matrix-vector product with C the transposed matrix. C subroutine dcopy(n,dx,incx,dy,incy) C Computes dy = dx. C subroutine m1t(x,y) C Computes x = M_1^{-T} * x. C subroutine m2t(x,y) C Computes x = M_2^{-T} * x. C C Noel M. Nachtigal C September 27, 1990 C X DOUBLE PRECISION B(*),X(*) C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Copy X to the local vector. C X CALL DCOPY (NROW,X,1,XTMP,1) C C Multiply by preconditioner matrix M_1^{-T}. C X IF (PRECON.EQ.1) CALL M1T (XTMP,B) C C Multiply by A. C X CALL ATMUX (NROW,XTMP,B,A,JA,IA) C C Multiply by preconditioner matrix M_2^{-T} C X IF (PRECON.EQ.1) CALL M2T (B,XTMP) C X RETURN X END C C********************************************************************** C X SUBROUTINE GETMAT C C Purpose: C This subroutine initializes the matrix by reading it from a user- C specified Boeing-Harwell format data file. Optionally, it will C also convert the matrix to a dense matrix and output the latter C to the ASCII file matrix.dat. Note that the Boeing-Harwell stuff C is stored in compressed sparse column format, while SPARSKIT uses C compressed sparse row format. C C External routines used: C subroutine readmt(nmax,nzmax,job,iounit,a,ja,ia,rhs,nrhs,guesol, C nrow,ncol,nnz,title,key,type,ierr) C SPARSKIT routine to read in a sparse Harwell-Boeing matrix. C subroutine csrcsc(n,job,ipos,a,ja,ia,ao,jao,iao) C Convert a CSR matrix to CSC format or vice-versa. C C Noel M. Nachtigal C September 27, 1990 C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X CHARACTER ANS*1, GUESOL*2, KEY*8, TITLE*72, TYPE*3 X INTEGER I, IERR, IOUNIT, J, J1, J2, JOB, K, NRHS, NNZ X DOUBLE PRECISION RHS C C Get the data file from the user. C X WRITE (6,'(A29,$)') 'Enter sparse data file name: ' X READ (5,'(A72)') TITLE C C Open the file. C X OPEN (10,FILE=TITLE) C C Call the SPARSKIT routine to read in the matrix. C X JOB = 2 X NRHS = 0 X IOUNIT = 10 X CALL READMT (NDIM,NZMAX,JOB,IOUNIT,LU,JLU,ILU,RHS,NRHS,GUESOL, X $ NROW,NCOL,NNZ,TITLE,KEY,TYPE,IERR) X CLOSE (10) C C Output the matrix parameters. C X WRITE (6,'(A7,I10)') 'NDIM :',NDIM X WRITE (6,'(A7,I10)') 'NZMAX :',NZMAX X WRITE (6,'(A7,I10)') 'JOB :',JOB X WRITE (6,'(A7,I10)') 'NRHS :',NRHS X WRITE (6,'(A7,A10)') 'GUESOL:',GUESOL X WRITE (6,'(A7,I10)') 'NROW :',NROW X WRITE (6,'(A7,I10)') 'NCOL :',NCOL X WRITE (6,'(A7,I10)') 'NNZ :',NNZ X WRITE (6,'(A7,A73)') 'TITLE :',TITLE X WRITE (6,'(A7,A11)') 'KEY :',KEY X WRITE (6,'(A7,A10)') 'TYPE :',TYPE X WRITE (6,'(A7,I10)') 'IERR :',IERR C C Check for errors. C X IF (IERR.NE.0) STOP C C Convert the matrix to CSR format (transpose it). C X CALL CSRCSC (NROW,1,1,LU,JLU,ILU,A,JA,IA) C C Order the indices in ascending order. C X DO 10 I = 1, NROW X K = IA(I) X J = IA(I+1) - K X CALL SORTJA (J,A(K),JA(K)) X K = ILU(I) X J = ILU(I+1) - K X CALL SORTJA (J,LU(K),JLU(K)) X 10 CONTINUE C C Does the user want a dense matrix file? C C WRITE (6,'(A38,$)') 'Produce dense matrix data file (Y/N)? ' C READ (5,'(A1)') ANS X ANS = 'N' C C Make the dense matrix file. The matrix is output by columns. C X IF ((ANS.EQ.'Y').OR.(ANS.EQ.'y')) THEN X OPEN (12,FILE='matrix.dat') X WRITE (12,'(I20)') NROW X DO 50 I = 1, NROW X J1 = 1 X DO 30 K = ILU(I), ILU(I+1)-1 X J2 = JLU(K) X DO 20 J = J1, J2-1 X WRITE (12,'(E25.17)') 0.0 X 20 CONTINUE X WRITE (12,'(E25.17)') LU(K) X J1 = J2 + 1 X 30 CONTINUE X DO 40 J = J1, NCOL X WRITE (12,'(E25.17)') 0.0 X 40 CONTINUE X 50 CONTINUE X CLOSE (12) X END IF C X RETURN X END C C********************************************************************** C X DOUBLE PRECISION FUNCTION GETNRM () C C Purpose: C This function returns an estimate for the norm of the matrix. The C current estimate is based on Gershgorin disks. C C External routines used: C subroutine cnrms(nrow,nrm,a,ja,ia,diag) C SPARSKIT routine to compute column norms. C integer function idamax(n,x,incx) C Computes the index of the largest entry in magnitude. C subroutine rnrms(nrow,nrm,a,ja,ia,diag) C SPARSKIT routine to compute row norms. C C Noel M. Nachtigal C September 27, 1990 C X INTEGER IDAMAX X EXTERNAL IDAMAX X INTRINSIC MIN C X INCLUDE 'dimblk.inc' X INCLUDE 'csr.inc' C C Local variables. C X INTEGER IDX, NRM X DOUBLE PRECISION DIAG(NDIM-1), DTMP1, DTMP2 C X NRM = 1 X CALL CNRMS (NROW,NRM,A,JA,IA,DIAG) X IDX = IDAMAX(NROW,DIAG,1) X DTMP1 = DIAG(IDX) X CALL RNRMS (NROW,NRM,A,JA,IA,DIAG) X IDX = IDAMAX(NROW,DIAG,1) X DTMP2 = DIAG(IDX) X GETNRM = MIN(DTMP1,DTMP2) C X RETURN X END C C********************************************************************** C X SUBROUTINE SORTJA (N,AROW,JROW) C C Purpose: C This is a SPARSKIT subroutine to sort the column indices of a row C of a matrix stored in compressed sparse row (CSR) format in order C (increasing). The routine is given the array of elements of the C row in AROW, with the corresponding column indices in JROW. Both C arrays are of length N; the elements are given so that they can C be kept in correspondence with the elements in AROW. C The routine uses Heapsort to carry out the sorting; the code is C copied verbatim from Numerical Recipes (minus the bugs). C C Parameters: C N = the length of the row (input). C AROW = the row to be sorted (input/output). C JROW = the array of matching column indices (input/output). C C Noel M. Nachtigal C October 28, 1990 C X INTRINSIC ABS C X INTEGER N, JROW(*) X DOUBLE PRECISION AROW(*) C C Local variables. C X INTEGER I, J, JTMP, K, L X DOUBLE PRECISION DTMP C X IF (N.LE.1) RETURN C X L = N / 2 + 1 X K = N X 10 IF (L.GT.1) THEN X L = L - 1 X DTMP = AROW(L) X JTMP = JROW(L) X ELSE X DTMP = AROW(K) X JTMP = JROW(K) X AROW(K) = AROW(1) X JROW(K) = JROW(1) X K = K - 1 X IF (K.LE.1) THEN X AROW(1) = DTMP X JROW(1) = JTMP X RETURN X END IF X END IF X I = L X J = L + L X 20 IF (J.LE.K) THEN X IF (J.LT.K) THEN X IF (JROW(J).LT.JROW(J+1)) J = J + 1 X END IF X IF (JTMP.LT.JROW(J)) THEN X AROW(I) = AROW(J) X JROW(I) = JROW(J) X I = J X J = J + J X ELSE X J = K + 1 X END IF X GO TO 20 X END IF X AROW(I) = DTMP X JROW(I) = JTMP X GO TO 10 C X END C C********************************************************************** X SHAR_EOF chmod 0600 dble/csr/dcsr.f || echo 'restore of dble/csr/dcsr.f failed' Wc_c="`wc -c < 'dble/csr/dcsr.f'`" test 10293 -eq "$Wc_c" || echo 'dble/csr/dcsr.f: original size 10293, current size' "$Wc_c" fi # ============= dble/csr/ilut.f ============== if test -f 'dble/csr/ilut.f' -a X"$1" != X"-c"; then echo 'x - skipping dble/csr/ilut.f (File already exists)' else echo 'x - extracting dble/csr/ilut.f (Text)' sed 's/^X//' << 'SHAR_EOF' > 'dble/csr/ilut.f' && C********************************************************************** C C Copyright (C) 1991 Roland W. Freund and Noel M. Nachtigal C All rights reserved. C C This code is part of a copyrighted package. For details, see the C file `cpyrit.doc' in the current directory. C C ***************************************************************** C ANY USE OF THIS CODE CONSTITUES ACCEPTANCE OF THE TERMS OF THE C COPYRIGHT NOTICE C ***************************************************************** C C********************************************************************** C C These are the setup routines for the ILUT preconditioners applied C to sparse matrices. C C There are two routines: C ILUT - carries out the ILUT decomposition; C SORT - used by ILUT. C C Noel M. Nachtigal C November 23, 1990 C C********************************************************************** C X SUBROUTINE ILUT (NROW,A,JA,IA,IDA,LU,JLU,ILU,TOL,FILL,ITMP,DN,DR, X $ DS) C C Purpose: C This subroutine computes the Incomplete LU decomposition of A C with a variant of the dual thresholding strategy. C C Parameters: C NROW = the dimension of A (input). C A = the matrix A in compressed sparse row