c    -*- fortran -*-
c    author:  travis e. oliphant, 2004
c
c    computational tools for sparse matrices in <_t>
c
c    add two matrices in compressed sparse column format
c        this assumes that the rowa indices are sorted for each column
c        and indexes are 0-based

      subroutine <_c>cscadd(n,a,rowa,ptra,nnzamax,b,rowb,ptrb,nnzbmax,
     $     c,rowc,ptrc,nnzcmax,ierr)
      <_t> a(0:nnzamax-1), b(0:nnzbmax-1)
      <_t> c(0:nnzcmax-1), val
      integer n, rowa(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
      integer ptra(0:n), ptrb(0:n), ptrc(0:n), nnzamax, nnzbmax
      integer ia, ib, ra, rb, j, nnzcmax
      integer ierr, nnzc
      
      ierr=0
      nnzc=0
      ia=ptra(0)
      ib=ptrb(0)

c  is there a need to initialize the output arrays?  assume no for now.
c      do j = 0, n
c         ptrc(j) = 0
c      end do
c      do j = 0, nnzcmax-1
c         rowc(j) = 0
c         c(j) = 0.0
c      end do

      do j = 0, n-1
         iamax = ptra(j+1)
         ibmax = ptrb(j+1)
         do while ((ia.lt.iamax).and.(ib.lt.ibmax))
            ra = rowa(ia)
            rb = rowb(ib)
            if (ra.eq.rb) then
c     this is a common element
               val = a(ia) + b(ib)
               ia = ia + 1
               ib = ib + 1
               if (val.eq.0.0) goto 10
               if (nnzc.ge.nnzcmax) goto 999
               c(nnzc) = val
               rowc(nnzc) = ra
            else if (ra .lt. rb) then
c     a has this but not b
               val = a(ia)
               ia = ia + 1 
               if (val.eq.0.0) goto 10 
               if (nnzc.ge.nnzcmax) goto 999
               c(nnzc) = val
               rowc(nnzc) = ra
            else
c     b has this but not a
               val = b(ib)
               ib = ib + 1
               if (val.eq.0.0) goto 10
               if (nnzc.ge.nnzcmax) goto 999
               c(nnzc) = val
               rowc(nnzc) = rb
            end if
            ptrc(j+1) = ptrc(j+1) + 1
            nnzc = nnzc + 1
 10         continue
         end do
         
         if (ia .eq. iamax) then
c     all finished with a for this column just copy the rest from b
            do while (ib .lt. ibmax)
               val = b(ib)
               rb = rowb(ib)
               ib = ib+1
               if (val.ne.0.0) then
                  if (nnzc.ge.nnzcmax) goto 999
                  c(nnzc) = val
                  rowc(nnzc) = rb
                  ptrc(j+1) = ptrc(j+1) + 1
                  nnzc = nnzc + 1
               end if
            end do
         else if (ib .eq. ibmax) then
c     all finished with b for this column just copy the rest from a
            do while (ia .lt. iamax)
               val = a(ia)
               ra = rowa(ia)
               ia = ia + 1
               if (val.ne.0.0) then
                  if (nnzc.ge.nnzcmax) goto 999
                  c(nnzc) = val
                  rowc(nnzc) = ra
                  ptrc(j+1) = ptrc(j+1) + 1
                  nnzc = nnzc + 1
               end if
            end do
         end if 
      end do

c   successful completion (fix the ptr array)
      cumsum = 0
      do k = 1, n
         cumsum = cumsum + ptrc(k)
         ptrc(k) = cumsum         
      end do
      return
      
      
 999  continue
      ierr = 1
      return
      end subroutine <_c>cscadd

c     element-by-element multiplication
c     can have at most the minimum of nnzamax and nnzbmax

      subroutine <_c>cscmul(n,a,rowa,ptra,nnzamax,b,rowb,ptrb,nnzbmax,
     $     c,rowc,ptrc,nnzcmax,ierr)
      <_t> a(0:nnzamax-1), b(0:nnzbmax-1)
      <_t> c(0:nnzcmax-1), val
      integer n, rowa(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
      integer ptra(0:n), ptrb(0:n), ptrc(0:n), nnzamax, nnzbmax
      integer ia, ib, ra, rb, j, nnzcmax
      integer ierr, k, cumsum
      
      ierr=0
      nnzc=0
      ia=ptra(0)
      ib=ptrb(0)
      do j = 0, n-1
         iamax = ptra(j+1)
         ibmax = ptrb(j+1)
         do while ((ia.lt.iamax).and.(ib.lt.ibmax))
            ra = rowa(ia)
            rb = rowb(ib)
            if (ra.eq.rb) then
c     this is a common element
               val = a(ia)*b(ib)
               ia = ia + 1
               ib = ib + 1
               if (val.eq.0.0) goto 10
               if (nnzc.ge.nnzcmax) goto 999
               c(nnzc) = val
               rowc(nnzc) = ra
               ptrc(j+1) = ptrc(j+1) + 1
               nnzc = nnzc + 1
            else if (ra .lt. rb) then
c     a has this but not b
               ia = ia + 1
            else
c     b has this but not a
               ib = ib + 1
            end if
 10         continue
         end do
      end do
      
c   successful completion (fix the ptr array)
      cumsum = 0
      do k = 1, n
         cumsum = cumsum + ptrc(k)
         ptrc(k) = cumsum         
      end do

      return
      
 999  continue
      ierr = 1
      return
      end subroutine <_c>cscmul


c  matrix-vector multiplication

      subroutine <_c>cscmux(a,rowa,ptra,nnzamax,ncol,x,mrow,y)
      <_t> a(0:nnzamax-1), x(0:ncol-1), y(0:mrow-1)
      integer rowa(0:nnzamax-1), ptra(0:ncol)
      integer nnzamax, mrow, ncol
      integer i, j, ia, ra
      
      do i = 0, mrow-1
         y(i) = 0.0
      end do

      do j = 0, ncol-1
         do ia = ptra(j), ptra(j+1)-1
            ra = rowa(ia)
            y(ra) = y(ra) + a(ia)*x(j)
         end do
      end do

      return
      end subroutine <_c>cscmux


      subroutine <_c>csrmux(a,cola,ptra,nnzamax,ncol,x,mrow,y)
      <_t> a(0:nnzamax-1),x(0:ncol-1),y(0:mrow-1)
      integer cola(0:nnzamax-1), ptra(0:mrow)
      integer nnzamax, mrow, ncol
      integer i, ja, ca
      
      do i = 0, mrow-1
         y(i) = 0.0
         do ja = ptra(i), ptra(i+1)-1
            ca = cola(ja)
            y(i) = y(i) + a(ja)*x(ca)           
         end do
      end do

      return
      end subroutine <_c>csrmux


c     matrix multiplication
c             c = a * b
c
c     where c is mxn csc matrix 
c           a is mxk csc matrix
c           b is kxn csr matrix
c
c     irow and kcol give position to start 
c     nnzc gives current last element in output array (initialize to 0)
c     intended to recall so that calculation can continue 
c     if memory ran-out at last attempt



      subroutine <_c>cscmucsr(m,k,n, a,rowa,ptra,nnzamax, b,colb,ptrb,
     $     nnzbmax, c,rowc,ptrc,nnzcmax, irow,kcol,ierr)

      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1)
      integer rowa(0:nnzamax-1), colb(0:nnzbmax-1), rowc(0:nnzcmax-1)
      integer ptra(0:k), ptrb(0:k), ptrc(0:n)
      integer nnzamax, nnzbmax, nnzcmax
      integer m, k, n, irow, kcol, ierr, cumsum

      integer kk, ii, jjb, jb, cb, ia, nnzc
      <_t> val


      nnzc = ierr
      ierr = 0
      do kk = kcol, n-1
         do ii = irow, m-1
            if (nnzc.ge.nnzcmax) goto 999
            irow = 0
            val = 0.0
c   loop through the column array of b using the ptrb array
c        so that we know which row we are on.
            do jjb = 0, k-1
               do jb = ptrb(jjb), ptrb(jjb+1)-1
                  cb = colb(jb)
                  if (cb.eq.kk) then
c     see if aij is nonzero and if so add to val
                     do ia = ptra(jjb), ptra(jjb+1)-1
                        if (rowa(ia).eq.ii) then
                           val = val + a(ia)*b(jb)
                        end if
                     end do
                  end if
               end do
            end do
            if (val.ne.0.0) then
c     there is a non-zero value for this value of i and k
               c(nnzc) = val
               rowc(nnzc) = ii
               ptrc(kk+1) = ptrc(kk+1) + 1
               nnzc = nnzc+1               
            end if
         end do         
      end do

c   successful completion (fix the ptr array)
      cumsum = 0
      do kk = 1, n
         cumsum = cumsum + ptrc(kk)
         ptrc(kk) = cumsum         
      end do
      return

      return

 999  continue
      kcol = kk
      irow = ii
      ierr = nnzc

      return
      end subroutine <_c>cscmucsr



c     matrix multiplication
c             c = a * b
c
c     where c is mxn csc matrix 
c           a is mxk csr matrix
c           b is kxn csc matrix
c
c     irow and jcol give position to start 
c     bnum gives current size of output array
c     intended to recall so that calculation can continue 
c     where it left off if memory runs-out during calculation
c
      subroutine <_c>csrmucsc(m,n,a,cola,ptra,nnzamax,b,rowb,ptrb,
     $     nnzbmax,c,rowc,ptrc,nnzcmax,irow,kcol,ierr)

      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1)
      integer cola(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
      integer ptra(0:m), ptrb(0:n), ptrc(0:n)
      integer nnzamax, nnzbmax, nnzcmax
      integer m, n, irow, kcol, ierr, cumsum

      integer kk, ii, nnzc
      <_t> val

      nnzc = ierr
      ierr = 0
      do kk = kcol, n-1
         do ii = irow, m-1
            if (nnzc.ge.nnzcmax) goto 999
            irow = 0
            val = 0.0
c   loop through the non-zero rows of b
            do ib = ptrb(kk), ptrb(kk+1)-1
               rb = rowb(ib)
               do ja = ptra(ii), ptra(ii+1)-1
                  if (cola(ja).eq.rb) then
                     val = val + a(ja)*b(ib)
                  end if
               end do
            end do
            if (val.ne.0.0) then
c     there is a non-zero value for this value of i and k
               c(nnzc) = val
               rowc(nnzc) = ii
               ptrc(kk+1) = ptrc(kk+1) + 1
               nnzc = nnzc+1 
            end if
         end do
      end do

c   successful completion (fix the ptr array)
      cumsum = 0
      do kk = 1, n
         cumsum = cumsum + ptrc(kk)
         ptrc(kk) = cumsum         
      end do
      return


      return

 999  continue
      kcol = kk
      irow = ii
      ierr = nnzc

      return
      end subroutine <_c>csrmucsc


c     matrix-matrix multiplication
c
c          c = a * b    where  a, b, and c are 
c                       compressed sparse column matrices
c
c     or it computes
c
c          c = b * a    where  b, a, and c are 
c                       compressed sparse row matrices
c   
c     intentended for re-entry if nnzcmax is not big enough

      subroutine <_c>cscmucsc(m,k,n,a,rowa,ptra,nnzamax,b,rowb,ptrb,
     $     nnzbmax,c,rowc,ptrc,nnzcmax,irow,kcol,ierr)

      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1)
      integer rowa(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
      integer ptra(0:k), ptrb(0:n), ptrc(0:n)
      integer nnzamax, nnzbmax, nnzcmax
      integer m, k, n, irow, kcol, ierr

      integer kk, ii, jb, ia, ra, nnzc, cumsum
      <_t> val
      
      nnzc = ierr
      ierr = 0
      do kk = kcol, n-1
         do ii = irow, m-1
            if (nnzc.ge.nnzcmax) goto 999
c   reset the irow to 0
            irow = 0
            val = 0.0
c   loop through the row array of b for this column
            do jb = ptrb(kk), ptrb(kk+1)-1
               ra = rowb(jb)
               do ia = ptra(ra), ptra(ra+1)-1
c     see if aij is nonzero and if so add to val
                  if (rowa(ia).eq.ii) then
                        val = val + a(ia)*b(jb)
                  end if
               end do
            end do
            if (val.ne.0.0) then
c     there is a non-zero value for this value of i and k
               c(nnzc) = val
               rowc(nnzc) = ii
               ptrc(kk+1) = ptrc(kk+1) + 1
               nnzc = nnzc+1               
            end if
         end do         
      end do

c   successful completion (fix the ptr array)
      cumsum = 0
      do kk = 1, n
         cumsum = cumsum + ptrc(kk)
         ptrc(kk) = cumsum         
      end do
      return
      
 999  continue
      kcol = kk
      irow = ii
      ierr = nnzc
      
      return
      end subroutine <_c>cscmucsc
     

