*-------------------------------------------------------------- * Date created: December 4, 1995 * Date last modified: December 8, 1995 *-------------------------------------------------------------- * * Main program for solving large-scale linear algebraic systems * * A*X=B or AT*X=B * ======= ======== * * taking into account the sparsity of the A matrix and using the * Product Form of Inverse representation of the A**(-1). (AT is * the transpouse of the matrix A) * * This is a package of subroutines dedicated to compute the * Product Form of Inverse (PFI) with a preassigned pivot procedure * of the A**(-1), to solve corresponding large-scale systems of * linear equations, A*X = B or AT*X = B, exploiting sparsity * in all cases. * * The calling sequence is as follows: * * MAIN ---- RF02A * RF02C ---- MRM03A * MFTRAN --- MRM02A * MBTRAN --- MRM02A * * where: * * RF02A determine a permutation of rows and columns of the * matrix to the "bump and spike" structure. * * RF02C compute the Product Form of Inverse of the matrix * permuted to this form "bump and spike" * * MRM03A compress the arrays a and indl in order to create * more room for PFI generation. * * MFTRAN is for forward transformation of the RHS term of the * system to compute the solution of A*X = B. * X = A**(-1)*B = QP*(Tn*......(T1*B)). * * MBTRAN is for backward transformation of the RHS term of the * system to compute the solution of AT*X = B. * X = A**(-T)*B = B*(((QP)*Tn)*......T1). * where Ti (i=1,...,N) is the i-th eta-vector of the * inverse A**(-1). * * MRM02A is for insitu permutation of a vector according to a * permutation vector. *-------------------------------------------------------------- * Neculai Andrei *-------------------------------------------------------------- parameter (iain=500000, iainn=100000) real*8 a(iain) real*8 b(iainn) real*8 tolpiv, tolzer integer indl(iain),indc(iain) integer ipel(iainn),ipec(iainn) integer ilp(iainn),icp(iainn) integer nzl(iainn),nzc(iainn) integer ics(iainn),ikm(iainn),mc(iainn) * common /tol/ tolpiv,tolzer * * Tolerances: tolinv is for pivot rejection. * tolzer is for flush out small non-zeros. * tolpiv=1.0D0-10 tolzer=1.0D0-12 * no=2 infinv=1 mesaj=0 * Input file open(unit=1,file='sys.dat',status='old') * Output file open(unit=2,file='sys.lis',status='unknown') * write(2,1) 1 format(//5x,'Solving linear algebraic systems:') write(2,2) 2 format(5x,' A*X = B or AT*X = B') write(2,3) 3 format(5x,'---------------------------------'/) * *-------------------------------------------------------------- * Read the data associated to the system: * * n = Number of unknowns (rows or columns) of the system. * nz = Number of non-zero elements in the A matrix. * * a = Array containing the non-zero elements of the matrix * given in natural column order. * indl Array containing the indices of rows of the above * non-zero coefficients of the matrix from a array. * ipec Array containing the address in a and indl where the * columns are stored. ipec(k) is the address in a and * indl array where starts the column k. (k=1,2,...,n) * b Array containing the elements of the RHS term. Vector * b contains all the elements (zero and non-zero) in natural * order. * * ia = Space allocated in the memory to compute the inverse * of the matrix. If ia is to small than it will be * necessary to compress the a and indl arrays by means * of the MRM03A subroutine. (Please see RF02C routine.) *------------------------------------------------------------- * read(1,10)n,nz,ia 10 format(3i4) read(1,20)(a(i),i=1,nz) 20 format(10f8.2) read(1,31)(indl(i),i=1,nz) read(1,31)(ipec(i),i=1,n+1) 31 format(10i2) read(1,20)(b(i),i=1,n) * *------------------------------------------------------------- * write(2,4) 4 format(4x,'Data of the system:') write(2,40) n,nz,ia 40 format(/5x,'n =',i4,2x,' nz =',i4,2x,' ia =',i4) write(2,11) 11 format(/5x,'A: Non-zero elements in column order') write(2,32)(a(i),i=1,nz) 32 format(5f8.2) write(2,12) 12 format(/5x,'INDL: Row indices of the above non-zeros') write(2,30)(indl(i),i=1,nz) 30 format(20i3) write(2,27) 27 format(/5x,'IPEC: The starting address of columns ', 1 'in indl array') write(2,30)(ipec(i),i=1,n+1) write(2,14) 14 format(/5x,'B: Right-Hand-Side term of the system') write(2,32)(b(i),i=1,n) * *-------------------------------------------------------------- * On input in RF02A we must supply only the structure of the * matrix: n, nz, indl, ipec. * The matrix must be introduced in columns order. * call rf02a(n,nz,indl,ipec,indc,ipel,ilp,icp,nzl,nzc, 1 ics,ikm,mc,ns,mesaj,no,infinv) * write(2,400)n,nz 400 format(5x,'After RF02A ',' n= ',i3,3x,' nz= ',i3) write(2,301) 301 format(5x,'Area a:') write(2,32)(a(i),i=1,nz) * write(2,342) 342 format(5x,'Area indl:') write(2,343)(indl(i),i=1,nz) 343 format(15i4) * write(2,356) 356 format(5x,'Area indc:') write(2,343)(indc(i),i=1,nz) * write(2,344) 344 format(5x,'Area ipec:') write(2,343)(ipec(i),i=1,n+1) * write(2,345) 345 format(5x,'Area ipel:') write(2,343)(ipel(i),i=1,n+1) * write(2,346)(ilp(i),i=1,n) 346 format(5x,'ilp=',15i3) write(2,347)(icp(i),i=1,n) 347 format(5x,'icp=',15i3) write(2,348)(nzl(i),i=1,n) 348 format(5x,'nzl=',15i3) write(2,349)(nzc(i),i=1,n) 349 format(5x,'nzc=',15i3) write(2,315)(ics(i),i=1,n) 315 format(5x,'ics=',15i3) write(2,378)(ikm(i),i=1,n) 378 format(5x,'ikm=',15i3) write(2,316)(mc(i),i=1,n) 316 format(5x,'mc=',15i3) write(2,317)ns 317 format(5x,'ns=',i4) * if(mesaj.lt.0) stop * *-------------------------------------------------------------- * On input RF02C considers: n, nz, ia, a, indl, ipec, ilp, * icp, ics, nzc, and ns, as they are given by RF02A. * call rf02c(n,nz,ia,a,indl,ipec,ipel,ilp,icp,nzl,nzc, 1 ics,mc,nreta,ns,mesaj,no,infinv) * * write(2,300)n,nz 300 format(5x,'After RF02C ',' n= ',i3,3x,' nz= ',i3) write(2,301) write(2,32)(a(i),i=1,nz) * write(2,302) 302 format(5x,'Area indl:') write(2,303)(indl(i),i=1,nz) 303 format(15i4) * write(2,304) 304 format(5x,'Area ipec:') write(2,303)(ipec(i),i=1,n+1) * write(2,305) 305 format(5x,'Area ipel:') write(2,303)(ipel(i),i=1,n+1) * write(2,306)(ilp(i),i=1,n) 306 format(5x,'ilp=',15i3) write(2,307)(icp(i),i=1,n) 307 format(5x,'icp=',15i3) write(2,308)(nzl(i),i=1,n) 308 format(5x,'nzl=',15i3) write(2,309)(nzc(i),i=1,n) 309 format(5x,'nzc=',15i3) write(2,310)(ics(i),i=1,n) 310 format(5x,'ics=',15i3) write(2,311)(mc(i),i=1,n) 311 format(5x,'mc=',15i3) write(2,312)nreta,ns 312 format(5x,'nreta=',i4,4x,'ns=',i4) * if(mesaj.lt.0) stop *-------------------------------------------------------------- * MFTRAN computes: A**(-1)*B * call mftran(a,indl,ipec,mc,n,nz,nreta,icp,b) * *-------------------------------------------------------------- * write(2,55) 55 format(/5x,'Solution vector: A**(-1)*B') write(2,56)(b(i),i=1,n) 56 format(5x,e14.7) * *============================================================== * * A new RHS term is considered * read(1,20)(b(i),i=1,n) write(2,60) 60 format(/5x,'Solve the system AT*X=B (AT: transpose of A)') write(2,61) 61 format(/5x,'Right-Hand-Side term of the system AT*X=B') write(2,32)(b(i),i=1,n) * * MBTRAN computes: A**(-T)*B * call mbtran(a,indl,ipec,mc,n,nz,nreta,ilp,b) * *-------------------------------------------------------------- * write(2,57) 57 format(/5x,'Solution vector: A**(-T)*B') write(2,56)(b(i),i=1,n) * stop end *-------------------------------------------------------------- * RF02A == RF02A * ############################################################# C------------------------------------------------------------- C Date created: December 14, 1988. C Date last modified: December 4, 1995 C------------------------------------------------------------- C Purpose: C -------- C This subroutine belongs to a package of subroutines dedicated C to compute the Product Form of Inverse (PFI) with preassigned C pivot procedure of a matrix, to solve corresponding large- C scale systems of linear equations, exploiting the sparsity in C all cases. C It implements the preassigned pivot procedure P4 of Hellerman C and Rarick, generating a pivot sequence which minimize the C number of non-zero elements in eta-vectors, and the work done C in forming these vectors. C RF02A finds a permutation of rwos and columns of the matrix to C "bump and spike" form. These permutations are stored in arrays C ilp and icp as it is explained below. C C Argument List: C -------------- C call RF02A(n,nz,indl,ipec,indc,ipel,ilp,icp,nzl,nzc,ics,ikm, C mc,ns,mesaj,no,infinv) C Inputs: C ------- C n The order of the matrix. It is not altered by rf02a C nz The number of non-zero elements from the matrix. C It is not altered by rf02a. C indl An integer array of dimension nz containing the row C indices of the elements of the matrix given in column C order. It is not altered by rf02a. C ipec An integer array of dimension n+1 containing the address C where the columns start in array indl. ipec(k) contains C the starting address of the k-th column stored in array C indl. Here we need only the information on the structure C of the matrix. Important: The matrix must be introduced C in column order. It is not altered by rf02a. C no User supplied number of the output device on which C output is to be written. C infinv Parameter used to write the structure synopsis of the C matrix in "bump and spike" form. If infinv=0 no such C information is supplied. C C Outputs: C -------- C ilp Array of dimension n with the row pivot sequence. C icp Array of dimension n with the column pivot sequence. C These two arrays contain the sequence of pivots. C nzc Array of dimension n with the number of non-zeros in C each column of the matrix. C ics Array of dimension n with the indices of the spikes C columns. C ikm Array of dimension n with the number of spikes in each C bump. C mc Array of dimension n with the dimension of each bump. C ns Total number of spikes columns from the matrix. C mesaj Integer used to transmit the error diagnostics. C mesaj=-1 if there are zero rows. C mesaj=-2 if the matrix is structurally singular. C C Working arrays: C --------------- C indc Array of dimension nz. C ipel Array of dimension n+1 C nzl Array of dimension n. C-------------------------------------------------------------- * Neculai Andrei *-------------------------------------------------------------- C SUBROUTINE RF02A(N,NZ,INDL,IPEC,INDC,IPEL,ILP,ICP, 1 NZL,NZC,ICS,IKM,MC,NS,MESAJ,NO,INFINV) C C parameter (iain=500000, iainn=100000) INTEGER INDL(iain),INDC(iain) INTEGER IPEL(iainn),IPEC(iainn) INTEGER ILP(iainn), ICP(iainn) INTEGER NZL(iainn), NZC(iainn) INTEGER ICS(iainn), IKM(iainn), MC(iainn) CHARACTER*3 LC(2) C DATA LC/'LIN','COL'/ C C C---- ******************** BODY OF THE PROGRAM (RF02A). C---- ************************************************* C C1. ------------------------- Initialization. C MIU=N NIU=1 IS=0 DO I=1,N NZL(I)=0 ILP(I)=0 ICP(I)=0 END DO C C2. ------------------------- Structure matrix determination. C Compute: NZL, NZC, IPEL, INDC DO 1 I=1,N 1 NZC(I)=IPEC(I+1)-IPEC(I) DO 2 K=1,NZ I=INDL(K) 2 NZL(I)=NZL(I)+1 C K=1 IPEL(1)=K DO 3 J=1,N I=NZL(J) IF(I.EQ.0) GO TO 1000 K=K+I 3 IPEL(J+1)=K C C3.-------------------------- Sort non-zeros into rows order. C DO 5 I=1,N KP=IPEC(I) KU=IPEC(I+1)-1 DO 4 K=KP,KU J=INDL(K) L=IPEL(J) INDC(L)=I 4 IPEL(J)=IPEL(J)+1 5 CONTINUE K=1 IPEL(1)=K DO 6 J=1,N K=K+NZL(J) 6 IPEL(J+1)=K C C4. ------------------------- Backward Triangularization. C Determination of the singleton C columns. 7 DO 11 I=1,N IF(NZC(I).NE.1) GO TO 11 ICP(MIU)=I NZC(I)=0 KP=IPEC(I) KU=IPEC(I+1)-1 DO 8 K=KP,KU J=INDL(K) IF(NZL(J).GT.0) GO TO 9 8 CONTINUE 9 ILP(MIU)=J NZL(J)=-1 KP=IPEL(J) KU=IPEL(J+1)-1 DO 10 K=KP,KU J=INDC(K) II=2 C C5. ------------------------- Test to detect singularities. C IF(NZC(J).EQ.1) GO TO 1002 10 NZC(J)=NZC(J)-1 MIU=MIU-1 IF(MIU.EQ.0) GO TO 17 GO TO 7 11 CONTINUE C C6. ------------------------- Forward Triangularization. C Determination of the singleton C rows. 12 DO 16 I=1,N IF(NZL(I).NE.1) GO TO 16 ILP(NIU)=I NZL(I)=0 KP=IPEL(I) KU=IPEL(I+1)-1 DO 13 K=KP,KU J=INDC(K) IF(NZC(J).GT.0) GO TO 14 13 CONTINUE 14 ICP(NIU)=J NZC(J)=-1 KP=IPEC(J) KU=IPEC(J+1)-1 DO 15 K=KP,KU J=INDL(K) II=1 C C7. ------------------------- Test to detect singularities. C IF(NZL(J).EQ.1) GO TO 1002 15 NZL(J)=NZL(J)-1 NIU=NIU+1 IF(NIU.GT.MIU) GO TO 17 IF(IS.EQ.0) GO TO 12 GO TO 40 16 CONTINUE C C8. ------------------------- Test for Completion or start of C a bump. C 17 IF(NIU.GT.MIU) GO TO 68 C C9. ------------------------- Set KK to minimum of non-zero row C count (NZL). 18 KK=N DO 19 I=1,N 19 IF(NZL(I).GT.0.AND.NZL(I).LT.KK) KK=NZL(I) C C10. ------------------------ Determine the set of columns C available for pivoting. C DO 20 I=1,N IKM(I)=0 20 MC(I)=I II=NIU-1 DO 21 I=1,II J=IABS(ICP(I)) 21 MC(J)=0 II=N-MIU IF(II.GT.0) THEN DO 22 I=1,II K=N+1-I J=IABS(ICP(K)) 22 MC(J)=0 END IF IF(IS.EQ.0) GO TO 24 DO 23 I=1,IS J=ICS(I) 23 MC(J)=0 24 J=0 DO 25 I=1,N IF(MC(I).EQ.0) GO TO 25 J=J+1 IKM(J)=MC(I) 25 CONTINUE DO 26 I=1,J MC(I)=IKM(I) 26 IKM(I)=0 C C11. ------------------------ Compute the tally function. C 27 DO 29 I=1,J II=MC(I) KP=IPEC(II) KU=IPEC(II+1)-1 DO 28 K=KP,KU L=INDL(K) 28 IF(NZL(L).GT.0.AND.NZL(L).LE.KK) IKM(I)=IKM(I)+1 29 CONTINUE C C12. ------------------------ Spike decisions. C MAXI=0 DO 30 I=1,J 30 IF(IKM(I).GT.MAXI) MAXI=IKM(I) II=0 DO 31 I=1,J IF(IKM(I).NE.MAXI) GO TO 31 II=II+1 NS=MC(I) 31 CONTINUE IF(MAXI.EQ.1) GO TO 33 IF(II.EQ.1) GO TO 38 MAX=0 DO 32 I=1,J IF(IKM(I).NE.MAXI) GO TO 32 K=MC(I) IF(NZC(K).LE.MAX) GO TO 32 MAX=NZC(K) NS=K 32 CONTINUE GO TO 38 33 DO 34 I=1,J IF(IKM(I).EQ.1) GO TO 34 MC(I)=0 34 CONTINUE K=0 DO 35 I=1,J IF(MC(I).EQ.0) GO TO 35 K=K+1 IKM(K)=MC(I) 35 CONTINUE J=K DO 36 I=1,J MC(I)=IKM(I) 36 IKM(I)=0 K1=KK KK=N+1 DO 37 I=1,N 37 IF(NZL(I).GT.K1.AND.NZL(I).LE.KK) KK=NZL(I) IF(KK.EQ.N+1) KK=K1 GO TO 27 C C13. ------------------------ Spike array entry. C 38 IS=IS+1 ICS(IS)=NS NZC(NS)=0 KP=IPEC(NS) KU=IPEC(NS+1)-1 DO 39 K=KP,KU J=INDL(K) 39 NZL(J)=NZL(J)-1 C C14. ------------------------ Path decisions. Another spike C or forward triangularization. C 40 MIN=N DO 41 I=1,N IF(NZL(I).GT.0.AND.NZL(I).LE.MIN) MIN=NZL(I) 41 CONTINUE IF(MIN.GT.1) GO TO 18 II=0 DO 42 I=1,N IF(NZL(I).NE.MIN) GO TO 42 II=II+1 42 CONTINUE IF(II.EQ.1.OR.IS.EQ.0) GO TO 12 KK=1 DO 43 I=1,N IKM(I)=0 43 MC(I)=I II=NIU-1 DO 44 I=1,II J=IABS(ICP(I)) 44 MC(J)=0 II=N-MIU IF(II.GT.0) THEN DO 45 I=1,II K=N+1-I J=IABS(ICP(K)) 45 MC(J)=0 END IF IF(IS.EQ.0) GO TO 47 DO 46 I=1,IS J=ICS(I) 46 MC(J)=0 47 J=0 DO 48 I=1,N IF(MC(I).EQ.0) GO TO 48 J=J+1 IKM(J)=MC(I) 48 CONTINUE DO 49 I=1,J MC(I)=IKM(I) 49 IKM(I)=0 C C15. ------------------------ Compute the tally function. C 50 DO 52 I=1,J II=MC(I) KP=IPEC(II) KU=IPEC(II+1)-1 DO 51 K=KP,KU L=INDL(K) 51 IF(NZL(L).GT.0.AND.NZL(L).LE.KK) IKM(I)=IKM(I)+1 52 CONTINUE C C16. ------------------------ Row and column pivot selection. C MAXI=0 DO 53 I=1,J 53 IF(IKM(I).GT.MAXI) MAXI=IKM(I) II=0 DO 54 I=1,J IF(IKM(I).NE.MAXI) GO TO 54 II=II+1 NP=MC(I) IR=I 54 CONTINUE IF(MAXI.EQ.1) GO TO 56 IF(II.EQ.1) GO TO 61 MAX=0 DO 55 I=1,J IF(IKM(I).NE.MAXI) GO TO 55 K=MC(I) IF(NZC(K).LE.MAX) GO TO 55 MAX=NZC(K) NP=K IR=I 55 CONTINUE GO TO 61 56 DO 57 I=1,J IF(IKM(I).EQ.1) GO TO 57 MC(I)=0 57 CONTINUE K=0 DO 58 I=1,J IF(MC(I).EQ.0) GO TO 58 K=K+1 IKM(K)=MC(I) 58 CONTINUE J=K DO 59 I=1,J MC(I)=IKM(I) 59 IKM(I)=0 K1=KK KK=N+1 DO 60 I=1,N 60 IF(NZL(I).GT.K1.AND.NZL(I).LE.KK) KK=NZL(I) IF(KK.EQ.N+1) KK=K1 GO TO 50 61 IQ=IKM(IR) C C17. ------------------------ Enter row and column pivot into C the pivot sequence list. C ICP(NIU)=NP NZC(NP)=0 KP=IPEC(NP) KU=IPEC(NP+1)-1 DO 62 K=KP,KU J=INDL(K) IF(NZL(J).EQ.1) GO TO 63 62 CONTINUE 63 ILP(NIU)=J NZL(J)=0 DO 64 K=KP,KU J=INDL(K) 64 NZL(J)=NZL(J)-1 65 NIU=NIU+1 IF(NIU.GT.MIU) GO TO 68 IQ=IQ-1 IF(IQ.EQ.0) GO TO 40 C C18. ------------------------ Enter spike into the pivot C sequence list C II=0 DO 66 I=1,N IF(NZL(I).NE.0) GO TO 66 II=I GO TO 67 66 CONTINUE 67 IF(II.EQ.0) GO TO 40 ILP(NIU)=II ICP(NIU)=-ICS(IS) IS=IS-1 IF(IS.EQ.0) ILP(NIU)=-II NZL(II)=-1 GO TO 65 68 CONTINUE C C19. ---------------- Determination of: C Number of bumps - KK. C Total number of spikes - NS. C Number of spikes in each bump - IKM array. C Dimension of each bump - MC array. C The set of spikes columns - ICS array. C NS=0 IS=0 KK=0 DO 71 I=1,N KP=IPEC(I) KU=IPEC(I+1)-1 NZC(I)=KU-KP+1 ICS(I)=0 J=ICP(I) IF(J.GT.0) GO TO 71 NS=NS+1 ICS(NS)=-J IS=IS+1 IF(ILP(I).GT.0) GO TO 71 KK=KK+1 KP=IPEC(IABS(J)) KU=IPEC(1-J)-1 DO 70 II=1,I-1 DO 69 K=KP,KU IF(IABS(ILP(II)).NE.INDL(K)) GO TO 69 ILP(II)=-ILP(II) MC(KK)=I-II+1 IKM(KK)=IS IS=0 GO TO 71 69 CONTINUE 70 CONTINUE 71 CONTINUE C C20. ------------------------ Write structure synopsis of the C matrix. (if the case) C IF(INFINV.EQ.0) GO TO 72 WRITE(NO,101) 101 FORMAT(//5X,'I. STRUCTURE SYNOPSIS'/5X,24('=')) WRITE(NO,102) 102 FORMAT(/6X,'TOTAL :'/6X,7('-')) WRITE(NO,103) N 103 FORMAT(7X,'Dimension of the matrix -- =',I4) WRITE(NO,104) KK 104 FORMAT(7X,'Number of bumps ---------- =',I4) WRITE(NO,105) NS 105 FORMAT(7X,'Number of spikes --------- =',I4) C IF(KK.EQ.0) GO TO 72 C WRITE(NO,106) 106 FORMAT(/6X,'DISTRIBUTION :'/6X,14('-')) WRITE(NO,107) 107 FORMAT(7X,'BUMP NR. DIMENSION NR. OF SPIKES'/ 1 7X,8('-'),3X,9('-'),3X,13('-')) WRITE(NO,108) (I,MC(I),IKM(I),I=1,KK) 108 FORMAT(9X,I4,8X,I4,10X,I4) WRITE(NO,109) 109 FORMAT(6x,37('-')) C 72 CONTINUE C RETURN C C21. ------------------------ Failure exists. C 1000 WRITE(NO,1001)J 1001 FORMAT(' ***Err: Matrix row number',I5,' has no elements') MESAJ=-1 RETURN 1002 WRITE(NO,1003) LC(II),J 1003 FORMAT(' ***Err: Matrix is structurally singular: ', 1 A3,I5,' cannot be pivoted') MESAJ=-2 RETURN END C-------------------------------------------------------------- * RF02C == RF02C * ############################################################# C------------------------------------------------------------- C Date created: December 14, 1988 C Date last modified: December 9 ,1995 C------------------------------------------------------------- C Purpose: C -------- C This subroutine belongs to a package of subroutines dedicated C to compute the Product Form of Inverse (PFI) with preassigned C pivot procedure of a matrix, to solve corresponding large- C scale systems of linear equations, exploiting the sparsity in C all cases. C RF02C compute the eta-vectors of the matrix according to a C sequence of pivots determined by RF02A subroutine. The pivots C list is in arrays ilp (row pivots list) and icp (column pivots C list). During the calculation of eta-vectors, if the pivot C element equals zero, RF02C takes advantage to change the spikes C of a bump. C C Argument list: C -------------- C call RF02C(n,nz,ia,indl,ipec,ipel,ilp,icp,nzl,nzc,ics,mc, C nreta,nrs,mesaj,no,infinv) C Inputs: C ------- C n The order of the matrix. It is not altered by RF02C. C nz The number of non-zero elements from the matrix. C It is altered by RF02C. On output it is equal with the C number of non-zero elements from the PFI representation C of the matrix, i.e. the number of non-zeros from the C eta-vectors. C ia Integer set by the user in main program for dimensioning C the arrays a, indl and indc. C The matrix is factorized into a product of elementary C transformation matrices Ti. The matrix Ti differ from the C unit matrix I on only one column. This column, ussualy C named eta-vector, is stored in sparse form in arrays a C indl and ipec. The non-zero values of eta-vectors are C stored in array a, their row indices in array indl, and C the starting point of eta-vector k-th in ipec(k). The C most critical situation from the space viewpoint is with C arrays a and indl. These two arrays need occasional C compression to release space of columns which was C transformed in eta-vectors. This compression, actually C performed by the subroutine MRM03A, will not add a C significant overhead to the computational cost if it C happens rarely. In fact, the compression is invocated C only when there is no enough room to store the eta- C vectors in a and indl arrays. So, ia the size of a and C indl must be assigned to such a value to avoid overfrequent C compression. C a Array containing, in column order, the non-zero elements C of the matrix. On exit from RF02C it contains the eta- C vectors of the current matrix. C indl Array containing the row indices of the elements of the C matrix given in array a. On exit from RF02C it contains C the row indices of non-zero elements from the eta- C vectors. C ipec Array containing the address where the columns of the C current matrix starts in arrays a and indl. On exit from C RF02C it contains the starting address of the eta-vectors, C i.e. ipec(k) is the address of starting the k-th eta- C vector. C ilp Array with the rows pivot sequence. C icp Array with the columns pivot sequence. C nzc Array with the number of non-zeros in each column of the C matrix. C ics Array with the indices of the spikes columns. C mc Array with the dimension of the bumps. C nrs The total number of spikes columns from the matrix. C no User supplied number of the output device on which C output is to be written. C infinv Parameter used to write the PFI synopsis of the matrix. C If infinv=0 no such information is supplied. C C Outputs: C -------- C nz The total number of non-zero elements in PFI of the C matrix. C a Array with the non-zeros of the eta-vectors. C indl Array with the row indices of non-zeros from eta-vectors. C ipec Array with the address where the eta-vectors starts in C a and indl arrays. C nreta The number of eta-vectors. A matrix of order n has n C eta-vectors, but RF02C do not store the unit eta-vectors. C mesaj Integer used to transmit the error diagnostics. C mesaj=-3 if ia is too small. C mesaj=-4 if a too small pivot has been found in an C unspike column. C mesaj=-5 if a too small pivot has been found in a spike C column. C mesaj=-6 if no pivot has been found in an unspike column. C C Working arrays: C --------------- C ipel Array of dimension n+1. C nzl Array of dimension n. C-------------------------------------------------------------- * Neculai Andrei *-------------------------------------------------------------- C SUBROUTINE RF02C(N,NZ,IA,A,INDL,IPEC,IPEL,ILP,ICP, 1 NZL,NZC,ICS,MC,NRETA,NRS,MESAJ,NO,INFINV) C C parameter (iain=500000, iainn=100000) REAL*8 A(iain) INTEGER INDL(iain) INTEGER IPEC(iainn), IPEL(iainn) INTEGER ILP(iainn), ICP(iainn) INTEGER NZL(iainn), NZC(iainn) INTEGER ICS(iainn), MC(iainn) C REAL*8 PIV,AUX1,AUX2 REAL*8 TOLPIV,TOLZER C CHARACTER*1 AST C COMMON/TOL/TOLPIV,TOLZER C C C---- ******************** BODY OF THE PROGRAM (RF02C). C---- ************************************************* C C1. ------------------------- The Product Form of Inverse. C The eta-vectors computation. C IA1=IA+1 NZ1=NZ+1 DO 25 I=1,NZ IN=IA1-I IV=NZ1-I A(IN)=A(IV) 25 INDL(IN)=INDL(IV) C II=IA-NZ DO 30 I=1,N 30 IPEL(I)=IPEC(I)+II C AST=' ' C C2. ------------------------- Start of main loop of PFI. C -------------------------- C IS=1 JP=1 NUL=0 IAP=II+1 NRETA=0 C DO 160 L=1,N I=IABS(ILP(L)) J=IABS(ICP(L)) C LL=NRETA+1 IPEC(LL)=JP MC(LL)=I KP=IPEL(J) KU=NZC(J)+KP-1 IF(KP.EQ.KU.AND.A(KP).EQ.1.0D0) GO TO 151 JU=JP+(KU-KP) C C3. ------------------------- Compress the column structure C files if there is not enough room C for new entry. C IF(JU.LT.IAP) GO TO 35 IF(IAP-JP+NUL.LT.KU-KP+1) GO TO 1000 CALL MRM03A(A,INDL,IPEL,NZC,N) KP=IPEL(J) KU=NZC(J)+KP-1 IAP=IAP+NUL NUL=0 35 CONTINUE C C4. ------------------------- Place the column J at the end of C eta-vectors file. C IK=JP-KP DO 40 K=KP,KU II=IK+K A(II)=A(K) 40 INDL(II)=INDL(K) IF(J.EQ.ICS(IS)) GO TO 45 C PIV=0.0D0 DO 41 JJ=JP,JU IF(INDL(JJ).NE.I) GO TO 41 PIV=A(JJ) GO TO 42 41 CONTINUE GO TO 1008 C C5. ------------------------- Perform the stability test. C 42 IF(DABS(PIV).LT.TOLPIV) GO TO 1002 GO TO 145 C 45 CONTINUE C C C6. ------------------------- Determine the eligible spikes C from the current bump. C NES=0 NCS=0 C IF(ICP(L).LT.0.AND.ILP(L).LT.0) GO TO 51 DO 50 K=L+1,N IF(ICP(K).GT.0) GO TO 50 NES=NES+1 NZL(NES)=-ICP(K) IF(ILP(K).GT.0) GO TO 50 GO TO 51 50 CONTINUE C C7. ------------------------- Compress the columns structure C files for possible fill'ins. C 51 CONTINUE IF(JP+N.LT.IAP) GO TO 55 CALL MRM03A(A,INDL,IPEL,NZC,N) KP=IPEL(J) KU=NZC(J)+KP-1 IAP=IAP+NUL NUL=0 55 CONTINUE C C8. ------------------------- Updating the spike column with C active eta-vectors. C 56 DO 80 LL=1,NRETA IC=MC(LL) DO 60 JJ=JP,JU IF(INDL(JJ).NE.IC) GO TO 60 PIV=A(JJ) GO TO 65 60 CONTINUE GO TO 80 65 IP=IPEC(LL) IU=IPEC(LL+1)-1 NZU=0 DO 75 II=IP,IU DO 70 JJ=JP,JU IF(INDL(II).NE.INDL(JJ)) GO TO 70 AUX1=A(JJ) AUX2=A(II)*PIV A(JJ)=AUX1+AUX2 IF(DABS(A(JJ)).LT.TOLZER*DMAX1(DABS(AUX1),DABS(AUX2))) 1 A(JJ)=0.0D0 IF(INDL(JJ).EQ.IC) A(JJ)=A(II)*PIV GO TO 75 70 CONTINUE C C9. ------------------------- New element has been created. C NZU=NZU+1 JJ=JU+NZU IF(JJ.GT.IAP) GO TO 1000 A(JJ)=A(II)*PIV INDL(JJ)=INDL(II) 75 CONTINUE C JU=JU+NZU C 80 CONTINUE C C C10. ------------------------ Decision on the way. Change C eligible spikes or pivoting. C PIV=0.0D0 DO 85 JJ=JP,JU IF(INDL(JJ).NE.I) GO TO 85 PIV=A(JJ) GO TO 90 85 CONTINUE C 90 IF(PIV.EQ.0.0D0) GO TO 95 C C11. ------------------------ Perform the stability test. C IF(DABS(PIV).GT.TOLPIV) GO TO 106 C C12. ------------------------ Change two eligible spikes. C 95 CONTINUE NCS=NCS+1 IF(NCS.GT.NES) GO TO 1004 IJ=NZL(NCS) KP=IPEL(IJ) KU=NZC(IJ)+KP-1 II=0 DO 100 LL=1,L DO 100 K=KP,KU IF(INDL(K).NE.IABS(ILP(LL))) GO TO 100 II=1 GO TO 105 100 CONTINUE C 105 CONTINUE IF(II.NE.0) GO TO 110 GO TO 95 C 106 IS=IS+1 IF(NCS.EQ.0) GO TO 145 GO TO 120 C C13. ------------------------ Place the column IJ at the end C of the eta-vectors file. C 110 IK=JP-KP DO 115 K=KP,KU II=IK+K A(II)=A(K) 115 INDL(II)=INDL(K) JU=JP+KU-KP GO TO 56 C C14. ------------------------ Updating ICP and ICS arrays C according to the permutation C of spikes. C 120 CONTINUE MESAJ=1 DO 125 LL=L+1,N IF(IABS(ICP(LL)).EQ.IJ) GO TO 130 125 CONTINUE 130 ICP(L)=-IJ ICP(LL)=-J DO 135 LL=1,NRS IF(ICS(LL).EQ.IJ) GO TO 140 135 CONTINUE 140 ICS(LL)=J ICS(IS-1)=IJ J=IJ C C15. ------------------------ Pivoting. C 145 PIV=1.0D0/PIV DO 146 JJ=JP,JU IF(INDL(JJ).EQ.I)THEN A(JJ)=PIV ELSE A(JJ)=-A(JJ)*PIV ENDIF 146 CONTINUE C C16. ------------------------ Flush out small elements after C the pivoting. C II=JP-1 DO 150 JJ=JP,JU IF(DABS(A(JJ)).LT.TOLPIV.AND.INDL(JJ).NE.I) GO TO 150 II=II+1 A(II)=A(JJ) INDL(II)=INDL(JJ) 150 CONTINUE JU=II C JP=JU+1 NRETA=NRETA+1 NZE=JU 151 IPEL(J)=-1 NUL=NUL+NZC(J) C 160 CONTINUE C C17. ------------------------ End of the main loop of the PFI C generation. C ------------------------------- C DO 165 K=1,N I=IABS(ILP(K)) J=IABS(ICP(K)) NZL(I)=J NZC(J)=I 165 CONTINUE C C18. ------------------------ Hold permutation (QP)*T in ILP. C Hold permutation (QP) in ICP. C DO 166 K=1,N ILP(K)=NZL(K) 166 ICP(K)=NZC(K) C C19. ------------------------ Product Form of Inverse Synopsis. C (if the case) C IF(INFINV.EQ.0) GO TO 170 IF(MESAJ.EQ.1)AST='*' WRITE(NO,209) 209 FORMAT(//5X,'II. PFI SYNOPSIS'/5X,19('=')) WRITE(NO,210) 210 FORMAT(/16X,'NO. COLUMNS',3X,'NO. NON-ZEROS',3X, 1 'FILLS (%)'/16X,11('-'),3X,13('-'),3X,10('-')) PIV=100./(FLOAT(N)*FLOAT(N)) PIV=PIV*FLOAT(NZ) WRITE(NO,211)N,NZ,PIV 211 FORMAT(6X,'Matrix',7X,I4,11X,I5,7X,F8.2,'%') NZ=NZE PIV=100./(FLOAT(N)*FLOAT(N)) PIV=PIV*FLOAT(NZ) WRITE(NO,212) NRETA,NZ,PIV,AST 212 FORMAT(6X,'Inverse',6X,I4,11X,I5,7X,F8.2,'%',3X,A1/) C 170 CONTINUE NZ=NZE C C RETURN C C20. ------------------------ Failure exists. C C 1000 WRITE(NO,1001) L 1001 FORMAT(' ***Err in RF02C : IA is too small at',I5, 1 '-th eta-vactor') MESAJ=-3 GO TO 1006 1002 WRITE(NO,1003) L,I,J,PIV 1003 FORMAT(' ***Err in RF02C : At',I5,'-th eta-vector', 1 ' too small pivot found in '/19X,'row',I5, 1 ' and unspike column',I5/19X,'pivot value =', 1 E21.13) MESAJ=-4 GO TO 1006 1004 WRITE(NO,1005) L,I,J,PIV 1005 FORMAT(' ***Err in RF02C : At',I5,'-th eta-vector ', 1 'too small pivot found in '/19X,'row',I5, 1 ' and spike column',I5/19X,'pivot value =', 1 E21.13/19X,'No other elligible spike found') MESAJ=-5 GO TO 1006 1008 WRITE(NO,1009)L,I,J 1009 FORMAT(' ***Err in RF02C : At',I5,'-th eta-vector', 1 ' no pivot found in '/19X,'row',I5, 1 ' and unspike column',I5) MESAJ=-6 1006 WRITE(NO,1007) MESAJ 1007 FORMAT(19X,'MESAJR =',I4) NZ=NZE N=L C RETURN END C------------------------------------------------------------- * MFTRAN == MFTRAN * ############################################################# C-------------------------------------------------------------- C Date created: December 14, 1988. C Date last modified: December 8, 1995. C-------------------------------------------------------------- C SUBROUTINE MFTRAN(A,INDL,IPEC,IETA,N,NZ,NRETA,ICP,B) C parameter (iain=500000, iainn=100000) REAL*8 A(iain),B(iainn) INTEGER INDL(iain) INTEGER IPEC(iainn),IETA(iainn),ICP(iainn) C REAL*8 BI,A1,A2 C REAL*8 TOLPIV,TOLZER C COMMON/TOL/TOLPIV,TOLZER C C C--- Forward transformation: C C A**(-1) * B = QP(E ...(E *E *B))...) C IETA(NRETA) IETA(2) IETA(1) C C In input: A,INDL,IPEC,IETA,B. C In output: B. C C---- *************** BODY OF THE PROGRAM (MFTRAN). C ********************************************** C DO 20 I=1,NRETA IC=IETA(I) BI=B(IC) IF(BI.EQ.0.0D0)GO TO 20 JP=IPEC(I) JU=NZ IF(I.NE.NRETA) JU=IPEC(I+1)-1 C C1. ------------------------- Compute : E *B ----> B C IETA(I) C DO 10 J=JP,JU K=INDL(J) IF(K.EQ.IC) THEN KJ=J GO TO 10 ENDIF A1=B(K) A2=A(J)*BI B(K)=A1+A2 IF(DABS(B(K)).LT.TOLZER*DMAX1(DABS(A1),DABS(A2))) 1 B(K)=0.0D0 10 CONTINUE C B(IC)=A(KJ)*BI C 20 CONTINUE C DO 30 I=1,N IF(DABS(B(I)).LE.TOLZER) B(I)=0.0D0 30 CONTINUE C C2. ------------------------- Reorder the vector B according C to the permutation vector ICP. C CALL MRM02A(N,B,ICP) C RETURN END C C------------------------------------------------------------- * MBTRAN == MBTRAN * ############################################################# C-------------------------------------------------------------- C Date created: December 14, 1988 C Date last modified: December 8, 1995 C--------------------------------------------------------------- C SUBROUTINE MBTRAN(A,INDL,IPEC,IETA,N,NZ,NRETA,IP,B) C parameter (iain=500000, iainn=100000) REAL*8 A(iain),B(iainn) INTEGER INDL(iain) INTEGER IPEC(iainn),IETA(iainn),IP(iainn) C REAL*8 SN,SP,AUX C REAL*8 TOLPIV,TOLZER C COMMON/TOL/TOLPIV,TOLZER C C--- Backward transformation: C C A**(-T) * B = (...((B*QP)*E )...*E ) C IETA(NRETA) IETA(1) C C In input: A,INDL,IPEC,IETA,B,IP. C In output: B. B C C--- *********************** BODY OF THE PROGRAM (MBTRAN). C ***************************************************** C DO 20 I=1,NRETA K=NRETA-I+1 IC=IETA(K) JP=IPEC(K) JU=NZ IF(K.NE.NRETA) JU=IPEC(K+1)-1 SN=0.0D0 SP=0.0D0 C C1. ------------------------- Compute : B * E ---> B C IETA(K) C DO 10 J=JP,JU K=INDL(J) K=IP(K) AUX=B(K)*A(J) IF(AUX.LT.0.0D0) THEN SN=SN+AUX ELSE SP=SP+AUX ENDIF 10 CONTINUE K=IP(IC) C B(K)=SP+SN IF(DABS(B(K)).LT.TOLZER*DMAX1(DABS(SP),DABS(SN))) 1 B(K)=0.0D0 C 20 CONTINUE C DO 30 I=1,N IF(DABS(B(I)).LE.TOLZER) B(I)=0.0D0 30 CONTINUE C C2. ------------------------- Reorder the vector B according C to the permutation vector IP. C CALL MRM02A(N,B,IP) C RETURN END C C------------------------------------------------------------- * MRM02A == MRM02A * ############################################################# C-------------------------------------------------------------- C Date created: December 8, 1995 C Date last modified: December 8, 1995 C-------------------------------------------------------------- C SUBROUTINE MRM02A(N,B,IP) C parameter (iain=500000, iainn=100000) REAL*8 B(iainn), BM INTEGER IP(iainn) C I=1 1 IF(I.EQ.IP(I)) GO TO 3 BM=B(I) J=I K=IP(J) 2 B(J)=B(K) IP(J)=-IP(j) J=K K=IP(J) IF(K.NE.I) GO TO 2 B(J)=BM IP(J)=-IP(J) 3 I=I+1 IF(I.GT.N) GO TO 4 IF(IP(I).LT.0) GO TO 3 GO TO 1 4 DO I=1,N IF(IP(I).LT.0) IP(I)=-IP(I) END DO C RETURN END C------------------------------------------------------------- * MRM03A == MRM03A * ############################################################# C-------------------------------------------------------------- C Date created: December 14, 1988 C Date last modified : December 4, 1995 C-------------------------------------------------------------- C SUBROUTINE MRM03A(A,INDL,IPEC,NZC,N) C parameter (iain=500000, iainn=100000) REAL*8 A(iain) INTEGER INDL(iain) INTEGER IPEC(iainn),NZC(iainn) C KK=0 K=N C 5 IPECF=IPEC(K) NZCF=NZC(K) IF(IPECF) 1,2,3 1 KK=KK+NZCF IPEC(K)=0 GO TO 2 3 IF(KK.NE.0) THEN II=IPECF+NZCF-1 4 A(II+KK)=A(II) INDL(II+KK)=INDL(II) II=II-1 IF(II.GE.IPECF) GO TO 4 IPEC(K)=IPECF+KK ENDIF 2 K=K-1 IF(K.GT.0) GO TO 5 C RETURN END C------------------------------------------------------------- * * Last line of the package (1508) Neculai Andrei * REFERENCES * * N. Andrei, * Basic FORTRAN subroutines for large-scale sparse linear programming * problems using PFI factorization. * I.C.I. Technical Report, LSSO-012, April 1983, pp.1-23. * * H.M. Markowitz, * The Elimination Form of the Inverse and its applications to linear programming. * Management Science, Vol.3, 1957, pp.255-269. * * E. Hellerman, D.Rarick, * Reinversion with the preassigned pivot procedure. * Mathematical Programming, Vol.1, 1971, pp.195-216 * * M. Arioli, I.S. Duff, N.I.M. Gould, J.K. Reid, * Use of the P4 and P5 algorithms for in-core factorization of sparse matrices. * SIAM J. Sci. Stat. Comput., Vol.11, No. 5, September 1990, pp.913-927. * * A.M. Erisman, R.G. Grimes, J.G. Lewis, W.G. Poole, Jr., * A structurally stable modification of Hellerman-Rarick's P4 algorithm for * reordering unsymmetric sparse matrices. * SIAM J. Numer. Anal., Vol.22, No.2, April 1985, pp.369-385. * * I.S. Duff, * On permutation to block triangular form, * J. Inst. Maths. Appl., Vol.19, 1977, pp.339-342 * * I.S. Duff, * MA28 - A set of FORTRAN subroutines for sparse unsymmetric linear equations. * AERE Harwell, R8730, Computer Science and Systems Division, Building 8.9, * Didcot, Oxon, June 1977, pp.1-153 * * I.S. Duff, J.K. Reid, * An implementation of Tarjan's algorithm for the block triangularization of a matrix. * ACM Transactions on Mathematical Software, Vol.4, No.2, June 1978, pp.137-147. * * I.S. Duff, J.K. Reid, * ALGORITHM 529. Permutations to Block Triangular Form [F1]. * ACM Transactions on Mathematical Software, Vol.4, No.2, June 1978, pp.189-192 * * I.S. Duff, * On algorithms for obtaining a maximum transversal. * AERE Harwell, CSS 49, Computer Science and Systems Division, Building 8.9, * Didcot, Oxon, October 1978, pp.1-33. * * J.K. Reid, * A sparsity-exploiting variant of the Bartels-Golub decomposition for linear * programming bases. * Mathematical Programming, Vol.24, 1982, pp.55-69. * * J.K. Reid, * TREESOLVE - A FORTRAN package for solving large sets of linear finite-element * equations. * AERE Harwell, CSS 155, Computer Science and Systems Division, Building 8.9, * Didcot, Oxon, March 1984, pp.1-15 * * U.H. Suhl, L. Aittoniemi, * Computing sparse LU-factorizations for large-scale linear programming bases. * Arbeitspapier Nr. 58 / 87, September 1987, * Freie Universitaet Berlin, Fachbereich Wirtschaftswissenschaft, * Angewandte Informatik, * Garystrasse 21, D-1000 Berlin 33, ISBN 3-88398-058-7 * * T.D. Lin, R.S.H. Mah, * Hierarchical partition - A new optimal pivoting algorithm. * Mathematical Programming, Vol.12, 1977, pp.260-278. * * T.D. Lin, R.S.H. Mah, * Sparse computation system for process design and simulation. Part I: Data * structures and processing techniques. * AIChE Journal, Vol.24, No.5, September 1978, pp.830-839. * Part II: A performance evaluation based on the simulation of a natural * gas liquefaction process. * AIChE Journal, Vol.24, No.5, September 1978, pp.839-848. * * R.V. Helgason, J.L. Kennington, * Spike swapping in basis reinversion. * Naval Research Logistics Quarterly, Vol.27, No.4, December 1980, pp.697-701. * * R.V. Helgason, J.L. Kennington, * A note on splitting the bump in an elimination factorization. * Naval Research Logistics Quarterly, Vol.29, No.1, March 1982, pp.169178. * * J.J.H. Forrest, J.A. Tomlin, * Updated triangular factors of the basis to maintain sparsity in the * Product Form Simplex Method. * Mathematical Programming, Vol.2, 1972, pp.263-278. * * R.D. McBride, * A bump triangular dynamic factorization algorithm for the simplex method. * Mathematical Programming, Vol.18, 1980, pp.49-61. * * R.D. McBride, * A spike collective dynamic factorization algorithm for the simplex method. * Management Science, Vol.24, No.10, June 1878, pp.1031-1042 * * R.E. Marsten, * The design of the XMP linear programming library. * ACM Transactions on Mathematical Software, Vol.7, No.4, December 1981, pp.481-497. * * U. Schendel, * Sparse Matrices: Numerical aspects with applications for scientists and engineers. * Ellis Horwood Limited, 1989, Chichester, England * * R.P. Tewarson, * Sparse Matrices * Academic Press, 1973, New York. * * Last line of References ----------------------------------- Neculai Andrei