#include "wwm_functions.h"
! I5 is under the assumptions that we have a lot of memory
! and tries to minimize the number of MPI exchanges and
! to have the processors be as busy as possible
! We use memory ordered as AC(MNP,MSC,MDC)
! I5B is the same as I5. We use memory ordered as AC(MSC,MDC,MNP)
! so reordering at the beginning but less operations later on.
#undef DEBUG
!#define DEBUG
! This is for the reordering of ASPAR_pc and hopefully higher speed
! in the application of the preconditioner.
#undef REORDER_ASPAR_PC
#define REORDER_ASPAR_PC
! This is for the computation of ASPAR_block by a block algorithm
! with hopefully higher speed.
#undef ASPAR_B_COMPUTE_BLOCK
#define ASPAR_B_COMPUTE_BLOCK
! Either we use the SCHISM exchange routine or ours that exchanges only
! the ghost nodes and not the interface nodes.
#undef NO_SELFE_EXCH
#define NO_SELFE_EXCH
! Repeated CX/CY computations but less memory used.
#undef NO_MEMORY_CX_CY
#define NO_MEMORY_CX_CY
! For the SOR preconditioner, we can actually compute directly
! from the matrix since it is so simple.
#undef SOR_DIRECT
#define SOR_DIRECT
!
#undef SINGLE_LOOP_AMATRIX
#define SINGLE_LOOP_AMATRIX
!
! More complexity! Some options excludes other!
!
#if defined REORDER_ASPAR_PC && defined SOR_DIRECT
# undef REORDER_ASPAR_PC
#endif
!**********************************************************************
!* We have to think on how the system is solved. Many questions are *
!* mixed: the ordering of the nodes, the ghost nodes, the aspar array *
!* Here is a repository of the conclusions that have been reached *
!* *
!* Ordering 1> should be that way: We have two global nodes i and j. *
!* -- If i and j belong to a common local grid, then we select the *
!* grid of lowest color and decide whether ipgl(i) < ipgl(j) *
!* -- If i and j belong to two different grid then *
!* ---If Color(i) < Color(j) or reverse we decide by that *
!* ---If Color(i) = Color(j) we decide by i<j or not (but it *
!* does not matter to the solution) *
!* The functions WRITE_EXPLICIT_ORDERING does exactly that and *
!* provides an ordering that can be used. That is we start with the *
!* nodes of lowest color and index until we arrive at highest color *
!* *
!* The ASPAR is computed correctly only on 1:NP_RES but this can be *
!* extended by exchange routines. *
!* We Compute on the resident nodes only. This means loops over *
!* IP=1,NP_RES and backwards. This means that we do not have to do *
!* exchanges of ASPAR values. Only the resident nodes are sent. *
!* This is smaller and this is all that we ever need. *
!* *
!* WRONG APPROACHES: *
!* to use all nodes 1:MNP may look simpler but it forces to have the *
!* following property of the NNZ, IA, JA arrays. If two vertices *
!* i and j are adjacent in a grid G, then they are adjacent in ANY *
!* of the grid in which they are both contained. *
!* This property is actually not satisfied in general. *
!* We could extend the IA, JA arrays by *
!* adding some vertices but that looks quite hazardous idea and it *
!* it is actually not needed by the ILU0 preconditioner and other *
!* *
!* PROBLEM: *
!* There is an asymmetry in the construction of the ordering. *
!* We start from low colors and upwards. If we had started with *
!* high colors and gone downwards, then we get a different ordering *
!* (even if we take the opposite, because of the interface nodes) *
!* This requires the construction of many mappings. *
!* Our approach is actually to rebuild separate node sets. *
!* *
!* CHECKS: *
!* ---The sum of number of non-zero entries in Jstatus_L over all *
!* nodes should be equal to the sum of number of non-zero entries *
!* of Jstatus_U over all nodes. *
!* This is because number of upper diagonal entries should be *
!* equal to number of lower diagonal entries. *
!* ---We CANNOT have Jstatus_L(J)=1 and Jstatus_U(J)=1, i.e. a matrix *
!* entry cannot be both lower and upper. *
!* ---We have sum of Jstatus_L + sum J_status_U + np_global should *
!* be equal to NNZ_global *
!* *
!* So, procedure is as follows: *
!* ---compute ASPAR on 1,NP_RES nodes and no synchronization *
!* ---compute on IP=1,NP_RES for L solving *
!* ---export to grids of higher rank, the values on nodes 1,NP_RES *
!* only. The other ghost points have invalid values or are *
!* resident of other grids of lower rank. (at this stage, some *
!* ghost values are wrong but are not exported) *
!* ---export to grid of lower rank in order to correct their ghost *
!* values and get the value *
!* ---compute on IP=NP_RES,1,-1 for U solving *
!* ---export to grid of lower rank. *
!* Do everything similarly to L solve. *
!* *
!* The basic approach is that we compute at a node S if and only if *
!* it is a resident node. The twist come because some nodes are *
!* resident for TWO domains. This is why we have the CovLower *
!* We need to create disjoint domains for each node, so that *
!* we have sum sum(CovLower) = np_global *
!* *
!* At the end of the resolution of the system, we need to do the *
!* synchronization with respect to the unused nodes. *
!* Mystery?: When we apply the function, we do it on 1:NP_RES and *
!* then call synchronizer. So, this means we need to do the *
!* synchronization after the call to the preconditioner. *
!* But there may be space for improvements here. *
!* *
!* Description of specific exchange arrays: *
!* ---wwm_p2dsend_type/wwm_p2drecv_type *
!* The points of 1:NP_RES are sent to nodes that contained them *
!* length=1 *
!* ---wwmtot_p2dsend_type/wwmtot_p2drecv_type *
!* same as above but length=MSC*MDC *
!* ---blk_p2dsend_type/blk_p2drecv_type *
!* same as above but length=maxBlockLength for matrix exchanges *
!* ---wwmmat_p2dsend_type/wwmmat_p2drecv_type *
!* exchange of correct matrix elements, i.e. elements A(I,J) *
!* with I<=NP_RES *
!* length is 1. *
!* ---u2l_p2dsend_type/u2l_p2drecv_type *
!* upper 2 lower exchange arrays, depends on CovLower, so on *
!* the coloring chosen. length=maxBlockLength *
!* exchange are from upper to lower. *
!* ---sync_p2dsend_type/sync_p2drecv_type *
!* synchronize value, i.e. the CovLower=1 values are send to all *
!* nodes. Length is MSC*MDC *
!**********************************************************************
#if defined WWM_SOLVER && defined MPI_PARALL_GRID
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EXCHANGE_P4D_WWM(LocalColor, AC)
USE DATAPOOL, only : MSC, MDC, rkind, LocalColorInfo
USE DATAPOOL, only : wwm_nnbr_send, wwm_nnbr_recv
USE DATAPOOL, only : wwm_ListNeigh_send, wwm_ListNeigh_recv
USE DATAPOOL, only : wwmtot_p2dsend_type, wwmtot_p2drecv_type
USE DATAPOOL, only : wwm_p2dsend_type, wwm_p2drecv_type
USE DATAPOOL, only : wwm_p2dsend_rqst, wwm_p2drecv_rqst
USE DATAPOOL, only : wwm_p2dsend_stat, wwm_p2drecv_stat
USE DATAPOOL, only : ZERO, NP_RES, MNP
USE datapool, only : comm, ierr, myrank
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
real(rkind), intent(inout) :: AC(MSC,MDC,MNP)
integer iSync, iRank
integer IS, ID, IP
real(rkind) SumErr, Lerror
real(rkind) :: ACtest(MSC,MDC,MNP)
real(rkind) :: U1(MNP), U2(MNP)
CALL I5B_TOTAL_COHERENCY_ERROR_NPRES(MSC, AC, Lerror)
! Print *, 'NP_RES cohenrency error=', Lerror
#ifdef DEBUG
WRITE(740+myrank,*) 'I5B_EXCHANGE_P4D_WWM, begin, sum(AC)=', sum(AC)
FLUSH(740+myrank)
#endif
ACtest=AC
SumErr=ZERO
DO IS=1,MSC
DO ID=1,MDC
U1=AC(IS,ID,:)
U2=AC(IS,ID,:)
DO iSync=1,wwm_nnbr_send
iRank=wwm_ListNeigh_send(iSync)
CALL mpi_isend(U1, 1, wwm_p2dsend_type(iSync), iRank-1, 1020, comm, wwm_p2dsend_rqst(iSync), ierr)
END DO
DO iSync=1,wwm_nnbr_recv
iRank=wwm_ListNeigh_recv(iSync)
call mpi_irecv(U2,1,wwm_p2drecv_type(iSync),iRank-1,1020,comm,wwm_p2drecv_rqst(iSync),ierr)
END DO
IF (wwm_nnbr_send > 0) THEN
call mpi_waitall(wwm_nnbr_send, wwm_p2dsend_rqst, wwm_p2dsend_stat,ierr)
END IF
IF (wwm_nnbr_recv > 0) THEN
call mpi_waitall(wwm_nnbr_recv, wwm_p2drecv_rqst, wwm_p2drecv_stat,ierr)
END IF
DO IP=1,NP_RES
SumErr=SumErr + abs(U1(IP) - U2(IP))
END DO
END DO
END DO
#ifdef DEBUG
WRITE(740+myrank,*) 'Total SumErr=', SumErr
FLUSH(740+myrank)
#endif
DO iSync=1,wwm_nnbr_send
iRank=wwm_ListNeigh_send(iSync)
CALL mpi_isend(ACtest, 1, wwmtot_p2dsend_type(iSync), iRank-1, 1020, comm, wwm_p2dsend_rqst(iSync), ierr)
END DO
DO iSync=1,wwm_nnbr_recv
iRank=wwm_ListNeigh_recv(iSync)
call mpi_irecv(AC,1,wwmtot_p2drecv_type(iSync),iRank-1,1020,comm,wwm_p2drecv_rqst(iSync),ierr)
END DO
IF (wwm_nnbr_send > 0) THEN
call mpi_waitall(wwm_nnbr_send, wwm_p2dsend_rqst, wwm_p2dsend_stat,ierr)
END IF
IF (wwm_nnbr_recv > 0) THEN
call mpi_waitall(wwm_nnbr_recv, wwm_p2drecv_rqst, wwm_p2drecv_stat,ierr)
END IF
SumErr=ZERO
DO IS=1,MSC
DO ID=1,MDC
DO IP=1,NP_RES
SumErr=SumErr + abs(AC(IS,ID,IP) - ACtest(IS,ID,IP))
END DO
END DO
END DO
#ifdef DEBUG
WRITE(740+myrank,*) 'SumErr=', SumErr
FLUSH(740+myrank)
#endif
#ifdef DEBUG
WRITE(740+myrank,*) 'I5B_EXCHANGE_P4D_WWM, end, sum(AC)=', sum(AC)
FLUSH(740+myrank)
#endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EXCHANGE_SL_WWM(LocalColor, AC)
USE DATAPOOL, only : MSC, MDC, rkind, LocalColorInfo
USE DATAPOOL, only : wwm_nnbr_send_sl, wwm_nnbr_recv_sl
USE DATAPOOL, only : wwm_ListNeigh_send_sl, wwm_ListNeigh_recv_sl, wwm_ListNeigh_send
USE DATAPOOL, only : wwmsl_send_type, wwmsl_recv_type
USE DATAPOOL, only : wwmsl_send_rqst, wwmsl_recv_rqst
USE DATAPOOL, only : wwmsl_send_stat, wwmsl_recv_stat
USE DATAPOOL, only : ZERO, NP_RES, MNP
USE datapool, only : comm, ierr, myrank
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
real(rkind), intent(inout) :: AC(MSC,MDC,MNP)
integer iSync, iRank
DO iSync=1,wwm_nnbr_send_sl
iRank=wwm_ListNeigh_send(iSync)
CALL mpi_isend(AC, 1, wwmsl_send_type(iSync), iRank-1, 1020, comm, wwmsl_send_rqst(iSync), ierr)
END DO
DO iSync=1,wwm_nnbr_recv_sl
iRank=wwm_ListNeigh_recv_sl(iSync)
call mpi_irecv(AC,1,wwmsl_recv_type(iSync),iRank-1,1020,comm,wwmsl_recv_rqst(iSync),ierr)
END DO
IF (wwm_nnbr_send_sl > 0) THEN
call mpi_waitall(wwm_nnbr_send_sl, wwmsl_send_rqst, wwmsl_send_stat,ierr)
END IF
IF (wwm_nnbr_recv_sl > 0) THEN
call mpi_waitall(wwm_nnbr_recv_sl, wwmsl_recv_rqst, wwmsl_recv_stat,ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EXCHANGE_ASPAR(LocalColor, ASPAR_bl)
USE DATAPOOL, only: comm, ierr, myrank, ierr, wwm_nnbr_m_send, wwm_ListNeigh_m_send
use datapool, only: wwmmat_p2dsend_type, wwmmat_p2dsend_rqst, wwm_nnbr_m_recv, wwm_ListNeigh_m_recv
use datapool, only: wwmmat_p2drecv_type, wwmmat_p2drecv_rqst, LocalColorInfo, rkind
use datapool, only: wwmmat_p2drecv_stat, wwmmat_p2dsend_stat, MSC, MDC, NNZ
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
real(rkind), intent(inout) :: ASPAR_bl(MSC,MDC,NNZ)
integer I, iProc
do I=1,wwm_nnbr_m_send
iProc=wwm_ListNeigh_m_send(I)
call mpi_isend(ASPAR_bl,1,wwmmat_p2dsend_type(I),iProc,991,comm,wwmmat_p2dsend_rqst(i),ierr)
enddo
do I=1,wwm_nnbr_m_recv
iProc=wwm_ListNeigh_m_recv(I)
call mpi_irecv(ASPAR_bl,1,wwmmat_p2drecv_type(I),iProc,991,comm,wwmmat_p2drecv_rqst(i),ierr)
enddo
IF (wwm_nnbr_m_recv .gt. 0) THEN
call mpi_waitall(wwm_nnbr_m_recv,wwmmat_p2drecv_rqst,wwmmat_p2drecv_stat,ierr)
END IF
IF (wwm_nnbr_m_send .gt. 0) THEN
call mpi_waitall(wwm_nnbr_m_send,wwmmat_p2dsend_rqst,wwmmat_p2dsend_stat,ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE COLLECT_ALL_COVLOWER(LocalColor)
USE DATAPOOL
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
integer, allocatable :: rbuf_int(:)
integer len, iProc, IP, idx, sumMNP
sumMNP=sum(ListMNP)
allocate(LocalColor % ListCovLower(sumMNP), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 18')
IF (myrank == 0) THEN
idx=0
DO IP=1,MNP
idx=idx+1
LocalColor % ListCovLower(idx)=LocalColor % CovLower(IP)
END DO
DO iProc=2,nproc
len=ListMNP(iProc)
allocate(rbuf_int(len), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 19')
CALL MPI_RECV(rbuf_int,len,itype, iProc-1, 809, comm, istatus, ierr)
DO IP=1,len
idx=idx+1
LocalColor % ListCovLower(idx)=rbuf_int(IP)
END DO
deallocate(rbuf_int)
END DO
DO iProc=2,nproc
CALL MPI_SEND(LocalColor % ListCovLower,sumMNP,itype, iProc-1, 811, comm, ierr)
END DO
ELSE
CALL MPI_SEND(LocalColor % CovLower,MNP,itype, 0, 809, comm, ierr)
CALL MPI_RECV(LocalColor % ListCovLower,sumMNP,itype, 0, 811, comm, istatus, ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE CREATE_WWM_P2D_EXCH
USE DATAPOOL
implicit none
integer :: ListFirst(nproc)
integer :: ListCommon_recv(nproc)
integer :: ListCommon_send(nproc)
integer :: ListMapped(np_global)
integer :: ListMappedB(np_global)
integer, allocatable :: dspl_send(:), dspl_recv(:)
integer, allocatable :: dspl_send_tot(:), dspl_recv_tot(:)
integer IP, IP_glob, iProc, MNPloc, idx, NP_RESloc
integer iNeigh, IPmap, nbCommon
integer nbCommon_send, nbCommon_recv, idx_send, idx_recv
integer nbCommon_send_sl, nbCommon_recv_sl
integer sumNbCommon_send, sumNbCommon_recv
integer idxDspl_send, idxDspl_recv
ListFirst=0
DO iProc=2,nproc
ListFirst(iProc)=ListFirst(iProc-1) + ListMNP(iProc-1)
END DO
!
! First the arrays so that coordinates 1:NP_RES got send
! to all nodes 1:MNP of neighboring components
!
ListMapped=0
DO IP=1,MNP
IP_glob=iplg(IP)
ListMapped(IP_glob)=IP
END DO
wwm_nnbr_send=0
wwm_nnbr_recv=0
ListCommon_send=0
ListCommon_recv=0
DO iProc=1,nproc
IF (iPROC .ne. myrank+1) THEN
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMappedB=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
ListMappedB(IP_glob)=IP
END DO
!
nbCommon_recv=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
IF (ListMapped(IP_glob).gt.0) THEN
nbCommon_recv=nbCommon_recv+1
END IF
END DO
IF (nbCommon_recv .gt. 0) THEN
wwm_nnbr_recv=wwm_nnbr_recv+1
ListCommon_recv(iProc)=nbCommon_recv
END IF
!
nbCommon_send=0
DO IP=1,NP_RES
IP_glob=iplg(IP)
IF (ListMappedB(IP_glob).gt.0) THEN
nbCommon_send=nbCommon_send+1
END IF
END DO
IF (nbCommon_send .gt. 0) THEN
wwm_nnbr_send=wwm_nnbr_send+1
ListCommon_send(iProc)=nbCommon_send
END IF
END IF
END DO
allocate(wwm_ListNbCommon_send(wwm_nnbr_send), wwm_ListNbCommon_recv(wwm_nnbr_recv), wwm_ListNeigh_send(wwm_nnbr_send), wwm_ListNeigh_recv(wwm_nnbr_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 27')
idx_send=0
idx_recv=0
sumNbCommon_send=0
sumNbCommon_recv=0
DO iProc=1,nproc
IF (ListCommon_send(iProc) .gt. 0) THEN
idx_send=idx_send+1
wwm_ListNeigh_send(idx_send)=iProc
nbCommon=ListCommon_send(iProc)
wwm_ListNbCommon_send(idx_send)=nbCommon
sumNbCommon_send=sumNbCommon_send+nbCommon
END IF
IF (ListCommon_recv(iProc) .gt. 0) THEN
idx_recv=idx_recv+1
wwm_ListNeigh_recv(idx_recv)=iProc
nbCommon=ListCommon_recv(iProc)
wwm_ListNbCommon_recv(idx_recv)=nbCommon
sumNbCommon_recv=sumNbCommon_recv+nbCommon
END IF
END DO
allocate(wwm_ListDspl_send(sumNbCommon_send), wwm_ListDspl_recv(sumNbCommon_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 28')
wwm_nnbr=0
DO iProc=1,nproc
IF ((ListCommon_send(iProc).gt.0).or.(ListCommon_recv(iProc).gt.0)) THEN
wwm_nnbr=wwm_nnbr+1
END IF
END DO
allocate(wwm_ListNeigh(wwm_nnbr), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 29')
idx=0
DO iProc=1,nproc
IF ((ListCommon_send(iProc).gt.0).or.(ListCommon_recv(iProc).gt.0)) THEN
idx=idx+1
wwm_ListNeigh(idx)=iProc
END IF
END DO
allocate(wwm_p2dsend_rqst(wwm_nnbr_send), wwm_p2drecv_rqst(wwm_nnbr_recv), &
&wwm_p2dsend_stat(MPI_STATUS_SIZE,wwm_nnbr_send), wwm_p2drecv_stat(MPI_STATUS_SIZE,wwm_nnbr_recv), &
&wwm_p2dsend_type(wwm_nnbr_send), wwm_p2drecv_type(wwm_nnbr_recv), &
&wwmtot_p2dsend_type(wwm_nnbr_send), wwmtot_p2drecv_type(wwm_nnbr_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 30')
idxDspl_send=0
DO iNeigh=1,wwm_nnbr_send
iProc=wwm_ListNeigh_send(iNeigh)
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMappedB=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
ListMappedB(IP_glob)=IP
END DO
nbCommon=wwm_ListNbCommon_send(iNeigh)
allocate(dspl_send(nbCommon), dspl_send_tot(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 31')
idx=0
DO IP=1,NP_RES
IP_glob=iplg(IP)
IF (ListMappedB(IP_glob).gt.0) THEN
IPmap=ListMappedB(IP_glob)
idx=idx+1
dspl_send(idx)=IP-1
dspl_send_tot(idx)=MSC*MDC*(IP-1)
idxDspl_send=idxDspl_send+1
wwm_ListDspl_send(idxDspl_send)=IP
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,1,dspl_send,rtype,wwm_p2dsend_type(iNeigh), ierr)
call mpi_type_commit(wwm_p2dsend_type(iNeigh), ierr)
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_send_tot,rtype,wwmtot_p2dsend_type(iNeigh), ierr)
call mpi_type_commit(wwmtot_p2dsend_type(iNeigh), ierr)
deallocate(dspl_send, dspl_send_tot)
END DO
idxDspl_recv=0
DO iNeigh=1,wwm_nnbr_recv
iProc=wwm_ListNeigh_recv(iNeigh)
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
nbCommon=wwm_ListNbCommon_recv(iNeigh)
allocate(dspl_recv(nbCommon), dspl_recv_tot(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 32')
idx=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
IF (ListMapped(IP_glob).gt.0) THEN
IPmap=ListMapped(IP_glob)
idx=idx+1
dspl_recv(idx)=IPmap-1
dspl_recv_tot(idx)=MSC*MDC*(IPmap-1)
idxDspl_recv=idxDspl_recv+1
wwm_ListDspl_recv(idxDspl_recv)=IPmap
#ifdef DEBUG
WRITE(800+myrank,*) 'Recv IP=', IPmap, 'IPglob=', IP_glob
FLUSH(800+myrank)
#endif
END IF
END DO
!
call mpi_type_create_indexed_block(nbCommon,1,dspl_recv,rtype,wwm_p2drecv_type(iNeigh),ierr)
call mpi_type_commit(wwm_p2drecv_type(iNeigh), ierr)
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_recv_tot,rtype,wwmtot_p2drecv_type(iNeigh),ierr)
call mpi_type_commit(wwmtot_p2drecv_type(iNeigh), ierr)
deallocate(dspl_recv, dspl_recv_tot)
END DO
!
! First the arrays so that coordinates 1:NP_RES got send
! to all nodes NP_RES+1:MNP of neighboring components
! "sl" : "super local"
!
ListMapped=0
DO IP=NP_RES+1,MNP
IP_glob=iplg(IP)
ListMapped(IP_glob)=IP
END DO
wwm_nnbr_send_sl=0
wwm_nnbr_recv_sl=0
ListCommon_send=0
ListCommon_recv=0
DO iNeigh=1,wwm_nnbr
iProc=wwm_ListNeigh(iNeigh)
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMappedB=0
DO IP=NP_RESloc+1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
ListMappedB(IP_glob)=IP
END DO
!
nbCommon_recv_sl=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
IF (ListMapped(IP_glob).gt.0) THEN
nbCommon_recv_sl=nbCommon_recv_sl+1
END IF
END DO
IF (nbCommon_recv_sl .gt. 0) THEN
wwm_nnbr_recv_sl=wwm_nnbr_recv_sl+1
ListCommon_recv(iProc)=nbCommon_recv_sl
END IF
!
nbCommon_send_sl=0
DO IP=1,NP_RES
IP_glob=iplg(IP)
IF (ListMappedB(IP_glob).gt.0) THEN
nbCommon_send_sl=nbCommon_send_sl+1
END IF
END DO
IF (nbCommon_send_sl .gt. 0) THEN
wwm_nnbr_send_sl=wwm_nnbr_send_sl+1
ListCommon_send(iProc)=nbCommon_send_sl
END IF
END DO
allocate(wwm_ListNbCommon_send_sl(wwm_nnbr_send_sl), wwm_ListNbCommon_recv_sl(wwm_nnbr_recv_sl), wwm_ListNeigh_send_sl(wwm_nnbr_send_sl), wwm_ListNeigh_recv_sl(wwm_nnbr_recv_sl), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 33')
idx_send=0
idx_recv=0
DO iProc=1,nproc
nbCommon=ListCommon_send(iProc)
IF (nbCommon .gt. 0) THEN
idx_send=idx_send+1
wwm_ListNeigh_send_sl(idx_send)=iProc
wwm_ListNbCommon_send_sl(idx_send)=nbCommon
END IF
nbCommon=ListCommon_recv(iProc)
IF (nbCommon .gt. 0) THEN
idx_recv=idx_recv+1
wwm_ListNeigh_recv_sl(idx_recv)=iProc
wwm_ListNbCommon_recv_sl(idx_recv)=nbCommon
END IF
END DO
allocate(wwmsl_send_rqst(wwm_nnbr_send_sl), wwmsl_recv_rqst(wwm_nnbr_recv_sl), &
&wwmsl_send_stat(MPI_STATUS_SIZE,wwm_nnbr_send_sl), wwmsl_recv_stat(MPI_STATUS_SIZE,wwm_nnbr_recv_sl), &
&wwmsl_send_type(wwm_nnbr_send_sl), wwmsl_recv_type(wwm_nnbr_recv_sl), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 35')
idxDspl_send=0
DO iNeigh=1,wwm_nnbr_send_sl
iProc=wwm_ListNeigh_send_sl(iNeigh)
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMappedB=0
DO IP=NP_RESloc+1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
ListMappedB(IP_glob)=IP
END DO
nbCommon=wwm_ListNbCommon_send_sl(iNeigh)
allocate(dspl_send(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 36')
idx=0
DO IP=1,NP_RES
IP_glob=iplg(IP)
IF (ListMappedB(IP_glob).gt.0) THEN
IPmap=ListMappedB(IP_glob)
idx=idx+1
dspl_send(idx)=MSC*MDC*(IP-1)
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_send,rtype,wwmsl_send_type(iNeigh), ierr)
call mpi_type_commit(wwmsl_send_type(iNeigh), ierr)
deallocate(dspl_send)
END DO
idxDspl_recv=0
DO iNeigh=1,wwm_nnbr_recv_sl
iProc=wwm_ListNeigh_recv_sl(iNeigh)
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
nbCommon=wwm_ListNbCommon_recv_sl(iNeigh)
allocate(dspl_recv(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 37')
idx=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
IF (ListMapped(IP_glob).gt.0) THEN
IPmap=ListMapped(IP_glob)
idx=idx+1
dspl_recv(idx)=MSC*MDC*(IPmap-1)
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_recv,rtype,wwmsl_recv_type(iNeigh),ierr)
call mpi_type_commit(wwmsl_recv_type(iNeigh), ierr)
deallocate(dspl_recv)
END DO
END SUBROUTINE
!**********************************************************************
!* This subroutine creates array for exchanging matrix entries *
!* Process cannot get all the matrix entries right, and so they need *
!* to exchange data. *
!* Matrix entries A(i,j) with i <= NP_RES are correct and are *
!* exported *
!**********************************************************************
SUBROUTINE CREATE_WWM_MAT_P2D_EXCH
USE DATAPOOL
implicit none
integer :: ListFirstMNP(nproc), ListFirstNNZ(nproc)
integer :: ListCommon_send(nproc), ListCommon_recv(nproc)
integer :: ListMapped(np_global)
integer :: ListMappedB(np_global)
integer, allocatable :: dspl_send(:), dspl_recv(:)
integer nbCommon_send, nbCommon_recv
integer IAfirst
integer IP, JP, I, J, J2, IP_glob, JP_glob, iProc
integer MNPloc, NP_RESloc, JP_j
integer IPloc, JPloc, Jfound, idx
integer iNeigh, nbCommon
integer sumNbCommon_send, sumNbCommon_recv
integer idxDspl_send, idxDspl_recv
ListFirstNNZ=0
ListFirstMNP=0
DO iProc=2,nproc
ListFirstMNP(iProc)=ListFirstMNP(iProc-1) + ListMNP(iProc-1)
ListFirstNNZ(iProc)=ListFirstNNZ(iProc-1) + ListNNZ(iProc-1)
END DO
ListMapped=0
DO IP=1,MNP
IP_glob=iplg(IP)
ListMapped(IP_glob)=IP
END DO
wwm_nnbr_m_recv=0
wwm_nnbr_m_send=0
ListCommon_recv=0
DO I=1,wwm_nnbr_recv
iProc=wwm_ListNeigh_recv(I)
NP_RESloc=ListNP_RES(iProc)
MNPloc=ListMNP(iProc)
ListMappedB=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirstMNP(iProc))
ListMappedB(IP_glob)=IP
END DO
nbCommon_recv=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirstMNP(iProc))
IPloc=ListMapped(IP_glob)
IF (IPloc.gt.0) THEN
IAfirst=ListFirstMNP(iProc) + iProc-1
DO J=ListIA(IP+IAfirst),ListIA(IP+IAfirst+1)-1
JP=ListJA(J+ListFirstNNZ(iProc))
JP_glob=ListIPLG(JP+ListFirstMNP(iProc))
JPloc=ListMapped(JP_glob)
IF (JPloc.gt.0) THEN
JFOUND=-1
DO J2=IA(IPloc),IA(IPloc+1)-1
IF (JA(J2) == JPloc) THEN
JFOUND=J2
END IF
END DO
IF (JFOUND .gt. 0) THEN
nbCommon_recv=nbCommon_recv+1
END IF
END IF
END DO
END IF
END DO
IF (nbCommon_recv .gt. 0) THEN
wwm_nnbr_m_recv=wwm_nnbr_m_recv+1
ListCommon_recv(iProc)=nbCommon_recv
END IF
END DO
ListCommon_send=0
DO I=1,wwm_nnbr_send
iProc=wwm_ListNeigh_send(I)
NP_RESloc=ListNP_RES(iProc)
MNPloc=ListMNP(iProc)
ListMappedB=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirstMNP(iProc))
ListMappedB(IP_glob)=IP
END DO
nbCommon_send=0
DO IP=1,NP_RES
IP_glob=iplg(IP)
IPloc=ListMappedB(IP_glob)
IF (IPloc.gt.0) THEN
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
JP_glob=iplg(JP)
JPloc=ListMappedB(JP_glob)
IF (JPloc.gt.0) THEN
IAfirst=ListFirstMNP(iProc) + iProc-1
JFOUND=-1
DO J2=ListIA(IPloc+IAfirst),ListIA(IPloc+IAfirst+1)-1
JP_j=ListJA(J2+ListFirstNNZ(iProc))
IF (JP_j == JPloc) THEN
JFOUND=J2
END IF
END DO
IF (JFOUND .gt. 0) THEN
nbCommon_send=nbCommon_send+1
END IF
END IF
END DO
END IF
END DO
IF (nbCommon_send .gt. 0) THEN
wwm_nnbr_m_send=wwm_nnbr_m_send+1
ListCommon_send(iProc)=nbCommon_send
END IF
END DO
allocate(wwmmat_p2dsend_rqst(wwm_nnbr_m_send), wwmmat_p2drecv_rqst(wwm_nnbr_m_recv), &
&wwmmat_p2dsend_stat(MPI_STATUS_SIZE,wwm_nnbr_m_send), &
&wwmmat_p2drecv_stat(MPI_STATUS_SIZE,wwm_nnbr_m_recv), wwmmat_p2dsend_type(wwm_nnbr_m_send), &
&wwmmat_p2drecv_type(wwm_nnbr_m_recv), wwm_ListNbCommon_m_send(wwm_nnbr_m_send), &
&wwm_ListNbCommon_m_recv(wwm_nnbr_m_recv), wwm_ListNeigh_m_recv(wwm_nnbr_m_recv), &
&wwm_ListNeigh_m_send(wwm_nnbr_m_send), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 38')
idx=0
sumNbCommon_send=0
DO iProc=1,nproc
IF (ListCommon_send(iProc) .gt. 0) THEN
idx=idx+1
wwm_ListNeigh_m_send(idx)=iProc-1
nbCommon=ListCommon_send(iProc)
wwm_ListNbCommon_m_send(idx)=nbCommon
sumNbCommon_send=sumNbCommon_send+nbCommon
END IF
END DO
idx=0
sumNbCommon_recv=0
DO iProc=1,nproc
IF (ListCommon_recv(iProc) .gt. 0) THEN
idx=idx+1
wwm_ListNeigh_m_recv(idx)=iProc-1
nbCommon=ListCommon_recv(iProc)
wwm_ListNbCommon_m_recv(idx)=nbCommon
sumNbCommon_recv=sumNbCommon_recv+nbCommon
END IF
END DO
allocate(wwm_ListDspl_m_send(sumNbCommon_send), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 39')
idxDspl_send=0
DO iNeigh=1,wwm_nnbr_m_send
iProc=wwm_ListNeigh_m_send(iNeigh)+1
MNPloc=ListMNP(iProc)
ListMappedB=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirstMNP(iProc))
ListMappedB(IP_glob)=IP
END DO
nbCommon=wwm_ListNbCommon_m_send(iNeigh)
allocate(dspl_send(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 40')
idx=0
DO IP=1,NP_RES
IP_glob=iplg(IP)
IPloc=ListMappedB(IP_glob)
IF (IPloc.gt.0) THEN
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
JP_glob=iplg(JP)
JPloc=ListMappedB(JP_glob)
IF (JPloc .gt. 0) THEN
IAfirst=ListFirstMNP(iProc) + iProc-1
JFOUND=-1
DO J2=ListIA(IPloc+IAfirst),ListIA(IPloc+IAfirst+1)-1
JP_j=ListJA(J2+ListFirstNNZ(iProc))
IF (JP_j == JPloc) THEN
JFOUND=J2
END IF
END DO
IF (JFOUND .gt. 0) THEN
idxDspl_send=idxDspl_send+1
wwm_ListDspl_m_send(idxDspl_send)=J
idx=idx+1
dspl_send(idx)=(J-1)*MSC*MDC
END IF
END IF
END DO
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_send,rtype,wwmmat_p2dsend_type(iNeigh), ierr)
call mpi_type_commit(wwmmat_p2dsend_type(iNeigh), ierr)
deallocate(dspl_send)
END DO
allocate(wwm_ListDspl_m_recv(sumNbCommon_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 41')
idxDspl_recv=0
DO iNeigh=1,wwm_nnbr_m_recv
iProc=wwm_ListNeigh_m_recv(iNeigh)+1
nbCommon=wwm_ListNbCommon_m_recv(iNeigh)
allocate(dspl_recv(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 42')
NP_RESloc=ListNP_RES(iProc)
idx=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirstMNP(iProc))
IPloc=ListMapped(IP_glob)
IF (IPloc.gt.0) THEN
IAfirst=ListFirstMNP(iProc) + iProc-1
DO J=ListIA(IP+IAfirst),ListIA(IP+IAfirst+1)-1
JP=ListJA(J+ListFirstNNZ(iProc))
JP_glob=ListIPLG(JP+ListFirstMNP(iProc))
JPloc=ListMapped(JP_glob)
IF (JPloc.gt.0) THEN
JFOUND=-1
DO J2=IA(IPloc),IA(IPloc+1)-1
IF (JA(J2) == JPloc) THEN
JFOUND=J2
END IF
END DO
IF (JFOUND /= -1) THEN
idxDspl_recv=idxDspl_recv+1
wwm_ListDspl_m_recv(idxDspl_recv)=Jfound
idx=idx+1
dspl_recv(idx)=(Jfound-1)*MSC*MDC
END IF
END IF
END DO
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_recv,rtype,wwmmat_p2drecv_type(iNeigh), ierr)
call mpi_type_commit(wwmmat_p2drecv_type(iNeigh), ierr)
deallocate(dspl_recv)
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE BUILD_MULTICOLORING(AdjGraph, ListColor)
USE datapool, only : myrank, Graph
implicit none
type(Graph), intent(in) :: AdjGraph
integer, intent(out) :: ListColor(AdjGraph%nbVert)
integer, allocatable :: CurrColor(:)
integer MaxDeg, iVert, eVert, eColor, eDeg
integer idx, I, ChromaticNr, nbVert
integer, allocatable :: ListPosFirst(:)
integer eColorF, iVertFound, eAdjColor
integer nbUndef, MinDeg, eAdj, MinUndef, PosMin
integer istat
MaxDeg=AdjGraph % MaxDeg
nbVert=AdjGraph % nbVert
allocate(CurrColor(MaxDeg+1), ListPosFirst(nbVert), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 43')
ListColor=0
idx=0
DO iVert=1,nbVert
ListPosFirst(iVert)=idx
idx=idx+AdjGraph % ListDegree(iVert)
END DO
MinDeg=MaxDeg+3
PosMin=-1
DO iVert=1,nbVert
eDeg=AdjGraph % ListDegree(iVert)
IF (eDeg .lt. MinDeg) THEN
MinDeg=eDeg
PosMin=iVert
END IF
END DO
idx=ListPosFirst(PosMin)
DO I=0,MinDeg
IF (I.eq.0) THEN
eVert=PosMin
ELSE
eVert=AdjGraph % ListEdge(idx+I,2)
END IF
ListColor(eVert)=I+1
END DO
DO
MinUndef=nbVert
iVertFound=0
DO iVert=1,nbVert
IF (ListColor(iVert) == 0) THEN
idx=ListPosFirst(iVert)
eDeg=AdjGraph % ListDegree(iVert)
nbUndef=0
DO I=1,eDeg
eAdj=AdjGraph % ListEdge(idx+I,2)
eAdjColor=ListColor(eAdj)
IF (eAdjColor == 0) THEN
nbUndef=nbUndef+1
END IF
END DO
IF (nbUndef .lt. MinUndef) THEN
MinUndef=nbUndef
iVertFound=iVert
END IF
END IF
END DO
IF (iVertFound == 0) THEN
EXIT
END IF
eDeg=AdjGraph % ListDegree(iVertFound)
idx=ListPosFirst(iVertFound)
CurrColor=0
DO I=1,eDeg
eVert=AdjGraph % ListEdge(idx+I,2)
eColor=ListColor(eVert)
IF (eColor.gt.0) THEN
CurrColor(eColor)=1
END IF
END DO
eColorF=-1
DO I=1,MaxDeg+1
IF (eColorF == -1) THEN
IF (CurrColor(I) == 0) THEN
eColorF=I
END IF
END IF
END DO
ListColor(iVertFound)=eColorF
END DO
deallocate(ListPosFirst, CurrColor)
ChromaticNr=maxval(ListColor)
# ifdef DEBUG
WRITE(740+myrank,*) 'ChromaticNr=', ChromaticNr
DO iVert=1,nbVert
WRITE(740+myrank,*) 'iVert=', iVert, 'eColor=', ListColor(iVert)
END DO
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE DeallocateGraph(TheGraph)
USE DATAPOOL, only : Graph
implicit none
type(Graph), intent(inout) :: TheGraph
deallocate(TheGraph % ListDegree)
deallocate(TheGraph % ListEdge)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_BLOCK_FREQDIR(LocalColor, Nblock)
USE DATAPOOL, only : MNP, MDC, MSC, LocalColorInfo, rkind, stat
USE DATAPOOL, only : wwm_nnbr_send, wwm_nnbr_recv
USE DATAPOOL, only : wwm_ListNbCommon_send, wwm_ListNbCommon_recv
USE DATAPOOL, only : XP, YP, rtype, ierr, myrank, iplg
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
integer, intent(in) :: Nblock
integer Ntot, Hlen, Delta, iBlock, idx, ID, IS
integer lenBlock, maxBlockLength
integer istat
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'ENTERING INIT_BLOCK_FREQDIR'
FLUSH(STAT%FHNDL)
Ntot=MSC*MDC
Hlen=INT(MyREAL(Ntot)/Nblock)
Delta=Ntot - Hlen*Nblock
iBlock=1
idx=1
LocalColor % Nblock=Nblock
IF (Delta == 0) THEN
maxBlockLength=Hlen
ELSE
maxBlockLength=Hlen+1
ENDIF
allocate(LocalColor % ISindex(Nblock, maxBlockLength), LocalColor % IDindex(Nblock, maxBlockLength), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 44')
DO IS=1,MSC
DO ID=1,MDC
LocalColor % ISindex(iBlock, idx)=IS
LocalColor % IDindex(iBlock, idx)=ID
IF (iBlock <= Delta) THEN
lenBlock=Hlen+1
ELSE
lenBlock=Hlen
END IF
idx=idx+1
IF (idx > lenBlock) THEN
iBlock=iBlock+1
idx=1
ENDIF
END DO
END DO
allocate(LocalColor % BlockLength(Nblock), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 45')
DO iBlock=1,Nblock
IF (iBlock <= Delta) THEN
lenBlock=Hlen+1
ELSE
lenBlock=Hlen
END IF
LocalColor % BlockLength(iBlock)=lenBlock
END DO
LocalColor % maxBlockLength = maxBlockLength
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'FINISHING INIT_BLOCK_FREQDIR'
FLUSH(STAT%FHNDL)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_BLK_L2U_ARRAY(LocalColor)
USE DATAPOOL, only : MNP, MSC, MDC, LocalColorInfo, rkind, stat
USE DATAPOOL, only : wwm_nnbr_send, wwm_nnbr_recv
USE DATAPOOL, only : wwm_ListNbCommon_send, wwm_ListNbCommon_recv
USE DATAPOOL, only : wwm_ListDspl_send, wwm_ListDspl_recv
USE DATAPOOL, only : wwm_ListNeigh_send, wwm_ListNeigh_recv
USE DATAPOOL, only : XP, YP, rtype, ierr, myrank, iplg
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
integer, allocatable :: ListNeed(:), IdxRev(:)
integer nbNeedSend_blk, nbNeedRecv_blk
integer idx, IP
integer nbUpp_send, nbLow_recv, iUpp, iLow, iRank
integer maxBlockLength, idxSend, idxRecv
integer I, IC, nbCommon, eFirst
integer ListFirstCommon_send(wwm_nnbr_send)
integer ListFirstCommon_recv(wwm_nnbr_recv)
integer, allocatable :: dspl_send(:), dspl_recv(:)
integer istat
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'ENTERING INIT_BLK_L2U_ARRAY'
FLUSH(STAT%FHNDL)
maxBlockLength=LocalColor % maxBlockLength
ListFirstCommon_send=0
DO I=2,wwm_nnbr_send
ListFirstCommon_send(I)=ListFirstCommon_send(I-1)+wwm_ListNbCommon_send(I-1)
END DO
ListFirstCommon_recv=0
DO I=2,wwm_nnbr_recv
ListFirstCommon_recv(I)=ListFirstCommon_recv(I-1)+wwm_ListNbCommon_recv(I-1)
END DO
nbUpp_send=LocalColor % nbUpp_send
nbLow_recv=LocalColor % nbLow_recv
allocate(LocalColor % l2u_p2dsend_type(nbUpp_send), LocalColor % l2u_p2drecv_type(nbLow_recv), LocalColor % l2u_ListNeigh_send(nbUpp_send), LocalColor % l2u_ListNeigh_recv(nbLow_recv), ListNeed(MNP), IdxRev(MNP), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 46')
ListNeed=0
IdxRev=0
nbNeedSend_blk=0
DO iUpp=1,nbUpp_send
I=LocalColor % ListIdxUpper_send(iUpp)
iRank=wwm_ListNeigh_send(I)-1
LocalColor % l2u_ListNeigh_send(iUpp)=iRank
nbCommon=wwm_ListNbCommon_send(I)
eFirst=ListFirstCommon_send(I)
DO IC=1,nbCommon
idxSend=wwm_ListDspl_send(eFirst+IC)
IF (ListNeed(idxSend) .eq. 0) THEN
ListNeed(idxSend)=1
nbNeedSend_blk=nbNeedSend_blk+1
END IF
END DO
END DO
LocalColor % nbNeedSend_blk=nbNeedSend_blk
allocate(LocalColor % IdxSend_blk(nbNeedSend_blk), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 47')
idx=0
DO IP=1,MNP
IF (ListNeed(IP) .eq. 1) THEN
idx=idx+1
LocalColor % IdxSend_blk(idx)=IP
IdxRev(IP)=idx
END IF
END DO
DO iUpp=1,nbUpp_send
I=LocalColor % ListIdxUpper_send(iUpp)
nbCommon=wwm_ListNbCommon_send(I)
eFirst=ListFirstCommon_send(I)
ALLOCATE(dspl_send(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 48')
DO IC=1,nbCommon
idxSend=wwm_ListDspl_send(eFirst+IC)
dspl_send(IC)=(IdxRev(idxSend)-1)*maxBlockLength
END DO
call mpi_type_create_indexed_block(nbCommon,maxBlockLength,dspl_send,rtype,LocalColor % l2u_p2dsend_type(iUpp),ierr)
call mpi_type_commit(LocalColor % l2u_p2dsend_type(iUpp), ierr)
DEALLOCATE(dspl_send)
END DO
!
!
ListNeed=0
IdxRev=0
nbNeedRecv_blk=0
DO iLow=1,nbLow_recv
I=LocalColor % ListIdxLower_recv(iLow)
iRank=wwm_ListNeigh_recv(I)-1
LocalColor % l2u_ListNeigh_recv(iLow)=iRank
nbCommon=wwm_ListNbCommon_recv(I)
eFirst=ListFirstCommon_recv(I)
DO IC=1,nbCommon
idxRecv=wwm_ListDspl_recv(eFirst+IC)
IF (ListNeed(idxRecv) .eq. 0) THEN
ListNeed(idxRecv)=1
nbNeedRecv_blk=nbNeedRecv_blk+1
END IF
END DO
END DO
LocalColor % nbNeedRecv_blk=nbNeedRecv_blk
allocate(LocalColor % IdxRecv_blk(nbNeedRecv_blk), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 49')
idx=0
DO IP=1,MNP
IF (ListNeed(IP) .eq. 1) THEN
idx=idx+1
LocalColor % IdxRecv_blk(idx)=IP
IdxRev(IP)=idx
END IF
END DO
DO iLow=1,nbLow_recv
I=LocalColor % ListIdxLower_recv(iLow)
nbCommon=wwm_ListNbCommon_recv(I)
eFirst=ListFirstCommon_recv(I)
ALLOCATE(dspl_recv(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 50')
DO IC=1,nbCommon
idxRecv=wwm_ListDspl_recv(eFirst+IC)
dspl_recv(IC)=(IdxRev(idxRecv)-1)*maxBlockLength
END DO
call mpi_type_create_indexed_block(nbCommon,maxBlockLength,dspl_recv,rtype,LocalColor % l2u_p2drecv_type(iLow),ierr)
call mpi_type_commit(LocalColor % l2u_p2drecv_type(iLow), ierr)
DEALLOCATE(dspl_recv)
END DO
deallocate(ListNeed, IdxRev)
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'FINISHING INIT_BLK_L2U_ARRAY'
FLUSH(STAT%FHNDL)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE SYMM_INIT_COLORING(LocalColor, NbBlock)
USE DATAPOOL, only : LocalColorInfo, MNP, rkind, XP, YP, stat
USE DATAPOOL, only : DO_SYNC_UPP_2_LOW, DO_SYNC_LOW_2_UPP, DO_SYNC_FINAL
USE datapool, only : myrank, nproc, iplg, Graph
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
integer, intent(in) :: NbBlock
type(Graph) :: AdjGraph
integer :: ListColor(nproc)
integer :: ListColorWork(nproc)
integer istat
# ifdef DEBUG
integer TheRes
# endif
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'ENTERING SYMM_INIT_COLORING'
FLUSH(STAT%FHNDL)
# ifdef DEBUG
CALL COMPUTE_TOTAL_INDEX_SHIFT(TheRes)
WRITE(740+myrank,*) 'Total residual shift=', TheRes
# endif
CALL COLLECT_ALL_IA_JA
# ifdef DEBUG
WRITE(740+myrank,*) 'After COLLECT_ALL_IA_JA'
# endif
CALL CREATE_WWM_P2D_EXCH
# ifdef DEBUG
WRITE(740+myrank,*) 'After CREATE_WWM_P2D_EXCH'
# endif
# ifdef DEBUG
CALL CHECK_I5B_EXCHANGE(LocalColor)
# endif
CALL CREATE_WWM_MAT_P2D_EXCH
# ifdef DEBUG
WRITE(740+myrank,*) 'After CREATE_WWM_MAT_P2D_EXCH'
# endif
CALL SYMM_GRAPH_BUILD_ADJ(AdjGraph)
CALL BUILD_MULTICOLORING(AdjGraph, ListColor)
# ifdef DEBUG
WRITE(740+myrank,*) 'After BUILD_MULTICOLORING'
# endif
CALL DeallocateGraph(AdjGraph)
ListColorWork=-ListColor
allocate(LocalColor % ListColor(nproc), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 51')
LocalColor % ListColor=ListColorWork
# ifdef DEBUG
WRITE(740+myrank,*) 'Before INIT_LOW_2_UPP_ARRAYS'
# endif
CALL INIT_LOW_2_UPP_ARRAYS(LocalColor, ListColorWork)
!
# ifdef DEBUG
WRITE(740+myrank,*) 'Before CALL_BLOCK_FREQDIR'
# endif
CALL INIT_BLOCK_FREQDIR(LocalColor, NbBlock)
!
# ifdef DEBUG
WRITE(740+myrank,*) 'Before INIT_BLK_L2U_ARRAY'
# endif
CALL INIT_BLK_L2U_ARRAY(LocalColor)
!
# ifdef DEBUG
WRITE(740+myrank,*) 'Before COLLECT_ALL_COVLOWER'
# endif
CALL COLLECT_ALL_COVLOWER(LocalColor)
!
# ifdef DEBUG
WRITE(740+myrank,*) 'Before INIT_COVLOWER_ARRAY'
# endif
CALL INIT_COVLOWER_ARRAY(LocalColor)
!
# ifdef DEBUG
WRITE(740+myrank,*) 'Before DETERMINE_JSTATUS_L_U'
# endif
CALL DETERMINE_JSTATUS_L_U(LocalColor)
!
DO_SYNC_UPP_2_LOW=.TRUE.
DO_SYNC_LOW_2_UPP=.TRUE.
DO_SYNC_FINAL=.TRUE.
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'FINISHED WITH SYMM_INIT_COLORING'
FLUSH(STAT%FHNDL)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_LOW_2_UPP_ARRAYS(LocalColor, ListColor)
USE DATAPOOL, only : LocalColorInfo, MNP, rkind, stat
USE DATAPOOL, only : NNZ, IA, JA, NP_RES
USE DATAPOOL, only : wwm_ListNbCommon_send, wwm_ListNbCommon_recv
USE DATAPOOL, only : wwm_nnbr_send, wwm_nnbr_recv
USE DATAPOOL, only : wwm_p2drecv_type, wwm_p2dsend_type
USE DATAPOOL, only : wwm_ListNeigh_send, wwm_ListNeigh_recv
USE DATAPOOL, only : wwm_ListDspl_recv
USE datapool, only : myrank, nproc, comm, ierr, nbrrank_p
USE datapool, only : MPI_STATUS_SIZE
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
integer, intent(in) :: ListColor(nproc)
real(rkind) :: p2d_data_send(MNP)
real(rkind) :: CovLower(MNP), CovLower_meth2(MNP)
real(rkind) :: SumErr
integer eColor, fColor, I, iRank
integer nbUpp_send, nbLow_recv
integer iLow, iUpp
integer IC, eFirst, nbCommon, IPloc
integer ListFirstCommon_send(wwm_nnbr_send)
integer ListFirstCommon_recv(wwm_nnbr_recv)
integer istat
# ifdef DEBUG
integer IP
# endif
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'ENTERING INIT_LOW_2_UPP_ARRAYS'
FLUSH(STAT%FHNDL)
eColor=ListColor(myrank+1)
nbUpp_send=0
DO I=1,wwm_nnbr_send
iRank=wwm_ListNeigh_send(I)
fColor=ListColor(iRank)
# ifdef DEBUG
WRITE(740+myrank,*) 'I=', I, 'iRank=', iRank, 'fColor=', fColor
IF (fColor.eq.eColor) THEN
call wwm_abort('Major error in the code')
END IF
# endif
IF (fColor.gt.eColor) THEN
nbUpp_send=nbUpp_send + 1
ENDIF
END DO
nbLow_recv=0
DO I=1,wwm_nnbr_recv
iRank=wwm_ListNeigh_recv(I)
fColor=ListColor(iRank)
# ifdef DEBUG
IF (fColor.eq.eColor) THEN
call wwm_abort('Major error in the code')
END IF
# endif
IF (fColor.lt.eColor) THEN
nbLow_recv=nbLow_recv + 1
ENDIF
END DO
ListFirstCommon_send=0
DO I=2,wwm_nnbr_send
ListFirstCommon_send(I)=ListFirstCommon_send(I-1)+wwm_ListNbCommon_send(I-1)
END DO
ListFirstCommon_recv=0
DO I=2,wwm_nnbr_recv
ListFirstCommon_recv(I)=ListFirstCommon_recv(I-1)+wwm_ListNbCommon_recv(I-1)
END DO
# ifdef DEBUG
WRITE(740+myrank,*) 'SIC: nbLow_recv=', nbLow_recv, ' nbUpp_send=', nbUpp_send
# endif
LocalColor % nbUpp_send=nbUpp_send
LocalColor % nbLow_recv=nbLow_recv
allocate(LocalColor % ListIdxUpper_send(nbUpp_send), LocalColor % ListIdxLower_recv(nbLow_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 52')
iUpp=0
DO I=1,wwm_nnbr_send
iRank=wwm_ListNeigh_send(I)
fColor=ListColor(iRank)
IF (fColor.gt.eColor) THEN
iUpp=iUpp + 1
LocalColor % ListIdxUpper_send(iUpp)=I
ENDIF
END DO
iLow=0
DO I=1,wwm_nnbr_recv
iRank=wwm_ListNeigh_recv(I)
fColor=ListColor(iRank)
IF (fColor.lt.eColor) THEN
iLow=iLow + 1
LocalColor % ListIdxLower_recv(iLow)=I
ENDIF
END DO
allocate(LocalColor % Upp_s_rq(nbUpp_send), LocalColor % Upp_s_stat(MPI_STATUS_SIZE, nbUpp_send), LocalColor % Low_r_rq(nbLow_recv), LocalColor % Low_r_stat(MPI_STATUS_SIZE, nbLow_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 53')
p2d_data_send=0
CovLower=1
CovLower_meth2=1
DO iUpp=1,nbUpp_send
I=LocalColor % ListIdxUpper_send(iUpp)
iRank=wwm_ListNeigh_send(i)
call mpi_isend(p2d_data_send,1,wwm_p2dsend_type(i),iRank-1,13,comm,LocalColor % Upp_s_rq(iUpp),ierr)
END DO
DO iLow=1,nbLow_recv
I=LocalColor % ListIdxLower_recv(iLow)
iRank=wwm_ListNeigh_recv(i)
nbCommon=wwm_ListNbCommon_recv(I)
eFirst=ListFirstCommon_recv(I)
DO IC=1,nbCommon
IPloc=wwm_ListDspl_recv(eFirst+IC)
CovLower_meth2(IPloc)=0
END DO
call mpi_irecv(CovLower,1,wwm_p2drecv_type(i),iRank-1,13,comm,LocalColor % Low_r_rq(iLow),ierr)
END DO
IF (nbUpp_send > 0) THEN
call mpi_waitall(nbUpp_send, LocalColor%Upp_s_rq, LocalColor%Upp_s_stat,ierr)
END IF
IF (nbLow_recv > 0) THEN
call mpi_waitall(nbLow_recv, LocalColor%Low_r_rq, LocalColor%Low_r_stat,ierr)
END IF
allocate(LocalColor % CovLower(MNP), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 54')
LocalColor % CovLower=INT(CovLower)
SumErr=sum(abs(CovLower-CovLower_meth2))
# ifdef DEBUG
WRITE(740+myrank,*) 'SumErr(meth1/meth2) CovLower=', SumErr
WRITE(740+myrank,*) 'MNP=', MNP, ' sum(CovLower)=', sum(CovLower)
IF (SumErr .gt. 0) THEN
DO IP=1,MNP
WRITE(740+myrank,*) IP, CovLower(IP), CovLower_meth2(IP)
END DO
ENDIF
# endif
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'FINISHED WITH INIT_LOW_2_UPP_ARRAYS'
FLUSH(STAT%FHNDL)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_COVLOWER_ARRAY(LocalColor)
USE DATAPOOL
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
integer :: ListFirst(nproc)
integer :: ListCommon_recv(nproc)
integer :: ListCommon_send(nproc)
integer :: ListMapped0(np_global)
integer :: ListMapped1(np_global)
integer :: ListMapped0_B(np_global)
integer :: ListMapped1_B(np_global)
integer, allocatable :: dspl_send(:), dspl_recv(:)
integer IP, IP_glob, iProc, WeMatch, MNPloc, idx, NP_RESloc
integer iNeigh, IPmap, nbCommon
integer nbCommon_send, nbCommon_recv, idx_send, idx_recv
integer u2l_nnbr_send, u2l_nnbr_recv
integer sync_nnbr_send, sync_nnbr_recv
integer eCov, eColor, fColor
integer maxBlockLength
integer nbMap0, nbMap1
integer DoOper
integer, allocatable :: ListNeedSend(:), ListNeedRecv(:), IdxRev(:)
integer nbNeedSend_u2l, nbNeedRecv_u2l
integer lenMNP
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'ENTERING INIT_COVLOWER_ARRAY'
FLUSH(STAT%FHNDL)
ListFirst=0
DO iProc=2,nproc
ListFirst(iProc)=ListFirst(iProc-1) + ListMNP(iProc-1)
END DO
ListMapped0=0
ListMapped1=0
nbMap0=0
nbMap1=0
DO IP=1,MNP
IP_glob=iplg(IP)
eCov=LocalColor % CovLower(IP)
IF (eCov == 0) THEN
ListMapped0(IP_glob)=IP
nbMap0=nbMap0+1
END IF
IF (eCov == 1) THEN
ListMapped1(IP_glob)=IP
nbMap1=nbMap1+1
END IF
END DO
# ifdef DEBUG
WRITE(740+myrank,*) 'nbMap0=', nbMap0, ' nbMap1=', nbMap1
# endif
!
! First the Upper to lower (u2l) block arrays
!
u2l_nnbr_send=0
u2l_nnbr_recv=0
ListCommon_send=0
ListCommon_recv=0
eColor=LocalColor % ListColor(myrank+1)
allocate(ListNeedRecv(MNP), ListNeedSend(MNP), IdxRev(MNP), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 55')
ListNeedRecv=0
ListNeedSend=0
IdxRev=0
nbNeedSend_u2l=0
nbNeedRecv_u2l=0
# ifdef DEBUG
WRITE(740+myrank,*) 'U2L eColor=', eColor
# endif
DO iNeigh=1,wwm_nnbr
iProc=wwm_ListNeigh(iNeigh)
fColor=LocalColor % ListColor(iProc)
# ifdef DEBUG
WRITE(740+myrank,*) 'U2L iNeigh=', iNeigh, ' iProc=', iProc, 'fColor=', fColor
# endif
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
IF (fColor .ge. eColor) THEN
nbCommon_recv=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov .eq. 1) THEN
IPmap=ListMapped0(IP_glob)
WeMatch=0
IF (IPmap .gt. 0) THEN
WeMatch=1
ELSE
IPmap=ListMapped1(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap .gt. NP_RES)) THEN
WeMatch=1
END IF
END IF
IF (WeMatch .eq. 1) THEN
nbCommon_recv=nbCommon_recv+1
IF (ListNeedRecv(IPmap) .eq. 0) THEN
nbNeedRecv_u2l=nbNeedRecv_u2l+1
ListNeedRecv(IPmap)=1
END IF
END IF
END IF
END DO
IF (nbCommon_recv .gt. 0) THEN
u2l_nnbr_recv=u2l_nnbr_recv+1
ListCommon_recv(iProc)=nbCommon_recv
END IF
# ifdef DEBUG
WRITE(740+myrank,*) ' U2L nbCommon_recv=', nbCommon_recv
# endif
END IF
IF (fColor .le. eColor) THEN
nbCommon_send=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IPmap=ListMapped1(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap .le. NP_RES)) THEN
IF ((eCov .eq. 0).or.(IP.gt.NP_RESloc)) THEN
nbCommon_send=nbCommon_send+1
IF (ListNeedSend(IPmap) .eq. 0) THEN
nbNeedSend_u2l=nbNeedSend_u2l+1
ListNeedSend(IPmap)=1
END IF
END IF
END IF
END DO
IF (nbCommon_send .gt. 0) THEN
u2l_nnbr_send=u2l_nnbr_send+1
ListCommon_send(iProc)=nbCommon_send
END IF
# ifdef DEBUG
WRITE(740+myrank,*) ' U2L nbCommon_send=', nbCommon_send
# endif
END IF
END DO
# ifdef DEBUG
WRITE(740+myrank,*) 'WWM_P2D: u2l_nnbr_send=', u2l_nnbr_send
WRITE(740+myrank,*) 'WWM_P2D: u2l_nnbr_recv=', u2l_nnbr_recv
# endif
LocalColor % u2l_nnbr_send=u2l_nnbr_send
LocalColor % u2l_nnbr_recv=u2l_nnbr_recv
allocate(LocalColor % u2l_ListNbCommon_send(u2l_nnbr_send), LocalColor % u2l_ListNbCommon_recv(u2l_nnbr_recv), LocalColor % u2l_ListNeigh_send(u2l_nnbr_send), LocalColor % u2l_ListNeigh_recv(u2l_nnbr_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 56')
idx_send=0
idx_recv=0
DO iProc=1,nproc
IF (ListCommon_send(iProc) .gt. 0) THEN
idx_send=idx_send+1
LocalColor % u2l_ListNeigh_send(idx_send)=iProc-1
nbCommon=ListCommon_send(iProc)
LocalColor % u2l_ListNbCommon_send(idx_send)=nbCommon
END IF
IF (ListCommon_recv(iProc) .gt. 0) THEN
idx_recv=idx_recv+1
LocalColor % u2l_ListNeigh_recv(idx_recv)=iProc-1
nbCommon=ListCommon_recv(iProc)
LocalColor % u2l_ListNbCommon_recv(idx_recv)=nbCommon
END IF
END DO
LocalColor % nbNeedSend_u2l=nbNeedSend_u2l
LocalColor % nbNeedRecv_u2l=nbNeedRecv_u2l
allocate(LocalColor % IdxSend_u2l(nbNeedSend_u2l), LocalColor % IdxRecv_u2l(nbNeedRecv_u2l), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 57')
!
! Now creating the u2l exchange
!
# ifdef DEBUG
WRITE(740+myrank,*) 'WWM_P2D: wwm_nnbr=', wwm_nnbr
WRITE(740+myrank,*) 'WWM_P2D: wwm_ListNeigh built'
# endif
allocate(LocalColor % u2l_p2dsend_rqst(u2l_nnbr_send), LocalColor % u2l_p2drecv_rqst(u2l_nnbr_recv), &
&LocalColor % u2l_p2dsend_stat(MPI_STATUS_SIZE,u2l_nnbr_send), &
&LocalColor % u2l_p2drecv_stat(MPI_STATUS_SIZE,u2l_nnbr_recv), &
&LocalColor % u2l_p2dsend_type(u2l_nnbr_send), LocalColor % u2l_p2drecv_type(u2l_nnbr_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 58')
# ifdef DEBUG
WRITE(740+myrank,*) 'WWM_P2D: alloc done'
# endif
maxBlockLength=LocalColor % maxBlockLength
idx=0
DO IP=1,MNP
IF (ListNeedSend(IP) .eq. 1) THEN
idx=idx+1
LocalColor % IdxSend_u2l(idx)=IP
IdxRev(IP)=idx
END IF
END DO
DO iNeigh=1,u2l_nnbr_send
iProc=LocalColor % u2l_ListNeigh_send(iNeigh)+1
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMapped0_B=0
ListMapped1_B=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov == 0) THEN
ListMapped0_B(IP_glob)=IP
END IF
IF (eCov == 1) THEN
ListMapped1_B(IP_glob)=IP
END IF
END DO
nbCommon=LocalColor % u2l_ListNbCommon_send(iNeigh)
allocate(dspl_send(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 59')
idx=0
DO IP=1,NP_RES
IF (LocalColor % CovLower(IP) .eq. 1) THEN
IP_glob=iplg(IP)
IPmap=ListMapped0_B(IP_glob)
WeMatch=0
IF (IPmap .gt. 0) THEN
WeMatch=1
ELSE
IPmap=ListMapped1_B(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap .gt. NP_RESloc)) THEN
WeMatch=1
END IF
END IF
IF (WeMatch .eq. 1) THEN
idx=idx+1
dspl_send(idx)=maxBlockLength*(IdxRev(IP)-1)
END IF
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,maxBlockLength,dspl_send,rtype,LocalColor % u2l_p2dsend_type(iNeigh), ierr)
call mpi_type_commit(LocalColor % u2l_p2dsend_type(iNeigh), ierr)
deallocate(dspl_send)
END DO
IdxRev=0
idx=0
DO IP=1,MNP
IF (ListNeedRecv(IP) .eq. 1) THEN
idx=idx+1
LocalColor % IdxRecv_u2l(idx)=IP
IdxRev(IP)=idx
END IF
END DO
DO iNeigh=1,u2l_nnbr_recv
iProc=LocalColor % u2l_ListNeigh_recv(iNeigh)+1
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMapped0_B=0
ListMapped1_B=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov == 0) THEN
ListMapped0_B(IP_glob)=IP
END IF
IF (eCov == 1) THEN
ListMapped1_B(IP_glob)=IP
END IF
END DO
nbCommon=LocalColor % u2l_ListNbCommon_recv(iNeigh)
allocate(dspl_recv(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 60')
idx=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov == 1) THEN
IPmap=ListMapped0(IP_glob)
WeMatch=0
IF (IPmap .gt. 0) THEN
WeMatch=1
ELSE
IPmap=ListMapped1(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap .gt. NP_RES)) THEN
WeMatch=1
END IF
END IF
IF (WeMatch .eq. 1) THEN
idx=idx+1
dspl_recv(idx)=maxBlockLength*(IdxRev(IPmap)-1)
END IF
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,maxBlockLength,dspl_recv,rtype,LocalColor % u2l_p2drecv_type(iNeigh),ierr)
call mpi_type_commit(LocalColor % u2l_p2drecv_type(iNeigh), ierr)
deallocate(dspl_recv)
END DO
deallocate(ListNeedRecv, ListNeedSend, IdxRev)
!
! Now the synchronization arrays
!
sync_nnbr_send=0
sync_nnbr_recv=0
ListCommon_send=0
ListCommon_recv=0
# ifdef DEBUG
WRITE(740+myrank,*) 'wwm_nnbr=', wwm_nnbr
WRITE(740+myrank,*) 'sum(ListCovLower)=', sum(LocalColor % ListCovLower)
# endif
DO iNeigh=1,wwm_nnbr
iProc=wwm_ListNeigh(iNeigh)
fColor=LocalColor % ListColor(iProc)
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMapped0_B=0
ListMapped1_B=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov == 0) THEN
ListMapped0_B(IP_glob)=IP
END IF
IF (eCov == 1) THEN
ListMapped1_B(IP_glob)=IP
END IF
END DO
nbCommon_recv=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov .eq. 1) THEN
IF (ListMapped0(IP_glob) .gt. 0) THEN
nbCommon_recv=nbCommon_recv+1
ELSE
IPmap=ListMapped1(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap .gt. NP_RES)) THEN
nbCommon_recv=nbCommon_recv+1
ENDIF
END IF
END IF
END DO
IF (nbCommon_recv .gt. 0) THEN
sync_nnbr_recv=sync_nnbr_recv+1
ListCommon_recv(iProc)=nbCommon_recv
END IF
nbCommon_send=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IPmap=ListMapped1(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap .le. NP_RES)) THEN
IF ((eCov .eq. 0).or.(IP.gt.NP_RESloc)) THEN
nbCommon_send=nbCommon_send+1
END IF
END IF
END DO
IF (nbCommon_send .gt. 0) THEN
sync_nnbr_send=sync_nnbr_send+1
ListCommon_send(iProc)=nbCommon_send
END IF
# ifdef DEBUG
WRITE(740+myrank,*) ' nbCommon(send/recv)=', nbCommon_send, nbCommon_recv
# endif
END DO
# ifdef DEBUG
WRITE(740+myrank,*) 'WWM_P2D: sync_nnbr_send=', sync_nnbr_send
WRITE(740+myrank,*) 'WWM_P2D: sync_nnbr_recv=', sync_nnbr_recv
# endif
LocalColor % sync_nnbr_send=sync_nnbr_send
LocalColor % sync_nnbr_recv=sync_nnbr_recv
allocate(LocalColor % sync_ListNbCommon_send(sync_nnbr_send), LocalColor % sync_ListNbCommon_recv(sync_nnbr_recv), LocalColor % sync_ListNeigh_send(sync_nnbr_send), LocalColor % sync_ListNeigh_recv(sync_nnbr_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 61')
idx_send=0
idx_recv=0
DO iProc=1,nproc
IF (ListCommon_send(iProc) .gt. 0) THEN
idx_send=idx_send+1
LocalColor % sync_ListNeigh_send(idx_send)=iProc-1
nbCommon=ListCommon_send(iProc)
LocalColor % sync_ListNbCommon_send(idx_send)=nbCommon
END IF
IF (ListCommon_recv(iProc) .gt. 0) THEN
idx_recv=idx_recv+1
LocalColor % sync_ListNeigh_recv(idx_recv)=iProc-1
nbCommon=ListCommon_recv(iProc)
LocalColor % sync_ListNbCommon_recv(idx_recv)=nbCommon
END IF
END DO
!
! Now creating the sync exchange
!
allocate(LocalColor % sync_p2dsend_rqst(sync_nnbr_send), &
&LocalColor % sync_p2drecv_rqst(sync_nnbr_recv), &
&LocalColor % sync_p2dsend_stat(MPI_STATUS_SIZE,sync_nnbr_send), &
&LocalColor % sync_p2drecv_stat(MPI_STATUS_SIZE,sync_nnbr_recv), &
&LocalColor % sync_p2dsend_type(sync_nnbr_send), &
&LocalColor % sync_p2drecv_type(sync_nnbr_recv), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 62')
# ifdef DEBUG
WRITE(740+myrank,*) 'SYNC sync_nnbr_send=', sync_nnbr_send
# endif
DO iNeigh=1,sync_nnbr_send
iProc=LocalColor % sync_ListNeigh_send(iNeigh)+1
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMapped0_B=0
ListMapped1_B=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov == 0) THEN
ListMapped0_B(IP_glob)=IP
END IF
IF (eCov == 1) THEN
ListMapped1_B(IP_glob)=IP
END IF
END DO
nbCommon=LocalColor % sync_ListNbCommon_send(iNeigh)
# ifdef DEBUG
WRITE(740+myrank,*) ' SYNC iNeigh=', iNeigh, ' nbCommon=', nbCommon
# endif
allocate(dspl_send(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 63')
idx=0
DO IP=1,NP_RES
IF (LocalColor % CovLower(IP) .eq. 1) THEN
DoOper=0
IP_glob=iplg(IP)
IPmap=ListMapped0_B(IP_glob)
IF (IPmap .gt. 0) THEN
DoOper=1
ELSE
IPmap=ListMapped1_B(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap.gt.NP_RESloc)) THEN
DoOper=1
END IF
END IF
IF (DoOper == 1) THEN
idx=idx+1
dspl_send(idx)=MSC*MDC*(IP-1)
END IF
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_send,rtype,LocalColor % sync_p2dsend_type(iNeigh), ierr)
call mpi_type_commit(LocalColor % sync_p2dsend_type(iNeigh), ierr)
deallocate(dspl_send)
END DO
DO iNeigh=1,sync_nnbr_recv
iProc=LocalColor % sync_ListNeigh_recv(iNeigh)+1
MNPloc=ListMNP(iProc)
NP_RESloc=ListNP_RES(iProc)
ListMapped0_B=0
ListMapped1_B=0
DO IP=1,MNPloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov == 0) THEN
ListMapped0_B(IP_glob)=IP
END IF
IF (eCov == 1) THEN
ListMapped1_B(IP_glob)=IP
END IF
END DO
nbCommon=LocalColor%sync_ListNbCommon_recv(iNeigh)
allocate(dspl_recv(nbCommon), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 64')
idx=0
DO IP=1,NP_RESloc
IP_glob=ListIPLG(IP+ListFirst(iProc))
eCov=LocalColor % ListCovLower(IP+ListFirst(iProc))
IF (eCov == 1) THEN
IPmap=ListMapped0(IP_glob)
DoOper=0
IF (IPmap .gt. 0) THEN
DoOper=1
ELSE
IPmap=ListMapped1(IP_glob)
IF ((IPmap .gt. 0).and.(IPmap.gt.NP_RES)) THEN
DoOper=1
END IF
END IF
IF (DoOper == 1) THEN
idx=idx+1
dspl_recv(idx)=MSC*MDC*(IPmap-1)
END IF
END IF
END DO
call mpi_type_create_indexed_block(nbCommon,MSC*MDC,dspl_recv,rtype,LocalColor % sync_p2drecv_type(iNeigh),ierr)
call mpi_type_commit(LocalColor % sync_p2drecv_type(iNeigh), ierr)
deallocate(dspl_recv)
END DO
lenMNP=0
IF (LocalColor % nbNeedSend_blk > lenMNP) THEN
lenMNP=LocalColor % nbNeedSend_blk
END IF
IF (LocalColor % nbNeedRecv_blk > lenMNP) THEN
lenMNP=LocalColor % nbNeedRecv_blk
END IF
IF (LocalColor % nbNeedSend_u2l > lenMNP) THEN
lenMNP=LocalColor % nbNeedSend_u2l
END IF
IF (LocalColor % nbNeedRecv_u2l > lenMNP) THEN
lenMNP=LocalColor % nbNeedRecv_u2l
END IF
allocate(LocalColor % ACexch(maxBlockLength, lenMNP), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 65')
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'FINISHED WITH INIT_COVLOWER_ARRAY'
FLUSH(STAT%FHNDL)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE DETERMINE_JSTATUS_L_U(LocalColor)
USE DATAPOOL, only : LocalColorInfo, MNP, rkind
USE DATAPOOL, only : NNZ, IA, JA, NP_RES, I_DIAG
USE DATAPOOL, only : wwm_nnbr_send, wwm_nnbr_recv
USE datapool, only : myrank, nproc, comm, ierr, nbrrank_p
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
integer Jstatus_L(NNZ), Jstatus_U(NNZ)
integer IP, J, JP, DoOper
integer istat
# ifdef REORDER_ASPAR_PC
integer nb, idx
# endif
Jstatus_L=0
Jstatus_U=0
DO IP=1,NP_RES
IF (LocalColor%CovLower(IP) == 1) THEN
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
IF (JP /= IP) THEN
IF (JP .lt. IP) THEN
DoOper=1
ELSE
IF (LocalColor % CovLower(JP) == 0) THEN
DoOper=1
ELSE
DoOper=0
END IF
END IF
ELSE
DoOper=0
END IF
Jstatus_L(J)=DoOper
END DO
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
IF (JP /= IP) THEN
IF (JP .lt. IP) THEN
DoOper=0
ELSE
IF (LocalColor % CovLower(JP) == 0) THEN
DoOper=0
ELSE
DoOper=1
END IF
END IF
ELSE
DoOper=0
END IF
Jstatus_U(J)=DoOper
END DO
END IF
END DO
# ifdef DEBUG
DO IP=1,NP_RES
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
IF ((Jstatus_L(J).eq.1).and.(Jstatus_U(J).eq.1)) THEN
WRITE(myrank+919,*) 'MNP=', MNP, 'NP_RES=', NP_RES
WRITE(myrank+919,*) 'IP=', IP, ' JP=', JP
WRITE(myrank+919,*) 'IPcovLower=', LocalColor%CovLower(IP)
WRITE(myrank+919,*) 'JPcovLower=', LocalColor%CovLower(JP)
WRITE(myrank+919,*) 'We have major error'
CALL WWM_ABORT('Please panic and debug')
END IF
END DO
END DO
WRITE(740+myrank,*) 'sum(Jstatus_L)=', sum(Jstatus_L)
WRITE(740+myrank,*) 'sum(Jstatus_U)=', sum(Jstatus_U)
# endif
# ifdef REORDER_ASPAR_PC
allocate(LocalColor % IA_L(NP_RES+1), LocalColor % IA_U(NP_RES+1), LocalColor % JA_LU(NNZ), LocalColor % JmapR(NNZ), stat=istat)
IF (istat/=0) CALL WWM_ABORT('determine_jstatus_L_U, error 1')
LocalColor % JmapR=-1
LocalColor % IA_L(1)=1
idx=0
DO IP=1,NP_RES
nb=0
DO J=IA(IP),IA(IP+1)-1
IF (Jstatus_L(J).eq.1) THEN
JP=JA(J)
nb=nb+1
idx=idx+1
LocalColor % JmapR(J)=idx
LocalColor % JA_LU(idx)=JP
END IF
END DO
LocalColor % IA_L(IP+1)=LocalColor % IA_L(IP)+nb
END DO
LocalColor % IA_U(1)=LocalColor % IA_L(NP_RES+1)
DO IP=1,NP_RES
nb=0
DO J=IA(IP),IA(IP+1)-1
IF (Jstatus_U(J).eq.1) THEN
JP=JA(J)
nb=nb+1
idx=idx+1
LocalColor % JmapR(J)=idx
LocalColor % JA_LU(idx)=JP
END IF
END DO
nb=nb+1
idx=idx+1
J=I_DIAG(IP)
LocalColor % JmapR(J)=idx
LocalColor % JA_LU(idx)=JP
LocalColor % IA_U(IP+1)=LocalColor % IA_U(IP)+nb
END DO
# endif
allocate(LocalColor % Jstatus_L(NNZ), LocalColor % Jstatus_U(NNZ), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 66')
LocalColor % Jstatus_L=Jstatus_L
LocalColor % Jstatus_U=Jstatus_U
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5_RECV_ASPAR_PC(LocalColor, SolDat)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData
USE DATAPOOL, only : MNP, MSC, MDC, rkind
USE datapool, only : ierr, comm, rtype, istatus, nbrrank_p
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
real(rkind), allocatable :: ASPAR_rs(:)
integer idx, iNNZ, jNNZ, IS, ID, NNZ_l, siz
integer iProc, i, iRank
integer istat
DO iProc=1,LocalColor % nbLow_send
i=LocalColor % ListIdxUpper_send(iProc)
iRank=nbrrank_p(i)
NNZ_l=LocalColor % NNZ_len_r(iProc)
siz=NNZ_l*MSC*MDC
allocate(ASPAR_rs(NNZ_l*MSC*MDC), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 67')
CALL mpi_recv(ASPAR_rs,siz,rtype,iRank,45,comm,istatus,ierr)
idx=0
DO iNNZ=1,NNZ_l
jNNZ=LocalColor % NNZ_index_r(iProc,iNNZ)
DO IS=1,MSC
DO ID=1,MDC
idx=idx+1
SolDat % ASPAR_pc(jNNZ,IS,ID)=ASPAR_rs(idx)
END DO
END DO
END DO
deallocate(ASPAR_rs)
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5_SEND_ASPAR_PC(LocalColor, SolDat)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData
USE DATAPOOL, only : MSC, MDC, rkind
USE datapool, only : ierr, comm, rtype, nbrrank_p
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
real(rkind), allocatable :: ASPAR_rs(:)
integer idx, iNNZ, jNNZ, IS, ID, NNZ_l, siz
integer iProc, iRank, i
integer istat
DO iProc=1,LocalColor % nbUpp_send
i=LocalColor % ListIdxUpper_send(iProc)
iRank=nbrrank_p(i)
NNZ_l=LocalColor % NNZ_len_s(iProc)
siz=NNZ_l*MSC*MDC
allocate(ASPAR_rs(NNZ_l*MSC*MDC), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 68')
idx=0
DO iNNZ=1,NNZ_l
jNNZ=LocalColor % NNZ_index_s(iProc,iNNZ)
DO IS=1,MSC
DO ID=1,MDC
idx=idx+1
ASPAR_rs(idx)=SolDat % ASPAR_pc(jNNZ,IS,ID)
END DO
END DO
END DO
CALL mpi_isend(ASPAR_rs,siz,rtype,iRank,45,comm,LocalColor%Upp_s_rq(iProc),ierr)
deallocate(ASPAR_rs)
END DO
IF (LocalColor % nbUpp_send > 0) THEN
call mpi_waitall(LocalColor %nbUpp_send, LocalColor % Upp_s_rq, LocalColor % Upp_s_stat,ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5_CREATE_PRECOND_ILU0(LocalColor, SolDat)
USE DATAPOOL, only : MNP, NP_RES, LocalColorInfo, I5_SolutionData
USE DATAPOOL, only : MSC, MDC, IA, JA, I_DIAG, rkind
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
integer IP, JP, J, JP2, J2, J_FOUND, IS, ID
integer :: ListJ(MNP)
real(rkind) tl
SolDat%ASPAR_pc=SolDat%ASPAR_block
CALL I5_RECV_ASPAR_PC(LocalColor, SolDat)
DO IP=1,NP_RES
IF (LocalColor % CovLower(IP) == 1) THEN
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
ListJ(JP)=J
END DO
DO J=IA(IP),I_DIAG(IP)-1
JP=JA(J)
DO IS=1,MSC
DO ID=1,MDC
tl=SolDat%ASPAR_pc(J,IS,ID)*SolDat%ASPAR_pc(I_DIAG(JP),IS,ID)
DO J2=IA(JP),IA(JP+1)-1
JP2=JA(J2)
J_FOUND=ListJ(JP2)
IF (J_FOUND.gt.0) THEN ! Here is ILU0 approximation
SolDat%ASPAR_pc(J_FOUND,IS,ID)=SolDat%ASPAR_pc(J_FOUND,IS,ID) - tl*SolDat%ASPAR_pc(J2,IS,ID)
END IF
END DO
SolDat%ASPAR_pc(J,IS,ID)=tl
END DO
END DO
END DO
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
ListJ(JP)=0
END DO
J=I_DIAG(IP)
DO IS=1,MSC
DO ID=1,MDC
SolDat%ASPAR_pc(J,IS,ID)=1.0_rkind/SolDat%ASPAR_pc(J,IS,ID)
END DO
END DO
END IF
END DO
CALL I5_SEND_ASPAR_PC(LocalColor, SolDat)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_CREATE_PRECOND_ILU0(LocalColor, SolDat)
USE DATAPOOL, only : MNP, NP_RES, LocalColorInfo, I5_SolutionData
USE DATAPOOL, only : MSC, MDC, IA, JA, I_DIAG, rkind
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
CALL WWM_ABORT('Please program it')
END SUBROUTINE
!**********************************************************************
!* We assign the values only for CovLower(IP)=1 *
!* We could with some effort assign values for all with some effort *
!* but the values would not be used *
!**********************************************************************
# ifndef SOR_DIRECT
SUBROUTINE I5B_CREATE_PRECOND_SOR(LocalColor, SolDat)
USE DATAPOOL, only : MNP, MDC, IA, JA, I_DIAG, NP_RES, rkind, ONE
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData
USE datapool, only : exchange_p4d_wwm
# ifdef DEBUG
USE datapool, only : myrank
# endif
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
integer IP, ID, IS, JP, JR, J1, J, IPglob, JPglob
real(rkind) eVal, Lerror
# if defined DEBUG
real(rkind) :: eSum, eSumB
# endif
DO IP=1,NP_RES
J=I_DIAG(IP)
SolDat%AC4(:,:,IP)=ONE/SolDat % ASPAR_block(:,:,J)
END DO
!# ifdef NO_SELFE_EXCH
! CALL I5B_EXCHANGE_P4D_WWM(LocalColor, SolDat%AC4)
!# else
! CALL I5B_TOTAL_COHERENCY_ERROR_NPRES(SolDat%AC4, Lerror)
! Print *, 'AC4_1 NP_RES cohenrency error=', Lerror
CALL EXCHANGE_P4D_WWM(SolDat%AC4)
!# endif
DO IP=1,NP_RES
IF (LocalColor%CovLower(IP) == 1) THEN
# if defined REORDER_ASPAR_PC
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor%Jstatus_L(J) == 1) THEN
JP=JA(J)
JR=LocalColor%JmapR(J)
SolDat % ASPAR_pc(:,:,JR)=SolDat % ASPAR_block(:,:,J)*SolDat%AC4(:,:,JP)
END IF
ENDDO
J=LocalColor% IA_U(IP+1)-1
SolDat % ASPAR_pc(:,:,J)=SolDat%AC4(:,:,IP)
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor%Jstatus_U(J) == 1) THEN
JP=JA(J)
JR=LocalColor%JmapR(J)
SolDat % ASPAR_pc(:,:,JR)=SolDat % ASPAR_block(:,:,J)
END IF
END DO
# else
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor%Jstatus_L(J) == 1) THEN
JP=JA(J)
SolDat % ASPAR_pc(:,:,J)=SolDat % ASPAR_block(:,:,J)*SolDat%AC4(:,:,JP)
END IF
ENDDO
J=I_DIAG(IP)
SolDat % ASPAR_pc(:,:,J)=SolDat%AC4(:,:,IP)
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor%Jstatus_U(J) == 1) THEN
JP=JA(J)
SolDat % ASPAR_pc(:,:,J)=SolDat % ASPAR_block(:,:,J)
END IF
END DO
# endif
END IF
END DO
END SUBROUTINE
# endif
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_CREATE_PRECOND(LocalColor, SolDat, TheMethod)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData, I_DIAG, NP_RES, rkind
USE datapool, only : exchange_p4d_wwm, myrank
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
integer, intent(in) :: TheMethod
integer IP, J
! real(rkind) :: Lerror
# ifdef SOR_DIRECT
IF (TheMethod.ne.1) THEN
CALL WWM_ABORT('With SOR_DIRECT, only SOR is possible')
END IF
DO IP=1,NP_RES
J=I_DIAG(IP)
SolDat%AC4(:,:,IP)=SolDat % ASPAR_block(:,:,J)
END DO
!# ifdef NO_SELFE_EXCH
! CALL I5B_EXCHANGE_P4D_WWM(LocalColor, SolDat%AC4)
!# else
! CALL I5B_TOTAL_COHERENCY_ERROR_NPRES(SolDat%AC4, Lerror)
! Print *, 'AC4_2 NP_RES cohenrency error=', Lerror
CALL EXCHANGE_P4D_WWM(SolDat%AC4)
!# endif
DO IP=1,NP_RES
J=I_DIAG(IP)
SolDat % ASPAR_block(:,:,J)=SolDat%AC4(:,:,IP)
END DO
!check this write ...
!WRITE(myrank+640,*) 'MIN ASPAR_block=', minval(SolDat % ASPAR_block)
# else
IF (TheMethod == 1) THEN ! SOR
CALL I5B_CREATE_PRECOND_SOR(LocalColor, SolDat)
ELSE IF (TheMethod == 2) THEN ! ILU0
CALL I5B_CREATE_PRECOND_ILU0(LocalColor, SolDat)
ELSE
CALL WWM_ABORT('Wrong choice of preconditioner')
ENDIF
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EXCHANGE_P3_LOW_2_UPP_Send(LocalColor, AC, iBlock)
USE DATAPOOL, only : LocalColorInfo, MNP, MSC, MDC, rkind
USE datapool, only : comm, ierr, myrank
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
REAL(rkind), intent(in) :: AC(MSC, MDC, MNP)
INTEGER, intent(in) :: iBlock
integer iUpp, iRank, idx, lenBlock, maxBlockLength, IS, ID, IP, nbUpp_send
integer idxIP
lenBlock=LocalColor % BlockLength(iBlock)
maxBlockLength=LocalColor % maxBlockLength
DO idxIP=1,LocalColor % nbNeedSend_blk
IP = LocalColor % IdxSend_blk(idxIP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
LocalColor % ACexch(idx,idxIP)=AC(IS,ID,IP)
END DO
END DO
nbUpp_send=LocalColor % nbUpp_send
DO iUpp=1,nbUpp_send
iRank=LocalColor % l2u_ListNeigh_send(iUpp)
CALL mpi_isend(LocalColor % ACexch, 1, LocalColor % l2u_p2dsend_type(iUpp), iRank, 7, comm, LocalColor%Upp_s_rq(iUpp), ierr)
END DO
IF (nbUpp_send > 0) THEN
call mpi_waitall(nbUpp_send, LocalColor % Upp_s_rq, LocalColor % Upp_s_stat,ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EXCHANGE_P3_LOW_2_UPP_Recv(LocalColor, AC, iBlock)
USE DATAPOOL, only : LocalColorInfo, MNP, MSC, MDC, rkind
USE datapool, only : comm, ierr, myrank
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
REAL(rkind), intent(inout) :: AC(MSC, MDC, MNP)
INTEGER, intent(in) :: iBlock
integer iProc, iRank, idx, lenBlock, IS, ID, IP, nbLow_recv
integer idxIP
lenBlock=LocalColor % BlockLength(iBlock)
nbLow_recv=LocalColor % nbLow_recv
DO iProc=1,nbLow_recv
iRank=LocalColor % l2u_ListNeigh_recv(iProc)
call mpi_irecv(LocalColor % ACexch,1,LocalColor % l2u_p2drecv_type(iProc),iRank,7,comm,LocalColor % Low_r_rq(iProc),ierr)
END DO
IF (nbLow_recv > 0) THEN
call mpi_waitall(nbLow_recv, LocalColor%Low_r_rq, LocalColor%Low_r_stat,ierr)
END IF
DO idxIP=1,LocalColor % nbNeedRecv_blk
IP = LocalColor % IdxRecv_blk(idxIP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
AC(IS,ID,IP) = LocalColor % ACexch(idx,idxIP)
END DO
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EXCHANGE_P3_UPP_2_LOW_Send(LocalColor, AC, iBlock)
USE DATAPOOL, only : LocalColorInfo, MNP, MSC, MDC, rkind
USE datapool, only : comm, ierr, myrank
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
REAL(rkind), intent(in) :: AC(MSC, MDC, MNP)
INTEGER, intent(in) :: iBlock
integer iProc, iRank, idx, lenBlock, IS, ID, nbLow_send, IP
integer idxIP
lenBlock=LocalColor % BlockLength(iBlock)
DO idxIP=1,LocalColor % nbNeedSend_u2l
IP = LocalColor % IdxSend_u2l(idxIP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
LocalColor % ACexch(idx,idxIP)=AC(IS,ID,IP)
END DO
END DO
nbLow_send=LocalColor % u2l_nnbr_send
DO iProc=1,nbLow_send
iRank=LocalColor % u2l_ListNeigh_send(iProc)
call mpi_isend(LocalColor % ACexch,1,LocalColor%u2l_p2dsend_type(iProc),iRank,1151,comm,LocalColor%u2l_p2dsend_rqst(iProc),ierr)
END DO
IF (nbLow_send > 0) THEN
call mpi_waitall(nbLow_send, LocalColor%u2l_p2dsend_rqst, LocalColor%u2l_p2dsend_stat,ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EXCHANGE_P3_UPP_2_LOW_Recv(LocalColor, AC, iBlock)
USE DATAPOOL, only : LocalColorInfo, MNP, MSC, MDC, rkind
USE datapool, only : comm, ierr, myrank
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
REAL(rkind), intent(inout) :: AC(MSC, MDC, MNP)
INTEGER, intent(in) :: iBlock
integer iProc, iRank, idx, lenBlock
integer nbUpp_recv, IS, ID, IP
integer idxIP
lenBlock=LocalColor % BlockLength(iBlock)
nbUpp_recv=LocalColor % u2l_nnbr_recv
DO iProc=1,nbUpp_recv
iRank=LocalColor % u2l_ListNeigh_recv(iProc)
call mpi_irecv(LocalColor % ACexch,1,LocalColor%u2l_p2drecv_type(iProc),iRank,1151,comm,LocalColor % u2l_p2drecv_rqst(iProc),ierr)
END DO
IF (nbUpp_recv > 0) THEN
call mpi_waitall(nbUpp_recv, LocalColor%u2l_p2drecv_rqst, LocalColor%u2l_p2drecv_stat,ierr)
END IF
DO idxIP=1,LocalColor % nbNeedRecv_u2l
IP = LocalColor % IdxRecv_u2l(idxIP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
AC(IS,ID,IP) = LocalColor % ACexch(idx,idxIP)
END DO
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
# ifdef DEBUG
SUBROUTINE COMPUTE_TOTAL_INDEX_SHIFT(TheRes)
USE DATAPOOL, only : NP_RES, IA, JA
implicit none
integer, intent(out) :: TheRes
integer :: IP, JP, J
TheRes=0
DO IP=1,NP_RES
DO J=IA(IP),IA(IP+1)-1
JP=JA(J)
TheRes=TheRes+abs(IP-JP)
END DO
END DO
END SUBROUTINE
# endif
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_PARTIAL_SOLVE_L(LocalColor, SolDat, iBlock, ACret)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData, I_DIAG, ONE
USE DATAPOOL, only : IA, JA, MSC, MDC, MNP, rkind, NP_RES, THR, THR8
USE datapool, only : myrank
# ifdef DEBUG
USE datapool, only : iplg
# endif
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(in) :: SolDat
integer, intent(in) :: iBlock
real(rkind), intent(inout) :: ACret(MSC,MDC,MNP)
real(rkind) :: eCoeff
! real(rkind) :: hVal
integer IP, idx, ID, IS, J
integer lenBlock, JP, Jb
! real(rkind) :: ACtest(MSC, MDC,MNP)
lenBlock=LocalColor % BlockLength(iBlock)
! DO IP=1,NP_RES
! J=I_DIAG(IP)
! ACtest(:,:,IP)=ONE/SolDat % ASPAR_block(:,:,J)
! END DO
DO IP=1,NP_RES
IF (LocalColor % CovLower(IP) == 1) THEN
# if defined REORDER_ASPAR_PC
DO J=LocalColor % IA_L(IP),LocalColor % IA_L(IP+1)-1
JP=LocalColor % JA_LU(J)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
eCoeff=SolDat % ASPAR_pc(IS,ID,J)
ACret(IS,ID,IP)=ACret(IS,ID,IP) - eCoeff*ACret(IS,ID,JP)
END DO
END DO
# elif defined SOR_DIRECT
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor % Jstatus_L(J) .eq. 1) THEN
JP=JA(J)
Jb=I_DIAG(JP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
eCoeff=SolDat % ASPAR_block(IS,ID,J)/SolDat % ASPAR_block(IS,ID,Jb)
ACret(IS,ID,IP)=ACret(IS,ID,IP) - eCoeff*ACret(IS,ID,JP)
END DO
END IF
END DO
# else
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor % Jstatus_L(J) .eq. 1) THEN
JP=JA(J)
Jb=I_DIAG(JP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
eCoeff=SolDat % ASPAR_pc(IS,ID,J)
! This is the mystery: the two formulas for eCoeffB gives
! different results
! hVal=ONE/SolDat % ASPAR_block(IS,ID,Jb)
! IF (abs(hVal - ACtest(IS,ID,JP)) .gt. THR8) THEN
! WRITE(740+myrank,*) 'hVal=', hVal, 'ACtest=', ACtest(IS,ID,JP)
! WRITE(740+myrank,*) ' diff=', hVal - ACtest(IS,ID,JP)
! FLUSH(740+myrank)
! END IF
! eCoeffB=SolDat % ASPAR_block(IS,ID,J)*hVal
! eCoeffB=SolDat % ASPAR_block(IS,ID,J)/SolDat % ASPAR_block(IS,ID,Jb)
! IF (abs(eCoeff - eCoeffB) .gt. THR) THEN
! WRITE(740+myrank,*) '1J=', J, 'eCoeff=', eCoeff, 'eCoeffB=', eCoeffB
! WRITE(740+myrank,*) ' diff=', eCoeff - eCoeffB
! FLUSH(740+myrank)
! END IF
ACret(IS,ID,IP)=ACret(IS,ID,IP) - eCoeff*ACret(IS,ID,JP)
END DO
END IF
END DO
# endif
ENDIF
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_PARTIAL_SOLVE_U(LocalColor, SolDat, iBlock, ACret)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData
USE DATAPOOL, only : MNP, IA, JA, I_DIAG, MSC, MDC, rkind, NP_RES, ONE, THR
USE datapool, only : myrank, iplg
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(in) :: SolDat
integer, intent(in) :: iBlock
real(rkind), intent(inout) :: ACret(MSC,MDC,MNP)
real(rkind) :: eCoeff
integer lenBlock, IP, JP, idx, J, IS, ID
lenBlock=LocalColor % BlockLength(iBlock)
DO IP=NP_RES,1,-1
IF (LocalColor % CovLower(IP) == 1) THEN
# if defined REORDER_ASPAR_PC
DO J=LocalColor % IA_U(IP),LocalColor % IA_U(IP+1)-2
JP=LocalColor % JA_LU(J)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
eCoeff=SolDat % ASPAR_pc(IS,ID,J)
ACret(IS,ID,IP)=ACret(IS,ID,IP) - eCoeff*ACret(IS,ID,JP)
END DO
END DO
J=LocalColor % IA_U(IP+1)-1
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
ACret(IS,ID,IP)=ACret(IS,ID,IP)*SolDat % ASPAR_pc(IS,ID,J)
END DO
# elif defined SOR_DIRECT
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor % Jstatus_U(J) .eq. 1) THEN
JP=JA(J)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
eCoeff=SolDat % ASPAR_block(IS,ID,J)
ACret(IS,ID,IP)=ACret(IS,ID,IP) - eCoeff*ACret(IS,ID,JP)
END DO
END IF
END DO
J=I_DIAG(IP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
ACret(IS,ID,IP)=ACret(IS,ID,IP)/SolDat % ASPAR_block(IS,ID,J)
END DO
# else
DO J=IA(IP),IA(IP+1)-1
IF (LocalColor % Jstatus_U(J) .eq. 1) THEN
JP=JA(J)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
eCoeff=SolDat % ASPAR_pc(IS,ID,J)
! eCoeffB=SolDat % ASPAR_block(IS,ID,J)
! IF (abs(eCoeff - eCoeffB) .gt. THR) THEN
! WRITE(740+myrank,*) '2J=', J, 'eCoeff=', eCoeff, 'eCoeffB=', eCoeffB
! WRITE(740+myrank,*) ' diff=', eCoeff - eCoeffB
! FLUSH(740+myrank)
! END IF
ACret(IS,ID,IP)=ACret(IS,ID,IP) - eCoeff*ACret(IS,ID,JP)
END DO
END IF
END DO
J=I_DIAG(IP)
DO idx=1,lenBlock
IS=LocalColor % ISindex(iBlock, idx)
ID=LocalColor % IDindex(iBlock, idx)
! eCoeff=SolDat % ASPAR_pc(IS,ID,J)
! eCoeffB=ONE/SolDat % ASPAR_block(IS,ID,J)
! IF (abs(eCoeff - eCoeffB) .gt. THR) THEN
! WRITE(740+myrank,*) '3J=', J, 'eCoeff=', eCoeff, 'eCoeffB=', eCoeffB
! WRITE(740+myrank,*) ' diff=', eCoeff - eCoeffB
! FLUSH(740+myrank)
! END IF
ACret(IS,ID,IP)=ACret(IS,ID,IP)*SolDat % ASPAR_pc(IS,ID,J)
END DO
# endif
ENDIF
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_SYNC_SENDRECV(LocalColor, AC)
USE DATAPOOL, only : LocalColorInfo, MNP, MSC, MDC, rkind
USE datapool, only : comm, ierr, myrank
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
real(rkind), intent(inout) :: AC(MSC, MDC, MNP)
INTEGER :: iRank
integer iSync
integer nbSync_send, nbSync_recv
nbSync_recv=LocalColor%sync_nnbr_recv
nbSync_send=LocalColor%sync_nnbr_send
DO iSync=1,nbSync_send
iRank=LocalColor % sync_ListNeigh_send(iSync)
CALL mpi_isend(AC, 1, LocalColor%sync_p2dsend_type(iSync), iRank, 1009, comm, LocalColor%sync_p2dsend_rqst(iSync), ierr)
END DO
DO iSync=1,nbSync_recv
iRank=LocalColor % sync_ListNeigh_recv(iSync)
call mpi_irecv(AC,1,LocalColor%sync_p2drecv_type(iSync),iRank,1009,comm,LocalColor%sync_p2drecv_rqst(iSync),ierr)
END DO
IF (nbSync_send > 0) THEN
call mpi_waitall(nbSync_send, LocalColor%sync_p2dsend_rqst, LocalColor%sync_p2dsend_stat,ierr)
END IF
IF (nbSync_recv > 0) THEN
call mpi_waitall(nbSync_recv, LocalColor%sync_p2drecv_rqst, LocalColor%sync_p2drecv_stat,ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_APPLY_PRECOND(LocalColor, SolDat, ACret)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData
USE DATAPOOL, only : MSC, MDC, MNP, rkind
USE datapool, only : myrank
USE DATAPOOL, only : DO_SYNC_UPP_2_LOW, DO_SYNC_LOW_2_UPP, DO_SYNC_FINAL
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
type(I5_SolutionData), intent(in) :: SolDat
real(rkind), intent(inout) :: ACret(MSC, MDC, MNP)
integer iBlock
DO iBlock=1,LocalColor % Nblock
IF (DO_SYNC_LOW_2_UPP) THEN
CALL I5B_EXCHANGE_P3_LOW_2_UPP_Recv(LocalColor, ACret, iBlock)
END IF
CALL I5B_PARTIAL_SOLVE_L(LocalColor, SolDat, iBlock, ACret)
IF (DO_SYNC_LOW_2_UPP) THEN
CALL I5B_EXCHANGE_P3_LOW_2_UPP_Send(LocalColor, ACret, iBlock)
END IF
END DO
DO iBlock=1,LocalColor%Nblock
IF (DO_SYNC_UPP_2_LOW) THEN
CALL I5B_EXCHANGE_P3_UPP_2_LOW_Recv(LocalColor, ACret, iBlock)
END IF
CALL I5B_PARTIAL_SOLVE_U(LocalColor, SolDat, iBlock, ACret)
IF (DO_SYNC_UPP_2_LOW) THEN
CALL I5B_EXCHANGE_P3_UPP_2_LOW_Send(LocalColor, ACret, iBlock)
END IF
END DO
IF (DO_SYNC_FINAL) THEN
CALL I5B_SYNC_SENDRECV(LocalColor, ACret)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_APPLY_FCT(LocalColor, SolDat, ACin, ACret)
USE DATAPOOL, only : I5_SolutionData, IA, JA, NP_RES, MDC, MSC, MNP, rkind
USE DATAPOOL, only : LocalColorInfo
USE datapool, only : exchange_p4d_wwm, myrank
implicit none
integer IP, J, idx
type(LocalColorInfo), intent(in) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
REAL(rkind), intent(in) :: ACin(MSC, MDC, MNP)
REAL(rkind), intent(inout) :: ACret(MSC, MDC, MNP)
REAL(rkind) :: eSum(MSC,MDC)
#ifdef DEBUG
REAL(rkind) :: Lerror
REAL(rkind) :: ACtest1(MSC, MDC, MNP)
REAL(rkind) :: ACtest2(MSC, MDC, MNP)
#endif
DO IP=1,NP_RES
eSum=0
DO J=IA(IP),IA(IP+1)-1
idx=JA(J)
eSum=eSum + SolDat % ASPAR_block(:,:,J)*ACin(:,:,idx)
END DO
ACret(:,:,IP)=eSum
END DO
#ifdef DEBUG
ACtest1=ACret
ACtest2=ACret
CALL I5B_EXCHANGE_P4D_WWM(LocalColor, ACtest1)
CALL EXCHANGE_P4D_WWM(ACtest2)
!
CALL I5B_TOTAL_COHERENCY_ERROR_NPRES(ACret, Lerror)
WRITE(740+myrank,*) 'ACret NP_RES cohenrency error=', Lerror
CALL I5B_TOTAL_COHERENCY_ERROR_NPRES(ACtest1, Lerror)
WRITE(740+myrank,*) 'ACtest1 NP_RES cohenrency error=', Lerror
CALL I5B_TOTAL_COHERENCY_ERROR_NPRES(ACtest2, Lerror)
WRITE(740+myrank,*) 'ACtest2 NP_RES cohenrency error=', Lerror
!
CALL I5B_TOTAL_COHERENCY_ERROR(ACret, Lerror)
WRITE(740+myrank,*) 'ACret MNP coherency error=', Lerror
CALL I5B_TOTAL_COHERENCY_ERROR(ACtest1, Lerror)
WRITE(740+myrank,*) 'ACtest1 MNP coherency error=', Lerror
CALL I5B_TOTAL_COHERENCY_ERROR(ACtest2, Lerror)
WRITE(740+myrank,*) 'ACtest2 MNP coherency error=', Lerror
!
FLUSH(740+myrank)
#endif
# ifdef NO_SELFE_EXCH
CALL I5B_EXCHANGE_P4D_WWM(LocalColor, ACret)
# else
CALL EXCHANGE_P4D_WWM(ACret)
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE REPLACE_NAN_ZERO(LocalColor, LScal)
USE DATAPOOL, only : rkind, MSC, MDC, LocalColorInfo
implicit none
type(LocalColorInfo), intent(in) :: LocalColor
real(rkind), intent(inout) :: LScal(MSC,MDC)
integer IS, ID
DO IS=1,MSC
DO ID=1,MDC
IF (LScal(IS,ID) .ne. LScal(IS,ID)) THEN
LScal(IS,ID)=0
END IF
END DO
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_L2_LINF(ACw1, ACw2, Norm_L2, Norm_LINF)
USE DATAPOOL, only : rkind, MNP, MDC, MSC
USE DATAPOOL, only : nwild_loc_res, NP_RES
USE datapool, only : myrank, comm, ierr, nproc, istatus, rtype
implicit none
real(rkind), intent(in) :: ACw1(MSC, MDC, MNP)
real(rkind), intent(in) :: ACw2(MSC, MDC, MNP)
real(rkind), intent(inout) :: Norm_L2(MSC, MDC)
real(rkind), intent(inout) :: Norm_LINF(MSC, MDC)
real(rkind) :: LScal(MSC, MDC, 2)
real(rkind) :: RScal(MSC, MDC, 2)
integer IP, iProc, IS, ID
LScal=0
DO IP=1,NP_RES
LScal(:,:,1)=LScal(:,:,1) + nwild_loc_res(IP)*((ACw1(:,:,IP) - ACw2(:,:,IP))**2)
DO IS=1,MSC
DO ID=1,MDC
LScal(IS,ID,2)=max(LScal(IS,ID,2), abs(ACw1(IS,ID,IP) - ACw2(IS,ID,IP)))
END DO
END DO
END DO
IF (myrank == 0) THEN
DO iProc=2,nproc
CALL MPI_RECV(RScal,MSC*MDC*2,rtype, iProc-1, 19, comm, istatus, ierr)
LScal(:,:,1) = LScal(:,:,1) + RScal(:,:,1)
DO IS=1,MSC
DO ID=1,MDC
LScal(IS,ID,2)=max(LScal(IS,ID,2), RScal(IS,ID,2))
END DO
END DO
END DO
DO iProc=2,nproc
CALL MPI_SEND(LScal,MSC*MDC*2,rtype, iProc-1, 23, comm, ierr)
END DO
ELSE
CALL MPI_SEND(LScal,MSC*MDC*2,rtype, 0, 19, comm, ierr)
CALL MPI_RECV(LScal,MSC*MDC*2,rtype, 0, 23, comm, istatus, ierr)
END IF
Norm_L2=LScal(:,:,1)
Norm_LINF=LScal(:,:,2)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_SCALAR(ACw1, ACw2, LScal)
USE DATAPOOL, only : rkind, MNP, MDC, MSC
USE DATAPOOL, only : nwild_loc_res, NP_RES
USE datapool, only : myrank, comm, ierr, nproc, istatus, rtype
implicit none
real(rkind), intent(in) :: ACw1(MSC, MDC, MNP)
real(rkind), intent(in) :: ACw2(MSC, MDC, MNP)
real(rkind), intent(inout) :: LScal(MSC, MDC)
real(rkind) :: RScal(MSC, MDC)
integer IP, iProc
LScal=0
DO IP=1,NP_RES
LScal=LScal + nwild_loc_res(IP)*ACw1(:,:,IP)*ACw2(:,:,IP)
END DO
IF (myrank == 0) THEN
DO iProc=2,nproc
CALL MPI_RECV(RScal,MSC*MDC,rtype, iProc-1, 19, comm, istatus, ierr)
LScal = LScal + RScal
END DO
DO iProc=2,nproc
CALL MPI_SEND(LScal,MSC*MDC,rtype, iProc-1, 23, comm, ierr)
END DO
ELSE
CALL MPI_SEND(LScal,MSC*MDC,rtype, 0, 19, comm, ierr)
CALL MPI_RECV(LScal,MSC*MDC,rtype, 0, 23, comm, istatus, ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_SUM_MAX(ACw, LSum, LMax)
USE DATAPOOL, only : rkind, MNP, MDC, MSC
USE DATAPOOL, only : nwild_loc_res, NP_RES
USE datapool, only : myrank, comm, ierr, nproc, istatus, rtype
implicit none
real(rkind), intent(in) :: ACw(MSC, MDC, MNP)
real(rkind), intent(inout) :: LSum(MSC, MDC)
real(rkind), intent(inout) :: LMax(MSC, MDC)
real(rkind) :: RScal(MSC, MDC)
integer IP, iProc, IS, ID
LSum=0
DO IP=1,NP_RES
LSum=LSum + nwild_loc_res(IP)*ACw(:,:, IP)
END DO
DO IS=1,MSC
DO ID=1,MDC
LMax(IS,ID)=maxval(ACw(IS,ID,:))
END DO
END DO
IF (myrank == 0) THEN
DO iProc=2,nproc
CALL MPI_RECV(RScal,MSC*MDC,rtype, iProc-1, 53, comm, istatus, ierr)
LSum = LSum + RScal
END DO
DO iProc=2,nproc
CALL MPI_RECV(RScal,MSC*MDC,rtype, iProc-1, 59, comm, istatus, ierr)
DO IS=1,MSC
DO ID=1,MDC
IF (RScal(IS,ID) .gt. LMax(IS,ID)) THEN
LMax(IS,ID)=RScal(IS,ID)
END IF
END DO
END DO
END DO
DO iProc=2,nproc
CALL MPI_SEND(LSum,MSC*MDC,rtype, iProc-1, 197, comm, ierr)
END DO
DO iProc=2,nproc
CALL MPI_SEND(LMax,MSC*MDC,rtype, iProc-1, 199, comm, ierr)
END DO
ELSE
CALL MPI_SEND(LSum,MSC*MDC,rtype, 0, 53, comm, ierr)
CALL MPI_SEND(LMax,MSC*MDC,rtype, 0, 59, comm, ierr)
CALL MPI_RECV(LSum,MSC*MDC,rtype, 0, 197, comm, istatus, ierr)
CALL MPI_RECV(LMax,MSC*MDC,rtype, 0, 199, comm, istatus, ierr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
FUNCTION I5B_SUMTOT(ACw)
USE DATAPOOL, only : rkind, MNP, MDC, MSC
implicit none
real(rkind), intent(in) :: ACw(MSC, MDC, MNP)
real(rkind) :: LSum(MSC, MDC)
real(rkind) :: LMax(MSC, MDC)
real(rkind) :: I5B_SUMTOT
CALL I5B_SUM_MAX(ACw, LSum, LMax)
I5B_SUMTOT=sum(LSum)
END FUNCTION
!**********************************************************************
!* *
!**********************************************************************
! We use the notations of
! http://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method
! and we use K1=Id
! In this algorithm, the use of v_{i-1}, v_i can be replace to just "v"
! The same for x, r
!
SUBROUTINE I5B_BCGS_REORG_SOLVER(LocalColor, SolDat, nbIter, Norm_L2, Norm_LINF)
USE DATAPOOL, only : MSC, MDC, MNP, NP_RES, NNZ, AC2, WAE_SOLVERTHR
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData, rkind
USE DATAPOOL, only : PCmethod, STAT, myrank
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
integer, intent(inout) :: nbIter
REAL(rkind), intent(inout) :: Norm_L2(MSC,MDC)
REAL(rkind), intent(inout) :: Norm_LINF(MSC,MDC)
REAL(rkind) :: Rho(MSC,MDC)
REAL(rkind) :: Prov(MSC,MDC)
REAL(rkind) :: Alpha(MSC,MDC)
REAL(rkind) :: Beta(MSC,MDC)
REAL(rkind) :: Omega(MSC,MDC)
REAL(rkind) :: MaxError, CritVal
integer :: MaxIter = 30
integer IP
MaxError=WAE_SOLVERTHR
CALL I5B_APPLY_FCT(LocalColor, SolDat, AC2, SolDat % AC3)
SolDat % AC1=0 ! y
SolDat % AC3=SolDat % B_block - SolDat % AC3 ! r residual
SolDat % AC4=SolDat % AC3 ! hat{r_0} term
SolDat % AC5=0 ! v
SolDat % AC6=0 ! p
SolDat % AC7=0 ! t
Rho=1
Alpha=1
Omega=1
nbIter=0
! WRITE(740+myrank,*) 'Beginning solution'
! FLUSH(740+myrank)
DO
nbIter=nbIter+1
! WRITE(740+myrank,*) 'nbIter=', nbIter
! FLUSH(740+myrank)
! L1: Rhoi =(\hat{r}_0, r_{i-1}
CALL I5B_SCALAR(SolDat % AC4, SolDat % AC3, Prov)
! L2: Beta=(RhoI/Rho(I-1)) * (Alpha/Omega(i-1))
Beta=(Prov/Rho)*(Alpha/Omega)
CALL REPLACE_NAN_ZERO(LocalColor, Beta)
Rho=Prov
! L3: Pi = r(i-1) + Beta*(p(i-1) -omega(i-1)*v(i-1))
DO IP=1,MNP
SolDat%AC6(:,:,IP)=SolDat%AC3(:,:,IP) &
& + Beta(:,:)*SolDat%AC6(:,:,IP) &
& - Beta(:,:)*Omega(:,:)*SolDat%AC5(:,:,IP)
END DO
! L4 y=K^(-1) Pi
SolDat%AC1=SolDat%AC6
IF (PCmethod .gt. 0) THEN
CALL I5B_APPLY_PRECOND(LocalColor, SolDat, SolDat%AC1)
ENDIF
! L5 vi=Ay
CALL I5B_APPLY_FCT(LocalColor, SolDat, SolDat%AC1, SolDat%AC5)
! L6 Alpha=Rho/(hat(r)_0, v_i)
CALL I5B_SCALAR(SolDat % AC4, SolDat % AC5, Prov)
Alpha(:,:)=Rho(:,:)/Prov(:,:)
CALL REPLACE_NAN_ZERO(LocalColor, Alpha)
! L6.1 x(i)=x(i-1) + Alpha y
DO IP=1,MNP
AC2(:,:,IP)=AC2(:,:,IP) &
& + Alpha(:,:)*SolDat%AC1(:,:,IP)
END DO
! L7 s=r(i-1) - alpha v(i)
DO IP=1,MNP
SolDat%AC3(:,:,IP)=SolDat%AC3(:,:,IP) &
& - Alpha(:,:)*SolDat%AC5(:,:,IP)
END DO
! L8 z=K^(-1) s
SolDat%AC1=SolDat%AC3
IF (PCmethod .gt. 0) THEN
CALL I5B_APPLY_PRECOND(LocalColor, SolDat, SolDat%AC1)
END IF
! L9 t=Az
CALL I5B_APPLY_FCT(LocalColor, SolDat, SolDat%AC1, SolDat%AC7)
! L10 omega=(t,s)/(t,t)
CALL I5B_SCALAR(SolDat % AC7, SolDat % AC3, Omega)
CALL I5B_SCALAR(SolDat % AC7, SolDat % AC7, Prov)
Omega(:,:)=Omega(:,:)/Prov(:,:)
CALL REPLACE_NAN_ZERO(LocalColor, Omega)
! L11 x(i)=x(i-1) + Omega z
DO IP=1,MNP
AC2(:,:,IP)=AC2(:,:,IP) &
& + Omega(:,:)*SolDat%AC1(:,:,IP)
END DO
! L12 If x is accurate enough finish
CALL I5B_APPLY_FCT(LocalColor, SolDat, AC2, SolDat%AC1)
CALL I5B_L2_LINF(SolDat%AC1, SolDat%B_block, Norm_L2, Norm_LINF)
CritVal=maxval(Norm_L2)
! WRITE(740+myrank,*) 'CritVal=', CritVal
! FLUSH(740+myrank)
IF (CritVal .lt. MaxError) THEN
EXIT
ENDIF
IF (nbIter .gt. MaxIter) THEN
EXIT
ENDIF
! L13 r=s-omega t
DO IP=1,MNP
SolDat%AC3(:,:,IP)=SolDat%AC3(:,:,IP) &
& - Omega(:,:)*SolDat%AC7(:,:,IP)
END DO
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_ALLOCATE(SolDat)
USE DATAPOOL, only : I5_SolutionData, MNP, MSC, MDC, NNZ
implicit none
type(I5_SolutionData), intent(inout) :: SolDat
integer istat
allocate(SolDat % AC1(MSC,MDC,MNP), SolDat % AC3(MSC,MDC,MNP), SolDat % AC4(MSC,MDC,MNP), SolDat % AC5(MSC,MDC,MNP), SolDat % AC6(MSC,MDC,MNP), SolDat % AC7(MSC,MDC,MNP), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 75')
allocate(SolDat % ASPAR_block(MSC,MDC,NNZ), SolDat % B_block(MSC,MDC, MNP), stat=istat)
# ifndef SOR_DIRECT
allocate(SolDat % ASPAR_pc(MSC,MDC,NNZ), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 77')
# endif
IF (istat/=0) CALL WWM_ABORT('wwm_parall_solver, allocate error 78')
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE WWM_SOLVER_INIT
implicit none
CALL I5B_SOLVER_INIT
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_SOLVER_INIT
USE DATAPOOL
implicit none
NblockFreqDir = NB_BLOCK
CALL SYMM_INIT_COLORING(MainLocalColor, NblockFreqDir)
# ifdef DEBUG
WRITE(myrank+740,*) 'After SYMM_INIT_COLORING'
FLUSH(myrank+740)
# endif
CALL I5B_ALLOCATE(SolDat)
# ifdef DEBUG
WRITE(myrank+740,*) 'After I5B_ALLOCATE'
FLUSH(myrank+740)
# endif
IF (PCmethod .eq. 2) THEN
! CALL CREATE_ASPAR_EXCHANGE_ARRAY(LocalColor)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_FREE(SolDat)
USE DATAPOOL, only : I5_SolutionData
implicit none
type(I5_SolutionData), intent(inout) :: SolDat
deallocate(SolDat % AC1, SolDat % AC3, SolDat % AC4, SolDat % AC5, SolDat % AC6, SolDat % AC7)
deallocate(SolDat % ASPAR_block, SolDat % B_block, SolDat % ASPAR_pc)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5_SUM(AC, eSum)
USE DATAPOOL, only : MNP, MSC, MDC, rkind, ZERO
implicit none
real(rkind), intent(in) :: AC(MNP,MSC,MDC)
real(rkind), intent(out) :: eSum
integer IP,IS,ID
eSum=ZERO
DO IP=1,MNP
DO IS=1,MSC
DO ID=1,MDC
eSum=eSum + AC(IP,IS,ID)
END DO
END DO
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE WWM_SOLVER_EIMPS(LocalColor, SolDat)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData, stat
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'ENTERING WWM_SOLVER_EIMPS'
FLUSH(STAT%FHNDL)
CALL I5B_EIMPS(LocalColor, SolDat)
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'FINISHING WWM_SOLVER_EIMPS'
FLUSH(STAT%FHNDL)
END SUBROUTINE
!**********************************************************************
!*
!**********************************************************************
SUBROUTINE EIMPS_B_BLOCK(U,B)
USE DATAPOOL
IMPLICIT NONE
REAL(rkind), intent(out) :: B(MSC, MDC, MNP)
REAL(rkind), intent(in) :: U(MSC, MDC, MNP)
INTEGER :: IP, ID, IS
INTEGER :: IPGL1, IPREL
REAL(rkind) :: TRIA03
!
# ifdef DEBUG
WRITE(740+myrank,*) 'Begin of EIMPS_B_BLOCK'
# endif
DO IP=1,MNP
DO ID=1,MDC
B(:,ID,IP) = U(:,ID,IP) * IOBPD(ID,IP)*IOBWB(IP)*IOBDP(IP)*SI(IP)
ENDDO
END DO
IF (LBCWA .OR. LBCSP) THEN
DO IP = 1, IWBMNP
IF (LINHOM) THEN
IPrel=IP
ELSE
IPrel=1
ENDIF
IPGL1 = IWBNDLC(IP)
B(:,:,IPGL1) = WBAC(:,:,IPrel) * SI(IPGL1) ! Overwrite ...
END DO
END IF
# if defined DEBUG
WRITE(3000+myrank,*) 'sum(B )=', sum(B)
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE EIMPS_ASPAR_B_BLOCK_SOURCES_TOTAL(U, ASPAR, B)
USE DATAPOOL
IMPLICIT NONE
REAL(rkind), intent(inout) :: ASPAR(MSC, MDC, NNZ)
REAL(rkind), intent(inout) :: B(MSC, MDC, MNP)
REAL(rkind), intent(in) :: U(MSC, MDC, MNP)
INTEGER :: POS_TRICK(3,2)
REAL(rkind) :: FL11(MSC,MDC), FL12(MSC,MDC), FL21(MSC,MDC), FL22(MSC,MDC), FL31(MSC,MDC), FL32(MSC,MDC)
REAL(rkind):: CRFS(MSC,MDC,3), K1(MSC,MDC), KM(MSC,MDC,3), K(MSC,MDC,3), TRIA03
REAL(rkind):: GTEMP2, DELFL, USFM, FLHAB, LIMFAC
# ifndef NO_MEMORY_CX_CY
REAL(rkind) :: CX(MSC,MDC,MNP), CY(MSC,MDC,MNP)
# else
REAL(rkind) :: CXY(2,MSC,MDC,3)
REAL(rkind) :: DIFRU, USOC, WVC
# endif
# ifndef SINGLE_LOOP_AMATRIX
REAL(rkind) :: DELTAL(MSC,MDC,3,MNE)
REAL(rkind) :: KP(MSC,MDC,3,MNE), NM(MSC,MDC,MNE)
INTEGER :: POS
# else
REAL(rkind) :: DELTAL(MSC,MDC,3)
REAL(rkind) :: KP(MSC,MDC,3), NM(MSC,MDC)
# endif
INTEGER :: I1, I2, I3
INTEGER :: IP, ID, IS, IE
INTEGER :: I, IPGL1, IPrel
REAL(rkind) :: DTK(MSC,MDC), TMP3(MSC,MDC)
REAL(rkind) :: LAMBDA(2,MSC,MDC)
# ifdef DEBUG
WRITE(740+myrank,*) 'Begin of EIMPS_ASPAR_B_BLOCK'
# endif
POS_TRICK(1,1) = 2
POS_TRICK(1,2) = 3
POS_TRICK(2,1) = 3
POS_TRICK(2,2) = 1
POS_TRICK(3,1) = 1
POS_TRICK(3,2) = 2
# ifndef NO_MEMORY_CX_CY
CALL CADVXY_VECTOR(CX, CY)
# endif
!
! Calculate countour integral quantities ...
!
# ifdef DEBUG
WRITE(740+myrank,*) ' Before MNE loop'
# endif
ASPAR = 0.0_rkind ! Mass matrix ...
B = 0.0_rkind ! Right hand side ...
DO IE = 1, MNE
# ifndef NO_MEMORY_CX_CY
I1 = INE(1,IE)
I2 = INE(2,IE)
I3 = INE(3,IE)
LAMBDA(1,:,:) = ONESIXTH * (CX(:,:,I1) + CX(:,:,I2) + CX(:,:,I3))
LAMBDA(2,:,:) = ONESIXTH * (CY(:,:,I1) + CY(:,:,I2) + CY(:,:,I3))
K(:,:,1) = LAMBDA(1,:,:) * IEN(1,IE) + LAMBDA(2,:,:) * IEN(2,IE)
K(:,:,2) = LAMBDA(1,:,:) * IEN(3,IE) + LAMBDA(2,:,:) * IEN(4,IE)
K(:,:,3) = LAMBDA(1,:,:) * IEN(5,IE) + LAMBDA(2,:,:) * IEN(6,IE)
FL11(:,:) = CX(:,:,I2)*IEN(1,IE)+CY(:,:,I2)*IEN(2,IE)
FL12(:,:) = CX(:,:,I3)*IEN(1,IE)+CY(:,:,I3)*IEN(2,IE)
FL21(:,:) = CX(:,:,I3)*IEN(3,IE)+CY(:,:,I3)*IEN(4,IE)
FL22(:,:) = CX(:,:,I1)*IEN(3,IE)+CY(:,:,I1)*IEN(4,IE)
FL31(:,:) = CX(:,:,I1)*IEN(5,IE)+CY(:,:,I1)*IEN(6,IE)
FL32(:,:) = CX(:,:,I2)*IEN(5,IE)+CY(:,:,I2)*IEN(6,IE)
# else
DO I=1,3
IP = INE(I,IE)
DO IS=1,MSC
DO ID=1,MDC
IF (LSECU .OR. LSTCU) THEN
CXY(1,IS,ID,I) = CG(IS,IP)*COSTH(ID)+CURTXY(IP,1)
CXY(2,IS,ID,I) = CG(IS,IP)*SINTH(ID)+CURTXY(IP,2)
ELSE
CXY(1,IS,ID,I) = CG(IS,IP)*COSTH(ID)
CXY(2,IS,ID,I) = CG(IS,IP)*SINTH(ID)
END IF
IF (LSPHE) THEN
CXY(1,IS,ID,I) = CXY(1,IS,ID,I)*INVSPHTRANS(IP,1)
CXY(2,IS,ID,I) = CXY(2,IS,ID,I)*INVSPHTRANS(IP,2)
END IF
IF (LDIFR) THEN
CXY(1,IS,ID,I) = CXY(1,IS,ID,I)*DIFRM(IP)
CXY(2,IS,ID,I) = CXY(2,IS,ID,I)*DIFRM(IP)
IF (LSECU .OR. LSTCU) THEN
IF (IDIFFR .GT. 1) THEN
WVC = SPSIG(IS)/WK(IS,IP)
USOC = (COSTH(ID)*CURTXY(IP,1) + SINTH(ID)*CURTXY(IP,2))/WVC
DIFRU = ONE + USOC * (ONE - DIFRM(IP))
ELSE
DIFRU = DIFRM(IP)
END IF
CXY(1,IS,ID,I) = CXY(1,IS,ID,I) + DIFRU*CURTXY(IP,1)
CXY(2,IS,ID,I) = CXY(2,IS,ID,I) + DIFRU*CURTXY(IP,2)
END IF
END IF
END DO
END DO
END DO
LAMBDA(:,:,:) = ONESIXTH * (CXY(:,:,:,1) + CXY(:,:,:,2) + CXY(:,:,:,3))
K(:,:,1) = LAMBDA(1,:,:) * IEN(1,IE) + LAMBDA(2,:,:) * IEN(2,IE)
K(:,:,2) = LAMBDA(1,:,:) * IEN(3,IE) + LAMBDA(2,:,:) * IEN(4,IE)
K(:,:,3) = LAMBDA(1,:,:) * IEN(5,IE) + LAMBDA(2,:,:) * IEN(6,IE)
FL11(:,:) = CXY(1,:,:,2)*IEN(1,IE)+CXY(2,:,:,2)*IEN(2,IE)
FL12(:,:) = CXY(1,:,:,3)*IEN(1,IE)+CXY(2,:,:,3)*IEN(2,IE)
FL21(:,:) = CXY(1,:,:,3)*IEN(3,IE)+CXY(2,:,:,3)*IEN(4,IE)
FL22(:,:) = CXY(1,:,:,1)*IEN(3,IE)+CXY(2,:,:,1)*IEN(4,IE)
FL31(:,:) = CXY(1,:,:,1)*IEN(5,IE)+CXY(2,:,:,1)*IEN(6,IE)
FL32(:,:) = CXY(1,:,:,2)*IEN(5,IE)+CXY(2,:,:,2)*IEN(6,IE)
# endif
CRFS(:,:,1) = - ONESIXTH * (TWO *FL31(:,:) + FL32(:,:) + FL21(:,:) + TWO * FL22(:,:) )
CRFS(:,:,2) = - ONESIXTH * (TWO *FL32(:,:) + TWO * FL11(:,:) + FL12(:,:) + FL31(:,:) )
CRFS(:,:,3) = - ONESIXTH * (TWO *FL12(:,:) + TWO * FL21(:,:) + FL22(:,:) + FL11(:,:) )
KM = MIN(ZERO,K)
# ifndef SINGLE_LOOP_AMATRIX
KP(:,:,:,IE) = MAX(ZERO,K)
DELTAL(:,:,:,IE) = CRFS(:,:,:) - KP(:,:,:,IE)
NM(:,:,IE)=ONE/MIN(-THR,KM(:,:,1) + KM(:,:,2) + KM(:,:,3))
# else
KP(:,:,:) = MAX(ZERO,K)
DELTAL(:,:,:) = CRFS(:,:,:)- KP(:,:,:)
NM(:,:)=ONE/MIN(-THR,KM(:,:,1) + KM(:,:,2) + KM(:,:,3))
TRIA03 = ONETHIRD * TRIA(IE)
DO I=1,3
IP=INE(I,IE)
I1=JA_IE(I,1,IE)
I2=JA_IE(I,2,IE)
I3=JA_IE(I,3,IE)
K1(:,:) = KP(:,:,I)
DO ID=1,MDC
DTK(:,ID) = K1(:,ID) * DT4A * IOBWB(IP) * IOBPD(ID,IP) * IOBDP(IP)
END DO
TMP3(:,:) = DTK(:,:) * NM(:,:)
ASPAR(:,:,I1) = TRIA03+DTK(:,:)- TMP3(:,:) * DELTAL(:,:,I ) + ASPAR(:,:,I1)
ASPAR(:,:,I2) = - TMP3(:,:) * DELTAL(:,:,POS_TRICK(I,1)) + ASPAR(:,:,I2)
ASPAR(:,:,I3) = - TMP3(:,:) * DELTAL(:,:,POS_TRICK(I,2)) + ASPAR(:,:,I3)
DO ID=1,MDC
B(:,ID,IP) = B(:,ID,IP) + U(:,ID,IP) * TRIA03 * IOBWB(IP) * IOBPD(ID,IP) * IOBDP(IP)
END DO
END DO
# endif
END DO
# ifndef SINGLE_LOOP_AMATRIX
J = 0 ! Counter ...
DO IP = 1, NP_RES
IF (IOBWB(IP) .EQ. 1 .AND. DEP(IP) .GT. DMIN) THEN
DO I = 1, CCON(IP)
J = J + 1
IE = IE_CELL(J)
POS = POS_CELL(J)
K1(:,:) = KP(:,:,POS,IE) ! Flux Jacobian
TRIA03 = ONETHIRD * TRIA(IE)
DO ID=1,MDC
DTK(:,ID) = K1(:,ID) * DT4A * IOBPD(ID,IP)
END DO
TMP3(:,:) = DTK(:,:) * NM(:,:,IE)
I1 = POSI(1,J) ! Position of the recent entry in the ASPAR matrix ... ASPAR is shown in fig. 42, p.122
I2 = POSI(2,J)
I3 = POSI(3,J)
ASPAR(:,:,I1) = TRIA03+DTK(:,:)- TMP3(:,:) * DELTAL(:,:,POS ,IE) + ASPAR(:,:,I1) ! Diagonal entry
ASPAR(:,:,I2) = - TMP3(:,:) * DELTAL(:,:,POS_TRICK(POS,1),IE) + ASPAR(:,:,I2) ! off diagonal entries ...
ASPAR(:,:,I3) = - TMP3(:,:) * DELTAL(:,:,POS_TRICK(POS,2),IE) + ASPAR(:,:,I3)
DO ID=1,MDC
B(:,ID,IP) = B(:,ID,IP) + IOBPD(ID,IP)*TRIA03 * U(:,ID,IP)
END DO
END DO
ELSE
DO I = 1, CCON(IP)
J = J + 1
IE = IE_CELL(J)
TRIA03 = ONETHIRD * TRIA(IE)
I1 = POSI(1,J) ! Position of the recent entry in the ASPAR matrix ... ASPAR is shown in fig. 42, p.122
ASPAR(:,:,I1) = TRIA03 + ASPAR(:,:,I1) ! Diagonal entry
B(:,:,IP) = ZERO
END DO
END IF
END DO
# endif
IF (LBCWA .OR. LBCSP) THEN
DO IP = 1, IWBMNP
IF (LINHOM) THEN
IPrel=IP
ELSE
IPrel=1
ENDIF
IPGL1 = IWBNDLC(IP)
ASPAR(:,:,I_DIAG(IPGL1)) = SI(IPGL1) ! Set boundary on the diagonal
B(:,:,IPGL1) = WBAC(:,:,IPrel) * SI(IPGL1)
END DO
END IF
IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0) THEN
DO IP = 1, NP_RES
ASPAR(:,:,I_DIAG(IP)) = ASPAR(:,:,I_DIAG(IP)) + IMATDAA(:,:,IP) * SI(IP) * DT4A !* IOBWB(IP) * IOBDP(IP) ! Add source term to the diagonal
B(:,:,IP) = B(:,:,IP) + IMATRAA(:,:,IP) * SI(IP) * DT4A !* IOBWB(IP) * IOBDP(IP) ! Add source term to the right hand side
END DO
ENDIF
! WRITE(*,*) SUM(IMATRAA), SUM(IMATDAA)
# if defined DEBUG
WRITE(3000+myrank,*) 'sum(ASPAR )=', sum(ASPAR)
WRITE(3000+myrank,*) 'sum(B )=', sum(B)
DO IS=1,MSC
WRITE(3000+myrank,*) 'IS, sum(ASPAR)=', IS, sum(ASPAR(IS,:,:))
END DO
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE I5B_EIMPS(LocalColor, SolDat)
USE DATAPOOL, only : LocalColorInfo, I5_SolutionData
USE DATAPOOL, only : rkind, MSC, MDC, AC2, MNP, NNZ
USE DATAPOOL, only : PCmethod, IOBPD, ZERO, STAT
USE datapool, only : myrank, exchange_p4d_wwm
implicit none
type(LocalColorInfo), intent(inout) :: LocalColor
type(I5_SolutionData), intent(inout) :: SolDat
integer nbIter
real(rkind) :: Norm_L2(MSC,MDC), Norm_LINF(MSC,MDC)
! real(rkind) :: Lerror
# ifndef ASPAR_B_COMPUTE_BLOCK
real(rkind) :: U(MNP), ASPAR(NNZ), B(MNP)
# endif
integer IS, ID, IP
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'ENTERING I5B_EIMPS'
FLUSH(STAT%FHNDL)
# ifdef DEBUG
WRITE(740+myrank,*) 'Begin I5B_EIMPS'
# endif
CALL EIMPS_ASPAR_B_BLOCK_SOURCES_TOTAL(AC2, SolDat%ASPAR_block, SolDat%B_block)
# ifdef DEBUG
WRITE(740+myrank,*) 'After ASPAR init'
# endif
!# ifdef NO_SELFE_EXCH
! CALL I5B_EXCHANGE_P4D_WWM(LocalColor, SolDat % B_block)
!# else
! CALL I5B_TOTAL_COHERENCY_ERROR_NPRES(SolDat%B_block, Lerror)
! Print *, 'B_block NP_RES cohenrency error=', Lerror
CALL EXCHANGE_P4D_WWM(SolDat % B_block)
!# endif
# ifdef DEBUG
WRITE(740+myrank,*) 'After EXCHANGE_P4D_WWM'
# endif
CALL I5B_EXCHANGE_ASPAR(LocalColor, SolDat%ASPAR_block)
# ifdef DEBUG
WRITE(740+myrank,*) 'After I5B_EXCHANGE_ASPAR'
# endif
CALL I5B_CREATE_PRECOND(LocalColor, SolDat, PCmethod)
# ifdef DEBUG
WRITE(740+myrank,*) 'After I5B_CREATE_PRECOND'
# endif
CALL I5B_BCGS_REORG_SOLVER(LocalColor, SolDat, nbIter, Norm_L2, Norm_LINF)
DO IP=1,MNP
DO ID=1,MDC
AC2(:,ID,IP)=MAX(ZERO, AC2(:,ID,IP))*MyREAL(IOBPD(ID,IP))
END DO
END DO
WRITE(STAT%FHNDL,*) 'nbIter=', nbIter, 'L2/LINF=', maxval(Norm_L2), maxval(Norm_LINF)
# ifdef DEBUG
IF (myrank == 1) THEN
Write(myrank+591,*) 'Clearing ENDING'
END IF
# endif
WRITE(STAT%FHNDL,'("+TRACE......",A)') 'FINISHING I5B_EIMPS'
FLUSH(STAT%FHNDL)
END SUBROUTINE
#endif