text:
module mod_solver
  use mpi_f08
  use mod_common_mpi, only: ierr
  use mod_types
  use, intrinsic :: iso_c_binding, only: C_PTR
  implicit none
  private
  public init_bc_rhs,init_matrix_3d,create_solver,setup_solver, &
         add_constant_to_diagonal,solve_helmholtz,finalize_matrix,finalize_solver, &
         hypre_solver, add_constant_to_boundary, &
         HYPRESolverSMG,HYPRESolverPFMG,HYPRESolverGMRES,HYPRESolverBiCGSTAB
#if defined(_FFT_X) || defined(_FFT_Y) || defined(_FFT_Z)
  public init_fft_reduction,init_n_2d_matrices,create_n_solvers,setup_n_solvers, &
         add_constant_to_n_diagonals,solve_n_helmholtz_2d, &
         finalize_n_matrices,finalize_n_solvers,init_comm_slab, &
         init_transpose_slab_uneven,transpose_slab,alltoallw, &
         init_n_3d_matrices,solve_n_helmholtz_3d
  type alltoallw
    integer            :: counts
    integer            :: disps
    type(MPI_DATATYPE) :: types
  end type alltoallw
#endif
!  ! parameterized derived types like below would make more sense but
!  ! they are very well supported yet
!  type alltoallw(n) ! Parameterized derived types are not working so well yet!
!    integer, len :: n
!    integer           , dimension(n) :: counts
!    integer           , dimension(n) :: disps
!    type(MPI_DATATYPE), dimension(n) :: types
!  end type alltoallw
  integer, parameter :: HYPRESolverSMG      = 1, &
                        HYPRESolverPFMG     = 2, &
                        HYPRESolverGMRES    = 3, &
                        HYPRESolverBiCGSTAB = 4
  type hypre_solver
    type(C_PTR) :: grid,stencil,precond,solver,mat,rhs,sol
    integer     :: stype,comm_hypre
  end type hypre_solver
  contains
  subroutine init_bc_rhs(cbc,bc,dl,is_bound,is_centered,lo,hi,periods, &
                         dx1,dx2,dy1,dy2,dz1,dz2,rhsx,rhsy,rhsz,bcx,bcy,bcz)
    !
    ! description
    !
    implicit none
    character(len=1)  , intent(in ), dimension(0:1,3) :: cbc
    real(rp)          , intent(in ), dimension(0:1,3) ::  bc
    real(rp)          , intent(in ), dimension(0:1,3) ::  dl
    logical           , intent(in ), dimension(0:1,3) ::  is_bound
    logical           , intent(in ), dimension(    3) ::  is_centered
    integer           , intent(in ), dimension(    3) :: lo,hi,periods
    real(rp)          , intent(in ), target, dimension(lo(1)-1:) :: dx1,dx2
    real(rp)          , intent(in ), target, dimension(lo(2)-1:) :: dy1,dy2
    real(rp)          , intent(in ), target, dimension(lo(3)-1:) :: dz1,dz2
    real(rp)          , intent(out), dimension(lo(2):,lo(3):,0:)    :: rhsx
    real(rp)          , intent(out), dimension(lo(1):,lo(3):,0:)    :: rhsy
    real(rp)          , intent(out), dimension(lo(1):,lo(2):,0:)    :: rhsz
    real(rp)          , intent(in ), dimension(lo(2):,lo(3):,0:), optional :: bcx
    real(rp)          , intent(in ), dimension(lo(1):,lo(3):,0:), optional :: bcy
    real(rp)          , intent(in ), dimension(lo(1):,lo(2):,0:), optional :: bcz
    integer, dimension(3) :: qqq
    real(rp), dimension(0:1,3) :: factor,sgn
    integer :: i,j,k,q,qq,idir,ib
    integer, dimension(3,3) :: eye
    real(rp) :: rhs
    !
    qqq(:) = 0
    where(.not.is_centered(:)) qqq(:) = 1
    factor(:,:) = 0._rp
    sgn(   :,:) = 0._rp
    do q=1,3
      do qq=0,1
        if(is_bound(qq,q)) then
          select case(cbc(qq,q))
          case('N')
            factor(qq,q) = 1._rp*dl(qq,q)*bc(qq,q)
            if(qq == 1) factor(qq,q) = -factor(qq,q)
            sgn(   qq,q) = 1._rp
          case('D')
            if(is_centered(q)) then
              factor(qq,q) = -2._rp*bc(qq,q)
              sgn(   qq,q) = -1._rp
            else
              factor(qq,q) = -1._rp*bc(qq,q)
              sgn(   qq,q) =  0._rp
            end if
          end select
        end if
      end do
    end do
    !
    eye(:,:) = 0
    do idir=1,3
      eye(idir,idir) = 1
    end do
    rhsx(:,:,:) = 0._rp
    rhsy(:,:,:) = 0._rp
    rhsz(:,:,:) = 0._rp
    do idir=1,3
      do k=lo(3),hi(3),max((hi(3)-lo(3))*eye(idir,3),1)
        do j=lo(2),hi(2),max((hi(2)-lo(1))*eye(idir,2),1)
          do i=lo(1),hi(1),max((hi(1)-lo(1))*eye(idir,1),1)
            if(periods(idir) == 0) then
              select case(idir)
              case(1)
                if(     i == lo(idir)) then
                  ib = 0
                else if(i == hi(idir)) then
                  ib = 1
                end if
                if(is_bound(ib,idir)) then
                  rhs = factor(ib,idir)/(dx1(i-(1-ib)+qqq(idir))*dx2(i))
                  if(present(bcx).and.bc(ib,idir)/=0._rp) rhs = rhs*bcx(j,k,ib)/bc(ib,idir)
                  rhsx(j,k,ib) = rhs
                end if
              case(2)
                if(     j == lo(idir)) then
                  ib = 0
                else if(j == hi(idir)) then
                  ib = 1
                end if
                if(is_bound(ib,idir)) then
                  rhs = factor(ib,idir)/(dy1(j-(1-ib)+qqq(idir))*dy2(j))
                  if(present(bcy).and.bc(ib,idir)/=0._rp) rhs = rhs*bcy(i,k,ib)/bc(ib,idir)
                  rhsy(i,k,ib) = rhs
                end if
              case(3)
                if(     k == lo(idir)) then
                  ib = 0
                else if(k == hi(idir)) then
                  ib = 1
                end if
                if(is_bound(ib,idir)) then
                  rhs = factor(ib,idir)/(dz1(k-(1-ib)+qqq(idir))*dz2(k))
                  if(present(bcz).and.bc(ib,idir)/=0._rp) rhs = rhs*bcz(i,j,ib)/bc(ib,idir)
                  rhsz(i,j,ib) = rhs
                end if
              end select
            end if
          end do
        end do
      end do
    end do
  end subroutine init_bc_rhs
  subroutine init_matrix_3d(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo,hi,periods, &
                            dx1,dx2,dy1,dy2,dz1,dz2,alpha,alpha_bc,asolver,lambda,is_bound_outflow)
    !
    ! description
    !
    implicit none
    integer, parameter :: nstencil = 7
    character(len=1)  , intent(in ), dimension(0:1,3) :: cbc
    real(rp)          , intent(in ), dimension(0:1,3) ::  bc
    real(rp)          , intent(in ), dimension(0:1,3) ::  dl
    logical           , intent(in ), dimension(0:1,3) ::  is_bound
    logical           , intent(in )                   ::  is_uniform_grid
    logical           , intent(in ), dimension(    3) ::  is_centered
    integer           , intent(in ), dimension(    3) :: lo,hi,periods
    real(rp)          , intent(in ), target, dimension(lo(1)-1:) :: dx1,dx2
    real(rp)          , intent(in ), target, dimension(lo(2)-1:) :: dy1,dy2
    real(rp)          , intent(in ), target, dimension(lo(3)-1:) :: dz1,dz2
    real(rp)          , intent(in )                              :: alpha
    real(rp)          , intent(in ), dimension(0:1,3)            :: alpha_bc
    type(hypre_solver), intent(out)                              :: asolver
    real(rp)          , intent(in ), optional, dimension(:)      :: lambda
    logical           , intent(in ), optional, dimension(0:1,3)  :: is_bound_outflow
    integer, dimension(3         ) :: qqq
    integer, dimension(3,nstencil) :: offsets
    real(rp), dimension(product(hi(:)-lo(:)+1)*nstencil) :: matvalues
    real(rp), dimension(0:1,3) :: factor,sgn
    type(C_PTR) :: grid,stencil,mat,rhs,sol
    integer :: i,j,k,q,qq
    real(rp) :: cc,cxm,cxp,cym,cyp,czm,czp
    integer            :: comm_hypre
    logical, dimension(0:1,3) :: is_bound_outflow_aux
    !
    comm_hypre = MPI_COMM_WORLD%MPI_VAL
    qqq(:) = 0
    where(.not.is_centered(:)) qqq(:) = 1
    factor(:,:) = 0._rp
    sgn(   :,:) = 0._rp
    do q=1,3
      do qq=0,1
        if(is_bound(qq,q)) then
          select case(cbc(qq,q))
          case('N')
            factor(qq,q) = 1._rp*dl(qq,q)*bc(qq,q)
            if(qq == 1) factor(qq,q) = -factor(qq,q)
            sgn(   qq,q) = 1._rp
          case('D')
            if(is_centered(q)) then
              factor(qq,q) = -2._rp*bc(qq,q)
              sgn(   qq,q) = -1._rp
            else
              factor(qq,q) = -1._rp*bc(qq,q)
              sgn(   qq,q) =  0._rp
            end if
          end select
        end if
      end do
    end do
    if(present(is_bound_outflow)) then
      is_bound_outflow_aux(:,:) = is_bound_outflow(:,:)
    else
      is_bound_outflow_aux(:,:) = .false.
    end if
    !
    ! create 3D grid object
    !
    call HYPRE_StructGridCreate(comm_hypre,3,grid,ierr)
    call HYPRE_StructGridSetPeriodic(grid,periods,ierr)
    call HYPRE_StructGridSetExtents(grid,lo,hi,ierr)
    call HYPRE_StructGridAssemble(grid,ierr)
    !
    ! setup the finite-difference stencil
    !
    call HYPRE_StructStencilCreate(3,nstencil,stencil,ierr)
    offsets = reshape([ 0, 0, 0, &
                       -1, 0, 0, &
                        1, 0, 0, &
                        0,-1, 0, &
                        0, 1, 0, &
                        0, 0,-1, &
                        0, 0, 1 ],shape(offsets))
    do q=1,nstencil
      call HYPRE_StructStencilSetElement(stencil,q-1,offsets(:,q),ierr)
    end do
    !
    ! create coefficient matrix, and solution & right-hand-side vectors
    !
    call HYPRE_StructMatrixCreate(comm_hypre,grid,stencil,mat,ierr)
    if(is_uniform_grid) call HYPRE_StructMatrixSetSymmetric(mat,1,ierr)
    call HYPRE_StructMatrixInitialize(mat,ierr)
    call HYPRE_StructVectorCreate(comm_hypre,grid,sol,ierr)
    call HYPRE_StructVectorInitialize(sol,ierr)
    call HYPRE_StructVectorCreate(comm_hypre,grid,rhs,ierr)
    call HYPRE_StructVectorInitialize(rhs,ierr)
    q = 0
    do k=lo(3),hi(3)
      do j=lo(2),hi(2)
        do i=lo(1),hi(1)
          q = q + 1
          cxm = 1._rp/(dx1(i-1+qqq(1))*dx2(i))
          cxp = 1._rp/(dx1(i  +qqq(1))*dx2(i))
          cym = 1._rp/(dy1(j-1+qqq(2))*dy2(j))
          cyp = 1._rp/(dy1(j  +qqq(2))*dy2(j))
          czm = 1._rp/(dz1(k-1+qqq(3))*dz2(k))
          czp = 1._rp/(dz1(k  +qqq(3))*dz2(k))
#ifdef _FFT_X
          cxm = 0._rp
          cxp = 0._rp
          qq  = i - lo(1) + 1
#elif  _FFT_Y
          cym = 0._rp
          cyp = 0._rp
          qq  = j - lo(2) + 1
#elif  _FFT_Z
          czm = 0._rp
          czp = 0._rp
          qq  = k - lo(3) + 1
#endif
          cc  = -(cxm+cxp+cym+cyp+czm+czp) + alpha
#if defined(_FFT_X) || defined(_FFT_Y) || defined(_FFT_Z)
          cc  = cc + lambda(qq)
#endif
          if(periods(1) == 0) then
            if(is_bound(0,1).and.i == lo(1)) then
              cc = cc + sgn(0,1)*cxm + alpha_bc(0,1)
              cxm = 0._rp
              if( is_bound_outflow_aux(0,1) ) then
                cc = cc + sgn(0,1)*cxp
                cxp = 0._rp
              end if
            end if
            if(is_bound(1,1).and.i == hi(1)) then
              cc = cc + sgn(1,1)*cxp + alpha_bc(1,1)
              cxp = 0._rp
              if( is_bound_outflow_aux(1,1) ) then
                cc = cc + sgn(1,1)*cxm
                cxm = 0._rp
              end if
            end if
          end if
          if(periods(2) == 0) then
            if(is_bound(0,2).and.j == lo(2)) then
              cc = cc + sgn(0,2)*cym + alpha_bc(0,2)
              cym = 0._rp
              if( is_bound_outflow_aux(0,2) ) then
                cc = cc + sgn(0,2)*cyp
                cyp = 0._rp
              end if
            end if
            if(is_bound(1,2).and.j == hi(2)) then
              cc = cc + sgn(1,2)*cyp + alpha_bc(1,2)
              cyp = 0._rp
              if( is_bound_outflow_aux(1,2) ) then
                cc = cc + sgn(1,2)*cym
                cym = 0._rp
              end if
            end if
          end if
          if(periods(3) == 0) then
            if(is_bound(0,3).and.k == lo(3)) then
              cc = cc + sgn(0,3)*czm + alpha_bc(0,3)
              czm = 0._rp
              if( is_bound_outflow_aux(0,3) ) then
                cc = cc + sgn(0,3)*czp
                czp = 0._rp
              end if
            end if
            if(is_bound(1,3).and.k == hi(3)) then
              cc = cc + sgn(1,3)*czp + alpha_bc(1,3)
              czp = 0._rp
              if( is_bound_outflow_aux(1,3) ) then
                cc = cc + sgn(1,3)*czm
                czm = 0._rp
              end if
            end if
          end if
          qq = (q-1)*nstencil
          matvalues(qq+1) = cc
          matvalues(qq+2) = cxm
          matvalues(qq+3) = cxp
          matvalues(qq+4) = cym
          matvalues(qq+5) = cyp
          matvalues(qq+6) = czm
          matvalues(qq+7) = czp
        end do
      end do
    end do
    call HYPRE_StructMatrixSetBoxValues(mat,lo,hi,nstencil, &
                                        [0,1,2,3,4,5,6],matvalues,ierr)
    call HYPRE_StructMatrixAssemble(mat,ierr)
    asolver%grid       = grid
    asolver%stencil    = stencil
    asolver%mat        = mat
    asolver%rhs        = rhs
    asolver%sol        = sol
    asolver%comm_hypre = comm_hypre
  end subroutine init_matrix_3d
  subroutine create_solver(maxiter,maxerror,stype,asolver)
    implicit none
    integer           ,         intent(   in) :: maxiter
    real(rp)          ,         intent(   in) :: maxerror
    integer           ,         intent(   in) :: stype
    type(hypre_solver), target, intent(inout) :: asolver
    type(C_PTR) :: solver,precond
    integer :: precond_id
    !
    ! setup solver
    !
    ! note: this part was taken from the Paris Simulator code
    !       freely available under a GPL license
    !       http://www.ida.upmc.fr/~zaleski/paris
    !
    if      ( stype == HYPRESolverSMG ) then
      call HYPRE_StructSMGCreate(asolver%comm_hypre,solver,ierr)
      call HYPRE_StructSMGSetMaxIter(solver,maxiter,ierr)
      call HYPRE_StructSMGSetTol(solver,maxerror,ierr)
      call hypre_structSMGsetLogging(solver,1,ierr)
      call HYPRE_StructSMGSetPrintLevel(solver,1,ierr)
    else if ( stype == HYPRESolverPFMG ) then
      call HYPRE_StructPFMGCreate(asolver%comm_hypre,solver,ierr)
      call HYPRE_StructPFMGSetMaxIter(solver,maxiter,ierr)
      call HYPRE_StructPFMGSetTol(solver,maxerror,ierr)
      call HYPRE_structPFMGsetLogging(solver,1,ierr)
      call HYPRE_StructPFMGSetPrintLevel(solver,1,ierr)
      call HYPRE_StructPFMGSetRelChange(solver,1,ierr)
      ! Relaxiation Method: 2 is the fastest if symm matrix
      ! 0: Jacobi
      ! 1: Weighted Jacobi (default)
      ! 2: Red/Black Gauss-Seidel (symmetric: RB pre- and post-relaxation)
      ! 3: Red/Black Gauss-Seidel (nonsymmetric: RB pre- and post-relaxation)
      call HYPRE_StructPFMGSetRelaxType(solver,1,ierr)
      call HYPRE_StructPFMGSetNumPreRelax(solver,1,ierr)
      call HYPRE_StructPFMGSetNumPostRelax(solver,1,ierr)
    else if ( stype == HYPRESolverGMRES .or. &
              stype == HYPRESolverBiCGSTAB   ) then
      if      ( stype == HYPRESolverGMRES ) then
        call HYPRE_StructGMRESCreate(asolver%comm_hypre,solver,ierr)
        call HYPRE_StructGMRESSetMaxIter(solver,maxiter,ierr)
        call HYPRE_StructGMRESSetTol(solver,maxerror,ierr)
        !call HYPRE_StructGMRESSetLogging(solver, 1 ,ierr)
      else if ( stype == HYPRESolverBiCGSTAB ) then
        call HYPRE_StructBiCGSTABCreate(asolver%comm_hypre,solver,ierr)
        call HYPRE_StructBiCGSTABSetMaxIter(solver,maxiter,ierr)
        call HYPRE_StructBiCGSTABSetTol(solver,maxerror,ierr)
      end if
      ! Use PFMG as preconditioner
      call HYPRE_StructPFMGCreate(asolver%comm_hypre,precond,ierr)
      call HYPRE_StructPFMGSetMaxIter(precond,10,ierr)
      call HYPRE_StructPFMGSetTol(precond,0._rp,ierr)
      call HYPRE_StructPFMGSetZeroGuess(precond,ierr)
      call HYPRE_StructPFMGSetRelChange(precond,1,ierr)
      call HYPRE_StructPFMGSetRelaxType(precond,2,ierr)
      precond_id = 1   ! Set PFMG as preconditioner
      if      ( stype == HYPRESolverGMRES ) then
        call HYPRE_StructGMRESSetPrecond(solver,precond_id,precond,ierr)
      else if ( stype == HYPRESolverBiCGSTAB ) then
        call HYPRE_StructBiCGSTABSetPrecond(solver,precond_id,precond,ierr)
      end if
      asolver%precond = precond
    end if
    asolver%solver  = solver
    asolver%stype   = stype
  end subroutine create_solver
  subroutine setup_solver(asolver)
    implicit none
    type(hypre_solver), target, intent(inout) :: asolver
    type(C_PTR), pointer :: solver,mat,rhs,sol
    integer    , pointer :: stype
    !
    solver => asolver%solver
    stype  => asolver%stype
    mat    => asolver%mat
    rhs    => asolver%rhs
    sol    => asolver%sol
    !
    ! setup solver
    !
    ! note: this part was taken from the Paris Simulator code
    !       freely available under a GPL license
    !       http://www.ida.upmc.fr/~zaleski/paris
    !
    if      ( stype == HYPRESolverSMG ) then
      call HYPRE_StructSMGSetup(solver,mat,rhs,sol,ierr)
    else if ( stype == HYPRESolverPFMG ) then
      call HYPRE_StructPFMGSetup(solver,mat,rhs,sol,ierr)
    else if ( stype == HYPRESolverGMRES .or. &
             stype == HYPRESolverBiCGSTAB   ) then
      if      ( stype == HYPRESolverGMRES ) then
        call HYPRE_StructGMRESSetup(solver,mat,rhs,sol,ierr)
      else if ( stype == HYPRESolverBiCGSTAB ) then
        call HYPRE_StructBiCGSTABSetup(solver,mat,rhs,sol,ierr)
      end if
    end if
  end subroutine setup_solver
  subroutine add_constant_to_diagonal(lo,hi,alpha,mat)
    implicit none
    integer    , intent(in   ), dimension(3) :: lo,hi
    real(rp)   , intent(in   ) :: alpha
    type(C_PTR), intent(inout) :: mat
    real(rp), dimension(product(hi(:)-lo(:)+1)) :: matvalues
    integer :: i,j,k,q
    q = 0
    do k=lo(3),hi(3)
      do j=lo(2),hi(2)
        do i=lo(1),hi(1)
          q=q+1
          matvalues(q) = alpha
        end do
      end do
    end do
    call HYPRE_StructMatrixAddToBoxValues(mat,lo,hi,1,[0],matvalues,ierr)
  end subroutine add_constant_to_diagonal
  subroutine solve_helmholtz(asolver,lo,hi,p,po)
    implicit none
    type(hypre_solver), target, intent(in   )               :: asolver
    integer           ,         intent(in   ), dimension(3) :: lo,hi
    real(rp)          ,         intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: p,po
    type(C_PTR), pointer :: solver,mat,rhs,sol
    integer    , pointer :: stype
    solver  => asolver%solver
    mat     => asolver%mat
    rhs     => asolver%rhs
    sol     => asolver%sol
    stype   => asolver%stype
    !
    call HYPRE_StructVectorSetBoxValues(rhs,lo,hi,p(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)),ierr)
    call HYPRE_StructVectorAssemble(rhs,ierr)
    !
    ! create soluction vector
    !
    call HYPRE_StructVectorSetBoxValues(sol,lo,hi,po(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)),ierr)
    call HYPRE_StructVectorAssemble(sol,ierr)
    !
    ! setup solver, and solve
    !
    ! note: this part was based on the the Paris Simulator code
    !       freely available under a GPL license; see:
    !       http://www.ida.upmc.fr/~zaleski/paris/
    !
    if      ( stype == HYPRESolverSMG ) then
      call HYPRE_StructSMGSolve(solver,mat,rhs,sol,ierr)
      !call HYPRE_StructSMGGetNumIterations(solver,num_iterations,ierr)
    else if ( stype == HYPRESolverPFMG ) then
      call HYPRE_StructPFMGSolve(solver,mat,rhs,sol,ierr)
      !call HYPRE_StructPFMGGetNumIteration(solver,num_iterations,ierr)
    else if ( stype == HYPRESolverGMRES ) then
      call HYPRE_StructGMRESSolve(solver,mat,rhs,sol,ierr)
      !call HYPRE_StructGMRESGetNumIteratio(solver, num_iterations,ierr)
    else if ( stype == HYPRESolverBiCGSTAB ) then
      call HYPRE_StructBiCGSTABSolve(solver,mat,rhs,sol,ierr)
      !call HYPRE_StructBiCGSTABGetNumItera(solver, num_iterations,ierr)
    end if ! stype
    !
    ! end of part based on the Paris Simulator code
    !
    ! fecth results
    !
    call HYPRE_StructVectorGetBoxValues(sol,lo,hi,p(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)),ierr)
    !$OMP WORKSHARE
    po(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) = p(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3))
    !$OMP END WORKSHARE
  end subroutine solve_helmholtz
  subroutine finalize_solver(asolver)
    implicit none
    type(hypre_solver), target, intent(in) :: asolver
    type(C_PTR), pointer :: precond,solver
    integer    , pointer :: stype
    !
    precond => asolver%precond
    solver  => asolver%solver
    stype   => asolver%stype
    !
    ! note: this part was based on the the Paris Simulator code
    !       freely available under a GPL license; see:
    !       http://www.ida.upmc.fr/~zaleski/paris/
    !
    if      ( stype == HYPRESolverSMG ) then
      call HYPRE_StructSMGDestroy(solver,ierr)
    else if ( stype == HYPRESolverPFMG ) then
      call HYPRE_StructPFMGDestroy(solver,ierr)
    else if ( stype == HYPRESolverGMRES ) then
      call HYPRE_StructGMRESDestroy(solver,ierr)
      call HYPRE_StructPFMGDestroy(precond,ierr)
    else if ( stype == HYPRESolverBiCGSTAB ) then
      call HYPRE_StructBiCGSTABDestroy(solver,ierr)
      call HYPRE_StructPFMGDestroy(precond,ierr)
    end if
  end subroutine finalize_solver
  subroutine finalize_matrix(asolver)
    implicit none
    type(hypre_solver), target, intent(in) :: asolver
    type(C_PTR), pointer :: grid,stencil,mat,rhs,sol
    !
    grid    => asolver%grid
    stencil => asolver%stencil
    mat     => asolver%mat
    rhs     => asolver%rhs
    sol     => asolver%sol
    !
    call HYPRE_StructGridDestroy(grid,ierr)
    call HYPRE_StructStencilDestroy(stencil,ierr)
    call HYPRE_StructMatrixDestroy(mat,ierr)
    call HYPRE_StructVectorDestroy(rhs,ierr)
    call HYPRE_StructVectorDestroy(sol,ierr)
  end subroutine finalize_matrix
  subroutine add_constant_to_boundary(lo,hi,is_bound,alpha_bc,mat)
    implicit none
    integer    , intent(in ), dimension(    3) :: lo,hi
    logical    , intent(in ), dimension(0:1,3) :: is_bound
    real(rp)   , intent(in ), dimension(0:1,3) :: alpha_bc
    type(C_PTR), intent(inout) :: mat
    real(rp), dimension(product(hi(:)-lo(:)+1)) :: matvalues
    integer :: i,j,k,q
    real(rp) :: cc
    q = 0
    do k=lo(3),hi(3)
      do j=lo(2),hi(2)
        do i=lo(1),hi(1)
          q=q+1
          cc = 0._rp
          if(is_bound(0,1).and.i == lo(1)) then
            cc = cc + alpha_bc(0,1)
          end if
          if(is_bound(1,1).and.i == hi(1)) then
            cc = cc + alpha_bc(1,1)
          end if
          if(is_bound(0,2).and.j == lo(2)) then
            cc = cc + alpha_bc(0,2)
          end if
          if(is_bound(1,2).and.j == hi(2)) then
            cc = cc + alpha_bc(1,2)
          end if
          if(is_bound(0,3).and.k == lo(3)) then
            cc = cc + alpha_bc(0,3)
          end if
          if(is_bound(1,3).and.k == hi(3)) then
            cc = cc + alpha_bc(1,3)
          end if
          matvalues(q) = cc
        end do
      end do
    end do
    call HYPRE_StructMatrixAddToBoxValues(mat,lo,hi,1,[0],matvalues,ierr)
  end subroutine add_constant_to_boundary
  subroutine add_constant_to_n_boundaries(n,lo,hi,is_bound,alpha_bc,mat)
    implicit none
    integer    , intent(in   )               :: n
    integer    , intent(in   ), dimension(2) :: lo,hi
    logical    , intent(in   ), dimension(0:1,2) :: is_bound
    real(rp)   , intent(in   ), dimension(0:1,2) :: alpha_bc
    type(C_PTR), intent(inout), dimension(:) :: mat
    integer :: q
    logical , dimension(0:1,3) :: is_bound_aux
    real(rp), dimension(0:1,3) :: alpha_bc_aux
    alpha_bc_aux(:,1:2) = alpha_bc(:,:)
    alpha_bc_aux(:,  3) = 0._rp
    is_bound_aux(:,1:2) = is_bound(:,:)
    is_bound_aux(:,  3) = .false.
    do q=1,n
      call add_constant_to_boundary([lo(1),lo(2),1],[hi(1),hi(2),1],is_bound_aux,alpha_bc_aux,mat(q))
    end do
  end subroutine add_constant_to_n_boundaries
  subroutine add_constant_to_n_3d_boundaries(nslices,lo_sp,hi_sp,is_bound,alpha_bc,mat)
    implicit none
    integer, intent(in)                       :: nslices
    integer, intent(in), dimension(3,nslices) :: lo_sp,hi_sp
    logical    , intent(in   ), dimension(0:1,3) :: is_bound
    real(rp)   , intent(in   ), dimension(0:1,3) :: alpha_bc
    type(C_PTR), intent(inout), dimension(:) :: mat
    integer :: q
    do q=1,nslices
      call add_constant_to_boundary(lo_sp(:,q),hi_sp(:,q),is_bound,alpha_bc,mat(q))
    end do
  end subroutine add_constant_to_n_3d_boundaries
#if defined(_FFT_X) || defined(_FFT_Y) || defined(_FFT_Z)
  subroutine init_matrix_2d(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo,hi,periods, &
                            dl1_1,dl1_2,dl2_1,dl2_2,alpha,alpha_bc,comm,asolver,is_bound_outflow)
    !
    ! description
    !
    implicit none
    integer, parameter :: nstencil = 5
    character(len=1)  , intent(in ), dimension(0:1,2) :: cbc
    real(rp)          , intent(in ), dimension(0:1,2) ::  bc
    real(rp)          , intent(in ), dimension(0:1,2) ::  dl
    logical           , intent(in )                   ::  is_uniform_grid
    logical           , intent(in ), dimension(0:1,2) ::  is_bound
    logical           , intent(in ), dimension(    2) ::  is_centered
    integer           , intent(in ), dimension(    2) :: lo,hi,periods
    real(rp)          , intent(in ), target, dimension(lo(1)-1:) :: dl1_1,dl1_2
    real(rp)          , intent(in ), target, dimension(lo(2)-1:) :: dl2_1,dl2_2
    real(rp)          , intent(in )                              :: alpha
    real(rp)          , intent(in ), dimension(0:1,2)            :: alpha_bc
    type(MPI_COMM)    , intent(in )                   :: comm
    type(hypre_solver), intent(out)                              :: asolver
    logical           , intent(in ), optional, dimension(0:1,2)  :: is_bound_outflow
    integer, dimension(2         ) :: qqq
    integer, dimension(2,nstencil) :: offsets
    real(rp), dimension(product(hi(:)-lo(:)+1)*nstencil) :: matvalues
    real(rp), dimension(0:1,2) :: factor,sgn
    type(C_PTR) :: grid,stencil,mat,rhs,sol
    integer :: i1,i2,q,qq
    real(rp) :: cc,c1m,c1p,c2m,c2p
    integer            :: comm_hypre
    logical, dimension(0:1,2) :: is_bound_outflow_aux
    !
    comm_hypre = comm%MPI_VAL
    !
    qqq(:) = 0
    where(.not.is_centered(:)) qqq(:) = 1
    factor(:,:) = 0._rp
    sgn(   :,:) = 0._rp
    do q=1,2
      do qq=0,1
        if(is_bound(qq,q)) then
          select case(cbc(qq,q))
          case('N')
            factor(qq,q) = 1._rp*dl(qq,q)*bc(qq,q)
            if(qq == 1) factor(qq,q) = -factor(qq,q)
            sgn(   qq,q) = 1._rp
          case('D')
            if(is_centered(q)) then
              factor(qq,q) = -2._rp*bc(qq,q)
              sgn(   qq,q) = -1._rp
            else
              factor(qq,q) = -1._rp*bc(qq,q)
              sgn(   qq,q) =  0._rp
            end if
          end select
        end if
      end do
    end do
    if(present(is_bound_outflow)) then
      is_bound_outflow_aux(:,:) = is_bound_outflow(:,:)
    else
      is_bound_outflow_aux(:,:) = .false.
    end if
    !
    ! create 2D grid object
    !
    call HYPRE_StructGridCreate(comm_hypre,2,grid,ierr)
    call HYPRE_StructGridSetPeriodic(grid,periods,ierr)
    call HYPRE_StructGridSetExtents(grid,lo,hi,ierr)
    call HYPRE_StructGridAssemble(grid,ierr)
    !
    ! setup the finite-difference stencil
    !
    call HYPRE_StructStencilCreate(2,nstencil,stencil,ierr)
    offsets = reshape([ 0, 0, &
                       -1, 0, &
                        1, 0, &
                        0,-1, &
                        0, 1 ],shape(offsets))
    do q=1,nstencil
      call HYPRE_StructStencilSetElement(stencil,q-1,offsets(:,q),ierr)
    end do
    !
    ! create coefficient matrix, and solution & right-hand-side vectors
    !
    call HYPRE_StructMatrixCreate(comm_hypre,grid,stencil,mat,ierr)
    if(is_uniform_grid) call HYPRE_StructMatrixSetSymmetric(mat,1,ierr)
    call HYPRE_StructMatrixInitialize(mat,ierr)
    call HYPRE_StructVectorCreate(comm_hypre,grid,sol,ierr)
    call HYPRE_StructVectorInitialize(sol,ierr)
    call HYPRE_StructVectorCreate(comm_hypre,grid,rhs,ierr)
    call HYPRE_StructVectorInitialize(rhs,ierr)
    q = 0
    do i2=lo(2),hi(2)
      do i1=lo(1),hi(1)
        c1m = 1._rp/(dl1_1(i1-1+qqq(1))*dl1_2(i1))
        c1p = 1._rp/(dl1_1(i1  +qqq(1))*dl1_2(i1))
        c2m = 1._rp/(dl2_1(i2-1+qqq(2))*dl2_2(i2))
        c2p = 1._rp/(dl2_1(i2  +qqq(2))*dl2_2(i2))
        cc  = -(c1m+c1p+c2m+c2p) + alpha
        if(periods(1) == 0) then
          if(is_bound(0,1).and.i1 == lo(1)) then
            cc = cc + sgn(0,1)*c1m + alpha_bc(0,1)
            c1m = 0._rp
            if( is_bound_outflow_aux(0,1) ) then
              cc = cc + sgn(0,1)*c1p
              c1p = 0._rp
            end if
          end if
          if(is_bound(1,1).and.i1 == hi(1)) then
            cc = cc + sgn(1,1)*c1p + alpha_bc(1,1)
            c1p = 0._rp
            if( is_bound_outflow_aux(1,1) ) then
              cc = cc + sgn(1,1)*c1m
              c1m = 0._rp
            end if
          end if
        end if
        if(periods(2) == 0) then
          if(is_bound(0,2).and.i2 == lo(2)) then
            cc = cc + sgn(0,2)*c2m + alpha_bc(0,2)
            c2m = 0._rp
            if( is_bound_outflow_aux(0,2) ) then
              cc = cc + sgn(0,2)*c2p
              c2p = 0._rp
            end if
          end if
          if(is_bound(1,2).and.i2 == hi(2)) then
            cc = cc + sgn(1,2)*c2p + alpha_bc(1,2)
            c2p = 0._rp
            if( is_bound_outflow_aux(1,2) ) then
              cc = cc + sgn(1,2)*c2m
              c2m = 0._rp
            end if
          end if
        end if
        q  = q + 1
        qq = (q-1)*nstencil
        matvalues(qq+1) = cc
        matvalues(qq+2) = c1m
        matvalues(qq+3) = c1p
        matvalues(qq+4) = c2m
        matvalues(qq+5) = c2p
      end do
    end do
    call HYPRE_StructMatrixSetBoxValues(mat,lo,hi,nstencil, &
                                        [0,1,2,3,4],matvalues,ierr)
    call HYPRE_StructMatrixAssemble(mat,ierr)
    asolver%grid       = grid
    asolver%stencil    = stencil
    asolver%mat        = mat
    asolver%rhs        = rhs
    asolver%sol        = sol
    asolver%comm_hypre = comm_hypre
  end subroutine init_matrix_2d
  subroutine solve_n_helmholtz_2d(asolver,lo_out,hi_out,nh,lo,hi,p,po)
    implicit none
    type(hypre_solver), target, intent(inout), dimension(:) :: asolver
    integer           ,         intent(in   )               :: lo_out,hi_out,nh
    integer           ,         intent(in   ), dimension(2) :: lo,hi
#ifdef _FFT_Z
    real(rp)          ,         intent(inout), dimension(lo(1)-nh:,lo(2)-nh:,lo_out-nh:) :: p,po
#elif  _FFT_Y
    real(rp)          ,         intent(inout), dimension(lo(1)-nh:,lo_out-nh:,lo(2)-nh:) :: p,po
#elif  _FFT_X
    real(rp)          ,         intent(inout), dimension(lo_out-nh:,lo(1)-nh:,lo(2)-nh:) :: p,po
#endif
    type(C_PTR), pointer :: solver,mat,rhs,sol
    integer    , pointer :: stype
    integer :: i_out,n
    do i_out=lo_out,hi_out
      n = i_out-lo_out+1
      solver  => asolver(n)%solver
      mat     => asolver(n)%mat
      rhs     => asolver(n)%rhs
      sol     => asolver(n)%sol
      stype   => asolver(n)%stype
      !
      ! setup solution and rhs vectors
      !
#ifdef _FFT_Z
      call HYPRE_StructVectorSetBoxValues(rhs,lo,hi, p(lo(1):hi(1),lo(2):hi(2),i_out),ierr)
      call HYPRE_StructVectorSetBoxValues(sol,lo,hi,po(lo(1):hi(1),lo(2):hi(2),i_out),ierr)
#elif  _FFT_Y
      call HYPRE_StructVectorSetBoxValues(rhs,lo,hi, p(lo(1):hi(1),i_out,lo(2):hi(2)),ierr)
      call HYPRE_StructVectorSetBoxValues(sol,lo,hi,po(lo(1):hi(1),i_out,lo(2):hi(2)),ierr)
#elif  _FFT_X
      call HYPRE_StructVectorSetBoxValues(rhs,lo,hi, p(i_out,lo(1):hi(1),lo(2):hi(2)),ierr)
      call HYPRE_StructVectorSetBoxValues(sol,lo,hi,po(i_out,lo(1):hi(1),lo(2):hi(2)),ierr)
#endif
      call HYPRE_StructVectorAssemble(rhs,ierr)
      call HYPRE_StructVectorAssemble(sol,ierr)
      !
      ! setup solver, and solve
      !
      ! note: this part was based on the the Paris Simulator code
      !       freely available under a GPL license; see:
      !       http://www.ida.upmc.fr/~zaleski/paris/
      !
      if      ( stype == HYPRESolverSMG ) then
        call HYPRE_StructSMGSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructSMGGetNumIterations(solver,num_iterations,ierr)
      else if ( stype == HYPRESolverPFMG ) then
        call HYPRE_StructPFMGSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructPFMGGetNumIteration(solver,num_iterations,ierr)
      else if ( stype == HYPRESolverGMRES ) then
        call HYPRE_StructGMRESSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructGMRESGetNumIteratio(solver, num_iterations,ierr)
      else if ( stype == HYPRESolverBiCGSTAB ) then
        call HYPRE_StructBiCGSTABSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructBiCGSTABGetNumItera(solver, num_iterations,ierr)
      end if ! stype
      !
      ! end of part based on the Paris Simulator code
      !
      ! fecth results
      !
#ifdef _FFT_Z
      call HYPRE_StructVectorGetBoxValues(sol,lo,hi,p(lo(1):hi(1),lo(2):hi(2),i_out),ierr)
      !$OMP WORKSHARE
      po(lo(1):hi(1),lo(2):hi(2),i_out) = p(lo(1):hi(1),lo(2):hi(2),i_out)
      !$OMP END WORKSHARE
#elif  _FFT_Y
      call HYPRE_StructVectorGetBoxValues(sol,lo,hi,p(lo(1):hi(1),i_out,lo(2):hi(2)),ierr)
      !$OMP WORKSHARE
      po(lo(1):hi(1),i_out,lo(2):hi(2)) = p(lo(1):hi(1),i_out,lo(2):hi(2))
      !$OMP END WORKSHARE
#elif  _FFT_X
      call HYPRE_StructVectorGetBoxValues(sol,lo,hi,p(i_out,lo(1):hi(1),lo(2):hi(2)),ierr)
      !$OMP WORKSHARE
      po(i_out,lo(1):hi(1),lo(2):hi(2)) = p(i_out,lo(1):hi(1),lo(2):hi(2))
      !$OMP END WORKSHARE
#endif
    end do
  end subroutine solve_n_helmholtz_2d
  subroutine solve_n_helmholtz_3d(asolver,nslices,lo_sp,hi_sp,nh,lo,hi,p,po)
    implicit none
    type(hypre_solver), target, intent(inout), dimension(:) :: asolver
    integer           ,         intent(in   )               :: nslices
    integer           ,         intent(in   ), dimension(3,nslices) :: lo_sp,hi_sp
    integer           ,         intent(in   )               :: nh
    integer           ,         intent(in   ), dimension(3) :: lo,hi
    real(rp)          ,         intent(inout), dimension(lo(1)-nh:,lo(2)-nh:,lo(3)-nh:) :: p,po
    type(C_PTR), pointer :: solver,mat,rhs,sol
    integer    , pointer :: stype
    integer, dimension(3) :: lo_s,hi_s
    integer :: q
    do q=1,nslices
      solver  => asolver(q)%solver
      mat     => asolver(q)%mat
      rhs     => asolver(q)%rhs
      sol     => asolver(q)%sol
      stype   => asolver(q)%stype
      lo_s(:) = lo_sp(:,q)
      hi_s(:) = hi_sp(:,q)
      !
      ! setup soluction and rhs vectors
      !
      call HYPRE_StructVectorSetBoxValues(rhs,lo_s,hi_s, p(lo_s(1):hi_s(1),lo_s(2):hi_s(2),lo_s(3):hi_s(3)),ierr)
      call HYPRE_StructVectorSetBoxValues(sol,lo_s,hi_s,po(lo_s(1):hi_s(1),lo_s(2):hi_s(2),lo_s(3):hi_s(3)),ierr)
      call HYPRE_StructVectorAssemble(rhs,ierr)
      call HYPRE_StructVectorAssemble(sol,ierr)
      !
      ! setup solver, and solve
      !
      ! note: this part was based on the the Paris Simulator code
      !       freely available under a GPL license; see:
      !       http://www.ida.upmc.fr/~zaleski/paris/
      !
      if      ( stype == HYPRESolverSMG ) then
        call HYPRE_StructSMGSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructSMGGetNumIterations(solver,num_iterations,ierr)
      else if ( stype == HYPRESolverPFMG ) then
        call HYPRE_StructPFMGSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructPFMGGetNumIteration(solver,num_iterations,ierr)
      else if ( stype == HYPRESolverGMRES ) then
        call HYPRE_StructGMRESSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructGMRESGetNumIteratio(solver, num_iterations,ierr)
      else if ( stype == HYPRESolverBiCGSTAB ) then
        call HYPRE_StructBiCGSTABSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructBiCGSTABGetNumItera(solver, num_iterations,ierr)
      end if ! stype
      !
      ! end of part based on the Paris Simulator code
      !
      ! fecth results
      !
      call HYPRE_StructVectorGetBoxValues(sol,lo_s,hi_s,p(lo_s(1):hi_s(1),lo_s(2):hi_s(2),lo_s(3):hi_s(3)),ierr)
      !$OMP WORKSHARE
      po(lo_s(1):hi_s(1),lo_s(2):hi_s(2),lo_s(3):hi_s(3)) = p(lo_s(1):hi_s(1),lo_s(2):hi_S(2),lo_s(3):hi_s(3))
      !$OMP END WORKSHARE
    end do
  end subroutine solve_n_helmholtz_3d
  subroutine init_n_2d_matrices(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo_out,hi_out,lo,hi,periods, &
                                dl1_1,dl1_2,dl2_1,dl2_2,alpha,alpha_bc,lambda,comm,asolver,is_bound_outflow)
    character(len=1)  , intent(in   ), dimension(0:1,2) :: cbc
    real(rp)          , intent(in   ), dimension(0:1,2) ::  bc
    real(rp)          , intent(in   ), dimension(0:1,2) ::  dl
    logical           , intent(in   )                   ::  is_uniform_grid
    logical           , intent(in   ), dimension(0:1,2) ::  is_bound
    logical           , intent(in   ), dimension(    2) ::  is_centered
    integer           , intent(in   )                   :: lo_out,hi_out
    integer           , intent(in   ), dimension(    2) :: lo,hi,periods
    real(rp)          , intent(in   ), target, dimension(lo(1)-1:) :: dl1_1,dl1_2
    real(rp)          , intent(in   ), target, dimension(lo(2)-1:) :: dl2_1,dl2_2
    real(rp)          , intent(in   )                              :: alpha
    real(rp)          , intent(in   ), dimension(0:1,2)            :: alpha_bc
    real(rp)          , intent(in   ), dimension(:) :: lambda
    type(MPI_COMM)    , intent(in   ), dimension(:) :: comm
    type(hypre_solver), intent(inout), dimension(:) :: asolver
    logical           , intent(in   ), optional, dimension(0:1,2)  :: is_bound_outflow
    type(hypre_solver) :: asolver_aux
    integer :: i_out,q
    logical, dimension(0:1,2) :: is_bound_outflow_aux
    is_bound_outflow_aux(:,:) = .false.
    if( present(is_bound_outflow) ) is_bound_outflow_aux(:,:) = is_bound_outflow(:,:)
    do i_out=lo_out,hi_out
      q = i_out-lo_out+1
      asolver_aux = asolver(q)
      call init_matrix_2d(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo,hi,periods, &
                          dl1_1,dl1_2,dl2_1,dl2_2,alpha,alpha_bc,comm(q),asolver_aux,is_bound_outflow_aux)
      call add_constant_to_diagonal([lo(1),lo(2),1],[hi(1),hi(2),1],lambda(q),asolver_aux%mat)
      asolver(q) = asolver_aux
    end do
  end subroutine init_n_2d_matrices
  subroutine init_n_3d_matrices(idir,nslice,cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo,periods, &
                                lo_sp,hi_sp,dl1_1,dl1_2,dl2_1,dl2_2,dl3_1,dl3_2,alpha,alpha_bc,lambda,asolver, &
                                is_bound_outflow)
    integer           , intent(in   )                   :: idir,nslice
    character(len=1)  , intent(in   ), dimension(0:1,3) :: cbc
    real(rp)          , intent(in   ), dimension(0:1,3) ::  bc
    real(rp)          , intent(in   ), dimension(0:1,3) ::  dl
    logical           , intent(in   )                   ::  is_uniform_grid
    logical           , intent(in   ), dimension(0:1,3) ::  is_bound
    logical           , intent(in   ), dimension(    3) ::  is_centered
    integer           , intent(in   ), dimension(    3) :: lo,periods
    integer           , intent(in   ), dimension(:,:  ) :: lo_sp,hi_sp
    real(rp)          , intent(in   ), target, dimension(lo(1)-1:) :: dl1_1,dl1_2
    real(rp)          , intent(in   ), target, dimension(lo(2)-1:) :: dl2_1,dl2_2
    real(rp)          , intent(in   ), target, dimension(lo(3)-1:) :: dl3_1,dl3_2
    real(rp)          , intent(in   )                              :: alpha
    real(rp)          , intent(in   ), dimension(0:1,3)            :: alpha_bc
    real(rp)          , intent(in   ), dimension(:) :: lambda
    type(hypre_solver), intent(inout), dimension(:) :: asolver
    logical           , intent(in   ), optional, dimension(0:1,3)  :: is_bound_outflow
    type(hypre_solver) :: asolver_aux
    integer :: q
    logical, dimension(0:1,3) :: is_bound_outflow_aux
    is_bound_outflow_aux(:,:) = .false.
    if( present(is_bound_outflow) ) is_bound_outflow_aux(:,:) = is_bound_outflow(:,:)
    do q=1,nslice
      asolver_aux = asolver(q)
      call init_matrix_3d(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo_sp(:,q),hi_sp(:,q),periods, &
                          dl1_1(lo_sp(1,q)-1:hi_sp(1,q)+1),dl1_2(lo_sp(1,q)-1:hi_sp(1,q)+1), &
                          dl2_1(lo_sp(2,q)-1:hi_sp(2,q)+1),dl2_2(lo_sp(2,q)-1:hi_sp(2,q)+1), &
                          dl3_1(lo_sp(3,q)-1:hi_sp(3,q)+1),dl3_2(lo_sp(3,q)-1:hi_sp(3,q)+1), &
                          alpha,alpha_bc, &
                          asolver_aux,lambda(lo_sp(idir,q)-lo(idir)+1:hi_sp(idir,q)-lo(idir)+1), &
                          is_bound_outflow_aux)
      asolver(q) = asolver_aux
    end do
  end subroutine init_n_3d_matrices
  subroutine create_n_solvers(n,maxiter,maxerror,stype,asolver)
    integer           , intent(   in) :: n
    integer           , intent(   in) :: maxiter
    real(rp)          , intent(   in) :: maxerror
    integer           , intent(   in) :: stype
    type(hypre_solver), intent(inout), dimension(:) :: asolver
    integer :: q
    do q=1,n
      call create_solver(maxiter,maxerror,stype,asolver(q))
    end do
  end subroutine create_n_solvers
  subroutine setup_n_solvers(n,asolver)
    integer           , intent(   in) :: n
    type(hypre_solver), intent(inout), dimension(:) :: asolver
    integer :: q
    do q=1,n
      call setup_solver(asolver(q))
    end do
  end subroutine setup_n_solvers
  subroutine add_constant_to_n_diagonals(n,lo,hi,alpha,mat)
    integer    , intent(in   )               :: n
    integer    , intent(in   ), dimension(2) :: lo,hi
    real(rp)   , intent(in   )               :: alpha
    type(C_PTR), intent(inout), dimension(:) :: mat
    integer :: q
    do q=1,n
      call add_constant_to_diagonal([lo(1),lo(2),1],[hi(1),hi(2),1],alpha,mat(q))
    end do
  end subroutine add_constant_to_n_diagonals
  subroutine finalize_n_matrices(n,asolver)
    integer           , intent(   in) :: n
    type(hypre_solver), intent(inout), dimension(:) :: asolver
    integer :: q
    do q=1,n
      call finalize_matrix(asolver(q))
    end do
  end subroutine finalize_n_matrices
  subroutine finalize_n_solvers(n,asolver)
    integer           , intent(   in) :: n
    type(hypre_solver), intent(inout), dimension(:) :: asolver
    integer :: q
    do q=1,n
      call finalize_solver(asolver(q))
    end do
  end subroutine finalize_n_solvers
  subroutine init_fft_reduction(idir,n,bc,is_centered,dl,arrplan,normfft,lambda)
    use mod_fft       , only: fftini,eigenvalues
    implicit none
    integer         , intent(in )                 :: idir
    integer         , intent(in ), dimension(3  ) :: n
    character(len=1), intent(in ), dimension(0:1) :: bc
    logical         , intent(in )                 :: is_centered
    real(rp)        , intent(in )                 :: dl
    type(C_PTR)     , intent(out), dimension(2  ) :: arrplan
    real(rp)        , intent(out)                 :: normfft
    real(rp)        , intent(out), dimension(:  ) :: lambda
    call fftini(idir,n,bc,is_centered,arrplan,normfft)
    call eigenvalues(n(idir),bc,is_centered,lambda)
    lambda(:) = lambda(:)/dl**2
  end subroutine init_fft_reduction
  subroutine solve_n_helmholtz_2d_old(asolver,maxiter,maxerror,lo_out,hi_out,lo,hi,lambda,alpha_old,p,po)
    implicit none
    type(hypre_solver), target, intent(inout)               :: asolver
    integer           ,         intent(in   )               :: maxiter
    real(rp)          ,         intent(in   )               :: maxerror
    integer           ,         intent(in   )               :: lo_out,hi_out
    integer           ,         intent(in   ), dimension(2) :: lo,hi
    real(rp)          ,         intent(inout)               :: alpha_old
    real(rp)          ,         intent(in   ), dimension(:) :: lambda
#ifdef _FFT_Z
    real(rp)          ,         intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo_out-1:) :: p,po
#elif  _FFT_Y
    real(rp)          ,         intent(inout), dimension(lo(1)-1:,lo_out-1:,lo(2)-1:) :: p,po
#elif  _FFT_X
    real(rp)          ,         intent(inout), dimension(lo_out-1:,lo(1)-1:,lo(2)-1:) :: p,po
#endif
    type(C_PTR), pointer :: solver,mat,rhs,sol
    integer    , pointer :: stype
    real(rp), dimension(product(hi(:)-lo(:)+1)) :: solvalues,rhsvalues
    integer :: i_out,ii,i1,i2,q
    real(rp) :: alpha
    solver  => asolver%solver
    mat     => asolver%mat
    rhs     => asolver%rhs
    sol     => asolver%sol
    stype   => asolver%stype
    do i_out=lo_out,hi_out
      q = 0
      do i2=lo(2),hi(2)
        do i1=lo(1),hi(1)
          q = q + 1
#ifdef _FFT_Z
          rhsvalues(q) = p( i1,i2,i_out)
          solvalues(q) = po(i1,i2,i_out)
#elif  _FFT_Y
          rhsvalues(q) = p( i1,i_out,i2)
          solvalues(q) = po(i1,i_out,i2)
#elif  _FFT_X
          rhsvalues(q) = p( i_out,i1,i2)
          solvalues(q) = po(i_out,i1,i2)
#endif
        end do
      end do
      ii = i_out-lo_out+1
      alpha = lambda(ii) - alpha_old
      call add_constant_to_diagonal([lo(1),lo(2),1],[hi(1),hi(2),1],alpha,mat)
      alpha_old = lambda(ii)
      !
      call HYPRE_StructVectorSetBoxValues(rhs,lo,hi, &
                                          rhsvalues,ierr)
      call HYPRE_StructVectorAssemble(rhs,ierr)
      !
      ! create soluction vector
      !
      call HYPRE_StructVectorSetBoxValues(sol,lo,hi, &
                                          solvalues,ierr)
      call HYPRE_StructVectorAssemble(sol,ierr)
      !
      ! setup solver, and solve
      !
      ! note: this part was based on the the Paris Simulator code
      !       freely available under a GPL license; see:
      !       http://www.ida.upmc.fr/~zaleski/paris/
      !
      call create_solver(maxiter,maxerror,stype,asolver)
      call setup_solver(asolver)
      if      ( stype == HYPRESolverSMG ) then
        call HYPRE_StructSMGSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructSMGGetNumIterations(solver,num_iterations,ierr)
      else if ( stype == HYPRESolverPFMG ) then
        call HYPRE_StructPFMGSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructPFMGGetNumIteration(solver,num_iterations,ierr)
      else if ( stype == HYPRESolverGMRES ) then
        call HYPRE_StructGMRESSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructGMRESGetNumIteratio(solver, num_iterations,ierr)
      else if ( stype == HYPRESolverBiCGSTAB ) then
        call HYPRE_StructBiCGSTABSolve(solver,mat,rhs,sol,ierr)
        !call HYPRE_StructBiCGSTABGetNumItera(solver, num_iterations,ierr)
      end if ! stype
      !
      ! end of part based on the Paris Simulator code
      !
      ! fecth results
      !
      call HYPRE_StructVectorGetBoxValues(sol,lo,hi,solvalues,ierr)
      q = 0
      do i2=lo(2),hi(2)
        do i1=lo(1),hi(1)
          q = q + 1
#ifdef _FFT_Z
          p( i1,i2,i_out) = solvalues(q)
          po(i1,i2,i_out) = p(i1,i2,i_out)
#elif  _FFT_Y
          p( i1,i_out,i2) = solvalues(q)
          po(i1,i_out,i2) = p(i1,i_out,i2)
#elif  _FFT_X
          p( i_out,i1,i2) = solvalues(q)
          po(i_out,i1,i2) = p(i_out,i1,i2)
#endif
        end do
      end do
      call finalize_solver(asolver)
    end do
  end subroutine solve_n_helmholtz_2d_old
  subroutine init_transpose_slab_uneven(idir,nhi,nho,dims,lo_p,lo_s,n_p,n_s,comm_block,transpose_params)
    implicit none
    integer, intent(in)                            :: idir,nhi,nho
    integer, intent(in), dimension(3)              :: dims,lo_p,lo_s,n_p,n_s ! n.b. lower bounds wrt [0,0,0]
    type(MPI_COMM)     , intent(in)                :: comm_block
    type(alltoallw), intent(out), dimension(:,:) :: transpose_params
    integer, dimension(3) :: n_i
    type(MPI_DATATYPE) :: type_sub_p,type_sub_s
    integer(MPI_ADDRESS_KIND) :: lb,ext_rp
    integer, dimension(3,3) :: eye
    integer :: i,j,k,ii,jj,kk,irank
    integer, dimension(3,product(dims)) :: lo_p_a,lo_s_a,n_p_a,n_s_a,n_i_a
    !
    eye(:,:) = 0
    do i=1,3
      eye(i,i) = 1
    end do
    call MPI_ALLGATHER(lo_p,3,MPI_INTEGER,lo_p_a,3,MPI_INTEGER,comm_block)
    call MPI_ALLGATHER(lo_s,3,MPI_INTEGER,lo_s_a,3,MPI_INTEGER,comm_block)
    call MPI_ALLGATHER(n_p ,3,MPI_INTEGER,n_p_a ,3,MPI_INTEGER,comm_block)
    call MPI_ALLGATHER(n_s ,3,MPI_INTEGER,n_s_a ,3,MPI_INTEGER,comm_block)
    transpose_params(:,:)%counts = 1
    call MPI_TYPE_GET_EXTENT(MPI_REAL_RP,lb,ext_rp)
    irank = 0
    do k=0,dims(3)-1
      do j=0,dims(2)-1
        do i=0,dims(1)-1
          ii = lo_s_a(1,irank+1)*eye(1,idir) + 1
          jj = lo_s_a(2,irank+1)*eye(2,idir) + 1
          kk = lo_s_a(3,irank+1)*eye(3,idir) + 1
          transpose_params(irank+1,1)%disps = extent_f(ii,jj,kk,n_p(:)+2*nhi)*ext_rp
          ii = lo_p_a(1,irank+1)*(1-eye(1,idir)) + 1
          jj = lo_p_a(2,irank+1)*(1-eye(2,idir)) + 1
          kk = lo_p_a(3,irank+1)*(1-eye(3,idir)) + 1
          transpose_params(irank+1,2)%disps = extent_f(ii,jj,kk,n_s(:)+2*nho)*ext_rp
          !
          n_i(:) = n_s_a(:,irank+1)*(eye(:,idir)) + n_p(:)*(1-eye(:,idir))
          call MPI_TYPE_CREATE_SUBARRAY(3,n_p(:)+2*nhi,n_i(:),[0,0,0],MPI_ORDER_FORTRAN,MPI_REAL_RP,type_sub_p)
          n_i(:) = n_s(:)*(eye(:,idir)) + n_p_a(:,irank+1)*(1-eye(:,idir))
          call MPI_TYPE_CREATE_SUBARRAY(3,n_s(:)+2*nho,n_i(:),[0,0,0],MPI_ORDER_FORTRAN,MPI_REAL_RP,type_sub_s)
          call MPI_TYPE_COMMIT(type_sub_p)
          call MPI_TYPE_COMMIT(type_sub_s)
          transpose_params(irank+1,1)%types = type_sub_p
          transpose_params(irank+1,2)%types = type_sub_s
          irank = irank + 1
        end do
      end do
    end do
  end subroutine init_transpose_slab_uneven
  subroutine init_transpose_slab(idir,nhi,nho,dims,n_p,n_s,transpose_params)
    implicit none
    integer, intent(in)                            :: idir,nhi,nho
    integer, intent(in), dimension(3)              :: dims,n_p,n_s
    type(alltoallw), intent(out), dimension(:,:)   :: transpose_params
    integer, dimension(3) :: n_i
    type(MPI_DATATYPE) :: type_sub_p,type_sub_s
    integer(MPI_ADDRESS_KIND) :: lb,ext_rp
    integer, dimension(3,3) :: eye
    integer :: i,j,k,ii,jj,kk,irank
    !
    eye(:,:) = 0
    do i=1,3
      eye(i,i) = 1
    end do
    select case(idir)
      case(1)
        n_i(:) = [n_s(1),n_p(2),n_p(3)]
      case(2)
        n_i(:) = [n_p(1),n_s(2),n_p(3)]
      case(3)
        n_i(:) = [n_p(1),n_p(2),n_s(3)]
    end select
    transpose_params(:,:)%counts = 1
    call MPI_TYPE_GET_EXTENT(MPI_REAL_RP,lb,ext_rp)
    irank = 0
    do k=0,dims(3)-1
      do j=0,dims(2)-1
        do i=0,dims(1)-1
          ii = irank*n_s(1)*eye(1,idir) + 1
          jj = irank*n_s(2)*eye(2,idir) + 1
          kk = irank*n_s(3)*eye(3,idir) + 1
          transpose_params(irank+1,1)%disps = extent_f(ii,jj,kk,n_p+2*nhi)*ext_rp
          ii = i*n_p(1)*(1-eye(1,idir)) + 1
          jj = j*n_p(2)*(1-eye(2,idir)) + 1
          kk = k*n_p(3)*(1-eye(3,idir)) + 1
          transpose_params(irank+1,2)%disps = extent_f(ii,jj,kk,n_s+2*nho)*ext_rp
          irank = irank + 1
        end do
      end do
    end do
    call MPI_TYPE_CREATE_SUBARRAY(3,n_p+2*nhi,n_i,[0,0,0],MPI_ORDER_FORTRAN,MPI_REAL_RP,type_sub_p)
    call MPI_TYPE_CREATE_SUBARRAY(3,n_s+2*nho,n_i,[0,0,0],MPI_ORDER_FORTRAN,MPI_REAL_RP,type_sub_s)
    call MPI_TYPE_COMMIT(type_sub_p)
    call MPI_TYPE_COMMIT(type_sub_s)
    transpose_params(:,1)%types = type_sub_p
    transpose_params(:,2)%types = type_sub_s
  end subroutine init_transpose_slab
  pure integer function extent_f(i,j,k,n) ! memory position of array with point i,j,k assuming Fortran ordering
     implicit none
     integer, intent(in) :: i,j,k,n(3)
     extent_f = (i-1) + (j-1)*n(1) + (k-1)*n(1)*n(2)
  end function extent_f
  subroutine transpose_slab(nhi,nho,t_param,comm_block,arr_in,arr_out)
    implicit none
    integer, intent(in)                            :: nhi,nho
    type(alltoallw), dimension(:,:), intent(in)    :: t_param
    type(MPI_COMM), intent(in)                     :: comm_block
    real(rp), intent(in ), dimension(1-nhi:,1-nhi:,1-nhi:) :: arr_in
    real(rp), intent(out), dimension(1-nho:,1-nho:,1-nho:) :: arr_out
    call MPI_ALLTOALLW(arr_in( 1,1,1),t_param(:,1)%counts,t_param(:,1)%disps,t_param(:,1)%types, &
                       arr_out(1,1,1),t_param(:,2)%counts,t_param(:,2)%disps,t_param(:,2)%types, &
                       comm_block)
  end subroutine transpose_slab
  subroutine init_comm_slab(lo_idir,hi_idir,lo_s_idir,hi_s_idir,myid,comms_slab)
    implicit none
    integer, intent(in) :: lo_idir,hi_idir,lo_s_idir,hi_s_idir,myid
    type(MPI_COMM), dimension(hi_s_idir-lo_s_idir+1), intent(out) :: comms_slab
    type(MPI_COMM) :: comm
    integer :: i,n,icolor
    do i=lo_idir,hi_idir
      n = i-lo_s_idir+1
      icolor = MPI_UNDEFINED
      if( i >= lo_s_idir .and. i <= hi_s_idir) icolor=i
      ! n.b. -- order of ranks could be set by e.g. memory layout (i-1 + (j-1)*n(1) + (k-1)*n(1)*n(2)*k)
      call MPI_COMM_SPLIT(MPI_COMM_WORLD,icolor,myid,comm)
      if(i >= lo_s_idir .and. i <= hi_s_idir) comms_slab(n) = comm
    end do
  end subroutine init_comm_slab
#endif
end module mod_solver

max_stars_repo_name: p-costa/SNaC

max_issues_repo_name: p-costa/SNaC

max_forks_repo_name: p-costa/SNaC


text:
PROGRAM AveGeom

  USE Util, ONLY : Startup, Cleanup
  IMPLICIT NONE

  ! read all the structures and average each FEP step
  CALL Startup(readCoords=.TRUE.)
  CALL MeanFepStructures()
  CALL Cleanup()

  CONTAINS

!*

  SUBROUTINE MeanFepStructures()

    USE Input, ONLY : trajectory
    USE Output, ONLY : WriteXYZ
    USE FileIO, ONLY : OpenFile, CloseFile
    IMPLICIT NONE
    INTEGER :: fepstep, timestep
    CHARACTER(50) :: title, fileName
    CHARACTER(5) :: elements(SIZE(trajectory,3))
    REAL(8) :: meanStructure(SIZE(trajectory,1),3,SIZE(trajectory,3))
    INTEGER, PARAMETER :: xyzUnit = 42

    elements(:) = "H"

    DO fepstep = 1, SIZE(trajectory,1)

      ! sum the atomic positions over all timesteps of this simulation
      DO timestep = 1, SIZE(trajectory,4)
        meanStructure(fepstep,:,:) = meanStructure(fepstep,:,:) + trajectory(fepstep,:,:,timestep)
      ENDDO

      meanStructure(fepstep,:,:) = meanStructure(fepstep,:,:) / SIZE(trajectory,4)

      WRITE(title,'(A,I4)') "MEAN STRUCTURE AT FEP STEP ",fepstep
      WRITE(fileName,'(A,I0.4,A)') "MeanStructure_", fepstep, ".xyz"

      CALL OpenFile(xyzUnit,TRIM(ADJUSTL(fileName)),"write")
      CALL WriteXYZ(title,elements,meanStructure(fepstep,:,:),xyzUnit)
      CALL CloseFile(xyzUnit)

    ENDDO

  END SUBROUTINE MeanFepStructures

END PROGRAM AveGeom

max_stars_repo_name: MJLMills/fepcat

max_issues_repo_name: MJLMills/fepcat

max_forks_repo_name: MJLMills/fepcat


!***********************************/
!	Name:内挿元となるドナー要素の検索・登録プログラム
!	Alias:Register4InterpolationMain
!	Description:ドナー要素の共有面積も計算する
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2017.12.19
!	Update:2017.12.19
!	Other:
!***********************************/
subroutine Register4InterpolationMain(LG,iGrid,Geom)
    use StructVar_Mod
!    use LoopVar_Mod
    use RegistVar_Mod4OMP
!$ use omp_lib
    implicit none
    type(LocalGrids), intent(inout) :: LG
    type(Geometry), intent(inout) :: Geom
    integer, intent(in) :: iGrid

    type(CountUpVertex) :: CUV
    integer :: iDonorGrid, iTimeCombination
    integer :: iTargetCell
    integer :: iKindOfDonor, iFindDonor
    integer :: iGlobalCellNumber
    integer :: iLoop
    type(PrivateVarOfRegist) :: PVR
    allocate(PVR%RG2L%CenterDistanceVector(3),PVR%RG2L%iCenter(3),PVR%RG2L%iCenter2(3))
    allocate(PVR%RIG%CellCenterCoords(3),PVR%RIG%GlobalGridNum(3))
    allocate(PVR%RSDILG%MWDIAG%AuxiliaryGridNum(3),PVR%RSDILG%SSDIG%CenterDistanceVector(3))
!内挿用のDonor格子を見つけて共有面積を保存する

!$omp parallel num_threads(CoreNumberOfCPU),shared(LG,iGrid,Geom),firstprivate(iTargetCell,iKindOfDonor,iFindDonor,CUV,iGlobalCellNumber,PVR,iLoop)
    do iTargetCell=1, LG%UG(iGrid)%GI%RealCells !Li格子の全実要素について

        call InputCUV(1,0,iTargetCell,LG%UG(iGrid),LG%MG(iGrid)%NSR,CUV,LG%MG(iGrid))
        !グローバル格子における位置判定+登録を行う
        call RegisterIntoGlobal(iGrid,LG%GG,iGlobalCellNumber,Geom,CUV,iLoop,PVR%RIG) !CutCell入れるなら消す？

        iFindDonor = 0
        do iKindOfDonor=1, 3
            if(iKindOfDonor <= 2 .or. iFindDonor == 0) then !上位格子+物体あり格子からの補間が終わった時点でまだ補間要素が見つかっていない時下位格子の検索を行う

                do iDonorGrid=1, LG%TotalGridsNumber !計算領域中の全格子について

                    if(LG%UG(iGrid)%Overlap(iTimeCombination+1,iDonorGrid) == iKindOfDonor) then !iKind=1で上位格子，iKind=2で物体ありの下位格子，iKind=3で下位格子
                        call RoughSearchDonorInLocalGrid(iTimeCombination,LG%MG(iDonorGrid),LG%UG(iDonorGrid),LG%MG(iGrid),CUV,iFindDonor)
                        call InputCUV(1,0,iTargetCell,LG%UG(iGrid),LG%MG(iGrid)%NSR,CUV,LG%MG(iGrid)) !Donor格子が変わるたびにその座標系に合わせる必要があるため，毎回グローバル座標系の値に戻す
                    end if

                end do

            else !上位格子+物体あり下位格子の検索が終わった時点で補間用要素が見つかっているとき = 積分不要
                LG%UG(iGrid)%GM%Interpolated(iTargetCell,1,1) = 1
            end if
        end do

    !この時点で要素と被る可能性のあるすべてのローカル格子との接触判定が終わっている
    !そのため，この時点でiFindDonorが0であるならばグローバル格子からの補間が必要となる
        if(LG%UG(iGrid)%GM%CellType(iTargetCell,1,1) == 99) then
            call RegisterGlobal2Local(iGlobalCellNumber,LG%MG(iGrid),iFindDonor,1,CUV,Geom,9999)
            iFindDonor = 1
        end if

        if(iFindDonor == 0) then
            call RegisterGlobal2Local(iGlobalCellNumber,LG%MG(iGrid),iFindDonor,0,CUV,Geom,9999)
        end if

    end do
return
end subroutine Register4InterpolationMain


text:
!======================================================================!
  SUBROUTINE CalcMn(n0, n1)
!----------------------------------------------------------------------!
!   Calculates time averaged velocity and velocity fluctuations.       !
!----------------------------------------------------------------------!
!------------------------------[Modules]-------------------------------!
  USE all_mod
  USE pro_mod
  USE les_mod
  USE rans_mod
!----------------------------------------------------------------------!
  IMPLICIT NONE
!-----------------------------[Parameters]-----------------------------!
  INTEGER :: n0, n1
!-------------------------------[Locals]-------------------------------!
  INTEGER :: c, n
!--------------------------------[CVS]---------------------------------!
!  $Id: CalcMn.f90,v 1.17 2008/12/05 13:22:06 IUS\mhadziabdic Exp $
!  $Source: /home/IUS/mhadziabdic/.CVSROOT/T-Rex/Process/CalcMn.f90,v $
!======================================================================!

  n=n1-n0

  if(n  > -1) then
    do c=-NbC,NC
!-----------------------!
!      mean values      !
!-----------------------!
      U % mean(c) = ( U % mean(c) * (1.*n) + U % n(c) ) / (1.*(n+1))
      V % mean(c) = ( V % mean(c) * (1.*n) + V % n(c) ) / (1.*(n+1))
      W % mean(c) = ( W % mean(c) * (1.*n) + W % n(c) ) / (1.*(n+1))
      P % mean(c) = ( P % mean(c) * (1.*n) + P % n(c) ) / (1.*(n+1))

!------------------------------!
!      fluctuating values      !
!------------------------------!
      uu % mean(c) = ( uu % mean(c)*(1.*n) + U % n(c) * U % n(c) ) &
                   / (1.*(n+1))
      vv % mean(c) = ( vv % mean(c)*(1.*n) + V % n(c) * V % n(c) ) &
                   / (1.*(n+1))
      ww % mean(c) = ( ww % mean(c)*(1.*n) + W % n(c) * W % n(c) ) &
                   / (1.*(n+1))

      uv % mean(c) = ( uv % mean(c)*(1.*n) + U % n(c) * V % n(c) ) &
                   / (1.*(n+1))
      uw % mean(c) = ( uw % mean(c)*(1.*n) + U % n(c) * W % n(c) ) &
                   / (1.*(n+1))
      vw % mean(c) = ( vw % mean(c)*(1.*n) + V % n(c) * W % n(c) ) &
                   / (1.*(n+1))

      VISt_mean(c) = ( VISt_mean(c)*(1.*n) + VISt(c) ) &
                   / (1.*(n+1))

      if(HOT==YES) then
        T % mean(c) = ( T % mean(c) * (1.*n) + T % n(c) ) / (1.*(n+1))
        TT % mean(c) = ( TT % mean(c)*(1.*n) + T % n(c) * T % n(c) ) &
                     / (1.*(n+1))
        uT % mean(c) = ( uT % mean(c)*(1.*n) + u % n(c) * T % n(c) ) &
                     / (1.*(n+1))
        vT % mean(c) = ( vT % mean(c)*(1.*n) + v % n(c) * T % n(c) ) &
                     / (1.*(n+1))
        wT % mean(c) = ( wT % mean(c)*(1.*n) + w % n(c) * T % n(c) ) &
                     / (1.*(n+1))
      end if
    end do
  end if

  RETURN

  END SUBROUTINE CalcMn

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
module retrieval
  use global_variables
  implicit none

contains

  function Resfun(data, theta) result(r)
    ! Residual function
    use forward_model, only : Beer
    implicit none
    type(obs), intent(in) :: data
    real(DP), dimension(:), intent(in) :: theta
    real(DP), dimension(data%nobs+size(theta)) :: r
    real(DP), dimension(data%nprior) :: prof
    prof = Redu2full(data%P, data%prior, data%air, theta)
    r(1:data%nobs) = data%t - Beer(data%cs, prof)
    r = [r(1:data%nobs)/data%noise, theta] ! add parameters to the end
  end function Resfun

  function Jacfun(data, theta) result(j)
    ! Jacobian function
    use math, only : Diag
    use forward_model, only : Beer_jac
    implicit none
    type(obs), intent(in) :: data
    real(DP), dimension(:), intent(in) :: theta
    real(DP), dimension(data%nobs+size(theta), size(theta)) :: j
    real(DP), dimension(data%nprior) :: prof
    j = 0.0
    prof = Redu2full(data%p, data%prior, data%air, theta)
    j(1:data%nobs,:) = 0.0001 * matmul(matmul(Beer_jac(data%cs, prof), Diag(data%air)), data%p)
    j(1:data%nobs,:) = j(1:data%nobs,:) / data%noise
  end function Jacfun

  subroutine Reduce_dim(c, d, p, q)
    ! Truncation of prior covariance c using
    ! d principal components so that p'p~=c
    use math, only : Svd, Diag
    implicit none
    real(DP), dimension(:,:), intent(in) :: c
    integer, intent(in) :: d
    real(DP), dimension(size(c,1), d), intent(out) :: p
    real(DP), dimension(size(c,1)), intent(out) :: q
    real(DP), dimension(size(c,2), size(c,2)) :: s, v, u
    integer :: nite
    call Svd(c, u, s, v, nite)
    q = Diag(s)
    p = matmul(u(:,1:d), Diag(sqrt(q(1:d))))
  end subroutine Reduce_dim

  function Redu2full(P, prior, air, theta) result(profile)
    ! New profile candidate from projection marix "P"
    ! and the parameters "theta"
    implicit none
    real(DP), dimension(:,:), intent(in) :: p
    real(DP), dimension(:), intent(in) :: theta, air, prior
    real(DP), dimension(size(prior)) :: profile
    profile = prior + matmul(p, theta)*air/1e9
  end function Redu2full

end module retrieval

max_stars_repo_name: tukiains/dimred-fortran

max_issues_repo_name: tukiains/dimred-fortran

max_forks_repo_name: tukiains/dimred-fortran


text:
Sue Greenwald was a member of the City Council.  She lost her reelection bid in the 2012 City Council Election.

 Background

Sue Greenwald served as Mayor Pro Tem from 2004 to 2005, and as Mayor from 2006 to 2008FootNote(The City Council votes to pick the mayor. However, in 1990 the Council enacted an ordinance that would hold them to vote for the council member who received the most votes. The council could move to override this ordinance and instead vote for whoever they wish, but has never done this.  Greenwald feared the council would repeal the ordinance to prevent her becoming mayor.).  She was the only member of the city council who voted against Target, though the majority (51.5%) in Davis voted in favor of Second Street Crossing coming to Davis. Greenwald was also the only council member to vote against Covell Village, where her position was favored by a majority (59.9%) of Davis voters. Greenwald won reelection in the 2008 City Council Election.

Sue has served on many boards and commissions during her twelve years on the council.   Among them are Board of Directors of the Capitol Corridor Joint Powers Authority, the Yolo County Transportation District,  the Lower Putah Creek Coordinating Committee, and the Yolo County Water Resources Association.  Sue has served on the League of California Cities Revenue and Taxation Policy Committee and Environmental Quality Policy Committee. She has served as council liaison to numerous city commissions and task forces, including the Finance and Budget Commission, Recreation and Parks Commission, Open Space and Habitat Commission, Senior Citizen Commission, and Bicycle Advisory Commission. Sue has served on the CityCounty 2 by 2, the DBA 2 by 2, the City School 2 by 2 and the MayorUCD Chancellor 1 by 1.

Sue lives in one of Davis oldest neighborhoods near campus with her husband, Mike Syvanen, who is a professor of Medical Microbiology and Immunology Medical Microbiology at UC Davis. Sue holds a B.A. and an M.A. from UC Berkeley, and has worked as a social science researcher.  Sue was motivated to run for City Council due to her passion for good city planning and for good  government.   Sue has an independent and analytical approach to problem solving.   She takes the time to research each issue thoroughly, knows how to reach out to the top outside experts and receive the best advice on the major issues.  Her circle of advisors include University economists, hydrologists, and environmental attorneys.

Living close to downtown, Sue is fortunate in that she can use her car very little and can rely heavily on walking.   People comment that they see Sue walking around town frequently.   Sue would like to bring more housing within walking distance to downtown so that more people can lead a life less reliant on the auto.

Sue enjoys Birds and Bird Watching birdwatching, photography, reading, Picnic in the Park, meeting with friends downtown and an occasional game of The Unofficial Davis Scrabble Club Scrabble in her spare time.

 Sue Greenwalds Candidate Statement

I have a passion for good city planning.   As your councilmember for the last 12 years and as your Mayor from 2006 to 2008,  I have worked hard to achieve neighborhood preservation, slow growth, environmentally sensitive city planning, genuine fiscal responsibility and a walkable, bikeable, vibrant downtown.  I have consistently voted to uphold these principles throughout three council terms, helping to ensure that Davis retains our special quality of life.

I have often been the lone voice for fiscal responsibility, back when the damage was still avoidable.  When necessary, I have not been afraid to stand firm and to stand alone on these important issues.   I didn’t prevail on every issue, but I succeeded in saving the City of Davis serious money.

I have taken the initiative and consistently voted to:
  Preserve and protect the quality of life in all of our neighborhoods.
  Bring more arts and entertainment downtown to create an environment worthy of a great University town.
  Contain urban sprawl while supporting quality infill and student housing    projects.
  Promote important environmental initiatives, protect habitat and agricultural land, and reestablish Davis as a leader in environmentally innovative planning.
  Support Davis’ tradition of excellent schools and respect for diversity and commitment to human rights.  I supported Measure C, our school parcel tax.
  Support hightech industry.
  Take the hard, necessary steps to bring our expenditures down.

If you go to the page http://suegreenwald.org/fordavis/WhereSueStands.aspx Where Sue Stands,  you can read about  concrete initiatives that I have personally  taken to achieve these goals in the past, and initiatives that I hope to take during the next four years.

We are facing difficult times with our budget crisis.  At the same time, we have a new city manager, our finance director is leaving and we have many new councilmembers and candidates.   We are losing our institutional memory.   I feel that I can make a unique contribution by bringing that history, that experience, and that institutional memory to the council.  A strong council has a balance of experience, stability, change, and viewpoints.

I hope you will consider supporting me for reelection.

 Sue Greenwald

 Controversy

 2010 Confrontation with Ruth Asmundson

On January 26, during a city council meeting, Sue Greenwald had a confrontation with then Mayor Ruth Asmundson that ended when Admusndon was temporarily hospitalized with an anxiety attack.  Like all city council meetings, it was televised.  A clip was posted on Youtube by David Greenwald of the The Vanguard; see http://www.youtube.com/watch?vIKudm65Dhek this video.  In the clip, Greenwald repeatedly accused Asmundson of lying about what said in closed session; Asmundson denied that Greenwald had made a motion concerning an MOU during the closed session, but Greenwald said that she had in fact made a motion and that Ruth was lying.  Greenwald asked Lamar Heystek whether she had made such a motion, and he confirmed that she had.  A majority of the council attempted to end the discussion, but Heystek and Greenwald voted to continue.  Asmundson blamed Greenwald for having caused her anxiety attacks in the past.  Then, as Asmundson was taken away by ambulance, there was further chaos in the foyer; according to the Davis Vanguard, Greenwald accused Asmundson of faking it and vainly attempted to continue political debate with several people.  The incident was covered by http://davisvanguard.org/index.php?optioncom_content&viewarticle&id3209:councildissintegratesintobickeringasmayorfallsill the Davis Vanguard (also http://davisvanguard.org/index.php?optioncom_content&viewarticle&id3211:whathappenedtuesdaynightafterthecameraswereoff&catid58:budgetfiscal&Itemid79&cpage150 here), http://www.sacbee.com/topstories/story/2497100.html#mi_rssTop%20Stories the Sacramento Bee , http://www.news10.net/news/local/story.aspx?storyid74164 News 10 TV, http://www.davisenterprise.com/ArchivedStories0/behavecouncilanddoyourjob/ the Davis Enterprise (including http://www.davisenterprise.com/ArchivedStories0/citycouncilsbickeringhitsanewlow/ Bob Dunning), and other media.  At first Greenwald said that Asmundson owed her an apology, but later Greenwald apologized in a city council meeting.

Greenwald had also had earlier confrontations with Asmundson during city council meetings, one of which was also posted http://www.youtube.com/watch?vfCbWEbQTIKQ to YouTube.

 2012 election mailer

On May 8, 2012, Davisites received an attack mailer citing this event, urging voters not to reward bad behavior.  The mailer stated that it was paid for by Plumbers and Pipefitters Local 447, IBEW Local 340, Operating Engineers Local Union 3 in Support of Stephen Souza and in Opposition to Sue Greenwald.  All four of the other candidates in the 2012 City Council Election denounced the attack mailer and disavowed any association with it.  Ruth Asmundson also denounced the attack mailer.

When the mailers first appeared, Sue Greenwald wrote the following in the comments on http://davisvanguard.org/index.php?optioncom_content&viewarticle&id5340:unionattackmailertargetsgreenwald&Itemid83 the Vanguard article concerning the mailer:

 A number of months ago, Jim Burchill of James Burchill & Associates approached me and told me that he was representing a building trades union, and that they really wanted the surface water project. He told me that the union was going to spend serious PAC money in our city election, and that I could be the beneficiary if I supported the project.

 I interpreted this as a threat as well as a bribe, i.e., that if I didnt support the project, they were going to engage in serious Karl Rove style attacks.

 I tried to explain to Jim that my job was to represent the best interests of citizens of Davis, and that although I was extremely sympathetic to the union members desire for jobs, our combined wastewater, surface water project costs, water rights purchase and related new infrastructure was still around $300 million, which was very high for a city of our size, and that I felt obligated to pursue less expensive approaches.

 I am afraid that this deplorable action by an outside PAC will intimidate many good, brave citizens from going into local Davis politics for a long time to come.

The Vanguard soon determined that http://davisvanguard.org/index.php?optioncom_content&viewarticle&id5341:politicaloperativejamesburchillappearstobebehindattackmailer&Itemid83 James Burchill did indeed seem to be behind the attack mailers.  The Vanguard states, Among Mr. Burchills clients are not only the building trade unions who would be involved in the construction of the surface water project, but also Angelo Tsakopoulos himself who as the owner of the company that owns Conaway Ranch, stands to make huge amounts from the water deal approved back in December of 2010.  The Davis Enterprise also http://www.davisenterprise.com/localnews/city/unioncampaignmailerattacksgreenwald/ covered the attack mailers, calling this an ugly turn in the Davis City Council campaign.  Bob Dunning also http://www.davisenterprise.com/opinion/dunning/nastymailerbackfiresineveryway/ denounced the mailer, calling it  the dirtiest campaign piece hed ever received, bar none.

On May 15, 2012, http://davisvanguard.org/index.php?optioncom_content&viewarticle&id5361:newplayeremergesongreenwaldattackmailer&catid50:elections&Itemid83 the Vanguard exposed an email that strongly implied Jon Li had a role in the attack mailer. A subsequent phone conversation between David Greenwald and Jon Li confirmed that Li had a small part in crafting the mailer, but, according to Li, his help was limited to a few suggestions and phone discussions with Burchill on the subject. Li ran an unsuccessful campaign for the 2010 City Council Election.

 2012 Stealing and blocking campaign signs

Another controversy developed in which Greenwalds husband Mike Syvanen was seen block campaign signs with signs for Sue Greenwald, and also taking campaign signs, on Covell Boulevard.  The story was covered in http://davis.patch.com/articles/photoscampaignsigncontroversystirsupdaviscitycouncilrace Davis Patch, http://davisvanguard.org/index.php?optioncom_content&viewarticle&id5401:accusationsofsignstealingandsignblockingemergethisweekend&catid50:elections&Itemid83 the Davis Vanguard, http://sacramento.cbslocal.com/2012/05/29/signshenanigansdistractfromdaviscitycouncilrace/ CBS 13 TV, and http://www.davisenterprise.com/localnews/city/greenwaldshusbandblockedothercandidatessigns/ the Davis Enterprise.  Greenwald said that Syvanen acted out of frustration because so many of her lawn signs had disappeared, that she had warned him to be careful lest political operatives exploit his actions, and that much of the evidence against Syvanen had been staged.

 Endorsements

Sue Greenwald was endorsed by the The Davis Enterprise for the 2012 City Council Election.

Sue http://media.www.californiaaggie.com/media/storage/paper981/news/2008/02/01/CampusNews/Councilmembers.Endorse.Primary.Candidates3182296.shtml endorsed Barack Obama for the February 2008 Primary Election and she wrote a http://media.www.californiaaggie.com/media/storage/paper981/news/2008/02/05/Opinion/Letters.To.The.Editor3188164.shtml letter to the editor of the California Aggie.

 Campaign Literature

 2012

 2008

20080612 05:53:56 nbsp Hi Sue I once waved hi to you on campus and you waved back!  You had an awesome 4th of July speech in 2007.  Go get em! Users/DavidLiu

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
SUBROUTINE AXPY(N1,N2,N3,localInterval,a,X,Y,Z)

  INTEGER(KIND=8), INTENT(IN)  :: N1, N2, N3, localInterval(6)
  REAL(KIND=8),    INTENT(IN)  :: a,X(N1*N2*N3)
  REAL(KIND=8),    INTENT(IN)  :: Y(N1*N2*N3)
  REAL(KIND=8),    INTENT(OUT) :: Z(N1*N2*N3)

  INTEGER(KIND=8) :: I, J, K, LI, LIK, LIJ, IS1, IE1
  INTEGER(KIND=8) :: xyPlane
  INTEGER(KIND=8) :: S1,E1,S2,E2,S3,E3

  S1 = localInterval(1)
  E1 = localInterval(2)
  S2 = localInterval(3)
  E2 = localInterval(4)
  S3 = localInterval(5)
  E3 = localInterval(6)

  yzPlane = (N2*N3)

  DO I = S1, E1
     !LI = (I-1)
     DO J = S2, E2
        LIJ = (J-1)*N1 + I
        DO K = S3, E3
           LIJK = (K-1)*yzPlane + LIJ
           Z(LIJK) = a*X(LIJK) + Y(LIJK)
        END DO
     END DO
  END DO

END SUBROUTINE AXPY

max_stars_repo_name: MTCam/OKLib

max_issues_repo_name: MTCam/OKLib

max_forks_repo_name: MTCam/OKLib


text:
	SUBROUTINE GPWOPT ( strmid, iret )
C************************************************************************
C* GPWOPT								*
C*									*
C* This subroutine displays a summary of the users options.		*
C*									*
C* GPWOPT ( STRMID, IRET )						*
C*									*
C* Input parameters:							*
C*	STRMID		CHAR*		Storm identifier		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -1 = exit			*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	 8/01	Copied from GPPOPT                      *
C* m.gamazaychikov/SAIC	03/04	Copied from GPKOPT                      *
C************************************************************************
	CHARACTER*(*)	strmid
C*
	LOGICAL		respnd
C------------------------------------------------------------------------
C
	iret = 0
	WRITE ( 6, 5000) strmid
C
	CALL IP_RESP ( respnd, ier )
	IF  ( respnd ) THEN
	    CALL TM_ACCP ( ier )
	    IF ( ier .eq. 2 ) iret = -1
	END IF
C
5000	FORMAT ( ' GPTCWW PARAMETERS:',//
     +           ' Strmid:              ', A )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak