module sparse_matmult_fly_sym_dp 
! sparse matrix multiply and its adjoint Ax=y and A'y=x
use DataModule, only : dp
use adj_mod_dp


implicit none
real    (dp),private,  dimension (:), pointer  :: bb
integer(kind=8) ,private                       :: nhxx,nhyy,naz,nax,nay,nhx,nhy,nz,nx,ny

contains

  subroutine smatmult_init_dp (bb_in,nhxx_in,nhyy_in,naz_in,nax_in,nay_in,nhx_in,nhy_in,nz_in,nx_in,ny_in)
    real (dp), dimension (:) , pointer :: bb_in 
    integer                            :: nhxx_in,nhyy_in,naz_in,nax_in,nay_in,nhx_in,nhy_in,nz_in,nx_in,ny_in

    bb => bb_in
    nhxx=nhxx_in
    nhyy=nhyy_in
    naz=naz_in
    nax=nax_in
    nay=nay_in
    nhx=nhx_in
    nhy=nhy_in
    nz=nz_in
    nx=nx_in
    ny=ny_in

  end subroutine

  function smatmult_lop_dp ( adj, add, x, y) result(stat)
    integer                 :: stat 
    logical,intent(in)      :: adj,add 
    real (dp),dimension(:)  :: x,y 

    call adjnull(       adj, add, x, y)

    call smatmult_lop3(adj,add,x,y )
    stat=0
  end function 

subroutine smatmult_lop3(adj,add,x,y)
    logical,intent(in)        :: adj,add 
    real (dp), dimension (:)  :: x 
    real (dp), dimension (:)  :: y 
    integer(kind=8)           :: irow,icol,iii 
    integer                   :: iaz,iax,iay,ihxx,ihyy,iz,ix,iy,ihx,ihy,ixx,iyy,izz



   do iy =1,ny
   do ix =1,nx
   do iz =1,nz
   do ihy=1,nhy
   do ihx=1,nhx
     irow=ihx+(ihy-1)*nhx         &
             +(iz -1)*nhx*nhy     &
             +(ix -1)*nhx*nhy*nz  &
             +(iy -1)*nhx*nhy*nz*nx
     irow=max(1,min(irow,nhx*nhy*nz*nx*ny))
     do iay =1,nay
       iyy=iy-(nay-1)/2+iay-1
       do iax =1,nax
         ixx=ix-(nax-1)/2+iax-1
         do iaz =1,naz
           izz=iz-(naz-1)+iaz-1
           do ihyy=1,nhyy
           do ihxx=1,nhxx
             icol=ihxx+(ihyy-1)*nhxx          &
                      +(izz -1)*nhxx*nhyy     &
                      +(ixx -1)*nhxx*nhyy*nz  &
                      +(iyy -1)*nhxx*nhyy*nz*nx
             icol=max(1,min(icol,nhx*nhy*nz*nx*ny))
             iii=ihxx+(ihyy-1)*nhxx                              &                      
                     +(iaz -1)*nhxx*nhyy                         &
                     +(iax -1)*nhxx*nhyy*naz                     &
                     +(iay -1)*nhxx*nhyy*naz*nax                 &
                     +(ihx -1)*nhxx*nhyy*naz*nax*nay             &
                     +(ihy -1)*nhxx*nhyy*naz*nax*nay*nhx         &
                     +(iz  -1)*nhxx*nhyy*naz*nax*nay*nhx*nhy     &
                     +(ix  -1)*nhxx*nhyy*naz*nax*nay*nhx*nhy*nz  &
                     +(iy  -1)*nhxx*nhyy*naz*nax*nay*nhx*nhy*nz*nx
             if ( adj) then
               if( irow == icol) then
                 x(icol)=x(icol)+bb(iii)*y(irow)
               else
                 x(icol)=x(icol)+bb(iii)*y(irow)
                 x(irow)=x(irow)+bb(iii)*y(icol)
               endif
             else
               if( irow == icol) then
                 y(irow)=y(irow)+bb(iii)*x(icol)
               else
                 y(irow)=y(irow)+bb(iii)*x(icol)
                 y(icol)=y(icol)+bb(iii)*x(irow)
               endif
             endif
           enddo
           enddo
         enddo
       enddo
     enddo
   enddo
   enddo
   enddo
   enddo
   enddo
  end subroutine

end module      
