LGSOLV - a set of FORTRAN-functions for Large Linear System solving. ======================================= Authors: P.G.Akishin, A.P.Sapoznikov. LCTA,JINR, Dubna,Russia. E-mail: akishin@main1.jinr.dubna.su sapoznikov@main1.jinr.dubna.su If you have to solve Large Linear Systems, then a program described below will be useful for you. If you haven't - you simply will get information about a typical Tool for such problems. The main problem for programmers working witn large matrices is how to keep them in computer's memory. Being written in a disk-storage, matrix involves a lot of disk exchanges, because all Matrix-Algorithms use it "row-by-row" and "column-by-column" at the same time. Here we're proposing an original modification of known Gauss method for pivot element elimination that uses internal rows-transposition-vector and works with matrix strictly "col-by-col". Our algorithm factorizes the source matrix so that created factor-matrix may be used several times for quick solution of Linear System with any Right Parts. Your Main Program ought to do only two things : 1. to prepare functions calculating matrix elements and system right parts; 2. to reserve 2 buffers in memory, the larger they are, the less will be a number of exchanges. So, let us start. The system A*X = C is being solved, where A is a square matrix [ NDIM : NDIM ], X,C - are NDIM-vectors. In practice it's often necessary to solve the whole family of NC systems: A * X(k) = C(k) k=1,2,...,NC with the same matrix A. Historycally FORTRAN was being used for similar problems, that's way our algorithm is, of course, FORTRAN-written. Here's our Operation Set for Large Linear System solving and results getting: INTEGER FUNCTION LGMINP( lbuf, mpf, ndim, elmatr ) INTEGER FUNCTION LGMINV( lbuf, mpf, ndim, elmatr ) INTEGER FUNCTION LGSOLV( nc, elvect ) REAL*8 FUNCTION G1SOL ( row, nrp ) REAL*8 FUNCTION G1RP ( row, nrp ) REAL*8 FUNCTION G1EMAT( row, col ) REAL*8 FUNCTION G1FMAT( row, col ) INTEGER lbuf,mpf,ndim,nc,row,col,nrp REAL*8 elmatr,elvect ndim - System Dimension nc - number of System Right Parts lbuf - size of 2 working buffers declared in main program as common/buf1/b1(lbuf) common/buf2/b2(lbuf) real*8 b1,b2 You are giving to SOLVER as many space, as you can, but no less than ndim. mpf - matrix packing factor (4 or 8). All the calculations use Double Precision but only mpf high bytes of each word will be written on disk. Default=8. Usage mpf=4 requires less working disk-storage but decreases the result's accuracy. elmatr - user function giving a Source Matrix element elvect - user function giving an element of Right Part LGMINP inputs a Source Matrix A invoking elmatr(i,j) ndim**2 times in a following sequence: (1,1),(2,1),...,(ndim,1),(1,2),(2,2),...,(ndim,2)... i.e. "col-by-col", and writes it to Scratch-file using logical number 1. Than LGMINP prepares two triangle matrices, L and R, getting A = ~L * R, and writes them to Scratch-file using logical number 2. Both files are created in the current user's working directory automatically and are deleted when user program ends. Buffers b1,b2 are used for the disk-exchanges. LGSOLV inputs all the System Right Parts C invoking elvect(i,k) ndim*nc times in a following sequence: (1,1),(2,1),...,(ndim,1), ..., (ndim-1,nc),(ndim,nc) and adds them to File-2. Than LGSOLV computes the System Solutions X = ~R * ( L * C ) and adds them to File-1. G1SOL - these 4 functions are simply extractors of information G1RP that SOLVER's files hold. G1SOL gives row-th element G1EMAT of nrp-th solution. G1RP gives row-th element of nrp-th G1FMAT Right Part. row,col=1,2,...,ndim; nrp=1,2,...,nc; G1EMAT gives (row,col)-th element of Source Matrix. G1FMAT - the same for factor-matrix (we don't imagine whom it will be useful, it was made "for symmetry"). LGMINV realises a special problem : Matrix Inversion. Having inputing the matrix, it just starts to solve system A * X = E and writes all ndim solutions to File-1 instead of Source Matrix. In this case G1EMAT gives elements of Inverse Matrix but G1SOL,G1RP - are nonsensible. WARNING: as we hold information "col-by-col", be careful during extraction : if you extract "row-by-row" then it takes a disk-exchange for almost every extracted element ! Return Codes All the operations don't print any message or stopped. REAL*8 entries simply return 0.0 if their parameters are wrong. INTEGER entries return codes are : 0 - all right 1 - Source Matrix is singular 2 - user buffer's size lbuf < ndim An Example In your application program there ought be the following : C C ..... C PARAMETER(lbuf=...) COMMON /buf1/ b1(lbuf) COMMON /buf2/ b2(lbuf) EXTERNAL fun1,fun2 REAL*8 b1,b2,fun1,fun2,G1SOL,G1EMAT,G1RP C ..... C 1. Input a System Dimension and Matrix into LGSOLV. C here ndim = , C mpf = C ner = LGMINP( lbuf, mpf, ndim, fun1 ) C C ..... C 2. System Solving for NC given Right Parts : C ner = LGSOLV( nc, fun2 ) C C ..... C 3. Results getting : i,j=1,2,...,ndim; k=1,2,...,nc C (i-th element of k-th solution, k-th Right Part C or j-th Source Matrix column ) : x = G1SOL( i, k ) c = G1RP ( i, k ) a = G1EMAT( i, j ) C ..... END REAL*8 FUNCTION fun1( i, j ) INTEGER*2 i,j C This function gives A(i,j) - element of System Matrix fun1 = ... RETURN END REAL*8 FUNCTION fun2( i, k ) INTEGER*2 i,k C This function gives C(i,k) - i-th element of C k-th System Right Part Vector fun2 = ... RETURN END The Demo Test Our Demo Test is prepared for IBM PC. When you run DEMOTEST.EXE file, you'll be asked about current lbuf,mpf,ndim,nc values and matrix type. Test solves the system A*X=C choosing C corresponing to known solution X(i)=i. During the solving process you can see the dynamic of source and factor matrices usage. After process' finish elapsed time, exchanges's counter and errors will be shown. Errors are the difference between computed and known solutions. Using the 286-processor we didn't want to compile in LARGE mode. Because of that you can't use lbuf > 8000 In the other hand, it seems that 8000 will be quite enough for all cases. You can to experiment with lbuf working with Demo Test. If you are going to do it, don't forget to experiment with Matrix Packing Factor also. You'll have seen, how the errors depend on MPF. You imagine, of course, that the sensible NDIM for test is about 50-100. NDIM=400 requires about 250 sec. for 386 processor. LGSOLV Packing List 1. DEMOTEST.EXE - run file for IBM PC 2. DEMOTEST.FOR - Source FORTRAN text of DEMOTEST Program 3. LGSOLV.OBJ - compiled by MS-FORTRAN 5.00 for 286 processor 4. README.DOC - this document If you want to join LGSOLV.OBJ to your own main program, add the following empty subroutine: subroutine vis0(m1,m2,job) integer*2 m1,m2 character*(*) title entry vist(title) entry vis(m1,m2,job) return end because LGSOLV Demo Version invokes these 3 subroutines for solving process' visualisation.