Skip to content
Snippets Groups Projects
Commit 48a0df89 authored by Joerg Behrens's avatar Joerg Behrens
Browse files

Fortran interface for yaxt (incomplete) + test

parent 2a11b8f5
No related branches found
No related tags found
No related merge requests found
......@@ -62,5 +62,6 @@ AC_CONFIG_FILES([tests/test_redist_p2p_parallel_run \
tests/test_ut_run \
tests/test_perf_run \
tests/test_sort_run \
tests/test_yaxt_run \
util/serialrun],[chmod a+x "$ac_file"])
AC_OUTPUT([Makefile src/Makefile tests/Makefile])
......@@ -45,5 +45,7 @@ libyaxt_a_SOURCES = \
sort_common.h \
xt_sort.c \
xt_sort.h \
yaxt_f2c.c \
yaxt.f90 \
yaxt.h
MODULE yaxt
!
! Fortran interface to yaxt implementation
!
USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_char, c_null_char, c_int, c_ptr, c_null_ptr, c_loc
IMPLICIT NONE
PRIVATE
PUBLIC :: xt_abort, xt_idxlist, xt_idxlist_delete, xt_idxvec_new, &
& xt_xmap, xt_xmap_new, xt_xmap_delete, &
& xt_redist, xt_redist_p2p_off_new, xt_redist_p2p_new, xt_redist_delete, &
& xt_redist_s_exchange1, xt_idxstripes_new, c_loc
INTEGER, PARAMETER, PUBLIC :: xt_count_kind = C_INT ! \todo this should come out of configure
TYPE, BIND(C), PUBLIC :: xt_stripe
INTEGER(C_INT) :: start
INTEGER(xt_count_kind) :: nstrides
INTEGER(C_INT) :: stride
END TYPE xt_stripe
TYPE, BIND(C) :: xt_idxlist
PRIVATE
TYPE(c_ptr) :: cptr = c_null_ptr
END TYPE xt_idxlist
TYPE, BIND(C) :: xt_xmap
PRIVATE
TYPE(c_ptr) :: cptr = c_null_ptr
END TYPE xt_xmap
TYPE, BIND(C) :: xt_redist
PRIVATE
TYPE(c_ptr) :: cptr = c_null_ptr
END TYPE xt_redist
INTERFACE
SUBROUTINE xt_abort_cmsl_f(comm_f, msg, source, line) BIND(C, name='xt_abort_cmsl_f')
IMPORT:: c_char, c_int
INTEGER, INTENT(in):: comm_f
CHARACTER(C_CHAR), DIMENSION(*), INTENT(in) :: msg
CHARACTER(C_CHAR), DIMENSION(*), INTENT(in) :: source
INTEGER(C_INT), VALUE, INTENT(in) :: line
END SUBROUTINE xt_abort_cmsl_f
END INTERFACE
INTERFACE
SUBROUTINE xt_abort_msl_f(msg, source, line) BIND(C, name='xt_abort_msl_f')
IMPORT:: c_char, c_int
CHARACTER(C_CHAR), DIMENSION(*), INTENT(in) :: msg
CHARACTER(C_CHAR), DIMENSION(*), INTENT(in) :: source
INTEGER(C_INT), VALUE, INTENT(in) :: line
END SUBROUTINE xt_abort_msl_f
END INTERFACE
INTERFACE xt_abort
MODULE PROCEDURE xt_abort_cmsl
MODULE PROCEDURE xt_abort_msl
END INTERFACE xt_abort
INTERFACE
SUBROUTINE xt_ut_destroy_transposition(handle) BIND(C, name='xt_ut_destroy_transposition')
IMPORT:: C_INT
IMPLICIT NONE
INTEGER(C_INT), INTENT(in) :: handle
END SUBROUTINE xt_ut_destroy_transposition
END INTERFACE
INTERFACE
SUBROUTINE xt_initialize(default_comm) BIND(C, name='xt_initialize_f')
IMPORT:: C_INT
IMPLICIT NONE
INTEGER, INTENT(in) :: default_comm
END SUBROUTINE xt_initialize
END INTERFACE
INTERFACE
FUNCTION xt_idxvec_new(idxvec, num_indices) BIND(C, name='xt_idxvec_new_f') RESULT(res)
IMPORT :: xt_idxlist
IMPLICIT NONE
INTEGER, INTENT(in) :: idxvec(*)
INTEGER, VALUE, INTENT(in) :: num_indices
TYPE(xt_idxlist) :: res
END FUNCTION xt_idxvec_new
END INTERFACE
INTERFACE
SUBROUTINE xt_idxlist_delete(idxlist) BIND(C, name='xt_idxlist_delete_f')
IMPORT :: xt_idxlist
IMPLICIT NONE
TYPE(xt_idxlist), INTENT(inout) :: idxlist
END SUBROUTINE xt_idxlist_delete
END INTERFACE
INTERFACE
FUNCTION xt_xmap_all2all_new(src_idxlist, dst_idxlist, comm) BIND(C, name='xt_xmap_all2all_new_f') RESULT(res)
IMPORT :: xt_idxlist, xt_xmap
IMPLICIT NONE
TYPE(xt_idxlist), INTENT(in) :: src_idxlist
TYPE(xt_idxlist), INTENT(in) :: dst_idxlist
INTEGER, INTENT(in) :: comm
TYPE(xt_xmap) :: res
END FUNCTION xt_xmap_all2all_new
END INTERFACE
INTERFACE
SUBROUTINE xt_xmap_delete(xmap) BIND(C, name='xt_xmap_delete_f')
IMPORT :: xt_xmap
IMPLICIT NONE
TYPE(xt_xmap), INTENT(inout) :: xmap
END SUBROUTINE xt_xmap_delete
END INTERFACE
INTERFACE
FUNCTION xt_redist_p2p_off_new(xmap, src_offsets, dst_offsets, comm) BIND(C, name='xt_redist_p2p_off_new_f') RESULT(res)
IMPORT :: xt_xmap, xt_redist, c_int
IMPLICIT NONE
TYPE(xt_xmap), INTENT(in) :: xmap
INTEGER(c_int), INTENT(in) :: src_offsets(*)
INTEGER(c_int), INTENT(in) :: dst_offsets(*)
INTEGER, INTENT(in) :: comm
TYPE(xt_redist) :: res
END FUNCTION xt_redist_p2p_off_new
END INTERFACE
INTERFACE
FUNCTION xt_redist_p2p_new(xmap, comm) BIND(C, name='xt_redist_p2p_new_f') RESULT(res)
IMPORT:: xt_xmap, xt_redist, c_int
IMPLICIT NONE
TYPE(xt_xmap), INTENT(in) :: xmap
INTEGER, INTENT(in) :: comm
TYPE(xt_redist) :: res
END FUNCTION xt_redist_p2p_new
END INTERFACE
INTERFACE
SUBROUTINE xt_redist_delete(redist) BIND(C, name='xt_redist_delete_f')
IMPORT :: xt_redist
IMPLICIT NONE
TYPE(xt_redist), INTENT(inout) :: redist
END SUBROUTINE xt_redist_delete
END INTERFACE
INTERFACE
SUBROUTINE xt_redist_s_exchange1(redist, src_data_cptr, dst_data_cptr) BIND(C, name='xt_redist_s_exchange1_f')
IMPORT:: xt_redist, c_ptr
TYPE(xt_redist), INTENT(in) :: redist
TYPE(c_ptr), INTENT(in) :: src_data_cptr
TYPE(c_ptr), INTENT(in) :: dst_data_cptr
END SUBROUTINE xt_redist_s_exchange1
END INTERFACE
INTERFACE
FUNCTION xt_idxstripes_new(stripes, num_stripes) BIND(C, name='xt_idxstripes_new_f') RESULT(res)
IMPORT:: xt_idxlist, xt_stripe, xt_count_kind
IMPLICIT NONE
TYPE(xt_stripe), INTENT(in) :: stripes(*)
INTEGER(xt_count_kind), VALUE, INTENT(in) :: num_stripes
TYPE(xt_idxlist) :: res
END FUNCTION xt_idxstripes_new
END INTERFACE
CONTAINS
SUBROUTINE xt_abort_cmsl(comm_f, msg, source, line)
INTEGER, INTENT(in):: comm_f
CHARACTER(len=*), INTENT(in) :: msg
CHARACTER(len=*), INTENT(in) :: source
INTEGER, INTENT(in) :: line
CALL xt_abort_cmsl_f(comm_f, TRIM(msg)//c_null_char, source, line)
END SUBROUTINE Xt_abort_cmsl
SUBROUTINE xt_abort_msl(msg, source, line)
CHARACTER(len=*), INTENT(in) :: msg
CHARACTER(len=*), INTENT(in) :: source
INTEGER, INTENT(in) :: line
CALL xt_abort_msl_f(TRIM(msg)//c_null_char, source, line)
END SUBROUTINE Xt_abort_msl
END MODULE yaxt
#include "core/core.h"
#include "xt_idxlist.h"
#include "xt_idxvec.h"
#include "xt_xmap.h"
#include "xt_xmap_all2all.h"
#include "xt_idxstripes.h"
#include "xt_redist_p2p.h"
struct xt_idxlist_f {
struct xt_idxlist *cptr;
};
struct xt_xmap_f {
struct Xt_xmap *cptr;
};
struct xt_redist_f {
struct xt_redist *cptr;
};
void xt_initialize_f(MPI_Fint *comm_f) {
MPI_Comm comm_c;
comm_c = MPI_Comm_f2c(*comm_f);
Xt_initialize(comm_c);
}
void xt_abort_cmsl_f(MPI_Fint *comm_f, const char *msg, const char *source, int line) {
MPI_Comm comm_c = MPI_Comm_f2c(*comm_f);
Xt_abort(comm_c, msg, source, line);
}
void xt_abort_msl_f(const char *msg, const char *source, int line) {
MPI_Comm comm_c = SymPrefix(default_comm);
Xt_abort(comm_c, msg, source, line);
}
struct xt_idxlist_f xt_idxvec_new_f(Xt_idx *idxvec, Xt_count num_indices) {
struct xt_idxlist_f idxlist_f;
idxlist_f.cptr = xt_idxvec_new(idxvec, num_indices);
return idxlist_f;
}
void xt_idxlist_delete_f(struct xt_idxlist_f *idxlist_f) {
xt_idxlist_delete(idxlist_f->cptr);
}
struct xt_xmap_f xt_xmap_all2all_new_f(struct xt_idxlist_f *src_idxlist_f, struct xt_idxlist_f *dst_idxlist_f, MPI_Fint *comm_f) {
MPI_Comm comm_c = MPI_Comm_f2c(*comm_f);
struct xt_xmap_f xmap_f;
xmap_f.cptr = xt_xmap_all2all_new(src_idxlist_f->cptr, dst_idxlist_f->cptr, comm_c);
return xmap_f;
}
void xt_xmap_delete_f(struct xt_xmap_f *xmap_f) {
xt_xmap_delete(xmap_f->cptr);
}
struct xt_redist_f xt_redist_p2p_off_new_f(struct xt_xmap_f *xmap_f, int *src_offsets, int *dst_offsets, MPI_Fint *datatype_f) {
MPI_Datatype datatype_c = MPI_Type_f2c(*datatype_f);
struct xt_redist_f redist_f;
redist_f.cptr = xt_redist_p2p_off_new(xmap_f->cptr, src_offsets, dst_offsets, datatype_c);
return redist_f;
}
struct xt_redist_f xt_redist_p2p_new_f(struct xt_xmap_f *xmap_f, MPI_Fint *datatype_f) {
MPI_Datatype datatype_c = MPI_Type_f2c(*datatype_f);
struct xt_redist_f redist_f;
redist_f.cptr = xt_redist_p2p_new(xmap_f->cptr, datatype_c);
return redist_f;
}
void xt_redist_delete_f(struct xt_redist_f *redist_f) {
xt_redist_delete(redist_f->cptr);
}
void xt_redist_s_exchange1_f(struct xt_redist_f *redist_f, void **src_data, void **dst_data) {
xt_redist_s_exchange(redist_f->cptr, src_data, 1, dst_data, 1);
}
struct xt_idxlist_f xt_idxstripes_new_f(struct Xt_stripe const *stripes, Xt_count num_stripes) {
struct xt_idxlist_f idxlist_f;
idxlist_f.cptr = xt_idxstripes_new(stripes, num_stripes);
return idxlist_f;
}
......@@ -10,7 +10,8 @@ noinst_PROGRAMS = \
test_xmap_parallel \
test_ut \
test_perf \
test_sort
test_sort \
test_yaxt
test_idxvec_SOURCES = test_idxvec.c tests.h test_idxlist.h
......@@ -26,6 +27,7 @@ test_ut_SOURCES = test_ut.f90
test_ut_LDADD = $(LDADD) $(MPI_FC_LIB)
test_perf_SOURCES = test_perf.f90
test_sort_SOURCES = test_sort.c tests.h
test_yaxt_SOURCES = test_yaxt.f90
AM_CFLAGS = -I$(top_srcdir)/src
AM_FCFLAGS = -I$(top_srcdir)/src -I../src $(MPI_FC_INCLUDE)
......@@ -43,6 +45,8 @@ TESTS = \
test_xmap_parallel_run \
test_ut_run \
test_perf_run \
test_sort_run
test_sort_run \
test_yaxt_run
AUTOMAKE_OPTIONS = color-tests
AUTOMAKE_OPTIONS = color-tests
\ No newline at end of file
PROGRAM test_yaxt
USE mpi
USE yaxt, ONLY: xt_abort, xt_idxlist, xt_idxlist_delete, xt_idxvec_new, &
& xt_xmap, xt_xmap_all2all_new, xt_xmap_delete, &
& xt_redist, xt_redist_p2p_new, xt_redist_delete, xt_redist_p2p_off_new, xt_redist_s_exchange1, &
& c_loc
IMPLICIT NONE
INTEGER, PARAMETER :: g_ie = 8, g_je = 4! global extents including halos
LOGICAL, PARAMETER :: verbose = .FALSE.
INTEGER, PARAMETER :: nlev = 3
INTEGER, PARAMETER :: undef_int = HUGE(undef_int)/2 - 1
INTEGER, PARAMETER :: undef_index = -1
INTEGER, PARAMETER :: nhalo = 1 ! 1dim. halo border size
INTEGER :: ie, je ! local extents, including halos
INTEGER :: p_ioff, p_joff ! offsets within global domain
INTEGER :: nprocx, nprocy ! process space extents
INTEGER :: nprocs ! == nprocx*nprocy
INTEGER :: mype, mypx, mypy ! process rank, process coords within (0:, 0:) process space
LOGICAL :: lroot ! true only for proc 0
INTEGER :: g_id(g_ie, g_je) ! global id
INTEGER :: g_tpex(g_ie, g_je) ! global "tripolar-like" toy bounds exchange
TYPE(xt_xmap) :: xmap_tpex
TYPE(xt_redist) :: redist_tpex
TYPE(xt_redist) :: redist_surf_tpex
INTEGER, ALLOCATABLE :: loc_id(:,:), loc_tpex(:,:)
INTEGER, ALLOCATABLE :: fval(:,:), gval(:,:)
INTEGER, ALLOCATABLE :: gval3d(:,:,:)
INTEGER, ALLOCATABLE :: id_pos(:,:), pos3d_surf(:,:)
! mpi & decomposition & allocate mem:
CALL init_all
! full global index space:
CALL id_map(g_id)
! local window of global index space:
CALL get_window(g_id, loc_id)
! define bounds exchange for full global index space
CALL def_exchange(g_id, g_tpex)
! local window of global bounds exchange:
CALL get_window(g_tpex, loc_tpex)
! template: loc_id -> loc_tpex
CALL gen_template(loc_id, loc_tpex, xmap_tpex) ! todo rename template to xmap
! transposition: loc_id:data -> loc_tpex:data
CALL gen_trans(xmap_tpex, MPI_INTEGER, MPI_INTEGER, redist_tpex)
! test 2d-to-2d transposition:
fval = loc_id
CALL exchange_int(redist_tpex, fval, gval)
CALL icmp('2d to 2d check', gval, loc_tpex)
! define positions of surface elements within (i,k,j) array
CALL id_map(id_pos)
CALL id_map(pos3d_surf)
CALL gen_pos3d_surf(pos3d_surf)
! generate surface transposition:
CALL gen_off_trans(xmap_tpex, MPI_INTEGER, id_pos(:,:)-1, MPI_INTEGER, pos3d_surf(:,:)-1, redist_surf_tpex)
! 2d to surface boundsexchange:
gval3d = -1
CALL exchange_int(redist_surf_tpex, fval, gval3d)
CALL icmp('surface check', gval3d(:,1,:), loc_tpex)
! check sub surface:
CALL icmp('sub surface check', gval3d(:,2,:), loc_tpex*0-1)
! cleanup:
CALL xt_xmap_delete(xmap_tpex)
CALL xt_redist_delete(redist_tpex)
CALL xt_redist_delete(redist_surf_tpex)
CALL exit_all
CONTAINS
SUBROUTINE gen_pos3d_surf(pos)
INTEGER, INTENT(inout) :: pos(:,:)
! positions for zero based arrays (ECHAM grid point dim order)
! old pos = i + j*ie
! new pos = i + k*ie + j*ie*nlev
INTEGER :: ii,jj, i,j,k, p,q
k = 0 ! surface
DO jj=1,je
DO ii=1,ie
p = pos(ii,jj) - 1 ! shift to 0-based index
j = p/ie
i = MOD(p,ie)
q = i + k*ie + j*ie*nlev
pos(ii,jj) = q + 1 ! shift to 1-based index
ENDDO
ENDDO
END SUBROUTINE gen_pos3d_surf
SUBROUTINE icmp(label, f,g)
CHARACTER(len=*), PARAMETER :: context = 'test_ut::icmp: '
CHARACTER(len=*), INTENT(in) :: label
INTEGER, INTENT(in) :: f(:,:)
INTEGER, INTENT(in) :: g(:,:)
INTEGER :: i, j, n1, n2
n1 = SIZE(f,1)
n2 = SIZE(f,2)
IF (SIZE(g,1) /= n1 .OR. SIZE(g,2) /= n2) CALL xt_abort(context//'internal error', __FILE__, __LINE__)
DO j = 1, n2
DO i = 1, n1
IF (f(i,j) /= g(i,j)) THEN
WRITE(0,*) context,'test failed: i, j, f(i,j), g(i,j) =',i, j, f(i,j), g(i,j)
CALL xt_abort(context//'test failed', __FILE__, __LINE__)
ENDIF
ENDDO
ENDDO
IF (verbose) WRITE(0,*) mype,':',context//label//' passed'
END SUBROUTINE icmp
SUBROUTINE init_all
CHARACTER(len=*), PARAMETER :: context = 'init_all: '
INTEGER :: ierror
CALL MPI_INIT(ierror)
IF (ierror /= MPI_SUCCESS) CALL xt_abort(context//'MPI_INIT failed', __FILE__, __LINE__)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierror)
IF (ierror /= MPI_SUCCESS) CALL xt_abort(context//'MPI_COMM_SIZE failed', __FILE__, __LINE__)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mype, ierror)
IF (ierror /= MPI_SUCCESS) CALL xt_abort(context//'MPI_COMM_RANK failed', __FILE__, __LINE__)
IF (mype==0) THEN
lroot = .true.
ELSE
lroot = .FALSE.
ENDIF
CALL factorize(nprocs, nprocx, nprocy)
IF (verbose .AND. lroot) WRITE(0,*) 'nprocx, nprocy=',nprocx, nprocy
mypy = mype / nprocx
mypx = MOD(mype, nprocx)
!CALL ut_init(decomp_size=30, comm_tmpl_size=30, comm_size=30, &
! & debug_lvl=0, mode=ut_mode_dt_p2p, debug_unit=0)
CALL deco
ALLOCATE(fval(ie,je), gval(ie,je))
ALLOCATE(loc_id(ie,je), loc_tpex(ie,je))
ALLOCATE(id_pos(ie,je), gval3d(ie,nlev,je), pos3d_surf(ie,je))
fval = undef_int
gval = undef_int
loc_id = undef_int
loc_tpex = undef_int
id_pos = undef_int
gval3d = undef_int
pos3d_surf = undef_int
END SUBROUTINE init_all
SUBROUTINE exit_all
CHARACTER(len=*), PARAMETER :: context = 'exit_all: '
INTEGER :: ierror
CALL MPI_FINALIZE(ierror)
IF (ierror /= MPI_SUCCESS) CALL xt_abort(context//'MPI_FINALIZE failed', __FILE__, __LINE__)
END SUBROUTINE exit_all
SUBROUTINE id_map(map)
INTEGER, INTENT(out) :: map(:,:)
INTEGER :: i,j,p
p = 0
DO j = 1, SIZE(map,2)
DO i = 1, SIZE(map,1)
p = p + 1
map(i,j) = p
ENDDO
ENDDO
END SUBROUTINE id_map
SUBROUTINE exchange_int(redist, f, g)
TYPE(xt_redist), INTENT(in) :: redist
INTEGER, TARGET, INTENT(in) :: f(*)
INTEGER, TARGET, VOLATILE, INTENT(out) :: g(*)
CALL xt_redist_s_exchange1(redist, c_loc(f), c_loc(g));
END SUBROUTINE exchange_int
SUBROUTINE gen_trans(xmap, send_dt, recv_dt, redist)
TYPE(xt_xmap), INTENT(in) :: xmap
INTEGER,INTENT(in) :: send_dt, recv_dt
TYPE(xt_redist),INTENT(out) :: redist
INTEGER :: dt
IF (send_dt /= recv_dt) CALL xt_abort('gen_trans: (send_dt /= recv_dt) unsupported', __FILE__, __LINE__)
dt = send_dt
redist = xt_redist_p2p_new(xmap, dt)
!CALL ut_init_transposition(itemp, dt, itrans)
END SUBROUTINE gen_trans
SUBROUTINE gen_off_trans(xmap, send_dt, send_off, recv_dt, recv_off, redist)
TYPE(xt_xmap), INTENT(in) :: xmap
INTEGER,INTENT(in) :: send_dt, recv_dt
INTEGER,INTENT(in) :: send_off(:,:), recv_off(:,:)
TYPE(xt_redist),INTENT(out) :: redist
!INTEGER :: send_offsets(SIZE(send_off)), recv_offsets(SIZE(recv_off))
!send_offsets = RESHAPE(send_off, (/SIZE(send_off)/) )
!recv_offsets = RESHAPE(recv_off, (/SIZE(recv_off)/) )
IF (recv_dt /= send_dt) THEN
CALL xt_abort('(datatype_in /= datatype_out) not supported', __FILE__, __LINE__)
ENDIF
redist = xt_redist_p2p_off_new(xmap, send_off, recv_off, send_dt);
!CALL ut_init_transposition(itemp, send_offsets, recv_offsets, send_dt, recv_dt, itrans)
END SUBROUTINE gen_off_trans
SUBROUTINE get_window(gval, win)
INTEGER, INTENT(in) :: gval(:,:)
INTEGER, INTENT(out) :: win(:,:)
INTEGER :: i, j, ig, jg
DO j = 1, je
jg = p_joff + j
DO i = 1, ie
ig = p_ioff + i
win(i,j) = gval(ig,jg)
ENDDO
ENDDO
END SUBROUTINE get_window
SUBROUTINE gen_template(local_src_idx, local_dst_idx, xmap)
INTEGER, INTENT(in) :: local_src_idx(:,:)
INTEGER, INTENT(in) :: local_dst_idx(:,:)
TYPE(xt_xmap), INTENT(out) :: xmap
TYPE(xt_idxlist) :: src_idxlist, dst_idxlist
src_idxlist = xt_idxvec_new(local_src_idx, g_ie*g_je)
dst_idxlist = xt_idxvec_new(local_dst_idx, g_ie*g_je)
xmap = xt_xmap_all2all_new(src_idxlist, dst_idxlist, MPI_COMM_WORLD)
CALL xt_idxlist_delete(src_idxlist)
CALL xt_idxlist_delete(dst_idxlist)
END SUBROUTINE gen_template
SUBROUTINE def_exchange(id_in, id_out)
INTEGER, INTENT(in) :: id_in(:,:)
INTEGER, INTENT(out) :: id_out(:,:)
LOGICAL, PARAMETER :: increased_north_halo = .FALSE.
LOGICAL, PARAMETER :: with_north_halo = .true.
INTEGER :: i, j
INTEGER :: g_core_is, g_core_ie, g_core_js, g_core_je
INTEGER :: north_halo
! global core domain:
g_core_is = nhalo + 1
g_core_ie = g_ie-nhalo
g_core_js = nhalo + 1
g_core_je = g_je-nhalo
! global tripolar boundsexchange:
id_out = undef_index
id_out(g_core_is:g_core_ie, g_core_js:g_core_je) = id_in(g_core_is:g_core_ie, g_core_js:g_core_je)
IF (with_north_halo) THEN
! north inversion, (maybe with increased north halo)
IF (increased_north_halo) THEN
north_halo = nhalo+1
ELSE
north_halo = nhalo
ENDIF
IF (2*north_halo > g_core_je) CALL xt_abort('def_exchange: grid too small (or halo too large) for tripolar north exchange', &
& __FILE__, __LINE__)
DO j = 1, north_halo
DO i = g_core_is, g_core_ie
id_out(i,j) = id_out(g_core_ie + (g_core_is-i), 2*north_halo + (1-j))
ENDDO
ENDDO
ELSE
DO j = 1, nhalo
DO i = nhalo+1, g_ie-nhalo
id_out(i,j) = id_in(i,j)
ENDDO
ENDDO
ENDIF
! south: no change
DO j = g_core_je+1, g_je
DO i = nhalo+1, g_ie-nhalo
id_out(i,j) = id_in(i,j)
ENDDO
ENDDO
! PBC
DO j = 1, g_je
DO i = 1, nhalo
id_out(g_core_is-i,j) = id_out(g_core_ie+(1-i),j)
ENDDO
DO i = 1, nhalo
id_out(g_core_ie+i,j) = id_out(nhalo+i,j)
ENDDO
ENDDO
CALL check_g_idx(id_out)
END SUBROUTINE def_exchange
SUBROUTINE check_g_idx(gidx)
INTEGER,INTENT(in) :: gidx(:,:)
INTEGER :: i,j
DO j=1,g_je
DO i=1,g_ie
IF (gidx(i,j) == undef_index) THEN
CALL xt_abort('check_g_idx: check failed', __FILE__, __LINE__)
ENDIF
ENDDO
ENDDO
END SUBROUTINE check_g_idx
SUBROUTINE deco
INTEGER :: cx0(0:nprocx-1), cxn(0:nprocx-1)
INTEGER :: cy0(0:nprocy-1), cyn(0:nprocy-1)
CALL regular_deco(g_ie-2*nhalo, cx0, cxn)
CALL regular_deco(g_je-2*nhalo, cy0, cyn)
! process local deco variables:
ie = cxn(mypx) + 2*nhalo
je = cyn(mypy) + 2*nhalo
p_ioff = cx0(mypx)
p_joff = cy0(mypy)
END SUBROUTINE deco
SUBROUTINE regular_deco(g_cn, c0, cn)
INTEGER, INTENT(in) :: g_cn
INTEGER, INTENT(out) :: c0(0:), cn(0:)
! convention: process space coords start at 0, grid point coords start at 1
integer :: tn
INTEGER :: d, m
INTEGER :: it
tn = SIZE(c0)
IF (tn<0) CALL xt_abort('(tn<0)', __FILE__, __LINE__)
IF (tn>g_cn) CALL xt_abort('regular_deco: too many task for such a core region', __FILE__, __LINE__)
d = g_cn/tn
m = MOD(g_cn, tn)
DO it = 0, m-1
cn(it) = d + 1
ENDDO
DO it = m, tn-1
cn(it) = d
ENDDO
c0(0)=0
DO it = 1, tn-1
c0(it) = c0(it-1) + cn(it-1)
ENDDO
IF (c0(tn-1)+cn(tn-1) /= g_cn) CALL xt_abort('regular_deco: internal error 1', __FILE__, __LINE__)
END SUBROUTINE regular_deco
SUBROUTINE factorize(c, a, b)
INTEGER, INTENT(in) :: c
INTEGER, INTENT(out) :: a, b ! c = a*b
INTEGER :: x0, i
IF (c<1) CALL xt_abort('factorize: invalid process space', __FILE__, __LINE__)
IF (c <= 3 .OR. c == 5 .OR. c == 7) THEN
a = c
b = 1
RETURN
ENDIF
! simple approach, we try to be near c = (2*x) * x
x0 = INT(SQRT(0.5*c) + 0.5)
a = 2*x0
f_loop: DO i = a, 1, -1
IF (MOD(c,a) == 0) THEN
b = c/a
EXIT f_loop
ENDIF
a = a - 1
ENDDO f_loop
END SUBROUTINE factorize
END PROGRAM test_yaxt
#! @SHELL@
@MPI_LAUNCH@ -n 1 ./test_yaxt
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment