module weiop_cig
  !! compute CIGs
  !!        offset domain (weihcig)
  !! ray parameter domain (weipcig)

  use wei_util
  use wei_pro
  use wei_par
  use wei_flags
  use wei_tds
  !---------------
  implicit none

  type(axis) :: akhx, akhy
  type(axis) :: aphx, aphy

  integer(kind=dbl),  private :: planF !! forward FFT
  integer(kind=dbl),  private :: planI !! inverse FFT
  real,               private :: scale

  real,               private :: pcigeps
  integer,            private :: jhx,jhy !! half nhx,nhy - even # of samples

  !! temp FFT arrays
  complex, dimension(:,:,:), allocatable, private :: hht
  complex, dimension(:,:),       pointer, private :: sftF, sftI !! shared

  !! mapping arrays
  integer, dimension(:,:),       pointer, private :: xilk, yilk !! index to lower kh
  logical, dimension(:,:),       pointer, private :: xobf, yobf !! off-bounds flag
  real,    dimension(:,:),       pointer, private :: xwui, ywui !! weight to upper kh

  !! tmp arrays
  complex, dimension(:,:),       pointer, private :: hxtmp, hytmp
  complex, dimension(:,:),       pointer, private :: pxtmp, pytmp

  !! tridiagonal matrices
  complex, dimension(:,:),       pointer, private :: xdiag, xoffd
  complex, dimension(:,:),       pointer, private :: ydiag, yoffd
contains
  !----------------------------------------------------------------  
  !! h-domain CIG
  integer function weiop_hcig_init() result(st)
    character(len=128), parameter :: unit='weiop_hcig_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif
    !! do nothing
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_hcig_init
  !----------------------------------------------------------------
  integer function weiop_hcig(wfld,iws,ith) result(st)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in) ::iws,ith
#ifdef DEBUG
    character(len=128), parameter :: unit='weiop_hcig'
    if(wpro%dbg>9) call in(unit)
#endif
    !! do nothing
#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
    st=weiOK
  end function weiop_hcig
  !----------------------------------------------------------------  
  !----------------------------------------------------------------  
  !! p-domain CIG
  integer function weiop_pcig_init() result(st)
    integer, dimension(3) :: ni
    logical               :: fftmeasure
    integer               :: fftsetup
    integer               :: er
    integer               :: ipx,ipy, jpx,jpy, ihx,ihy
    real                  :: arg
    integer :: mpierr
    character(len=128), parameter :: unit='weiop_pcig_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    !! Note: the offset axes should have an even number of samples
    jhx = wpro%ahx%n/2
    jhy = wpro%ahy%n/2

    if(wpar%impi==0) call from_param("pcigeps",pcigeps,1.0)
    if(wpar%nmpi >0) call MPI_BCAST(pcigeps, 1,MPI_REAL,0,MPI_COMM_WORLD,mpierr)

    !! ----- init FFTs -----
    if(wpar%impi==0) call from_param("fftmeasure",fftmeasure,.false.)
    if(wpar%nmpi >0) call MPI_BCAST(fftmeasure, 1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr)
    if(fftmeasure) then
       fftsetup = FFTW_MEASURE
    else
       fftsetup = FFTW_ESTIMATE
    end if

    ni = (/ wpro%ahx%n, wpro%ahy%n, 1 /)
    scale = 1./sqrt(1.*product(ni))
    !! ----- init FFT plans -----
    call fftwnd_f77_create_plan(planF,3,ni,FFTW_FORWARD  &
    ,         fftsetup + FFTW_THREADSAFE + FFTW_IN_PLACE)
    call fftwnd_f77_create_plan(planI,3,ni,FFTW_BACKWARD &
    ,         fftsetup + FFTW_THREADSAFE + FFTW_IN_PLACE)
    !! ----- FFT tmp array -----
    call weimem_add(8,ni(1)*ni(2)*wpar%nomp,unit)
    allocate(hht(     ni(1),ni(2),wpar%nomp),stat=er)
    if(er/=0) call seperr('cant allocate hht')
    !! ----- init FFT shift -----
    call weimem_add(8,wpro%ahx%n*wpro%ahy%n,unit)
    allocate(sftF(    wpro%ahx%n,wpro%ahy%n), stat=er)
    if(er /=0) call seperr("cant allocate sftF")
    call weimem_add(8,wpro%ahx%n*wpro%ahy%n,unit)
    allocate(sftI(    wpro%ahx%n,wpro%ahy%n), stat=er)
    if(er /=0) call seperr("cant allocate sftI")

    do ihx=1,wpro%ahx%n
       do ihy=1,wpro%ahy%n
          arg = 2*wpro%pi &
          * ( (ihx-1.)/wpro%ahx%n * wpro%ahx%o/wpro%ahx%d &
          +   (ihy-1.)/wpro%ahy%n * wpro%ahy%o/wpro%ahy%d )

          sftF(ihx,ihy) = cmplx(cos(arg),-(+1)*sin(arg))
          sftI(ihx,ihy) = cmplx(cos(arg),-(-1)*sin(arg))
       end do
    end do

    !! ----- tmp arrays ------
    call weimem_add(8,wpro%ahx%n*wpar%nomp,unit)
    allocate(hxtmp(   wpro%ahx%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate hxtmp") !! hxtmp(nhx,nomp)

    call weimem_add(8,wpro%ahy%n*wpar%nomp,unit)
    allocate(hytmp(   wpro%ahy%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate hytmp") !! hytmp(nhy,nomp)

    call weimem_add(8,wpro%apx%n*wpar%nomp,unit)
    allocate(pxtmp(   wpro%apx%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate pxtmp") !! pxtmp(npx,nomp)

    call weimem_add(8,wpro%apy%n*wpar%nomp,unit)
    allocate(pytmp(   wpro%apy%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate pytmp") !! pytmp(npy,nomp)

    !! ------ stretch arrays ------
    call weimem_add(1,wpro%apx%n*wpar%nomp,unit)
    allocate(xilk(    wpro%apx%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate xilk") !! xilk(npx,nomp)

    call weimem_add(1,wpro%apx%n*wpar%nomp,unit)
    allocate(xobf(    wpro%apx%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate xobf") !! xobf(npx,nomp)

    call weimem_add(4,wpro%apx%n*wpar%nomp,unit)
    allocate(xwui(    wpro%apx%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate xwui") !! xwui(npx,nomp)


    call weimem_add(1,wpro%apy%n*wpar%nomp,unit)
    allocate(yilk(    wpro%apy%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate yilk") !! yilk(npy,nomp)

    call weimem_add(1,wpro%apy%n*wpar%nomp,unit)
    allocate(yobf(    wpro%apy%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate yobf") !! yobf(npy,nomp)

    call weimem_add(4,wpro%apy%n*wpar%nomp,unit)
    allocate(ywui(    wpro%apy%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate ywui") !! ywui(npy,nomp)

    !! ----- tridiagonal solver -----
    call weimem_add(8,wpro%apx%n*wpar%nomp,unit)
    allocate(xdiag(   wpro%apx%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate xdiag") !! xdiag(npx,nomp)

    call weimem_add(8,wpro%apx%n*wpar%nomp,unit)
    allocate(xoffd(   wpro%apx%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate xoffd") !! xoffd(npx,nomp)

    call weimem_add(8,wpro%apy%n*wpar%nomp,unit)
    allocate(ydiag(   wpro%apy%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate ydiag") !! ydiag(npy,nomp)

    call weimem_add(8,wpro%apy%n*wpar%nomp,unit)
    allocate(yoffd(   wpro%apy%n,wpar%nomp), stat=er)
    if(er /=0) call seperr("cant allocate yoffd") !! yoffd(npy,nomp)

    !! ----- init tridiagonal solver -----
    call weitds_init(wpro%apx%n,wpro%apy%n)

    !! ----- init kh and ph axes -----
    !! offset x axes (khx,phx)
    akhx%n = wpro%ahx%n
    akhx%o = -   wpro%pi /  wpro%ahx%d
    akhx%d = 2 * wpro%pi / (wpro%ahx%d * wpro%ahx%n)
    if(wpro%ahx%n ==1) then
       akhx%o=0.
       akhx%d=1.
    end if
    aphx = wpro%apx

    !! offset y axes (khy,phy)
    akhy%n = wpro%ahy%n
    akhy%o = -   wpro%pi /  wpro%ahy%d
    akhy%d = 2 * wpro%pi / (wpro%ahy%d * wpro%ahy%n)
    if(wpro%ahy%n ==1) then
       akhy%o=0.
       akhy%d=1.
    end if
    aphy = wpro%apy

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_pcig_init
  !----------------------------------------------------------------
  integer function weiop_pcig(wfld,iws,ith) result(st)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in) ::iws,ith
    character(len=128), parameter :: unit='weipcig'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    !! interpolator map (A)
    call weipcig_map(iws,ith, aphx,akhx, xobf,xilk,xwui)
    call weipcig_map(iws,ith, aphy,akhy, yobf,yilk,ywui)

    if(wpro%adj) then !! MIGRATION
       call weipcig_fft(wfld,ith) !!            kh = Ft h
       call weipcig_h2p(wfld,ith) !! ph = At Mt kh
    else             !! MODELING
       call weipcig_p2h(wfld,ith) !!         kh = A  M  ph
       call weipcig_ift(wfld,ith) !!  h = F  kh
    end if

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
    st=weiOK
  end function weiop_pcig
  !----------------------------------------------------------------  
  !! ph to kh conversion
  subroutine weipcig_p2h(wfld,ith)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in) :: ith
    integer             :: imx,imy, ihx,ihy, ipx,ipy
    integer             :: il,iu
    real                :: wl,wu 
    character(len=128), parameter :: unit='weipcig_p2h'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    xp2h: if(wpro%apx%n >1) then
       !!setup tridiagonal solver
       call weipcig_i3d(ith,wpro%apx%n,xdiag,xoffd,xobf,xilk,xwui)

       do imx=1,wpro%amx%n
          do imy=1,wpro%amy%n
             do ihy=1,wpro%ahy%n
                hxtmp(:,ith) = 0.
                pxtmp(:,ith) = wfld(imx,imy,:,ihy,ith)

                !! solve (At A + e^2 Rt R ) ph = ph 
                call weitdsx(ith,xdiag,xoffd,pxtmp)
                do ipx=1,wpro%apx%n
                   if(xobf(ipx,ith)) cycle
                   il=xilk(ipx,ith); iu= il +  1
                   wu=xwui(ipx,ith); wl= 1. - wu 

                   !! kh = A ph
                   hxtmp(il,ith) = hxtmp(il,ith) + pxtmp(ipx,ith) * wl
                   hxtmp(iu,ith) = hxtmp(iu,ith) + pxtmp(ipx,ith) * wu
                end do

                wfld(imx,imy,:,ihy,ith) = hxtmp(:,ith)
             end do
          end do
       end do
    end if xp2h

    yh2p: if(wpro%apy%n >1) then
       !!setup tridiagonal solver
       call weipcig_i3d(ith,wpro%apy%n,ydiag,yoffd,yobf,yilk,ywui)

       do imx=1,wpro%amx%n
          do imy=1,wpro%amy%n
             do ihx=1,wpro%ahx%n
                hytmp(:,ith) = 0.
                pytmp(:,ith) = wfld(imx,imy,ihx,:,ith)

                !! solve (At A + e^2 Rt R ) ph = ph
                call weitdsy(ith,ydiag,yoffd,pytmp)
                do ipy=1,wpro%apy%n
                   if(yobf(ipy,ith)) cycle
                   il=yilk(ipy,ith); iu= il +  1
                   wu=ywui(ipy,ith); wl= 1. - wu 

                   !! kh = A ph
                   hytmp(il,ith) = hytmp(il,ith) + pytmp(ipy,ith) * wl
                   hytmp(iu,ith) = hytmp(iu,ith) + pytmp(ipy,ith) * wu
                end do

                wfld(imx,imy,ihx,:,ith) = pytmp(:,ith)
             end do
          end do
       end do

    end if yh2p

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_p2h
  !----------------------------------------------------------------
  !! kh to ph conversion
  subroutine weipcig_h2p(wfld,ith)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in) :: ith
    integer             :: imx,imy, ihx,ihy, ipx,ipy
    integer             :: il,iu
    real                :: wl,wu
    character(len=128), parameter :: unit='weipcig_h2p'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    xh2p: if(wpro%apx%n >1) then
       !!setup tridiagonal solver
       call weipcig_i3d(ith,wpro%apx%n,xdiag,xoffd,xobf,xilk,xwui)

       do imx=1,wpro%amx%n
          do imy=1,wpro%amy%n
             do ihy=1,wpro%ahy%n
                hxtmp(:,ith) = wfld(imx,imy,:,ihy,ith)
                pxtmp(:,ith) = 0.

                do ipx=1,wpro%apx%n
                   if(xobf(ipx,ith)) cycle
                   il=xilk(ipx,ith); iu= il +  1
                   wu=xwui(ipx,ith); wl= 1. - wu 

                   !! ph = At kh
                   pxtmp(ipx,ith) = pxtmp(ipx,ith) + hxtmp(il,ith) * wl
                   pxtmp(ipx,ith) = pxtmp(ipx,ith) + hxtmp(iu,ith) * wu
                end do
                !! solve (At A + e^2 Rt R ) ph = At kh 
                call weitdsx(ith,xdiag,xoffd,pxtmp)

                wfld(imx,imy,:,ihy,ith) = pxtmp(:,ith)
             end do
          end do
       end do
    end if xh2p

    yh2p: if(wpro%apy%n >1) then
       !!setup tridiagonal solver
       call weipcig_i3d(ith,wpro%apy%n,ydiag,yoffd,yobf,yilk,ywui)

       do imx=1,wpro%amx%n
          do imy=1,wpro%amy%n
             do ihx=1,wpro%ahx%n
                hytmp(:,ith) = wfld(imx,imy,ihx,:,ith)
                pytmp(:,ith) = 0.

                do ipy=1,wpro%apy%n
                   if(yobf(ipy,ith)) cycle
                   il=yilk(ipy,ith); iu= il +  1
                   wu=ywui(ipy,ith); wl= 1. - wu 

                   !! ph = At kh
                   pytmp(ipy,ith) = pytmp(ipy,ith) + hytmp(il,ith) * wl
                   pytmp(ipy,ith) = pytmp(ipy,ith) + hytmp(iu,ith) * wu
                end do
                !! solve (At A + e^2 Rt R ) ph = At kh 
                call weitdsy(ith,ydiag,yoffd,pytmp)

                wfld(imx,imy,ihx,:,ith) = pytmp(:,ith)
             end do
          end do
       end do
    end if yh2p
#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_h2p
  !----------------------------------------------------------------
  !! PULL stretch mapping - definition of A and At
  !! (emphasize migration over modeling) 
  !! loop over ph; map every ph into kh
  subroutine weipcig_map( iws,ith, aph,akh, obf,ilk,wui)
    integer,    intent(in) ::iws,ith
    type(axis), intent(in) :: aph, akh
    logical, dimension(:,:), pointer :: obf
    integer, dimension(:,:), pointer :: ilk
    real,    dimension(:,:), pointer :: wui
    integer :: ip
    real    :: p2k,w,ph
    character(len=128), parameter :: unit='weipcig_map'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    w  = wpro%bw%o + (iws-1)*wpro%bw%d

    p2k=0.
    obf(:,ith)=.false.
    do ip=1,aph%n; ph = aph%o + (ip-1) * aph%d
       p2k = ph * w  !! kh = f(ph) = ph * w

       !! ILK = Index to Lower kh
       ilk(ip,ith) = floor(1 + (p2k - akh%o)/akh%d)

       !! WUI = Weight to Upper Index in the kh array
       wui(ip,ith) = p2k - (akh%o + (ilk(ip,ith)-1)*akh%d)
       wui(ip,ith) = wui(ip,ith)/akh%d

    end do
    where(ilk(:,ith)<1 .or. ilk(:,ith)>akh%n-1)
       !! OBF = Off-Bounds Flag
       obf(:,ith) = .true.
    end where

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_map
  !----------------------------------------------------------------
  !! init tridiagonal matrix
  subroutine weipcig_i3d(ith,n,diag,offd,obf,ilk,wui)
    integer, intent(in) :: n, ith
    complex, dimension(:,:), pointer :: diag, offd
    logical, dimension(:,:), pointer :: obf
    integer, dimension(:,:), pointer :: ilk
    real,    dimension(:,:), pointer :: wui
    integer :: i, il,iu
    real    ::    wl,wu
    character(len=128), parameter :: unit='weipcig_i3d'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    diag(:,ith) = 0.
    offd(:,ith) = 0.
    do i=1,n
       if(obf(i,ith)) cycle
       il=ilk(i,ith); iu= il +  1
       wu=wui(i,ith); wl= 1. - wu

       !! (At A)
       diag(i,ith) = wl*wl + wu*wu
    end do
    !! (At A + e^2 Rt R)
    diag(:,ith) = diag(:,ith) + 2 * pcigeps*pcigeps
    offd(:,ith) =             -     pcigeps*pcigeps
#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_i3d
  !----------------------------------------------------------------
  subroutine weipcig_fft(wfld,ith)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in) :: ith
    integer :: imx,imy!!, ihx,ihy
    character(len=128), parameter :: unit='weipcig_fft'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif
    !! FFT on hx and hy
    do imx=1,wpro%amx%n
       do imy=1,wpro%amy%n
          hht(:,:,ith) = wfld(imx,imy,:,:,ith)
          call fftwnd_f77_one(planF,hht(:,:,ith),0)
          wfld(imx,imy,:,:,ith) = hht(:,:,ith) * scale
       end do
    end do

    call weipcig_flp(wfld,ith)
    call weipcig_sft(wfld,ith,sftF)

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_fft
  !----------------------------------------------------------------
  subroutine weipcig_ift(wfld,ith)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in) ::ith
    integer :: imx,imy!!, ihx,ihy
    character(len=128), parameter :: unit='weipcig_ift'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    call weipcig_sft(wfld,ith,sftI)
    call weipcig_flp(wfld,ith)

    !! IFT on hx and hy
    do imx=1,wpro%amx%n
       do imy=1,wpro%amy%n
          hht(:,:,ith) = wfld(imx,imy,:,:,ith)
          call fftwnd_f77_one(planI,hht(:,:,ith),0)
          wfld(imx,imy,:,:,ith) = hht(:,:,ith) * scale
       end do
    end do
#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_ift
  !----------------------------------------------------------------
  !! origin shift
  subroutine weipcig_sft(wfld,ith,sft)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    complex, dimension(:,:),       pointer :: sft
    integer,                    intent(in) ::ith
    integer :: imx,imy
    character(len=128), parameter :: unit='weipcig_sft'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif
    do imx=1,wpro%amx%n
       do imy=1,wpro%amy%n
          wfld(imx,imy,:,:,ith) = &
          wfld(imx,imy,:,:,ith) * sft
       end do
    end do
#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_sft
  !----------------------------------------------------------------
  subroutine weipcig_flp(wfld,ith)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in)::ith
    integer :: imx,imy, ihx,ihy
    character(len=128), parameter :: unit='weipcig_flp'
#ifdef DEBUG
    if(wpro%dbg>9) call in(unit)
#endif

    do imx=1,wpro%amx%n
       do imy=1,wpro%amy%n

          if(wpro%ahx%n>1) then
             do ihy=1,wpro%ahy%n
                hxtmp(    1:    jhx,ith) = wfld(imx,imy,jhx+1:jhx+jhx,ihy,ith)
                hxtmp(jhx+1:jhx+jhx,ith) = wfld(imx,imy,    1:    jhx,ihy,ith)
                wfld(imx,imy,:,ihy,ith)  = hxtmp(:,ith)
             end do
          end if

          if(wpro%ahy%n>1) then
             do ihx=1,wpro%ahx%n
                hytmp(    1:    jhy,ith) = wfld(imx,imy,ihx,jhy+1:jhy+jhy,ith)
                hytmp(jhy+1:jhy+jhy,ith) = wfld(imx,imy,ihx,    1:    jhy,ith)
                wfld(imx,imy,ihx,:,ith)  = hytmp(:,ith)
             end do
          end if

       end do
    end do

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
  end subroutine weipcig_flp
  !----------------------------------------------------------------
end module weiop_cig

!! ---------------------------------------------------
!!
!!        A        ph  =  kh
!!  . . . . . . .   .      .
!!  . . . . . . .   .      .
!!  . . . . . . .   .      .
!!  . . . . . . .   .      .
!!                  .
!!                  .
!!                  .
!!           A             ph =    kh
!! LSQ:  (At A + e^2 Rt R) ph = At kh
!! 
!!
!! MIGRATION: forward pCIG transformation
!!       ph = [(At A + e^2 Rt R)^-1]  At  kh
!!          = Mt  At  kh
!!
!!  MODELING: adjoint pCIG transformation
!!       kh = A   M   ph
!!          = A   [(At A + e^2 Rt R)^-1]  ph  
!!
!! ---------------------------------------------------
!!
!! weipcig_map - defines the pull mapping (stretch) A
!!
!!
!! weipcig_h2p - pCIG forward transformation
!!               * compute       ph = At kh
!!               * apply weight  ph = Mt ph
!!
!! weipcig_p2h - pCIG adjoint transformation
!!               * apply weight  ph = M  ph
!!               * compute       kh = A  ph 
!!
!! ---------------------------------------------------
