MODULE sparse_matrix

! module managing sparse matrices with a symmetrical profile

! Fortran interface :

!  INTEGER(C_LONG) :: i,j,n,noper,nlist,list(*),tj(*)
!  REAL(C_DOUBLE)  :: val,vector(*),result(*),a(nlist,nlist),tval(*)

!  CALL sparse_matrix_init(matrix,n[,mode,eps,noper])
!  CALL sparse_matrix_add(matrix,i,j,val)
!  CALL sparse_matrix_add_vect(matrix,i,n,tj,tval)
!  CALL sparse_matrix_add_matrix(m1,m2)
!  CALL sparse_matrix_omp_add_matrix(m1,m2)
!  CALL sparse_matrix_put(matrix,i,j,val)
!  CALL sparse_matrix_get(matrix,i,j,val)
!  CALL sparse_matrix_reset(matrix)
!  CALL sparse_matrix_optimize(matrix)
!  CALL sparse_matrix_lu_factorization(matrix)
!  CALL sparse_matrix_omp_lu_factorization(matrix)
!  CALL sparse_matrix_lu_vector(matrix,vector)
!  CALL sparse_matrix_vector(matrix,vector,result)
!  CALL sparse_matrix_lu_extract(matrix,nlist,list,a)

! C interface :

!  void*matrix;
!  int i,j,noper,nlist,*list;
!  double val,eps;
!  double*vector,*result,*a;

!  matrix=sparse_matrix_init(matrix,n,mode,eps,noper);
!  sparse_matrix_add(matrix,i,j,val);
!  sparse_matrix_add_vect(matrix,i,n,tj,tval);
!  sparse_matrix_add_matrix(m1,m2);
!  sparse_matrix_omp_add_matrix(m1,m2);
!  sparse_matrix_put(matrix,i,j,val);
!  sparse_matrix_reset(matrix);
!  sparse_matrix_optimize(matrix);
!  sparse_matrix_lu_factorization(matrix);
!  sparse_matrix_omp_lu_factorization(matrix);
!  sparse_matrix_lu_vector(matrix,vector);
!  sparse_matrix_vector(matrix,vector,result);
!  sparse_matrix_lu_extract(matrix,nlist,list,a);

USE, intrinsic :: iso_c_binding

IMPLICIT NONE

PRIVATE

PUBLIC :: sparse_matrix_type
PUBLIC :: fast_type ! for stupid compilers only
PUBLIC :: sparse_matrix_init,sparse_matrix_put,sparse_matrix_add,sparse_matrix_reset,sparse_matrix_print,&
          sparse_matrix_optimize,sparse_matrix_lu_factorization,sparse_matrix_lu_vector,sparse_matrix_lu_extract,&
          sparse_matrix_vector,sparse_matrix_get,sparse_matrix_add_matrix,&
          sparse_matrix_omp_lu_factorization, sparse_matrix_omp_add_matrix
          
INTEGER,PARAMETER :: mterm=4 ! default number of terms when putting a first value in an empty row

TYPE row_type                                ! description of a matrix row
  INTEGER(C_LONG) :: nterm=0                  ! number of values registered in the row
  INTEGER(C_LONG) :: degree=0                 ! the degree of the row (number of terme on the right of the diagonal)
  INTEGER(C_LONG),ALLOCATABLE :: colind(:)    ! column indexes
  REAL(C_DOUBLE) ,ALLOCATABLE :: values(:)    ! values
END TYPE
  
TYPE fast_type                               ! the structure used to build up a sparse matrix when
  TYPE(row_type),ALLOCATABLE :: row(:)       ! its structure is changing
END TYPE

TYPE sparse_matrix_type
  INTEGER(C_LONG)             :: state=-1     ! status of the matrix
  INTEGER(C_LONG)             :: n=0
  INTEGER(C_LONG)             :: nterm        ! number of terms
  LOGICAL                     :: lsym=.FALSE. ! symmetric flag
  INTEGER(C_LONG)             :: lorder=0     ! ordering mode : 0 (normal), 1 (no ordering), -1 (imposed ordering)
  REAL(C_DOUBLE) ,ALLOCATABLE :: values(:)    ! jacobian terms outside the current level
  INTEGER(C_LONG),ALLOCATABLE :: order(:)     ! best ordering        new_vector(i)=old_vector(order (i))  <=> iold = order (inew)
  INTEGER(C_LONG),ALLOCATABLE :: permut(:)    ! permut(order(i))=i   old_vector(i)=new_vector(permut(i))  <=> inew = permut(iold)
  INTEGER(C_LONG),ALLOCATABLE :: rowptr(:)    ! starting of rows/colowns (the profile is symmetrical)
  INTEGER(C_LONG),ALLOCATABLE :: colind(:)    ! column/row indexes
  INTEGER(C_LONG),ALLOCATABLE :: diagonal(:)  ! position of diagonal terms
  INTEGER(C_LONG),ALLOCATABLE :: lo(:)        ! list of operations for LU factorization
  TYPE(fast_type)             :: fast         ! the form used when modif > 0
  REAL(C_DOUBLE)              :: small=0      ! values less than that are not added
  INTEGER(C_LONG)             :: nomax=0      ! maximum number of registered operations
END TYPE

CONTAINS

SUBROUTINE sparse_matrix_init(matrix,n,mode,small,nomax)
  TYPE(sparse_matrix_type),INTENT(out) :: matrix
  INTEGER(C_LONG)          ,INTENT(in)  :: n
  INTEGER(C_LONG)          ,INTENT(in),OPTIONAL :: mode  ! 1 => symmetrical profile
  REAL(C_DOUBLE)          ,INTENT(in),OPTIONAL :: small ! smallest value taken into account by sparse_matrix_add
  INTEGER(C_LONG)          ,INTENT(in),OPTIONAL :: nomax ! maximum number of registered operations
  matrix%n=n
  matrix%state=1
  ALLOCATE(matrix%diagonal(n),matrix%rowptr(n+1),matrix%order(n),matrix%permut(n))
  matrix%diagonal=0
  matrix%rowptr=0
  ALLOCATE(matrix%fast%row(n))
  IF(PRESENT(small)) matrix%small=small
  IF(PRESENT(nomax)) matrix%nomax=nomax
  IF(PRESENT(mode)) THEN
    IF(mode == 1) matrix%lsym=.TRUE.
  ENDIF
END SUBROUTINE

SUBROUTINE sparse_matrix_reset(matrix)
  TYPE(sparse_matrix_type),INTENT(inout) :: matrix
  INTEGER(C_LONG) :: i
  IF(matrix%state == 1) THEN ! only the fast part is used
    DO i=1,matrix%n
      IF(matrix%fast%row(i)%nterm > 0) &
      matrix%fast%row(i)%values=0
    ENDDO
  ELSE
    matrix%state =0 ! a normal matrix and not a LU
    matrix%values=0
  ENDIF
END SUBROUTINE

SUBROUTINE sparse_matrix_put(matrix,irow,jcol,val)

  TYPE(sparse_matrix_type),INTENT(inout),target :: matrix
  INTEGER(C_LONG)          ,INTENT(in)           :: irow,jcol
  REAL(C_DOUBLE)          ,INTENT(in)           :: val
  
  TYPE(row_type),POINTER :: row
  INTEGER(C_LONG)         :: k,j
  
  IF(matrix%state == 1) THEN ! the fast structure is used
  
    row => matrix%fast%row(irow)
    
    DO j=1,row%nterm ! terms are not ordered in the fast structure
      IF(row%colind(j)==jcol) THEN
        row%values(j)=val
        RETURN
      ENDIF
    ENDDO
    
    CALL addterm(row,jcol,val)
    
    ! storing also the symmetrical term to keep a symmetrical profile
    
    IF(matrix%lsym) THEN
      IF(jcol /= irow) THEN
        row => matrix%fast%row(jcol)
        CALL addterm(row,irow,0._c_double)
      ENDIF
    ENDIF
    
    RETURN
    
  ELSE
  
    k=locate(matrix,irow,jcol)
    
    IF(k == 0) THEN ! A change in the structure is observed => the fast structure must be used
    
      WRITE(*,*) 'sparse_matrix_put : dangerous change of matrix state'
      matrix%state=1
      CALL create_fast(matrix)
      row => matrix%fast%row(irow)
      CALL addterm(row,jcol,val)
    
      ! storing also the symmetrical term to keep a symmetrical profile
    
      IF(matrix%lsym) THEN
        IF(jcol /= irow) THEN
          row => matrix%fast%row(jcol)
          CALL addterm(row,irow,0._c_double)
        ENDIF
      ENDIF
      
      RETURN
      
    ELSE
    
      matrix%values(k)=val
      
    ENDIF
    
  ENDIF
  
END SUBROUTINE

SUBROUTINE sparse_matrix_add_vect(matrix,irow,nval,jcol,val)

  TYPE(sparse_matrix_type),INTENT(inout),target :: matrix
  INTEGER(C_LONG)          ,INTENT(in)          :: irow,nval,jcol(nval)
  REAL(C_DOUBLE)          ,INTENT(in)           :: val(nval)
  
  TYPE(row_type),POINTER :: row
  INTEGER(C_LONG)         :: i,j,k,l,kbeg,kend
  
  IF(nval <= 0) RETURN
  
  IF(matrix%state == 1) THEN ! the fast structure is used
  
    row => matrix%fast%row(irow)
  
    loop1 : DO l=1,nval
    
      IF(ABS(val(l)) <= matrix%small) CYCLE
    
      DO j=1,row%nterm
        IF(row%colind(j)==jcol(l)) THEN
          row%values(j)=row%values(j)+val(l)
          CYCLE loop1
        ENDIF
      ENDDO
    
      CALL addterm(row,jcol(l),val(l))
    
      ! storing also a term in the symmetrical part to keep a symmetrical profile
    
      IF(matrix%lsym) THEN
        IF(jcol(l) /= irow) THEN
          CALL addterm(matrix%fast%row(jcol(l)),irow,0._c_double)
        ENDIF
      ENDIF
      
    ENDDO loop1
    
  ELSE IF(matrix%lorder == 1) THEN
  
    ! no ordering => taking advantage of the ordering of jcol
    
    k=matrix%rowptr(irow)
    kend=matrix%rowptr(irow+1)-1
    l=1
    
    DO l=1,nval
      IF(ABS(val(l)) <= matrix%small) CYCLE
      DO WHILE(jcol(l) /= matrix%colind(k))
        k=k+1
        IF(k > kend) GOTO 10
      ENDDO
      !$OMP ATOMIC
      matrix%values(k)=matrix%values(k)+val(l)
    ENDDO
    
  ELSE
  
    i=matrix%permut(irow)
    
    loop2 : DO l=1,nval
    
      IF(ABS(val(l)) <= matrix%small) CYCLE
    
      j=matrix%permut(jcol(l))
       
      kbeg=matrix%rowptr(i)
      kend=matrix%rowptr(i+1)-1
      
      DO
        IF(kend < kbeg) GOTO 10
        k=(kend+kbeg)/2
        IF(matrix%colind(k) == j) THEN
          !$OMP ATOMIC
          matrix%values(k)=matrix%values(k)+val(l)
          CYCLE loop2
        ELSEIF(matrix%colind(k) > j) THEN
          kend=k-1
        ELSEIF(matrix%colind(k) < j) THEN
          kbeg=k+1
        ENDIF
      ENDDO
      
    ENDDO loop2
      
  ENDIF  
    
  RETURN
  
10 CONTINUE  
   WRITE(*,*) 'sparse_matrix_add_vect : the position ',jcol(l),' is not foreseen in the row ',irow
   STOP

END SUBROUTINE

SUBROUTINE sparse_matrix_add(matrix,irow,jcol,val)

  TYPE(sparse_matrix_type),INTENT(inout),target :: matrix
  INTEGER(C_LONG)          ,INTENT(in)           :: irow,jcol
  REAL(C_DOUBLE)          ,INTENT(in)           :: val
  
  TYPE(row_type),POINTER :: row
  INTEGER(C_LONG)         :: k,j
  
  IF(abs(val) < matrix%small) RETURN
  
  IF(matrix%state == 1) THEN ! the fast structure is used
  
    row => matrix%fast%row(irow)
    
    DO j=1,row%nterm ! in fast, terms are not ordered
      IF(row%colind(j)==jcol) THEN
        !$OMP ATOMIC
        row%values(j)=row%values(j)+val
        RETURN
      ENDIF
    ENDDO
    
    CALL addterm(row,jcol,val)
    
    ! storing also a term in the symmetrical part to keep a symmetrical profile
    
    IF(matrix%lsym) THEN
      IF(jcol /= irow) THEN
        row => matrix%fast%row(jcol)
        CALL addterm(row,irow,0._c_double)
      ENDIF
    ENDIF
    
    RETURN

  ELSE
  
    k=locate(matrix,irow,jcol)
    
    IF(k == 0) THEN ! A change in the structure is observed
    
      WRITE(*,*) 'sparse_matrix_add : dangerous change of matrix state'
    
      matrix%state=1
      CALL create_fast(matrix) ! copying the existing structure into fast
      row => matrix%fast%row(irow)
      CALL addterm(row,jcol,val)
      IF(matrix%lsym) THEN
        IF(jcol /= irow) THEN
          row => matrix%fast%row(jcol)
          CALL addterm(row,irow,0._c_double)
        ENDIF
      ENDIF
      RETURN
      
    ELSE
    
      !$OMP ATOMIC
      matrix%values(k)=matrix%values(k)+val ! the normal case
      
    ENDIF
    
  ENDIF
  
END SUBROUTINE

SUBROUTINE sparse_matrix_get(matrix,irow,jcol,val)

  TYPE(sparse_matrix_type),INTENT(inout),target :: matrix
  INTEGER(C_LONG)          ,INTENT(in)           :: irow,jcol
  REAL(C_DOUBLE)          ,INTENT(out)          :: val
  
  TYPE(row_type),POINTER :: row
  INTEGER(C_LONG)         :: k,j
  
  IF(matrix%state == 1) THEN ! the fast structure is used
  
    row => matrix%fast%row(irow)
    
    DO j=1,row%nterm ! in fast, terms are not ordered
      IF(row%colind(j)==jcol) THEN
        val=row%values(j)
        RETURN
      ENDIF
    ENDDO
    
    val=0
    
  ELSE IF(irow == jcol) THEN
  
    IF(matrix%lorder == 1) THEN
      val=matrix%values(matrix%diagonal(irow))
    ELSE
      val=matrix%values(matrix%diagonal(matrix%permut(irow)))
    ENDIF
    
  ELSE
  
    k=locate(matrix,irow,jcol)
    
    IF(k <= 0) THEN ! term not found
      val=0
    ELSE
      val=matrix%values(k)
    ENDIF
    
  ENDIF
  
END SUBROUTINE

FUNCTION locate(matrix,ii,jj) RESULT(k) ! locating a term in the "normal" matrix (not the fast part)
  TYPE(sparse_matrix_type),INTENT(in) :: matrix
  INTEGER(C_LONG)          ,INTENT(in) :: ii,jj
  INTEGER(C_LONG)                      :: k
  INTEGER(C_LONG)                      :: kbeg,kend,i,j
  IF(matrix%lorder == 1) THEN
    i=ii
    j=jj
  ELSE
    i=matrix%permut(ii)
    IF(ii == jj) THEN
      k=matrix%diagonal(i)
      RETURN
    ENDIF
    j=matrix%permut(jj)
  ENDIF
  kbeg=matrix%rowptr(i)
  kend=matrix%rowptr(i+1)-1
  DO                                  ! dichotomic search terms being ordered
    IF(kend < kbeg) THEN
      k=0
      RETURN
    ENDIF
    k=(kend+kbeg)/2
    IF(matrix%colind(k) == j) THEN
      RETURN
    ELSEIF(matrix%colind(k) > j) THEN
      kend=k-1
    ELSEIF(matrix%colind(k) < j) THEN
      kbeg=k+1
    ENDIF
  ENDDO
END FUNCTION
    
SUBROUTINE create_fast(matrix) ! copying the "normal" matrix into fast because the structure of the matrix is changing

  TYPE(sparse_matrix_type),INTENT(inout),target :: matrix
  
  INTEGER(C_LONG)         :: i,kbeg,kend,nterm,ii
  TYPE(row_type),POINTER :: row
  
  IF(ALLOCATED(matrix%fast%row)) DEALLOCATE(matrix%fast%row)
  ALLOCATE(matrix%fast%row(matrix%n))
  
  DO i=1,matrix%n
    row => matrix%fast%row(i)
    ii=matrix%permut(i) ! the right row in the normal format
    kbeg=matrix%rowptr(ii)
    kend=matrix%rowptr(ii+1)-1
    nterm=kend-kbeg+1
    row%nterm=nterm
    IF(nterm > 0) THEN
      ALLOCATE(row%colind(nterm))
      ALLOCATE(row%values(nterm))
      row%colind(1:nterm)=matrix%order(matrix%colind(kbeg:kend)) ! reordering the indexes
      row%values(1:nterm)=matrix%values(kbeg:kend)
    ENDIF
  ENDDO
  
END SUBROUTINE

SUBROUTINE sparse_matrix_lu_factorization(matrix)

  TYPE(sparse_matrix_type),INTENT(inout) :: matrix
  
  IF(matrix%state == 0) THEN
    ! the matrix is ready for LU decomposition
    IF(ALLOCATED(matrix%lo)) THEN
      CALL lufactorizationsparse(matrix%n,matrix%lo,matrix%values)
    ELSE
      CALL lufactsparse(matrix%n,matrix%rowptr,matrix%colind,matrix%values)
    ENDIF
    matrix%state=3 ! LU done 
  ELSE
    WRITE(*,*) 'sparse_matrix_lu_factorization : the matrix state is not equal to 0'
    CALL abort()
  ENDIF
  
END SUBROUTINE

SUBROUTINE sparse_matrix_omp_lu_factorization(matrix)

  TYPE(sparse_matrix_type),INTENT(inout) :: matrix
  
  IF(matrix%state == 0) THEN
    ! the matrix is ready for LU decomposition
    IF(ALLOCATED(matrix%lo)) THEN
      CALL lufactorizationsparse(matrix%n,matrix%lo,matrix%values)
    ELSE
      CALL omplufactsparse(matrix%n,matrix%rowptr,matrix%colind,matrix%values)
    ENDIF
    matrix%state=3 ! LU done 
  ELSE
    WRITE(*,*) 'sparse_matrix_lu_factorization : the matrix state is not equal to 0'
    CALL abort()
  ENDIF
  
END SUBROUTINE

SUBROUTINE sparse_matrix_lu_vector(matrix,vector)
  TYPE(sparse_matrix_type),INTENT(in)    :: matrix
  REAL(C_DOUBLE)          ,INTENT(inout) :: vector(*)
  IF(matrix%lorder /= 1) vector(1:matrix%n)=vector(matrix%order) ! reordering the values
  IF(matrix%state == 3) THEN
    CALL luvectorsparse(matrix%n,matrix%rowptr,matrix%colind,matrix%diagonal,matrix%values,vector)
  ELSE
    STOP 'sparse_matrix_lu_vector : the matrix state is not equal to 3 (LU form)'
  ENDIF
  IF(matrix%lorder /= 1) vector(1:matrix%n)=vector(matrix%permut) ! restoring the initial order
END SUBROUTINE

SUBROUTINE sparse_matrix_vector(matrix,vector,result)
  TYPE(sparse_matrix_type),INTENT(in),target    :: matrix
  REAL(C_DOUBLE)          ,INTENT(in)    :: vector(*)
  REAL(C_DOUBLE)          ,INTENT(out)   :: result(*)
  INTEGER(C_LONG) :: i,j,jj
  TYPE(row_type),POINTER :: row
  result(1:matrix%n)=0
  IF(matrix%state == 1) THEN
    !$OMP PARALLEL DO PRIVATE(i,row,j)
    DO i=1,matrix%n
      row => matrix%fast%row(i)
      DO j=1,row%nterm
        result(i)=result(i)+row%values(j)*vector(row%colind(j))
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO
  ELSE
    !$OMP PARALLEL DO PRIVATE(i,j,jj)
    DO i=1,matrix%n
      DO j=matrix%rowptr(i),matrix%rowptr(i+1)-1
        jj=matrix%order(matrix%colind(j))
        result(i)=result(i)+matrix%values(j)*vector(jj)
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO
    result(1:matrix%n)=result(matrix%permut)
  ENDIF
END SUBROUTINE

SUBROUTINE sparse_matrix_lu_extract(matrix,m,list,a)
  TYPE(sparse_matrix_type),INTENT(in)    :: matrix
  INTEGER(C_LONG)          ,INTENT(in)    :: m       ! number of retained unknowns
  INTEGER(C_LONG)          ,INTENT(in)    :: list(m) ! list of retained unknowns
  REAL(C_DOUBLE)          ,INTENT(out)   :: a(m,m)  ! extracted inversed matrix
  INTEGER(C_LONG)  :: l(m) 
  IF(matrix%state == 3) THEN
    l=matrix%permut(list)
    CALL invert(matrix%n,matrix%rowptr,matrix%diagonal,matrix%colind,matrix%values,m,l,a)
  ELSE
    STOP 'sparse_matrix_lu_extract : the matrix state is not equal to 3'
  ENDIF
END SUBROUTINE

SUBROUTINE sparse_matrix_optimize(matrix,noper,iopt)

  ! minimum degree algorithm for a matrix with a symmetrical profile

  ! this is the simplest algorithm which does not manage optimizations
  ! like supernodes, multi elimination, approximation of minimum degree
  ! or parallel techniques. Anyway, it is enough for my needs because I
  ! just want to prepare a correct ordering of unknowns for a transient
  ! code : the same ordering will be used many times so I can sacrify some
  ! CPU time to compute it. And matrices I want to optimize are sparse and
  ! relatively small (n <= 2000)

  ! this routine starts from a compressed profile (matrix%fast) and computes
  ! - a permutation of indexes (permut(n) and the reverse function order(n))
  ! - a new compressed profile, using the new ordering and taking into account
  !   the fill in terms during a LU factorization
  ! - The new matrix after permutation is stored into matrix%colind matrix%rowptr and matrix%values
  ! - in addition, the work to compute LU factorization is estimated (matrix%no) and the list
  !   of operations is registered into matrix%lo
  ! - if the matrix is already in its optimized form, then noper is return equal to 0, else it is
  !   returned equal to matrix%no

  TYPE(sparse_matrix_type),INTENT(inout),target :: matrix             ! current matrix
  INTEGER(C_LONG)          ,INTENT(out)          :: noper              ! number of operations for LU factorization (0 => already optimized)
  INTEGER(C_LONG)          ,INTENT(in) ,OPTIONAL :: iopt               ! =1 => no ordering, -1=use a pre ordering

  TYPE(row_type),POINTER :: row(:)
  INTEGER(C_LONG) :: inew,iold,jnew,jold,k,knew,kold,l,lold,best,mindeg,npos,i,nterm,n,nlist
  INTEGER(C_LONG) :: list(matrix%n),ipos(matrix%n),colind(matrix%n)
  REAL(C_DOUBLE) :: values(matrix%n)
  INTEGER(C_LONG),POINTER :: order(:),permut(:)

  noper=0
  IF(matrix%state .NE. 1) RETURN ! already optimized
  
  IF(PRESENT(iopt)) matrix%lorder=iopt
  
  ! default ordering
  
  n=matrix%n
  
  order  => matrix%order
  permut => matrix%permut
  
  IF(matrix%lorder /= -1) THEN ! already ordered if =-1      
    DO i=1,n
      order (i)=i
      permut(i)=i
    ENDDO
  ENDIF
             
  row=>matrix%fast%row
  
  row(:)%degree=row(:)%nterm        ! to initialize the degree of rows

  noper=0

  DO inew=1,n                       ! loop on rows in the new order
    iold=order(inew)                ! old index of the current row
    mindeg=row(iold)%degree         ! degree of that row
    IF(matrix%lorder /= 1) THEN     ! do not change the ordering if iopt == 1
      best=inew                     ! for the moment, the current row is the best one
      DO jnew=inew+1,n              ! looking for the minimum degree 
        jold=order(jnew)
        IF(row(jold)%degree < mindeg) THEN
          best=jnew                 ! jnew is better
          mindeg=row(jold)%degree
        ENDIF
      ENDDO
      IF(best /= inew) THEN         ! permuting best and inew if necessary
        order(inew)=order(best)
        order(best)=iold
        permut(order(inew))=inew
        permut(order(best))=best
        iold=order(inew)
      ENDIF
    ENDIF
    ! computing fill in and updating the degree of touched rows
    npos=0
    DO k=1,row(iold)%nterm
      kold=row(iold)%colind(k)
      knew=permut(kold)
      IF(knew <= inew) CYCLE
      npos=npos+1
      ipos(npos)=kold               ! term which propagates
    ENDDO
    DO l=1,npos                     ! rows where fill in terms appear (here we use the symmetric profile)
      lold=ipos(l)
      DO k=1,npos
        colind(ipos(k))=0           ! unmark terms which may appear
      ENDDO
      nlist=0
      DO k=1,row(lold)%nterm
        kold=row(lold)%colind(k)
        knew=permut(kold)
        IF(knew <= inew) CYCLE
        colind(kold)=1              ! mark terms which are already in the row
        nlist=nlist+1
      ENDDO
      DO k=1,npos
        kold=ipos(k)
        IF(colind(kold) == 0) THEN
          CALL addterm(row(lold),kold,0._c_double)
          nlist=nlist+1
        ENDIF
      ENDDO
      row(lold)%degree=nlist
    ENDDO
    noper=noper+2+npos+npos*(npos+1)
  ENDDO

  ! building up the new profile

  nterm=SUM(row(:)%nterm)
  IF(ALLOCATED(matrix%colind)) DEALLOCATE(matrix%colind)
  IF(ALLOCATED(matrix%values)) DEALLOCATE(matrix%values)
  ALLOCATE(matrix%colind(nterm),matrix%values(nterm))
  
  k=0
  DO inew=1,n
  
    matrix%rowptr(inew)=k+1
    iold=order(inew)
    nterm=row(iold)%nterm
    
    DO l=1,nterm
      colind(l)=permut(row(iold)%colind(l))
      values(l)=row(iold)%values(l)
    ENDDO
    
    ! ordering to get new column indexes in increasing order
    
    DO l=1,nterm ; list(l)=l ; ENDDO
    
    CALL ordering(nterm,colind,list)
    
    matrix%diagonal(inew)=0
   
    DO l=1,nterm
      k=k+1
      matrix%colind(k)=colind(list(l))
      matrix%values(k)=values(list(l))
      IF(matrix%colind(k) == inew) matrix%diagonal(inew)=k ! registering the position of the diagonal
    ENDDO 
    
    IF(matrix%diagonal(inew) == 0) THEN
      WRITE(*,*) 'sparse_matrix_optimize : the row ',iold,' has no diagonal term '
      STOP
    ENDIF
     
    DO l=matrix%rowptr(inew),k-1
      IF(matrix%colind(l) >= matrix%colind(l+1)) THEN
        write(*,*) 'line ',inew,' invalid ',matrix%colind(l),matrix%colind(l+1)
        write(*,*) 'iold ',iold
        write(*,*) colind(1:nterm)
        STOP
      ENDIF
    ENDDO

  ENDDO
  
  matrix%rowptr(n+1)=k+1
  matrix%state=0               ! the matrix is optimized now
  DEALLOCATE(matrix%fast%row)  ! the fast form is not used anymore
  
  ! registering operations for further LU decompositions
  
  IF(ALLOCATED(matrix%lo)) DEALLOCATE(matrix%lo)
  IF(noper <= matrix%nomax) THEN
    ALLOCATE(matrix%lo(noper))
    CALL luregistersparse(matrix%n,matrix%rowptr,matrix%colind,matrix%diagonal,noper,matrix%lo)
  ENDIF

END SUBROUTINE          

SUBROUTINE addterm(row,col,val)

  TYPE(row_type) ,INTENT(inout) :: row
  INTEGER(C_LONG),INTENT(in)    :: col
  REAL(C_DOUBLE) ,INTENT(in)    :: val
  
  INTEGER(C_LONG),ALLOCATABLE :: colind(:)
  REAL(C_DOUBLE) ,ALLOCATABLE :: values(:)
  INTEGER(C_LONG)             :: k,ifind
  
  IF(row%nterm == 0) THEN
    ALLOCATE(row%colind(mterm),row%values(mterm))
  ELSE IF(row%nterm == SIZE(row%colind)) THEN
    ALLOCATE(colind(2*row%nterm),values(2*row%nterm))
    colind(1:row%nterm)=row%colind(1:row%nterm)
    values(1:row%nterm)=row%values(1:row%nterm)
    CALL move_alloc(colind,row%colind)
    CALL move_alloc(values,row%values)
  ENDIF
  
  ifind=0
  DO k=1,row%nterm
    IF(row%colind(k) == col) THEN
      ifind=1
      EXIT
    ENDIF
  ENDDO
  
  IF(ifind == 0) THEN
    row%nterm=row%nterm+1
    row%colind(row%nterm)=col
    row%values(row%nterm)=val
  ELSE
    row%values(row%nterm)=row%values(row%nterm)+val
  ENDIF
  
END SUBROUTINE

SUBROUTINE sparse_matrix_print(matrix)
  TYPE(sparse_matrix_type),INTENT(in) :: matrix
  INTEGER(C_LONG) :: i,ii,k
  WRITE(*,*) 'matrix state : ',matrix%state
  SELECT CASE(matrix%state)
  CASE(0,3)
    DO ii=1,matrix%n
      i=matrix%permut(ii)
      WRITE(*,*) 'row ',ii,' nterm=',matrix%rowptr(i+1)-matrix%rowptr(i)
      DO k=matrix%rowptr(i),matrix%rowptr(i+1)-1
        WRITE(*,*) '   column ',matrix%order(matrix%colind(k)),' value ',matrix%values(k)
      ENDDO
    ENDDO
  CASE(1)
    DO i=1,matrix%n
      WRITE(*,*) 'row ',i,' nterm=',matrix%fast%row(i)%nterm
      DO k=1,matrix%fast%row(i)%nterm
        WRITE(*,*) '   column ',matrix%fast%row(i)%colind(k),' value ',matrix%fast%row(i)%values(k)
      ENDDO
    ENDDO
  END SELECT
END SUBROUTINE

SUBROUTINE luregistersparse(n,rowptr,colind,diag,no,lo)
  INTEGER(C_LONG),INTENT(in)    :: n
  INTEGER(C_LONG),INTENT(in)    :: rowptr(n+1)
  INTEGER(C_LONG),INTENT(in)    :: colind(*)
  INTEGER(C_LONG),INTENT(in)    :: diag(n)
  INTEGER(C_LONG),INTENT(inout) :: no
  INTEGER(C_LONG),INTENT(out)   :: lo(*)
  INTEGER(C_LONG) :: i,j,k,l,li,listi(n),listk(n),listkk(n),m,coef,noper
  noper=no
  no=0
  DO k=1,n
    no=no+1 ; lo(no)=diag(k) ! diagonal
    m=0
    DO l=diag(k)+1,rowptr(k+1)-1
      m=m+1
      listk(m)=colind(l)
      listkk(m)=l
    ENDDO
    no=no+1 ; lo(no)=m
    DO l=1,m
      no=no+1 ; lo(no)=listkk(l)
    ENDDO
    DO li=1,m
      i=listk(li)
      coef=0
      DO l=1,m
        j=listk(l)
        listi(j)=0
      ENDDO
      DO l=rowptr(i),rowptr(i+1)-1
        j=colind(l)
        IF(j == k) THEN
          coef=coef+1
          no=no+1 ; lo(no)=l
        ELSE IF(j > k) THEN
          listi(j)=l
        ENDIF
      ENDDO
      IF(coef /= 1) THEN
        write(*,*) 'cycle k=',k,' i=',i
        STOP
      ENDIF
      DO l=1,m
        j=listk(l)
        IF(listi(j) == 0) THEN
          write(*,*) 'problem row ',i,' column',j,' k=',k
          STOP
        ENDIF
        no=no+1 ; lo(no)=listi(j)
      ENDDO
    ENDDO
  ENDDO
  IF(no > noper) THEN
    write(*,*) 'luregistersparse : too low dimension for lo ',no,noper
    STOP
  ENDIF
END SUBROUTINE

SUBROUTINE omplufactsparse(n,rowptr,colind,mat) ! factorization without using the registered operations
  INTEGER(C_LONG),INTENT(in)    :: n
  INTEGER(C_LONG),INTENT(in)    :: rowptr(n+1)
  INTEGER(C_LONG),INTENT(in)    :: colind(*)
  REAL(C_DOUBLE),INTENT(inout) :: mat(*)
  REAL(C_DOUBLE) :: diag,coef
  INTEGER(C_LONG) :: k,kdiag,li,i,idiag,kpos,ipos
  DO k=1,n
    ! looking for the diagonal term
    DO kdiag=rowptr(k),rowptr(k+1)-1
      IF(colind(kdiag)==k) EXIT
    ENDDO
    diag=1/mat(kdiag)
    mat(kdiag)=diag
    !$OMP PARALLEL DO PRIVATE(li,i,idiag,coef,kpos,ipos)
    DO li=kdiag+1,rowptr(k+1)-1
      i=colind(li)
      DO idiag=rowptr(i),rowptr(i+1)-1
        IF(colind(idiag)==k) EXIT
      ENDDO
      coef=-mat(idiag)*diag
      mat(idiag)=-coef
      kpos=kdiag+1
      ipos=idiag+1
      DO
        IF(colind(kpos) < colind(ipos)) THEN
          kpos=kpos+1
          IF(kpos ==  rowptr(k+1)) EXIT
        ELSE IF(colind(kpos) == colind(ipos)) THEN
          mat(ipos)=mat(ipos)+coef*mat(kpos)
          kpos=kpos+1
          IF(kpos ==  rowptr(k+1)) EXIT
          ipos=ipos+1
        ELSE
          ipos=ipos+1
        ENDIF
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO
  ENDDO
END SUBROUTINE      

SUBROUTINE lufactsparse(n,rowptr,colind,mat) ! factorization without using the registered operations
  INTEGER(C_LONG),INTENT(in)    :: n
  INTEGER(C_LONG),INTENT(in)    :: rowptr(n+1)
  INTEGER(C_LONG),INTENT(in)    :: colind(*)
  REAL(C_DOUBLE),INTENT(inout) :: mat(*)
  REAL(C_DOUBLE) :: diag,coef
  INTEGER(C_LONG) :: k,kdiag,li,i,idiag,kpos,ipos
  DO k=1,n
    ! looking for the diagonal term
    DO kdiag=rowptr(k),rowptr(k+1)-1
      IF(colind(kdiag)==k) EXIT
    ENDDO
    diag=1/mat(kdiag)
    mat(kdiag)=diag
    DO li=kdiag+1,rowptr(k+1)-1
      i=colind(li)
      DO idiag=rowptr(i),rowptr(i+1)-1
        IF(colind(idiag)==k) EXIT
      ENDDO
      coef=-mat(idiag)*diag
      mat(idiag)=-coef
      kpos=kdiag+1
      ipos=idiag+1
      DO
        IF(colind(kpos) < colind(ipos)) THEN
          kpos=kpos+1
          IF(kpos ==  rowptr(k+1)) EXIT
        ELSE IF(colind(kpos) == colind(ipos)) THEN
          mat(ipos)=mat(ipos)+coef*mat(kpos)
          kpos=kpos+1
          IF(kpos ==  rowptr(k+1)) EXIT
          ipos=ipos+1
        ELSE
          ipos=ipos+1
        ENDIF
      ENDDO
    ENDDO
  ENDDO
END SUBROUTINE      

SUBROUTINE lufactorizationsparse(n,lo,mat) ! factorization using the registered operations
  INTEGER(C_LONG),INTENT(in)    :: n
  INTEGER(C_LONG),INTENT(in)    :: lo(*)
  REAL(C_DOUBLE),INTENT(inout) :: mat(*)
  INTEGER(C_LONG) :: i,j,k,l,no,nij
  REAL(C_DOUBLE) :: row(n),diag,coef
  no=0
  DO k=1,n
    no=no+1 ; l=lo(no) ; 
    IF(ABS(mat(l)) < 1.d-100) THEN 
      diag=0 
    ELSE 
      diag=1/mat(l) 
    ENDIF 
    mat(l)=diag
    no=no+1 ; nij=lo(no)
    DO j=1,nij
      no=no+1 ; row(j)=mat(lo(no))
    ENDDO
    DO i=1,nij
      no=no+1 ; l=lo(no) ; coef=-mat(l)*diag ; mat(l)=-coef
      DO j=1,nij
        no=no+1 ; l=lo(no) ; mat(l)=mat(l)+coef*row(j)
      ENDDO
    ENDDO
  ENDDO
END SUBROUTINE

SUBROUTINE luvectorsparse(n,rowptr,colind,diag,mat,vect)
  INTEGER(C_LONG),INTENT(in)    :: n
  INTEGER(C_LONG),INTENT(in)    :: rowptr(n+1)
  INTEGER(C_LONG),INTENT(in)    :: colind(*)
  INTEGER(C_LONG),INTENT(in)    :: diag(n)
  REAL(C_DOUBLE),INTENT(in)    :: mat(*)      
  REAL(C_DOUBLE),INTENT(inout) :: vect(n)
  INTEGER(C_LONG) :: k,l
  ! lower triangular matrix
  DO k=1,n
    DO l=rowptr(k),diag(k)-1
      vect(k)=vect(k)-mat(l)*vect(colind(l))
    ENDDO
  ENDDO
  ! upper triangular matrix
  DO k=n,1,-1
    DO l=diag(k)+1,rowptr(k+1)-1
      vect(k)=vect(k)-mat(l)*vect(colind(l))
    ENDDO
    vect(k)=vect(k)*mat(diag(k))
  ENDDO
END SUBROUTINE
  
RECURSIVE SUBROUTINE ordering(n,values,list)
  INTEGER(C_LONG),INTENT(in)    :: n
  INTEGER(C_LONG),INTENT(in)    :: values(*)
  INTEGER(C_LONG),INTENT(inout) :: list(n) ! must be initialized
  INTEGER(C_LONG) :: n1,n2,l(n)
  IF(n <= 1) RETURN
  n1=n/2
  n2=n-n/2
  l(:)=list
  CALL ordering(n1,values,l(1:n1))
  CALL ordering(n2,values,l(n1+1:n2))
  CALL concat(values,n1,l(1:n1),n2,l(n1+1:n2),list)
END SUBROUTINE

SUBROUTINE concat(values,n1,list1,n2,list2,list)
  INTEGER(C_LONG),INTENT(in)  :: values(*)
  INTEGER(C_LONG),INTENT(in)  :: n1,list1(n1),n2,list2(n2)
  INTEGER(C_LONG),INTENT(out) :: list(n1+n2)
  INTEGER(C_LONG) :: i1,i2,i,n
  i1=1
  i2=1
  n=0
  DO
    IF(i1 > n1 .AND. i2 > n2) RETURN
    IF(i1 > n1) THEN
      DO i=i2,n2
        n=n+1
        list(n)=list2(i)
      ENDDO
      EXIT
    ELSE IF(i2 > n2) THEN
      DO i=i1,n1
        n=n+1
        list(n)=list1(i)
      ENDDO
      EXIT
    ENDIF
    IF(values(list1(i1)) <= values(list2(i2))) THEN
      n=n+1
      list(n)=list1(i1)
      i1=i1+1
    ELSE
      n=n+1
      list(n)=list2(i2)
      i2=i2+1
    ENDIF
  ENDDO
END SUBROUTINE

SUBROUTINE invert(n,rowptr,diagonal,colind,matrix,m,list,a)

  ! Computing a limited part of the invert matrix knowing
  ! the LU factorization. Of course m << n
  
  INTEGER(C_LONG),INTENT(in)  :: n           ! matrix size
  INTEGER(C_LONG),INTENT(in)  :: rowptr(n+1) ! row pointers
  INTEGER(C_LONG),INTENT(in)  :: diagonal(n) ! position of diagonal terms
  INTEGER(C_LONG),INTENT(in)  :: colind(*)   ! column indexes
  REAL(C_DOUBLE) ,INTENT(in)  :: matrix(*)   ! LU factorization
  INTEGER(C_LONG),INTENT(in)  :: m           ! number of retained indexes
  INTEGER(C_LONG),INTENT(in)  :: list(m)     ! list of retained indexes
  REAL(C_DOUBLE) ,INTENT(out) :: a(m,m)      ! part of the invert matrix
  
  INTEGER(C_LONG) :: i,j,k,kj
  REAL(C_DOUBLE) :: b(n,m)
  
  b(:,:)=0
  DO k=1,m
    b(list(k),k)=1
  ENDDO
  
  DO k=1,m
    DO i=list(k)+1,n
      DO kj=rowptr(i),rowptr(i+1)-1
        j=colind(kj)
        IF(j < list(k)) CYCLE
        IF(j >= i) CYCLE
        b(i,k)=b(i,k)-matrix(kj)*b(j,k)
      ENDDO
    ENDDO
  ENDDO
  
  DO k=1,m
    DO i=n,1,-1
      DO kj=rowptr(i),rowptr(i+1)-1
        j=colind(kj)
        IF(j <= i) CYCLE
        b(i,k)=b(i,k)-matrix(kj)*b(j,k)
      ENDDO
      b(i,k)=b(i,k)*matrix(diagonal(i))
    ENDDO
  ENDDO
  
  DO i=1,m
    DO j=1,m
      a(i,j) = b(list(i),j)
    ENDDO
  ENDDO
  
END SUBROUTINE

SUBROUTINE sparse_matrix_omp_add_matrix(m1,m2) ! computing m1 <- m1+m2

  TYPE(sparse_matrix_type),INTENT(inout),TARGET :: m1
  TYPE(sparse_matrix_type),INTENT(in)   ,TARGET :: m2
  
  TYPE(row_type),POINTER :: r1,r2
  INTEGER(C_LONG)        :: i,j,k,l
  
  IF(m2%state == -1) RETURN
  IF(m1%n /= m2%n)   STOP "the matrices have not the same size"
  
  SELECT CASE(10*m1%state+m2%state)
  
  CASE(0)
    !$OMP PARALLEL DO PRIVATE(i,j,k,l)
    DO i=1,m2%n
      j=m2%permut(i)
      DO k=m2%rowptr(i),m2%rowptr(i+1)-1
        l=m2%permut(m2%colind(k))
        CALL sparse_matrix_add(m1,j,l,m2%values(k))
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO
  CASE(1)
    !$OMP PARALLEL DO PRIVATE(i,j,k)
    DO i=1,m2%n
      r2 => m2%fast%row(i)
      DO k=1,r2%nterm
        CALL sparse_matrix_add(m1,i,r2%colind(k),r2%values(k))
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO
  CASE(10)
    !$OMP PARALLEL DO PRIVATE(i,j,k,l,r1)
    DO i=1,m2%n
      r1 => m1%fast%row(m2%permut(i))
      main2 : DO k=m2%rowptr(i),m2%rowptr(i+1)-1
        j=m2%permut(m2%colind(k))
        DO l=1,r1%nterm
          IF(r1%colind(l) == j) THEN
            r1%values(l)=r1%values(l)+m2%values(k)
            CYCLE main2
          ENDIF
        ENDDO
        CALL addterm(r1,j,m2%values(k))
      ENDDO main2
    ENDDO
    !$OMP END PARALLEL DO
  CASE(11)
    !$OMP PARALLEL DO PRIVATE(i,k,l,r1,r2)
    DO i=1,m2%n
      r1 => m1%fast%row(i)
      r2 => m2%fast%row(i)
      main1 : DO k=1,r2%nterm
        j=r2%colind(k)
        DO l=1,r1%nterm
          IF(r1%colind(l) == j) THEN
            r1%values(l)=r1%values(l)+r2%values(k)
            CYCLE main1
          ENDIF
        ENDDO
        CALL addterm(r1,j,r2%values(k))
      ENDDO main1
    ENDDO
    !$OMP END PARALLEL DO
  CASE DEFAULT
    STOP 'sparse_madrix_add : invalid matrix status '
  END SELECT
  
END SUBROUTINE

SUBROUTINE sparse_matrix_add_matrix(m1,m2) ! computing m1 <- m1+m2

  TYPE(sparse_matrix_type),INTENT(inout),TARGET :: m1
  TYPE(sparse_matrix_type),INTENT(in)   ,TARGET :: m2
  
  TYPE(row_type),POINTER :: r1,r2
  INTEGER(C_LONG)        :: i,j,k,l
  
  IF(m2%state == -1) RETURN
  IF(m1%n /= m2%n)   STOP "the matrices have not the same size"
  
  SELECT CASE(10*m1%state+m2%state)
  
  CASE(0)
    DO i=1,m2%n
      j=m2%permut(i)
      DO k=m2%rowptr(i),m2%rowptr(i+1)-1
        l=m2%permut(m2%colind(k))
        CALL sparse_matrix_add(m1,j,l,m2%values(k))
      ENDDO
    ENDDO
  CASE(1)
    DO i=1,m2%n
      r2 => m2%fast%row(i)
      DO k=1,r2%nterm
        CALL sparse_matrix_add(m1,i,r2%colind(k),r2%values(k))
      ENDDO
    ENDDO
  CASE(10)
    DO i=1,m2%n
      r1 => m1%fast%row(m2%permut(i))
      main2 : DO k=m2%rowptr(i),m2%rowptr(i+1)-1
        j=m2%permut(m2%colind(k))
        DO l=1,r1%nterm
          IF(r1%colind(l) == j) THEN
            r1%values(l)=r1%values(l)+m2%values(k)
            CYCLE main2
          ENDIF
        ENDDO
        CALL addterm(r1,j,m2%values(k))
      ENDDO main2
    ENDDO
  CASE(11)
    DO i=1,m2%n
      r1 => m1%fast%row(i)
      r2 => m2%fast%row(i)
      main1 : DO k=1,r2%nterm
        j=r2%colind(k)
        DO l=1,r1%nterm
          IF(r1%colind(l) == j) THEN
            r1%values(l)=r1%values(l)+r2%values(k)
            CYCLE main1
          ENDIF
        ENDDO
        CALL addterm(r1,j,r2%values(k))
      ENDDO main1
    ENDDO
  CASE DEFAULT
    STOP 'sparse_madrix_add : invalid matrix status '
  END SELECT
  
END SUBROUTINE

!! portable C Interface

SUBROUTINE sparse_matrix_init_C(m,n,mode,eps,noper) BIND(C,name="sparse_matrix_init")
  TYPE(C_PTR)   ,INTENT(out)      :: m
  INTEGER(C_LONG),INTENT(in),VALUE :: n,mode,noper
  REAL(C_DOUBLE),INTENT(in),VALUE :: eps
  TYPE(sparse_matrix_type),POINTER :: matrix
  ALLOCATE(matrix)
  CALL sparse_matrix_init(matrix,n,mode,eps,noper)
  m=C_LOC(matrix)
END SUBROUTINE

SUBROUTINE sparse_matrix_add_C(m,i,j,val) BIND(C,name="sparse_matrix_add")
  TYPE(C_PTR),VALUE :: m
  INTEGER(C_LONG),INTENT(in),VALUE :: i,j
  REAL(C_DOUBLE),INTENT(in),VALUE :: val
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_add(matrix,i+1,j+1,val)
END SUBROUTINE

SUBROUTINE sparse_matrix_add_vect_C(m,i,n,j,val) BIND(C,name="sparse_matrix_add_vect")
  TYPE(C_PTR),VALUE :: m
  INTEGER(C_LONG),INTENT(in),VALUE :: i,n
  INTEGER(C_LONG),INTENT(in)       :: j(n)
  REAL(C_DOUBLE) ,INTENT(in)       :: val(n)
  INTEGER(C_LONG) :: jj(n)
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  jj=j+1
  CALL sparse_matrix_add_vect(matrix,i+1,n,jj,val)
END SUBROUTINE

SUBROUTINE sparse_matrix_put_C(m,i,j,val) BIND(C,name="sparse_matrix_put")
  TYPE(C_PTR),VALUE :: m
  INTEGER(C_LONG),INTENT(in),VALUE :: i,j
  REAL(C_DOUBLE),INTENT(in),VALUE :: val
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_put(matrix,i+1,j+1,val)
END SUBROUTINE

SUBROUTINE sparse_matrix_reset_C(m) BIND(C,name="sparse_matrix_reset")
  TYPE(C_PTR),VALUE :: m
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_reset(matrix)
END SUBROUTINE

SUBROUTINE sparse_matrix_optimize_C(m,noper,iopt) BIND(C,name="sparse_matrix_optimize")
  TYPE(C_PTR),VALUE :: m
  INTEGER(C_LONG),VALUE,INTENT(in) :: iopt ! 0=normal , 1=no ordering, -1 already ordered
  INTEGER(C_LONG),INTENT(out) :: noper
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_optimize(matrix,noper,iopt)
END SUBROUTINE

SUBROUTINE sparse_matrix_lu_factorization_C(m) BIND(C,name="sparse_matrix_lu_factorization")
  TYPE(C_PTR),VALUE :: m
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_lu_factorization(matrix)
END SUBROUTINE

SUBROUTINE sparse_matrix_omp_lu_factorization_C(m) BIND(C,name="sparse_matrix_omp_lu_factorization")
  TYPE(C_PTR),VALUE :: m
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_omp_lu_factorization(matrix)
END SUBROUTINE

SUBROUTINE sparse_matrix_lu_vector_C(m,vector) BIND(C,name="sparse_matrix_lu_vector")
  TYPE(C_PTR),VALUE,INTENT(in) :: m
  REAL(C_DOUBLE),INTENT(inout) :: vector(*)
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_lu_vector(matrix,vector)
END SUBROUTINE

SUBROUTINE sparse_matrix_vector_C(m,vector,result) BIND(C,name="sparse_matrix_vector")
  TYPE(C_PTR),VALUE,INTENT(in)  :: m
  REAL(C_DOUBLE)   ,INTENT(in)  :: vector(*)
  REAL(C_DOUBLE)   ,INTENT(out) :: result(*)
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_vector(matrix,vector,result)
END SUBROUTINE

SUBROUTINE sparse_matrix_lu_extract_C(m,nlist,list,a) BIND(C,name="sparse_matrix_lu_extract")
  TYPE   (C_PTR) ,VALUE,INTENT(in)  :: m
  INTEGER(C_LONG),VALUE,INTENT(in)  :: nlist
  INTEGER(C_LONG),      INTENT(in)  :: list(nlist)
  REAL(C_DOUBLE) ,      INTENT(out) :: a(nlist,nlist)
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_lu_extract(matrix,nlist,list,a)
END SUBROUTINE

SUBROUTINE sparse_matrix_print_C(m) BIND(C,name="sparse_matrix_print")
  TYPE(C_PTR),VALUE,INTENT(in)  :: m
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_print(matrix)
END SUBROUTINE

FUNCTION sparse_matrix_get_C(m,i,j) RESULT(val) BIND(C,name="sparse_matrix_get")
  TYPE   (C_PTR) ,VALUE,INTENT(in)  :: m
  INTEGER(C_LONG),VALUE,INTENT(in)  :: i,j
  REAL(C_DOUBLE)                    :: val
  TYPE(sparse_matrix_type),POINTER  :: matrix
  CALL c_f_pointer(m,matrix)
  CALL sparse_matrix_get(matrix,i+1,j+1,val)
END FUNCTION

SUBROUTINE sparse_matrix_free_C(m) BIND(C,name="sparse_matrix_free")
  TYPE   (C_PTR),VALUE,INTENT(in)  :: m
  TYPE(sparse_matrix_type),POINTER :: matrix
  CALL c_f_pointer(m,matrix)
  DEALLOCATE(matrix)
END SUBROUTINE

SUBROUTINE sparse_matrix_copy_C(m2,m1) BIND(C,name="sparse_matrix_copy")
  TYPE   (C_PTR),VALUE,INTENT(in)  :: m1
  TYPE   (C_PTR)      ,INTENT(out) :: m2
  TYPE(sparse_matrix_type),POINTER :: mat1,mat2
  CALL c_f_pointer(m1,mat1)
  ALLOCATE(mat2)
  mat2=mat1
  m2=C_LOC(mat2)
END SUBROUTINE

SUBROUTINE sparse_matrix_add_matrix_C(m1,m2) BIND(C,name="sparse_matrix_add_matrix")
  TYPE   (C_PTR),VALUE :: m1
  TYPE   (C_PTR),VALUE :: m2
  TYPE(sparse_matrix_type),POINTER :: mat1,mat2
  CALL c_f_pointer(m1,mat1)
  CALL c_f_pointer(m2,mat2)
  CALL sparse_matrix_add_matrix(mat1,mat2)
END SUBROUTINE

SUBROUTINE sparse_matrix_omp_add_matrix_C(m1,m2) BIND(C,name="sparse_matrix_omp_add_matrix")
  TYPE   (C_PTR),VALUE :: m1
  TYPE   (C_PTR),VALUE :: m2
  TYPE(sparse_matrix_type),POINTER :: mat1,mat2
  CALL c_f_pointer(m1,mat1)
  CALL c_f_pointer(m2,mat2)
  CALL sparse_matrix_omp_add_matrix(mat1,mat2)
END SUBROUTINE

SUBROUTINE sparse_matrix_get_order_C(matrix,order) BIND(C,name="sparse_matrix_get_order")
  TYPE   (C_PTR),VALUE :: matrix
  TYPE   (C_PTR),INTENT(out) :: order
  TYPE(sparse_matrix_type),POINTER :: m
  INTEGER(C_LONG),POINTER :: o(:)
  CALL c_f_pointer(matrix,m)
  ALLOCATE(o(m%n))
  o=m%order-1
  order=C_LOC(o)
END SUBROUTINE

SUBROUTINE sparse_matrix_get_permut_C(matrix,permut) BIND(C,name="sparse_matrix_get_permut")
  TYPE   (C_PTR),VALUE :: matrix
  TYPE   (C_PTR),INTENT(out) :: permut
  TYPE(sparse_matrix_type),POINTER :: m
  INTEGER(C_LONG),POINTER :: p(:)
  CALL c_f_pointer(matrix,m)
  ALLOCATE(p(m%n))
  p=m%permut-1
  permut=C_LOC(p)
END SUBROUTINE

SUBROUTINE sparse_matrix_change_order_C(matrix,lorder) BIND(C,name="sparse_matrix_change_order")
  TYPE   (C_PTR),VALUE :: matrix
  INTEGER(C_INT),VALUE :: lorder
  TYPE(sparse_matrix_type),POINTER :: m
  CALL c_f_pointer(matrix,m)
  m%lorder=lorder
END SUBROUTINE

END MODULE    
