module sparse_matmult_fly_sym 
! sparse matrix multiply and its adjoint Ax=y and A'y=x
use adj_mod

implicit none
real    ,private,  dimension (:), pointer  :: bb
integer(kind=8) ,private                   :: nhxx,nhyy,naz,nax,nay,nhx,nhy,nz,nx,ny
!nhxx-num. of off-diag. x-sub-offset
!nhyy-num. of off-diag. y-sub-offset
!naz -num. of off-diag in z -->depth
!nax -num. of off-diag in x
!nay -num. of off-diag in y
!nhx -num. model space x sub-offset
!nhy -num. model space y sub-offset
!nx  -num. model space x 
!ny  -num. model space y 
!nz  -num. model space z

contains

  subroutine smatmult_init (bb_in,nhxx_in,nhyy_in,naz_in,nax_in,nay_in,nhx_in,nhy_in,nz_in,nx_in,ny_in)
    real , 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 ( adj, add, x, y) result(stat)
  function smatmult_lop ( adj, add) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    call smatmult_lop3(adj,add)
    stat=0
  end function 

subroutine smatmult_lop3(adj,add)
    logical,intent(in)   :: adj,add 
    integer(kind=8)      :: irow,icol,iii 
    integer              :: iaz,iax,iay,ihxx,ihyy,iz,ix,iy,ihx,ihy,ixx,iyy,izz
    real,parameter       :: PI=3.14159265
if ( adj) then
   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))
             write(0,*)icol
             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( 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
           enddo
           enddo
         enddo
       enddo
     enddo
   enddo
   enddo
   enddo
   enddo
   enddo
else
do iy =1,ny
        write(0,*)iy
        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
                                        !if(modulo(irow,2).ne.0) write(0,*)irow
                                        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
                                                                                if(modulo(icol,10).ne.110)then
                                                                                        write(0,*)'iyy,ixx,izz,ihyy,ihxx,icol,row'
                                                                                        write(0,*)iyy,ixx,izz,ihyy,ihxx,icol
                                                                                        write(0,*)nax,naz
                                                                                 endif
                                                                                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( 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
                                                                        enddo
                                                                enddo
                                                        enddo
                                                enddo
                                        enddo
                                        STOP
                                enddo
                        enddo
                enddo
        enddo
enddo
endif

  end subroutine

  subroutine smatmult_close()
  end subroutine 

end module      
