*-------------------------------------------------------------- * Date created: April 4, 1983 * Date modified: November 10, 1995 (Updated) *-------------------------------------------------------------- * * Example 12 * A random generated system with known solution * ------------------------------------------------- * * Solving linear algebraic systems of equations * * Ax=b or ATx=b * ======================= * * by * LU factorization of sparse matrices. * * * Dr. Neculai Andrei * Research Institute for Informatics * Center for Advanced Modeling and Optimization * 8-10, Averescu Avenue, Bucharest 1, Romania * E-mail: nandrei@ici.ro * * * MAIN PROGRAM * ================ * Main program for a large-scale sparse linear system equations * Ax=b solving taking into account the sparsity of the A matrix * using LU factorization of the coefficient matrix with Markowitz's * pivot strategy. (AT is the transpouse of the matrix A.) * * * The calling sequence is as follows: * * * LU -----> RF01A * | ---> RM01A, * | ---> RM04A. * ----> RS01A * * where: * * LU Main program for solving linear algebraic sysyems of * equations Ax=b or ATx=b using LU factorization of the * matrix with Markowitz's pivot selection strategy. * * RF01A Subroutine for LU factorization of the matrix A. * * RS01A Subroutine for solving large-scale systems of linear * equations computing: A**(-1)*b or A**(-T)*b, for a * given vector b, using the LU factorization of the * matrix A (given by RF01A subroutine), exploiting the * sparsity in all cases. * * RM01A Subroutine for sorting the non-zeros of a sparse matrix * from arbitrary order to column order, but unordered * within each column. * * RM04A Subroutine for compressing the column / row array of U * factor (from the LU factorization of a sparse matrix) in * order to eliminate the spaces between columns / rows. * * * Neculai Andrei ****************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) double precision a(8000000), b(2000) double precision z(2000) integer indl(8000000), indc(8000000), ip(2000,2), iz(2000,8) integer lc(2,3) double precision mesaj logical transp double precision mat(1000,1000), x(2000) open(unit=2,file='sys.out',status='unknown') n=1000 ia=8000000 * Generate the random matrix do i=1,n do j=1,n mat(i,j)=random() end do end do * Generate the sparse matrix do i=1,n do j=1,n if(mat(i,j) .lt. 0.5d0) mat(i,j) = 0.d0 end do end do * Generate the solution vector. do i=1,n x(i) = float(i) end do * Generate the RHS term. do i=1,n b(i)=0.d0 do j=1,n b(i)=b(i)+mat(i,j)*x(j) end do end do * Generate INDL, INDC, A vectors k=1 do i=1,n do j=1,n if(mat(i,j) .ne. 0.d0) then indl(k)=i indc(k)=j a(k)=mat(i,j) k=k+1 end if end do end do nz=k-1 write(2,10)n, nz, ia 10 format(4x,'n=',i7,4x,'nz=',i7,4x,'ia=',i7) C--- LU factorization of the matrix with Markowitz pivot * strategy. mesaj=0.d0 call rf01a(a,indl,indc,n,nz,ia,ip,iz,z,mesaj,0.1d0) C--- Solve the system Ax=b. transp = .false. call rs01a(a,indl,indc,n,ia,ip,iz,z,b,mesaj,transp) C--- Write out the solution of the system. write(2,4) (i,b(i),i=1,n) 4 format(5x,i4,4x,e20.13) stop end *************************************************************** * REFERENCES * ============== * * 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. C-------------------------------------------------------------------------------- * * C Dr. Neculai Andrei C Research Institute for Informatics C 8-10, Averescu Avenue, Bucharest 1, Romania C E-mail: nandrei@ici.ro C-------------------------------------------------------------- * Last line LU.FOR *--------------------------------------------------------------- * subroutine rf01a(a,indl,indc,n,nz,ia,ip,iz,z,mesaj,cpv) * C0. Purpose. * ======== * This subroutine belongs to a package of subroutines dedicated * to LU factorization of a large-scale sparse matrix, solve the * corresponding systems of linear equations and update the LU * factorization when a column of the matrix is modified, * exploiting the sparsity in all cases. * * Its primary application is to be used for large-scale sparse * matrices inversion in LU representation of the matrices and * handling large-scale linear programming bases in simplex * algorithms. * * It implements a sparse variant of the Bartles-Golub algorithm * with the Markowitz pivot preassigned procedure. * * The RF01A subroutine produces a LU factorization of a given * sparse matrix, where L is an inferior-triangular matrix whose * inverse is held as a product of elementary matrices, and U is * a permutation of a superior-triangular matrix. * * For more details about the LU factorization of a large-scale * sparse matrix with Markowitz's pivotal strategy, please see; * * 1. J.K. Reid, FORTRAN subroutines for handling sparse linear * programming bases. Report AERE - R.8269, Harwell, 1976. * 2. J.K. Reid, A sparsity-exploiting variant of Bartles-Golub * decomposition for linear programming bases. Mathematical * Programming, vol. 24, (1982) pp.55-69. * 3. N. Andrei, C. Rasturnoiu, Sparse Matrices and their * Applications. (Matrice rare si aplicatiile lor). Editura * Tehnica, Bucuresti, 1983. *----------------------------------------------------------------- * C1. Argument list. * ============== * CALL RF01A(A,INDL,INDC,N,NZ,IA,IP,IZ,Z,MESAJ,CPV) * * Parameters: * =========== * A A real array of length IA. On entry to RF01A it must be * set to contain, in any order, the non-zeros of the matrix. * On exit from RF01A it contains the L and U factors of the * matrix. * INDL An integer array of length IA. On entry to RF01A, INDL(k) * must be set to contain the row number of the non-zero held * in A(k), k=1,2,...,NZ. * On exit from RF01A it contains the row number of the non- * zero from L factor of the matrix. * INDC An integer array of length IA. On entry to RF01A, INDC(k) * must be set to the column number of the non-zero held in * A(k), k=1,2,...,NZ. * On exit from RF01A it contains the column number of the * non-zero from U factor of the matrix. * N Integer, must be set by the user to the order of the matrix. * It is not altered by the RF01A subroutine. * NZ Integer, must be set to the number of non-zero elements in * the matrix. It is not altered by RF01A. * IA Integer, must be set to indicate the size of arrays A, * INDL and INDC. IA must exceed NZ by a margin sufficient to * avoid overfrequent compression of the A, INDL and INDC * arrays. It is not altered by RF01A. * IP Integer working array of dimension (N,2). * IZ Integer working array of dimension (N,8). * Z Real working array of length N. * MESAJ Real parameter used to transmit the error diagnostics and * information about the numerical stability of the * factorization. After a successful entry to RF01A, MESAJ * is positive and equal to the modulus of the largest element * in any of the reduced matrices. After an unsuccessful entry * a mesaje is issued and MESAJ is set negative to indicate * one of the following situations: * -1 N or NZ is not positive, * -2 LIN (COL) I has no elements, * -3 K-th element is out of the range with row I and * column J, * -4 Duplicate element in row I and column J, * -5 The following LIN (COL) are dependent, * -6 Ia is too small. * CPV Real parameter set in the interval (0,1] to control the * choice of pivots. The value 0.1 has been found to be * satisfactory in the test examples. * C2. Use of COMMON. * ============== * RF01A contains the following COMMON BLOCK: * COMMON /RF01C/mic,no,eps,ieml,iemu,ncp,ilin,icol * The following elements must be set by data in common: * mic real parameter to control the stability. Any elements * of the upper triangular factor that are less than MIC are * reset to zero value. Suggested value for MIC: 0.0. * NO Integer, index of the line printer. * EPS Real parameter set by the user to the relative accuracy * of the floating-point computations. * C3. Other Subroutines. * ================== * The following subroutines are called by RF01A: * RM01A to sort the non-zeros of a sparse matrix from arbitrary * order to column order, but unordered withing each column. * RM04A to compress an array of positive integers with the * columns or rows structure of the U and L factors. * * April, 1983 * Neculai Andrei *----------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) double precision a(ia), z(n), mic, mesaj integer indl(ia), indc(ia), iz(n,8), ip(n,2), lc(2,3) common /rf01c/mic,no,eps,ieml,iemu,ncp,ilin,icol data lc(1,1),lc(1,2),lc(1,3), lc(2,1),lc(2,2),lc(2,3)/ *1hL,1hI,1hN, 1hC,1hO,1hL/ if(cpv .gt. 1.d0) cpv=1.d0 if(cpv .lt. eps) cpv = eps if(n .lt. 1 .or. nz .lt. 1) go to 58 mesaj = 0.d0 do 1 i=1,n z(i) = 0.d0 do 1 j=1,5 1 iz(i,j) = 0 * * Flush out small entries, count elements in row and columns. * l=1 iemu=nz do 4 ifa=1,nz if(l .gt. iemu) go to 5 do 2 k=l,iemu if(dabs(a(k)) .le. mic) go to 3 i=indl(k) j=indc(k) mesaj=amax1(dabs(a(k)),mesaj) if(i .lt. 1 .or. i .gt. n) go to 62 if(j .lt. 1 .or. j .gt. n) go to 62 iz(i,1)=iz(i,1) + 1 2 iz(j,2) = iz(j,2) +1 go to 5 3 l=k a(l)=a(iemu) indl(l)=indl(iemu) indc(l)=indc(iemu) 4 iemu=iemu-1 5 ieml=0 ilin=iemu icol=ilin * * MCP is the maximum number of compresses permited. * mcp=max0(n/10,20) ncp=0 k=1 do 6 ir=1,n k=k+iz(ir,2) ip(ir,2)=k do 6 l=1,2 if(iz(ir,l) .le. 0) go to 60 6 continue * * Reorder the matrix by rows. * call rm01a(a,indc,indl,n,iemu,0,ip) kl=iemu do 8 ii=1,n ir=n+1-ii kp=ip(ir,1) do 7 k=kp,kl j=indc(k) if(iz(j,5) .eq. ir) go to 64 iz(j,5)=ir kr=ip(j,2)-1 ip(j,2)=kr 7 indl(kr)=ir 8 kl=kp-1 * * Set up linked lists of rows and columns with equal number * of non-zero elements. * do 9 l=1,2 do 9 i=1,n nz=iz(i,l) in=iz(nz,l+2) iz(nz,l+2)=i iz(i,l+6)=in iz(i,l+4)=0 9 if(in .ne. 0) iz(in,l+4)=i * * Start of the main elimination loop. * ====================== * do 53 ipv=1,n * * Find the pivot using the Markowitz strategy. * JCOST is cheapest Markowitz cost found so far in row IPP and * column JP. * jcost=n*n do 18 nz=1,n if(jcost .le. (nz-1)**2) go to 19 j=iz(nz,4) * * Determine columns with NZ non-zero elements. * do 13 ifa=1,n if(j .le. 0) go to 14 kp=ip(j,2) kl=kp+iz(j,2)-1 do 12 k=kp,kl i=indl(k) kcost=(nz-1)*(iz(i,1)-1) if(kcost .ge. jcost) go to 12 if(nz .eq. 1) go to 11 * * Find the largest element in row of potential pivot. * amax=0.d0 k1=ip(i,1) k2=iz(i,1)+k1-1 do 10 kk=k1,k2 amax=amax1(amax,dabs(a(kk))) 10 if(indc(kk) .eq. j) kj=kk * * Stability test. * if(dabs(a(kj)) .lt. amax*cpv) go to 12 11 jcost=kcost ipp=i jp=j if(jcost .le. (nz-1)**2) go to 19 12 continue 13 j=iz(j,8) * * Determine the rows with NZ non-zero elements. * 14 i=iz(nz,3) do 17 ifa=1,n if(i .le. 0) go to 18 amax=0.d0 kp=ip(i,1) kl=kp+iz(i,1)-1 * * Find the largest element in the row. * do 15 k=kp,kl 15 amax=amax1(dabs(a(k)), amax) au=amax*cpv do 16 k=kp,kl * * Stability test. * if(dabs(a(k)) .lt. au) go to 16 j=indc(k) kcost=(nz-1)*(iz(j,2)-1) if(kcost .ge. jcost) go to 16 jcost=kcost ipp=i jp=j if(jcost .le. (nz-1)**2) go to 19 16 continue 17 i=iz(i,7) 18 continue * * Remove rows and columns involved in elimination. * 19 kp=ip(jp,2) kl=iz(jp,2)+kp-1 do 24 l=1,2 do 23 k=kp,kl if(l .eq. 2) go to 20 i=indl(k) go to 21 20 i=indc(k) 21 il=iz(i,l+4) in=iz(i,l+6) if(il .eq. 0) go to 22 iz(il,l+6)=in go to 23 22 nz=iz(i,l) iz(nz,l+2)=in 23 if(in .gt. 0) iz(in,l+4)=il kp=ip(ipp,1) 24 kl=kp+iz(ipp,1)-1 * * Store the pivot. * iz(ipp,5)=-ipv iz(jp,6)=-ipv * Eliminate pivot row from the column array and find the pivot * in row array. * do 27 k=kp,kl j=indc(k) kpc=ip(j,2) iz(j,2)=iz(j,2)-1 klc=kpc+iz(j,2) do 25 kc=kpc,klc if(ipp .eq. indl(kc)) go to 26 25 continue 26 indl(kc)=indl(klc) indl(klc)=0 27 if(j .eq. jp) kr=k au=a(kr) a(kr)=a(kp) a(kp)=au indc(kr)=indc(kp) indc(kp)=jp * * Perform the elimination. * nzc=iz(jp,2) if(nzc .eq. 0) go to 48 do 47 nc=1,nzc kc=ip(jp,2)+nc-1 ir=indl(kc) kr=ip(ir,1) krl=kr+iz(ir,1)-1 do 28 knp=kr,krl if(jp .eq. indc(knp)) go to 29 28 continue 29 am=a(knp) a(knp)=a(kr) a(kr)=am indc(knp)=indc(kr) indc(kr)=jp am=-a(kr)/a(kp) * * Compress the row array. * if(ilin+iz(ir,1)+iz(ipp,1)+ieml .le. ia) go to 30 if(ncp.ge.mcp .or. iemu+iz(ir,1)+iz(ipp,1)+ieml.gt.ia) go to 69 * call rm04a(a,ip,n,ia,indc,iz,.true.) * kp=ip(ipp,1) kr=ip(ir,1) 30 krl=kr+iz(ir,1)-1 kq=kp+1 kpl=kp+iz(ipp,1)-1 * * Place the pivot row in z array. * if(kq .gt. kpl) go to 32 do 31 k=kq,kpl j=indc(k) 31 z(j)=a(k) 32 ip(ir,1)=ilin+1 indc(kr)=0 kr=kr+1 if(kr .gt. krl) go to 37 do 36 ks=kr,krl j=indc(ks) au=a(ks)+am*z(j) indc(ks)=0 * * Elimination of the very small elements from U array. * if(dabs(au) .le. mic) go to 33 mesaj=amax1(mesaj,dabs(au)) ilin=ilin+1 a(ilin)=au indc(ilin)=j go to 36 33 iemu=iemu-1 k=ip(j,2) kl=k+iz(j,2)-1 iz(j,2)=kl-k do 34 kk=k,kl if(indl(kk) .eq. ir) go to 35 34 continue 35 indl(kk)=indl(kl) indl(kl)=0 36 z(j)=0.d0 * * Scan pivot row for fills. * 37 if(kq .gt. kpl) go to 45 do 44 ks=kq,kpl j=indc(ks) au=am*z(j) if(dabs(au) .le. mic) go to 44 ilin=ilin+1 a(ilin)=au indc(ilin)=j iemu=iemu+1 * * Create new element in column array. * nz=iz(j,2) k=ip(j,2) kl=k+nz-1 * * If it is possible place new element at the end of the present * entry. * if(kl .ne. icol) go to 38 if(icol+ieml .ge. ia) go to 40 icol=icol+1 go to 39 38 if(indl(kl+1) .ne. 0) go to 40 39 indl(kl+1)=ir go to 43 40 if(icol+ieml+nz+1 .lt. ia) go to 41 if(ncp .ge. mcp .or. iemu+ieml+nz+1 .ge. ia) go to 69 * * Compress the column array. * call rm04a(a,ip(1,2),n,ia,indc,iz(1,2),.false.) k=ip(j,2) kl=k+nz-1 41 ip(j,2)=icol+1 do 42 kk=k,kl icol=icol+1 indl(icol)=indl(kk) 42 indl(kk)=0 icol=icol+1 indl(icol)=ir 43 mesaj=amax1(mesaj,dabs(au)) iz(j,2)=nz+1 44 z(j)=0.d0 45 iz(ir,1)=ilin+1-ip(ir,1) if(ieml+icol+1 .le. ia) go to 46 * * Compress the column array if it is necessary. * if(ncp .ge. mcp) go to 69 call rm04a(a,ip(1,2),n,ia,indc,iz(1,2),.false.) 46 k=ia-ieml ieml=ieml+1 a(k)=am indl(k)=ipp indc(k)=ir iemu=iemu-1 47 continue * * Insert the rows and the columns involved in elimination in * linked lists of equal numbers of non-zero elements. * 48 k1=ip(jp,2) k2=iz(jp,2)+k1-1 iz(jp,2)=0 do 53 l=1,2 if(k2 .lt. k1) go to 52 do 51 k=k1,k2 if(l .eq. 2) go to 49 ir=indl(k) indl(k)=0 go to 50 49 ir=indc(k) 50 nz=iz(ir,l) if(nz .le. 0) go to 66 in=iz(nz,l+2) iz(ir,l+6)=in iz(ir,l+4)=0 iz(nz,l+2)=ir 51 if(in .ne. 0) iz(in,l+4)=ir 52 k1=ip(ipp,1)+1 53 k2=iz(ipp,1)+k1-2 do 54 i=1,n j=-iz(i,5) iz(j,3)=i j=-iz(i,6) iz(j,4)=i 54 iz(i,2)=0 do 55 i=1,n kp=ip(i,1) kl=iz(i,1)+kp-1 do 55 k=kp,kl j=indc(k) 55 iz(j,2)=iz(j,2)+1 k=1 do 56 i=1,n k=k+iz(i,2) 56 ip(i,2)=k icol=k-1 do 57 ii=1,n i=iz(ii,3) kp=ip(i,1) kl=iz(i,1)+kp-1 do 57 k=kp,kl j=indc(k) kn=ip(j,2)-1 ip(j,2)=kn 57 indl(kn)=i go to 73 * * The following lines implements te failure exists. * 58 if(no .gt. 0) write(2,59) 59 format(/5x,'N or NZ is not positive') mesaj=-1.d0 go to 71 60 if(no .gt. 0) write(2,61) (lc(l,i),i=1,3),ir 61 format(/5x,3a1,i5,'has no elements') mesaj=-2.d0 go to 71 62 if(no .gt. 0) write(2,63) k,i,j 63 format(/5x,i4,'-th element is out of the range', * 'with row ',i4,' and column ',i4) mesaj=-3.d0 go to 71 64 if(no .gt. 0) write(2,65) ir,j 65 format(/5x,'Duplicate element in row ',i4,' and column ',i4) mesaj=-4.d0 go to 71 66 ipv=ipv+1 iz(ipv,1)=ir do 67 i=1,n ii=-iz(i,l+4) 67 if(ii .gt. 0) iz(ii,1)=i if(no .gt. 0) write(2,68) (lc(l,i),i=1,3),(iz(i,1),i=1,ipv) 68 format(/5x,'The following ',3a1,' are dependent:', * /5x,10(3x,i5)) mesaj=-5.d0 go to 71 69 if(no .gt. 0) write(2,70) 70 format(/5x,'IA is too small') mesaj=-6.d0 71 if(no .gt.0) write(2,72) 72 format(/5x,'Error in RF01A subroutine',/) 73 return end *------------------------------------------------------------- * Last line RF01A BLOCK DATA *------------------------------------------------------------- common/rf01c/mic,no,eps,ieml,iemu,ncp,ilin,icol double precision mic data mic, no, eps / 0.d0, 6, 1.0d-7/ end *-------------------------------------------------------------- * Last line block data *--------------------------------------------------------------- * subroutine rs01a(a,indl,indc,n,ia,ip,iz,z,b,mesaj,transp) * C0. Purpose. * ======== * This subroutine belongs to a package of subroutines dedicated * to solve large-scale systems of linear equations by computing * A**(-1)*B, or A**(-T)*B for a given vector B, using the LU * factorization of the matrix A, exploiting sparsity in all * cases. * It is recommended that RS01A subroutine to be used in direct * conjunction with RF01A subroutine, or any other equivalent. * For more details about the solving of a large-scale sparse * system of linear equations with the coefficient matrix LU * factorized, please see: * 1. J.K. Reid, FORTRAN subroutines for handling sparse linear * programming bases. Report AERE - R.8269, Harwell, 1976. * 2. J.K. Reid, A sparsity-exploiting variant of Bartles-Golub * decomposition for linear programming bases. Mathematical * Programming, vol. 24, (1982) pp.55-69. * 3. N. Andrei, C. Rasturnoiu, Sparse Matrices and their * Applications. (Matrice rare si aplicatiile lor). Editura * Tehnica, Bucuresti, 1983. *----------------------------------------------------------------- * C1. Argument list. * ============== * CALL RS01A(A,INDL,INDC,N,IA,IP,IZ,Z,B,MESAJ,TRANSP) * * A Real array of length IA. On entry to RS01A it must be set * to contain the non-zeros of the L and U factors. * INDL Integer array of length IA. On entry to RS01A, INDL(k) * must be set to contain the row number of the non-zero held * in A(k), k=1,2,...,IA. * It is not altered by RS01A. * INDC Integer array of length IA. On entry to RS01A, INDC(k) * must be set to contain the column number of the non-zero * held in A(k), k=1,2,...,IA. * It is not altered by RS01A. * N Integer, must be set to the order of the matrix. * It is not altered by RS01A. * IA Integer, must be set to indicate the size of A, INDL and * INDC arrays. * It is not altered by RS01A. * IP Integer working array of dimension (N,2). * IZ Integer working array of dimension (N,4). * Z Real working array of length N. * B Real array of length N. On entry to RS01A it contains the * right-hand-side term of the system. * On exit from RS01A it contains the solution: * A**(-1)*B if TRANSP=.FALSE. or * A**(-T)*B if TRANSP=.TRUE. . * MESAJ Integer used to transmit the error diagnostics. * A non-negative value on entry to RS01A determine the * running of the subroutine. In a string of subroutines * if a previous subroutine gives a negative value for * MESAJ, then RS01A is suppressed. * TRANSP Logical parameter which must be set by the user to: * .FALSE. if A**(-1)*B is required from RS01A or * .TRUE. if A**(-T)*B is required from RS01A. * C2. Use of common. * ============== * RS01A contains the following common block: * COMMON/RF01C/MIC,NO,EPS,IEML,IEMU,NCP,ILIN,ICOL * (Please see RF01A subroutine where all these parameters * are explained.) * * April, 1983 * Neculai Andrei *---------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) double precision a(ia),b(n),z(n),mic,mesaj integer indl(ia),indc(ia),iz(n,4),ip(n,2) logical transp common/rf01c/mic,no,eps,ieml,iemu,ncp,ilin,icol * if(mesaj .lt. 0.d0) go to 13 kl=ia-ieml+1 if(transp) go to 8 * * Compute A**(-1)*B where A is LU factorized. * * Partial updating of B by the inverse of L. * if(ieml .le. 0) go to 2 l=ia+1 do 1 kk=1,ieml k=l-kk i=indl(k) if(b(i) .eq. 0.d0) go to 1 j=indc(k) b(j)=b(j)+a(k)*b(i) 1 continue 2 do 3 i=1,n z(i)=b(i) 3 b(i)=0.d0 * * Partial updating of B by the inverse of U. * n1=n+1 do 7 ii=1,n i=n1-ii i=iz(i,3) am=z(i) kp=ip(i,1) if(kp .gt. 0) go to 5 kp=-kp ip(i,1)=kp nz=iz(i,1) k1=kp+nz-1 k2=kp+1 do 4 k=k2,k1 j=indc(k) 4 am=am-a(k)*b(j) 5 if(am .eq. 0.d0) go to 7 j=indc(kp) b(j)=am/a(kp) kp1=ip(j,2) k1=iz(j,2)+kp1-1 if(k1 .eq. kp1) go to 7 k2=kp1+1 do 6 k=k2,k1 i=indl(k) 6 ip(i,1)=-iabs(ip(i,1)) 7 continue go to 15 * * Compute A**(-T)*B where A is LU factorized. * * Partial updating of B by the inverse of U transpose. * 8 do 9 i=1,n z(i)=b(i) 9 b(i)=0.d0 do 11 ii=1,n i=iz(ii,4) am=z(i) if(am .eq. 0.d0) go to 11 j=iz(ii,3) kp=ip(j,1) am=am/a(kp) b(j)=am k1=iz(j,1)+kp-1 if(kp .eq. k1) go to 11 k2=kp+1 do 10 k=k2,k1 i=indc(k) 10 z(i)=z(i)-am*a(k) 11 continue * * Partial updating of B by inverse of L transpose. * if(kl .gt. ia) return do 12 k=kl,ia j=indc(k) if(b(j) .eq. 0.d0) go to 12 i=indl(k) b(i)=b(i)+a(k)*b(j) 12 continue go to 15 13 if(no.gt.0) write(2,14) 14 format(/5x,'Error in RS01A because the previous', *' subroutine have gave error') 15 return end * *------------------------------------------------------------- * Last line RS01A *--------------------------------------------------------------- * subroutine rm01a(a,indl,indc,n,nz,jdep,ipec) * C0. Purpose. * ======== * This subroutine sorts the non-zeros of a sparse matrix from * arbitrary order to column order, but unordered within each * column. * C1. Argument list. * ============== * A Real array of length NZ. On entry to RM01A it must be * set to contain in any order the non-zero elements of * matriz A. * On exit from RM01A the non-zero elements of the matrix * are reordered so that column 1 precedes column 2 which * precedes column 3, etc. However, the order within * columns is arbitrary. * INDL Integer array of length NZ. On entry to and exit from * RM01A the absolute value of INDL(k) is the row number * of the element in A(k). * INDC Integer array of length NZ. On entry to and exit from * RM01A INDC(k)+JDEP is the column number of the element * held in A(k). * N Integer, must be set by the user to the number of * matrix columns. It is not altered by RM01A. * Nz Integer, must be set by the user to the number of * non-zeros in A. It is not altered by RM01A. * JDEP Integer, must be set by the user to his required * displacement for column numbers, in range [0,32767]. * Normally zero will be suitable, but positive values * permit matrices with up to 65534 columns to be handeled. * JDEP is not altered by RM01A. * IPEC Integer array of length N. On exit from RM01A, IPEC(i) * contains the pinr to the start of column i, i=1,2,...,N. * * Note: This subroutine will be called by other library * subroutines, but not by the user directly. * * February, 1983 * Neculai Andrei *-------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer indl(nz), indc(nz), ipec(n) double precision a(nz) * nul=-jdep do 1 j=1,n 1 ipec(j)=0 do 2 k=1,nz j=indc(k)+jdep ipec(j)=ipec(j)+1 2 continue k=1 do 3 j=1,n kl=k+ipec(j) ipec(j)=k 3 k=kl do 5 i=1,nz j1=indc(i)+jdep if(j1 .eq. 0) go to 5 a1=a(i) i1=indl(i) indc(i)=nul do 4 j=1,nz loc=ipec(j1) ipec(j1)=ipec(j1)+1 a2=a(loc) i2=indl(loc) j2=indc(loc) a(loc)=a1 indl(loc)=i1 indc(loc)=nul if(j2 .eq. nul) go to 5 a1=a2 i1=i2 4 j1=j2+jdep 5 continue * l1=1 do 6 j=1,n l2=ipec(j) ipec(j)=l1 6 l1=l2 * n1=n-1 do 8 i=1,n1 ks=ipec(i+1)-ipec(i) do 7 k=1,ks i1=ipec(i)+k-1 7 indc(i1)=i 8 continue * ks=nz+1-ipec(n) do 9 k=1,ks i1=ipec(n)+k-1 9 indc(i1)=n * return end *---------------------------------------------------------------- * Last line RM01A *--------------------------------------------------------------- * subroutine rm04a(a,ip,n,ia,izl,iz,lincol) * C0. Purpose. * ======== * This subroutine compress the column / row arrays of U factor * (from LU factorization of a sparse matrix) in order to * eliminate the spaces between columns / rows. * The length of compressed array is placed in ILIN if LINCOL= * .true., or in ICOL if LINCOL=.false. * C1. Argument list. * ============== * A Real array of length IA. On entry to RM04A it contains * the elements of the U factor of the matrix. On exit * from RM04A it contains the non-zeros of the U factor, * forward compressed. * IP Integer array of length N. On entry to and exit from * RM04A, IP(i) contains the poinr to the start of column/ * row i, i=1,2,...,n. * N Integer, must be set to the order of the matrix. It is * not altered by RM04A. * IA Integer, must be set to indicate the size of arrays * A and IZL. It is not altered by RM04A. * IZL Integer working array of length IA. * LINCOL Logical parameter. * * April, 1983 * Neculai Andrei *----------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) double precision a(ia), mic integer ip(n), izl(ia), iz(n) logical lincol common/rf01c/mic,no,eps,ieml,iemu,ncp,ilin,icol * ncp=ncp+1 do 1 j=1,n nz=iz(j) if(nz .le. 0) go to 1 k=ip(j)+nz-1 iz(j)=izl(k) izl(k)=-j 1 continue * k1=0 ip1=0 k2=icol if(lincol) k2=ilin * do 3 k=1,k2 if(izl(k) .eq. 0) go to 3 k1=k1+1 if(lincol) a(k1)=a(k) if(izl(k) .ge. 0) go to 2 j=-izl(k) izl(k)=iz(j) ip(j)=ip1+1 iz(j)=k1-ip1 ip1=k1 2 izl(k1)=izl(k) 3 continue * if(lincol) ilin=k1 if(.not. lincol) icol=k1 * return end * * April, 1983 * Neculai Andrei *--------------------------------------------------------------- * Last line RM04A *------------------------ Last line LU package ----------------*