subroutine genmatrix(L,mm0,kr,mdw,mnt,k6,k7,
     &     n,nz,big,iwork,dwork,icomm,rcomm,bnd)
      logical bnd
      implicit double precision (a-h,o-z)
      dimension dwork(mdw),iwork(kr)
      dimension rcomm(*),icomm(*)

      k1 = 1
      k2 = k1+L
      k3 = k2+L
      k4 = k3+L*L
      call wpack(L,iwork(k1),needs,iwork(k4))
      k5 = k4+needs
      k6 = k5+MM0+1
      kr = kr - k6+1

      iguess = (mnt+needs+3)
      mypart = kr/iguess
      k7  = k6 +mypart*mnt
      k8  = k7 +mypart
      k9  = k8 +mypart*needs
      k10 = k9 +mypart
      k11 = k10 +mypart
c     j2  = mdw -mypart+1
      j2  = mdw
      nzmax = mypart*mnt
      if (nzmax.ge.j2) nzmax = j2


      call statematx(L,iwork(k1),iwork(k2),iwork(k3),iwork(k4),
     * iwork(k5),iwork(k6),iwork(k7),iwork(k8),iwork(k9),iwork(k10),
     * n,nz,big,mm0,mypart,nzmax,needs,dwork,icomm,rcomm,bnd)

      return
      end



      subroutine wpack(L,max,needs,ints)
      dimension max(L),ints(L)

      do i=1,L
	 ints(i) = -1
      end do

      needs = 1
      ntot = max(1)+1
      k1 = 1
      do while (2**k1.lt.ntot)
	  k1 = k1+1
      end do

      do j=2,L
	 k2 = 1
	 do while (2**k2.lt.max(j)+1)
	    k2 = k2+1
         end do
	 if (k1+k2.le.28) then
            ntot = ntot*(max(j)+1)
	    do while (2**k1.lt.ntot)
	       k1 = k1+1
            end do
	 else
	    ints(needs) = j-1
	    needs = needs+1
	    ntot = max(j)+1
	    k1 = k2
         endif
      end do
      ints(needs) = L

      if (needs.gt.1) then
	 do j = needs,2,-1
	    ints(j) = ints(j)-ints(j-1)
         end do
      end if

      return
      end


      

      subroutine statematx(L,max,istate,tran,ints,
     &   mhash,ja,ia,list,nhash,into,
     &   n,nz,big,mm0,nmax,nzmax,needs,a,icomm,rcomm,bnd)

      implicit double precision (a-h,o-z)
      integer tran
      logical norate,bnd
      dimension max(L),istate(L),tran(L*L),ints(needs),list(nmax*needs)
      dimension a(nzmax),ja(nzmax),ia(nmax)
      dimension mhash(0:mm0),nhash(nmax),into(nmax),rcomm(*),icomm(*)
      dimension prob(1000)

      Ls1 = 1
      Ls2 = 2
      norate = .false. 
      call rate(Ls1,Ls2,istate,into,L,r,icomm,rcomm)
      if (Ls1.eq.-1 .and. Ls2.eq.-1) norate = .true.
    
      icmax = nmax/L
      ist = 0
      N1=1
      BIG=0.0D0
      MARKER=1
      LIMIT=1
      IA(1)=1
      next = needs
      mm1 = mm0+1
      do i = 0,mm0
	 mhash(i) = -1
      end do
      do i = 1,nmax
         nhash(i) = -1 
      end do

 
      call getkey(L,max,istate,list,needs,ints,key) 
 
      JX=MOD(key,MM1)
      MHASH(JX)=1

  100 IF(MARKER.GT.LIMIT)GO TO 90
 
      N4=N1
      A(N4)=0.0D0
      ja(n4)=marker
    
      ist = ist+1
      nextst = list(ist)
      if (ints(1).eq.1) goto 10
      do i=1,ints(1)-1
         istate(i) = mod(nextst,max(i)+1)
         nextst = (nextst-istate(i))/(max(i)+1)
      end do
   10 istate(ints(1)) = nextst

      if (needs.gt.1) then
	 j2 = ints(1)
	 do j1 = 2,needs
            ist = ist+1
	    nextst = list(ist)
            if (ints(j1).eq.1) goto 12
	    do i = 1,ints(j1)-1
	       j2 = j2+1
	       istate(j2) = mod(nextst,max(j2)+1)
               nextst = (nextst-istate(j2))/(max(j2)+1)
            end do
   12       j2 = j2+1
	    istate(j2) = nextst
         end do
      end if

 
      if(marker.eq.100)write(*,910)marker
      if(marker.eq.200)write(*,911)marker
      if(marker.eq.500)write(*,911)marker
      if (marker.ge.100000) then
	 if (mod(marker,100000).eq.0)write(*,911)marker
	 else
	 if (marker.ge.10000) then
	    if (mod(marker,10000).eq.0)write(*,911)marker
	    else
	    if (marker.ge.1000) then
	       if (mod(marker,1000).eq.0)write(*,911)marker
            endif
         endif
      endif

      if ( norate ) then

         r = 1.0d0
         do i = 1,L
            into(i) = istate(i)
         end do
         call instant(l,into,miw,prob,ic,icomm,rcomm)
         call stappend(limit,ic,L,max,into,list,needs,ints,mm0,
     &        mhash,nhash,a,ja,ia,n1,n4,r,prob,nmax,nzmax,next)

      else 

         if (bnd) then
         do 70 i=1,L
            if (istate(i).eq.0) go to 70
            do 60 j=1,L
               if (istate(j).eq.max(j).or.tran((i-1)*l+j).eq.0) goto 60
               call rate(i,j,istate,into,l,r,icomm,rcomm)
               if(r.le.0.0d0)go to 60
               call instant(l,into,miw,prob,ic,icomm,rcomm)
               call stappend(limit,ic,L,max,into,list,needs,ints,mm0,
     &              mhash,nhash,a,ja,ia,n1,n4,r,prob,nmax,nzmax,next)
   60       continue
   70    continue
         else
         do 71 i=1,L
            if (istate(i).lt.0) go to 71
            do 61 j=1,L
               if (istate(j).gt.max(j).or.tran((i-1)*l+j).eq.0) goto 61
               call rate(i,j,istate,into,l,r,icomm,rcomm)
               if(r.le.0.0d0)go to 61
               call instant(l,into,miw,prob,ic,icomm,rcomm)
               call stappend(limit,ic,L,max,into,list,needs,ints,mm0,
     &              mhash,nhash,a,ja,ia,n1,n4,r,prob,nmax,nzmax,next)
   61       continue
   71    continue
         end if 

      end if 
 
   80 MARKER=MARKER+1
      IF(BIG.GT.A(N4))BIG=A(N4)
      N1=N1+1
      IF(N1.GT.NZMAX .OR. MARKER.GT.NMAX)GO TO 999
      IA(MARKER)=N1
      GO TO 100
   
   90 N=LIMIT
      IF(N.GE.100)WRITE(*,912)N
      NZ=N1-1
      RETURN

  999 print *,'   ERROR:  Insufficient memory allocated.'
      stop
 
  910 FORMAT(/,'  Currently at State:',i7)
  911 FORMAT('                     ',i7)
  912 FORMAT('          Last State:',i7)
 
      END
 

      subroutine stappend(limit,ic,L,max,into,list,needs,ints,mm0,
     &           mhash,nhash,a,ja,ia,n1,n4,r,prob,nmax,nzmax,next)
      implicit double precision (a-h,o-z)
      dimension max(L),ints(needs),list(nmax*needs)
      dimension a(nzmax),ja(nzmax),ia(nmax)
      dimension mhash(0:mm0),nhash(nmax),into(nmax)
      dimension prob(nmax)
      dimension ip(10)
      dimension label(10)
      data ip / 1, 1531, 97, 47, 23, 13, 11, 7, 5, 3 / 

      mm1 = mm0+1
      r1 = r
      icc=1
      iccn=0

   15 call getkey(L,max,into,label,needs,ints,key) 

      nx = mod(key,mm1)
      ir = mhash(nx)
      if (ir.eq.-1) then
         mhash(nx) = limit+1
         goto 32
      endif
      do k = 1,needs
         if (list( (ir-1)*needs+k ).ne.label(k)) goto 30
      end do
      goto 32
   30 if (nhash(ir).eq.-1) then
         nhash(ir) = limit+1
         ir = -1
         goto 32
      end if
      ir = nhash(ir)
      do k = 1,needs
         if (list( (ir-1)*needs+k ).ne.label(k)) goto 30
      end do
   32 if (ir.ne.-1) goto 40
 
      limit = limit+1
      do k = 1,needs
         next = next+1
         list(next) = label(k)
      end do
      ir = limit
 
   40 if (ic.ge.1) r = r1*prob(icc)
         a(n4)=a(n4)-r
         do k=n4,n1
            if (ja(k).eq.ir) then 
               a(k)=a(k)+r
               go to 50 
         endif
      end do 
      n1=n1+1
      a(n1)=r
      ja(n1)=ir
 
   50 if(icc.ge.ic)return
      icc=icc+1
      iccn=iccn+L
      do k=1,L
         into(k)=into(k+iccn)
      end do
      go to 15

  906 format(d10.4,2x,i7,4x,10i6)
      end 


      subroutine searchlist(L,max,istate,list,label,n,needs,ints,
     &                      mhash,nhash,mm0,ir)
      implicit double precision (a-h,o-z)
      dimension max(L),istate(L),list(needs*n),ints(needs)
      dimension mhash(0:mm0),nhash(n),label(needs)

      call getkey(L,max,istate,label,needs,ints,key)

      mm1 = mm0+1 
      nx = mod(key,mm1)
      ir = mhash(nx)
      if (ir.eq.-1)return
      do k = 1,needs
         if (list( (ir-1)*needs+k ).ne.label(k)) goto 10
      end do
      return
   10 if (nhash(ir).eq.-1) then
	 ir = -1
         return
      end if
      ir = nhash(ir)
      do k = 1,needs
         if (list( (ir-1)*needs+k ).ne.label(k)) goto 10
      end do

      return
      end 


      subroutine getkey(L,max,istate,label,needs,ints,key)
      implicit double precision (a-h,o-z)
      dimension max(L),istate(L),ints(needs),label(needs)
      dimension ip(10)
      data ip / 1, 1531, 97, 47, 23, 13, 11, 7, 5, 3 / 

      key = 0
      do k= ints(1),2,-1
         key = (key + istate(k))*(max(k-1)+1)
      end do
      key = key + istate(1)
      label(1) = key

      if (needs.gt.1) then
         last = ints(1)
         do j1 = 2,needs
            index = last+1
 	    last = last + ints(j1)
 	    label(j1) = 0
 	    if (ints(j1).eq.1) goto 5
 	    do j2 = last,index+1,-1
               label(j1) = (label(j1) + istate(j2))*(max(j2-1)+1)
            end do
    5       label(j1) = label(j1) + istate(index)
            key = key + label(j1)*ip(j1)
            if (key.lt.0) key = -key
         end do
      end if

      return
      end
    
      
      subroutine writeToFile(n,nz,aa,ja,ia,io_unit)
      implicit double precision (a-h,o-z)
      dimension aa(nz),ja(nz),ia(n+1)

      iend = 0
      next = 0
      write(io_unit,*)
      write(io_unit,*)n,n,nz,'        Number of rows, columns, nonzeros'
      do j = 1,n
	 istart = iend+1
	 iend = ia(j+1)-1
	 next = next+1
	 do k = istart,iend
 	    write(io_unit,901) j,ja(k),aa(k)
         end do
      end do
      write (io_unit,*) ' -1  -1  -1'
      write(io_unit,*)
  901 format(i6,4x,i6,8x,f24.14)
      return
      end