C  Main Fortran program for solving a linear symmetric indefinite sparse system
C  by the routine BLKFCLT.f
C
C  Authors: Valeria RUGGIERO, Silvia BONETTINI, Federica TINTI
C  Last version: July 6, 2005
C
        program example
        implicit none
        integer*4       nmax,max,maxsub
        parameter     (nmax=100,max=50000,maxsub=100)
        integer*4      neqns   ! Number of equations
	  integer*4      cachsz  ! Dimension of the cache memory
	  integer*4      nsup    ! Rows number of the positive definite block
C
	  integer*4      adjncy(max), xadj(nmax) ! row indices and column pointers of the lower triangular part of the matrix
        real*8         a(max)                  !nonzero entries of the lower triangular part of the matrix af
	  integer*4      adjf(max), xadjf(nmax)  !row indices and column pointers of the
	  real*8         af(max)                  !nonzero entries of the matrix af
	  integer*4      perm(nmax),invp(nmax)   !permutation vector and its inverse
		
        integer*4     lindx(maxsub), xlindx(nmax) !row indices for each supernode
	  integer*4     xlnz(nmax)                  !column pointer of the Choleski factor
	  real*8        lnz(max)                    !choleski factor
	  integer*4     iwork(7*nmax)            ! integer working array
	  real*8        tmpvec(nmax)             ! double precision working storage
        integer*4     temp(max)                ! temporaty vector fot the lower trianguler part of af
        integer*4     offset(max)
	  integer*4     snode(nmax)              !maps each column into the supernode containing it
	  integer*4     split(nmax)              !splitting of supernode so that they fit into cache
	  integer*4     colcnt(nmax)             !array containing the number of nonzeros in each column of the
					       !Cholesky factor
        integer*4     xsuper(nmax+1)           !supernode partitioning
	  integer*4     tmpsiz                   !size of working storage required by BLKFCT
	  integer*4     nofsub                   !upper boud for the number of nonzero entries of the Cholesky factor
	  integer*4     iwsiz                    !size of integer working storage
	  integer*4     nnzl                     !number of nonzero entries of the Cholesky factor
	  integer*4     nsub                     !number of subscripts
	  integer*4     nsuper                   !number of supernodes
	  integer*4     nnza                     !number of nonzero entries of af
	  integer*4     iflag                    !
        integer*4     ij,i

        real*8        diag(nmax)               !diagonal entries of af
        real*8        b(nmax)                  !right hand side
	  real*8        sol(nmax)                !solution of the system
        external mmpym,smxpym

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !This program computes the solution of the system af*sol=b where!
      !                                                               !
      !             3 1 1 1 1         3                               !
      !             1 3   4          -3                               !
      ! af=         1   3        b=   1                               !
      !             1 4               1                               !
      !             1                 1                               !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !All the matrices have to be stored in Compressed Column Storage form

      !!!!!!!!Store in a the nonzero entries of af, without the diagonal elements

        do i=1,5
           a(i) = 1.
        end do
        a(6)    = 4.
        a(7)    = 1.
        a(8)    = 1.
        a(9)    = 4.
        a(10)   = 1.
        a(11)   = 0.

      !!!!!!!!!!!!adjncy contains the row indices of a

        adjncy(1) = 2
        adjncy(2) = 3
        adjncy(3) = 4
        adjncy(4) = 5
        adjncy(5) = 1
        adjncy(6) = 4
        adjncy(7) = 1
        adjncy(8) = 1
        adjncy(9) = 2
        adjncy(10)= 1

       !!!!!!!!!!!!!xadj contains the column pointers of a
        xadj(1) = 1
        xadj(2) = 5
        xadj(3) = 7
        xadj(4) = 8
        xadj(5) = 10
        xadj(6) = 11

        neqns = 5 ! number of equations
        iwsiz = 4*neqns  ! integer working storage for the minimum degree reordering
        nnza=10 ! number of nonzero entries of a
        !!!Save the adjncy vector before the minimum degree reordering
         do i=1,nnza
          temp(i)=adjncy(i)
         enddo
	 !minimum degree reordering
        call ordmmd(neqns,xadj,temp,invp,perm,iwsiz,iwork,&
                   nofsub,iflag)


        iwsiz = 7*neqns+3 ! integer working storage for the factorization subroutines

	!!!!!!!!!!!!!!!Initialization and determination of the supernodes

        call sfinit(neqns,nnza,xadj,adjncy,perm,invp,colcnt,&
                   nnzl,nsub,nsuper,snode,xsuper,iwsiz,&
                   iwork,iflag)


       !!!!Symbolic factorization

        call symfct(neqns,nnza,xadj,adjncy,perm,invp,colcnt,&
                   nsuper,xsuper,snode,nsub,xlindx,lindx,&
                   xlnz,iwsiz,iwork,iflag)

        cachsz=16   ! cache size in Kb
        call bfinit(neqns,nsuper,xsuper,snode,xlindx,lindx,&
           cachsz,tmpsiz,split)
        do i=1,neqns
          xadjf(i)=xadj(i)+i-1
          adjf(xadjf(i))=i
              do ij=xadj(i),xadj(i+1)-1
               adjf(ij+i)=adjncy(ij)
              enddo
        enddo
        xadjf(neqns+1)=xadj(neqns+1)+neqns

      ! Computation of the full matrix af by inserting the diagonal elements
      ! The output vectors xadjf and adjf contain the row indices and the column pointers
      ! of af.

        diag(1)=3.
        diag(2)=3.
        diag(3)=3.
        diag(4)=0
        diag(5)=0
        do i=1,neqns
          af(xadjf(i))=diag(i)
              do ij=xadj(i),xadj(i+1)-1
               af(ij+i)=a(ij)
              enddo
        enddo


        !Insert the numerical values in the data structures

        call inpnv(neqns,xadjf,adjf,af,perm,invp,nsuper,&
                  xsuper,xlindx,lindx,xlnz,lnz,offset)
        nsup=3 ! number of rows of the positive definite part
	!!!!!!!!!!!!block factorization of af
        call blkfct(neqns,nsuper,xsuper,snode,split,&
                xlindx,lindx,xlnz,lnz,iwsiz,iwork,tmpsiz,&
                tmpvec,iflag,mmpym,smxpym,perm,nsup)
        b(1)=3
        b(2)=-3
        b(3)=1
        b(4)=1
        b(5)=1
	 !!!!!!!!!!!!! solution of the sysyem !!!!!!!!!!
         do i=1,neqns
         sol(i)=b(perm(i))
         enddo
        call blkslv(nsuper,xsuper,xlindx,lindx,xlnz,lnz,sol,perm,nsup)
         do i=1,neqns
          b(perm(i))=sol(i)
         enddo
        print*,'sol=',(b(i),i=1,neqns)
        end program example