program localtocsr
implicit none

! include "mkl.fi" !including these include statements gives syntax errors I cannot resolve
! include "mkl_spblas.fi" 
 integer(kind=4),parameter    :: nep=4,nelem=3,nnodes=6
 integer(kind=4)              :: i,j,intma_G(nelem,nep)
 double precision             :: M1(nep,nep),M2(nep,nep),M3(nep,nep),ME(nep,nep),MG(nnodes,nnodes)

 !Intel mkl specific
 integer                      :: job(8),m,n,lda,info,nzmax
 integer,allocatable          :: ja(:)
 integer,allocatable          :: ia(:)
 double precision,allocatable :: acsr(:)
 double precision,allocatable :: adns(:,:)

 integer,allocatable          :: ib(:),jb(:),ic(:),jc(:),iG(:),jG(:)
 double precision,allocatable :: bcsr(:),ccsr(:),Gcsr(:)
 integer                      :: iadd,jadd,idum,sizeint,request,sort,iaux
 double precision             :: csradd,ddum,beta
 pointer  (iadd_ptr,iadd(1)),(jadd_ptr,jadd(1)),(csradd_ptr,csradd(1)),(iaux_ptr,iaux(1))
 character*1                  ::trans

! #ifdef _IA32 !including these conditional statements again gives syntax errors
! integer*4 MKL_MALLOC
! integer*4 alloc_size
! #else
 integer*8 MKL_MALLOC
 integer*8 alloc_size     
! #endif

sizeint = 8 

!Initialize the local matrices
M1(1,:) = (/1.0d00,0.0d00,2.0d00,4.0d00/)
M1(2,:) = (/2.0d00,4.0d00,0.0d00,6.0d00/)
M1(3,:) = (/9.0d00,8.0d00,0.0d00,0.0d00/)
M1(4,:) = (/0.0d00,0.0d00,0.0d00,1.0d00/)

M2(1,:) = (/1.0d00,0.0d00,2.0d00,6.0d00/)
M2(2,:) = (/0.0d00,0.0d00,5.0d00,9.0d00/)
M2(3,:) = (/1.0d00,2.0d00,6.0d00,8.0d00/)
M2(4,:) = (/0.0d00,0.0d00,9.0d00,9.0d00/)

M3(1,:) = (/1.0d00,2.0d00,3.0d00,4.0d00/)
M3(2,:) = (/0.0d00,0.0d00,6.0d00,0.0d00/)
M3(3,:) = (/0.0d00,8.0d00,0.0d00,6.0d00/)
M3(4,:) = (/0.0d00,1.0d00,2.0d00,3.0d00/)

!Print matrices on screen
write(*,*)'M1='
do i = 1,nep
  write(*,'(4f6.3)')(M1(i,j), j=1,nep)
enddo

write(*,*)'M2='
do i = 1,nep
  write(*,'(4f6.3)')(M2(i,j), j=1,nep)
enddo

write(*,*)'M3='
do i = 1,nep
  write(*,'(4f6.3)')(M3(i,j), j=1,nep)
enddo

!Call dense to csr converter from intel mkl library
job(1) = 0 !convert from dense to csr
job(2) = 1 !one based indexing is used for dense matrix
job(3) = 1 !one based indexing to be used for the CSR format
job(4) = 2 !indicates adns is a whole matrix

!call subroutine findnz to find the number of nonzero entries in the matrix
m = nep  ! number of rows in dense matrix
n = nep  ! number of columns in dense matrix allocate(adns(m,n))
allocate(adns(m,n))
adns = M1
call findnz(adns,m,nzmax)
allocate(ia(m+1),ja(nzmax),acsr(nzmax))
job(5) = nzmax !? 

job(6) = 1 ! should be >0 to indicate that all the 3 CSR arrays must be generated
lda = max(1,m)
write(*,*)'calling mkl_ddnscsr'
call mkl_ddnscsr(job, m, n, adns, lda, acsr, ja, ia, info)

write(*,*)'***************************************'
write(*,*)'CSR matrix for M1'
write(*,*)'ia=',(ia(i),i=1,(m+1))
write(*,*)'ja=',(ja(i),i=1,nzmax)
write(*,*)'acsr=',(acsr(i),i=1,nzmax)
write(*,*)'info=',info
write(*,*)'***************************************'

if (info.ne.0) then
    write(*,*)'Error:dns to csr failed for M1'
    stop
endif


!DNS to CSR for M2
!call subroutine findnz to find the number of nonzero entries in the matrix
adns = M2
call findnz(adns,m,nzmax)
allocate(ib(m+1),jb(nzmax),bcsr(nzmax))
job(5) = nzmax !? 

job(6) = 1 ! should be >0 to indicate that all the 3 CSR arrays must be generated
lda = max(1,m)
write(*,*)'calling mkl_ddnscsr'
call mkl_ddnscsr(job, m, n, adns, lda, bcsr, jb, ib, info)

write(*,*)'***************************************'
write(*,*)'CSR matrix for M2'
write(*,*)'ib=',(ib(i),i=1,(m+1))
write(*,*)'jb=',(jb(i),i=1,nzmax)
write(*,*)'bcsr=',(bcsr(i),i=1,nzmax)
write(*,*)'info=',info
write(*,*)'***************************************'

if (info.ne.0) then
    write(*,*)'Error:dns to csr failed for M2'
    stop
endif


!DNS to CSR for M3
!call subroutine findnz to find the number of nonzero entries in the matrix
adns = M3
call findnz(adns,m,nzmax)
allocate(ic(m+1),jc(nzmax),ccsr(nzmax))
job(5) = nzmax !? 

job(6) = 1 ! should be >0 to indicate that all the 3 CSR arrays must be generated
lda = max(1,m)
write(*,*)'calling mkl_ddnscsr'
call mkl_ddnscsr(job, m, n, adns, lda, ccsr, jc, ic, info)

write(*,*)'***************************************'
write(*,*)'CSR matrix for M3'
write(*,*)'ic=',(ic(i),i=1,(m+1))
write(*,*)'jc=',(jc(i),i=1,nzmax)
write(*,*)'ccsr=',(ccsr(i),i=1,nzmax)
write(*,*)'info=',info
write(*,*)'***************************************'
     
if (info.ne.0) then
    write(*,*)'Error:dns to csr failed for M3'
    stop
endif


!Initialize intma_G
intma_G(1,:) = (/1,2,3,4/)
intma_G(2,:) = (/3,4,5,6/)
intma_G(3,:) = (/5,6,1,2/)

!Form the assembled global matrix just to confirm the results in the end
MG = 0.0d00
do i = 1,nelem

    if(i.eq.1) then
          ME = M1
    elseif(i.eq.2) then
          ME = M2
    elseif(i.eq.3) then
          ME = M3
    else
          write(*,*)'Limited Capability Program! Not intended for meshes &
                     with more than 3 elements! Here nelem = ',nelem
    endif

    MG(intma_G(i,:),intma_G(i,:)) = MG(intma_G(i,:),intma_G(i,:)) + ME

enddo

! Print assembled matrix
write(*,*)'MG = '
do i = 1,nnodes
write(*,'(6f6.3)')(MG(i,j),j = 1,nnodes)
enddo

!DNS to CSR for MG
deallocate(adns)
m = nnodes  ! number of rows in dense matrix
n = nnodes  ! number of columns in dense matrix allocate(adns(m,n))
allocate(adns(m,n))
adns = MG
call findnz(adns,m,nzmax)
allocate(iG(m+1),jG(nzmax),Gcsr(nzmax))
job(5) = nzmax ! 

write(*,*)'nz for MG=',nzmax 

job(6) = 1 ! should be >0 to indicate that all the 3 CSR arrays must be generated
lda = max(1,m)
write(*,*)'calling mkl_ddnscsr'
call mkl_ddnscsr(job, m, n, adns, lda, Gcsr, jG, iG, info)

write(*,*)'***************************************'
write(*,*)'CSR matrix for MG'
write(*,*)'iG=',(iG(i),i=1,(m+1))
write(*,*)'jG=',(jG(i),i=1,nzmax)
write(*,*)'Gcsr=',(Gcsr(i),i=1,nzmax)
write(*,*)'info=',info
write(*,*)'***************************************'

if (info.ne.0) then
    write(*,*)'Error:dns to csr failed for MG'
    stop
endif

! MATRIX ADDITION 1 START                                               
! Calculation of number of non-zero entries                                    
alloc_size = sizeint * (m+1)                                                   
iadd_ptr = mkl_malloc(alloc_size,128)                                          
if (iadd_ptr.eq.0) then
write(*,*)'cannot allocate pointer array of length',m+1
stop 1
 endif

trans = 'N'
request = 1
!sort,idum and ddum are not initialized 
beta = 1.0d00
m = nep
n = nep

call mkl_dcsradd(trans,request , sort, m, n, acsr, ja, ia, beta, bcsr, jb, ib, &
                 ddum, idum, iadd, idum, idum)

nzmax = iadd(m+1) - 1                       
write(*,*)'nzmax_add1 =',nzmax
alloc_size = sizeint * nzmax
jadd_ptr = mkl_malloc(alloc_size,128)
if (jadd_ptr.eq.0) then
    write(*,*)'cannot allocate pointer array of length',nzmax
    stop 1
endif

alloc_size = 8 * nzmax
csradd_ptr = mkl_malloc(alloc_size,128)
if (csradd_ptr.eq.0) then
write(*,*)'cannot allocate pointer array of length',nzmax
stop 1
endif

request = 2

call mkl_dcsradd(trans,request , sort, m, n, acsr, ja, ia, beta, bcsr, jb, ib, &
                 csradd, jadd, iadd, nzmax, info)

if (info.ne.0) then                       
    write(*,*)'den to csr failed'
    stop
endif

!mkl_free - Frees the aligned memory buffer allocated by mkl_malloc.
call mkl_free(iadd_ptr)
call mkl_free(jadd_ptr)
call mkl_free(csradd_ptr)

write(*,*)'***************************************'
write(*,*)'CSR matrix for Madd'
write(*,*)'iadd=',(iadd(i),i=1,(m+1))
write(*,*)'jadd=',(jadd(i),i=1,nzmax)
write(*,*)'csradd=',(csradd(i),i=1,nzmax)
write(*,*)'info=',info
write(*,*)'***************************************'

! *************************************************PROBLEM 1*
 write(*,*)'size of iadd=',size(iadd,1)                    !*   
 write(*,*)'size of jadd=',size(jadd,1)                    !*
 write(*,*)'size of csradd=',size(csradd,1)                !*
 ! **********************************************************

 
 ! MATRIX ADDITION 1 END

! MATRIX ADDITION 2 START
! add M3 (i.e. ic,jc,c) to global matrix (i.e. iadd,jadd,csradd) 

 alloc_size = sizeint * (m+1)                                                   
 iaux_ptr = mkl_malloc(alloc_size,128)                                          
 if (iaux_ptr.eq.0) then
     write(*,*)'cannot allocate pointer array of length',m+1
     stop 1
 endif

trans = 'N'
request = 1
sort = 0
idum = 0
ddum = 0 
beta = 1.0d00
m = nep
n = nep


!call mkl_dcsradd(trans,request , sort, m, n, ccsr, jc, ic, beta, csradd, jadd, iadd, &
!                 ddum, idum, iaux, idum, idum)
call mkl_dcsradd(trans,request , sort, m, n, ccsr, jc, ic, beta, csradd, jadd, iadd, &
                 ddum, idum, iadd, idum, idum)

! nzmax = iaux(m+1) - 1     
! write(*,*)'nzmax=',nzmax   
! write(*,*)'iaux(m+1)',iaux(m+1)
! write(*,*)'iaux=',iaux   
 
      nzmax = iadd(m+1) - 1     
      write(*,*)'nzmax=',nzmax   
      write(*,*)'iadd(m+1)',iadd(m+1)
      write(*,*)'iadd=',iadd  
alloc_size = sizeint * nzmax
jadd_ptr = mkl_malloc(alloc_size,128)
if (jadd_ptr.eq.0) then
    write(*,*)'cannot allocate pointer array of length',nzmax
    stop 1
endif

alloc_size = 8 * nzmax
csradd_ptr = mkl_malloc(alloc_size,128)
if (csradd_ptr.eq.0) then
    write(*,*)'cannot allocate pointer array of length',nzmax
    stop 1
endif

request = 2

call mkl_dcsradd(trans,request , sort, m, n, ccsr, jc, ic, beta, csradd, jadd, iadd, &
                 csradd, jadd, iadd, nzmax, info)

if (info.ne.0) then                       
    write(*,*)'den to csr failed'
    stop
endif

!mkl_free - Frees the aligned memory buffer allocated by mkl_malloc.
!      call mkl_free(iadd_ptr)
!      call mkl_free(jadd_ptr)
!      call mkl_free(csradd_ptr)

write(*,*)'***************************************'
write(*,*)'CSR matrix for Madd2'
write(*,*)'iadd=',(iadd(i),i=1,(m+1))
write(*,*)'jadd=',(jadd(i),i=1,nzmax)
write(*,*)'csradd=',(csradd(i),i=1,nzmax)
write(*,*)'info=',info
write(*,*)'***************************************'

! MATRIX ADDITION 2 END
end program 


subroutine findnz(mat,matsize,nonzero)
implicit none

integer(kind=4)    :: matsize,nonzero
double precision   :: mat(matsize,matsize)

!Local
integer(kind=4)    :: nzi,nzj

!Find the number of non-zero entries in matden
!NUMBER OF NON-ZERO ENTRIES = nonzero
nonzero = 0
do nzi = 1,matsize
do nzj = 1,matsize
       if (mat(nzi,nzj).ne.0.0d00 ) then
           nonzero = nonzero + 1 
       endif  

enddo ! loop nzj
enddo ! loop nzi
end subroutine findnz

