*--------------------------------------------------------------
* Date created:   April 4, 1983
* Date modified:  November 10, 1995  (Updated)
*--------------------------------------------------------------
*
*                           Example 10
*           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
      nz=n*n
      ia=8000000
      
* Generate the random matrix
      do i=1,n
        do j=1,n
          mat(i,j)=random()
        end do
      end do
      k=1
      do i=1,n
        do j=1,n
          indl(k)=i
          indc(k)=j
          a(k)=mat(i,j)
          k=k+1
        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    

      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 ----------------*              