c    -*- fortran -*-
c    convert csc matrix a to csr matrix b
c    convert csr matrix a to csc matrix b
c    transpose csc matrix
c    transpose csr matrix
c    all depends on how you interpret the results
         
      subroutine <_c>transp(m,n,a,rowa,ptra,nnzamax,b,colb,ptrb)
      <_t> a(0:nnzamax-1), b(0:nnzamax-1)
      integer rowa(0:nnzamax-1), colb(0:nnzamax-1)
      integer ptra(0:n), ptrb(0:m)
      integer j, bnum, cola, ia, ra

      ptrb(0) = 0
      bnum = 0
      do j = 0, m-1
         do cola = 0, n-1
            do ia = ptra(cola), ptra(cola+1)-1
               ra = rowa(ia)
               if (ra.eq.j) then
                  b(bnum) = a(ia)
                  colb(bnum) = cola
                  bnum = bnum+1
               end if
            end do
         end do
         ptrb(j+1) = bnum
      end do

      return
      end subroutine <_c>transp

c      returns the index into the data array
c      and the value
      subroutine <_c>cscgetel(a,rowa,ptra,nnzamax,n,row,col,ind,val)
      
      <_t> a(0:nnzamax-1), val
      integer rowa(0:nnzamax-1), ptra(0:n)
      integer nnzamax, row, col, ind

      integer ia

      ind = -1
      val = 0.0
      do ia = ptra(col), ptra(col+1)-1
         if (rowa(ia).eq.row) then
            ind = ia
            val = a(ia)
            goto 10
         end if
      end do

 10   continue
      return 

      end subroutine <_c>cscgetel

      
c      set an element into the data array (assumes there is room)
c
      subroutine <_c>cscsetel(a,rowa,ptra,nnzamax,n,row,col,val)

      <_t> a(0:nnzamax-1), val
      integer rowa(0:nnzamax-1), ptra(0:n)
      integer nnza, n, row, col, nnzamax

      integer ra, newia

      nnza = ptra(n)
      newia = ptra(col)
      do ia = newia, ptra(col+1)-1
         ra = rowa(ia)
         if (ra.eq.row) then
c     just replace and be done
            a(ia) = val
            goto 10
         else if (ra.gt.row) then 
            newia = ia
            goto 5
         end if
      end do

c     we are past where it should be stored
c         (assumes sorted) so make room for this new element
c     here so that it will run if this is the first entry in the
c         column
 5    continue
      ia = newia
      do iia = nnza, ia+1, -1
         a(iia) = a(iia-1)
         rowa(iia) = rowa(iia-1)
      end do
      a(ia) = val
      rowa(ia) = row
      do iia = col+1, n
         ptra(iia) = ptra(iia)+1
      end do

      return

 10   continue
      
      return
      end subroutine <_c>cscsetel

c     assumes sorted by column and then row 

      subroutine <_c>cootocsc(n,vals,row,col,nnz,a,rowa,ptra,nnzamax,
     $     ierr)
      
      <_t> vals(0:nnz-1), a(0:nnzamax-1)
      integer rowa(0:nnzamax-1), ptra(0:n)
      integer row(0:nnz-1), col(0:nnz-1)
      integer n, nnzamax, nnz, ierr

      <_t> val
      integer j, k, cumsum

      if (nnz.gt.nnzamax) goto 999
      ierr=0
      nnza=0
      do k = 0, n-1
         ptra(k) = 0
      end do
      do ia = 0, nnzamax-1
         a(ia) = 0.0
         rowa(ia) = 0
      end do      
      do k = 0, nnz-1
         j = col(k)
         val = vals(k)
         if (val.ne.0.0) then
            a(nnza) = val
            rowa(nnza) = row(k)
            ptra(j+1) = ptra(j+1) + 1
            nnza = nnza + 1
         end if
      end do

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

      return

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

      subroutine <_c>csctocoo(n, vals, row, col, a, rowa, ptra, nnzamax)
      <_t> vals(0:nnzamax-1), a(0:nnzamax-1)
      integer rowa(0:nnzamax-1), ptra(0:n)
      integer row(0:nnzamax-1), col(0:nnzamax-1)
      integer n, nnzamax

      integer i, j, k

      j = 0
      do k = 0, n-1
         do i = ptra(k), ptra(k+1)-1
            row(j) = rowa(i)
            col(j) = k
            vals(j) = a(i)
            j = j + 1
         end do
      end do

      return
      end subroutine <_c>csctocoo

c     intended for re-entry in case a is too small

      subroutine <_c>fulltocsc(m,n,fulla,a,rowa,ptra,nnzamax,
     $     irow,jcol,ierr)
      
      <_t> fulla(0:m-1,0:n-1), a(0:nnzamax-1)
      integer rowa(0:nnzamax-1), ptra(0:n)
      integer m, n, nnzamax, irow, jcol, ierr

      <_t> val
      integer nnza, i, j, k, cumsum

      ierr=0
      nnza = ierr
      do j = jcol, n-1
         do i = irow, m-1
            val = fulla(i,j)
            if (val.ne.0.0) then
               if (nnza.ge.nnzamax) goto 999
               a(nnza) = val
               rowa(nnza) = i
               ptra(j+1) = ptra(j+1) + 1
               nnza = nnza + 1
            end if
         end do
      end do      

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

 999  continue
      ierr=nnza
      irow = i
      jcol = j
      return
      end subroutine <_c>fulltocsc


      subroutine <_c>csctofull(m, n, fulla, a, rowa, ptra, nnzamax)

      <_t> fulla(0:m-1, 0:n-1), a(0:nnzamax-1)
      integer rowa(0:nnzamax-1), ptra(0:n)
      integer m, n, nnzamax

      integer ia, k
c      integer i, j

c      do j = 0, n-1
c         do i = 0, m-1
c            fulla(i,j) = 0.0
c         end do
c      end do

      do k = 0, n-1
         do ia = ptra(k), ptra(k+1)-1
            fulla(rowa(ia),k) = a(ia)
         end do
      end do
      
      return
      end subroutine <_c>csctofull


c  extract a sub-matrix from a sparse matrix
c     c = a(ibeg:iend, jbeg:jend)  (inclusive of end-points)
c
c     intended for re-entry if nnzcmax is not big enough 
c      irow and jcol should initially be 0 and 0

      subroutine <_c>cscextract(n,a,rowa,ptra,nnzamax,ibeg,iend,jbeg,
     $     jend, c, rowc, ptrc, nnzcmax, irow, jcol, ierr)

      <_t> a(0:nnzamax-1), c(0:nnzcmax-1)
      integer rowa(0:nnzamax-1), rowc(0:nnzamax-1)
      integer ptra(0:n), ptrc(0:jend-jbeg+1)
      integer ibeg, iend, jbeg, jend
      integer nnzamax, nnzcmax, irow, jcol, ierr

      integer nnza, nc, j, iabeg, ia, k, cumsum

      nc = jend-jbeg+1

      nnza = ptra(n)
      nnzc = ierr

      if (jcol .lt. jbeg) jcol = jbeg

      do j = jcol, jend
         jc = j - jbeg
c   look through row indices in this column, copy all those that are
c      in required range.
         iabeg = max(irow,ptra(j))
         do ia = iabeg, ptra(j+1)-1
            ra = rowa(ia)
            if ((ra.le.iend).and.(ra.ge.ibeg)) then
               if (nnzc.ge.nnzcmax) goto 999
               c(nnzc) = a(ia)
               rowc(nnzc) = ra-ibeg
               ptrc(jc+1) = ptrc(jc+1) + 1
               nnzc = nnzc + 1
            end if
         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


      return

 999  continue
      ierr = nnzc
      irow = ia
      jcol = j
      return
      end subroutine <_c>cscextract


      
      subroutine <_c>diatocsc(m,n,diags,numdia,diasize,offsets,
     $     a,rowa,ptra,nzmax, ierr)
      
      integer n, numdia, diasize, nzmax, ierr
      <_t> diags(0:numdia-1, 0:diasize-1), a(0:nzmax-1)
      <_t> val
      integer offsets(0:numdia-1), rowa(0:nzmax-1)
      integer ptra(0:n)

      integer col, nnza, jj, row, idiag, ia

      nnza = 0
      ierr = 0
      write (*,*) m, n, numdia, diasize, nzmax
      do col = 0, n-1
c   Loop through the offsets 
         do jj = 0, numdia-1
            row = col - offsets(jj)
            if ((row.ge.0).and.(row.lt.m)) then
               idiag = min(row, col)
               val = diags(jj, idiag)
               if (val.ne.0.0) then
                  if (nnza.ge.nzmax) goto 999
c   find place to insert this row (inserted in sorted order)
                  ia = ptra(col)
 10               if ((ia.lt.ptra(col+1)).and.(rowa(ia).lt.row)) then
                     ia = ia + 1
                     goto 10
                  end if
                  do iia = nnza, ia+1, -1
                     a(iia) = a(iia-1)
                     rowa(iia) = rowa(iia-1)
                  end do
                  a(ia) = val
                  rowa(ia) = row
                  nnza = nnza + 1
               end if
            end if
         end do
         ptra(col+1) = nnza
      end do

      return 

 999  continue
      ierr = 1
      return
      

      end subroutine <_c>diatocsc
