! $Id: ESMF_FieldSMMEx.F90,v 1.33 2011/07/11 21:52:34 svasquez Exp $ ! ! Earth System Modeling Framework ! Copyright 2002-2011, University Corporation for Atmospheric Research, ! Massachusetts Institute of Technology, Geophysical Fluid Dynamics ! Laboratory, University of Michigan, National Centers for Environmental ! Prediction, Los Alamos National Laboratory, Argonne National Laboratory, ! NASA Goddard Space Flight Center. ! Licensed under the University of Illinois-NCSA License. ! !============================================================================== program FieldSMMEx !------------------------------------------------------------------------- !ESMF_MULTI_PROC_EXAMPLE String used by test script to count examples. !============================================================================== ! ! !PROGRAM: ESMF_FieldSMMEx - Field SMMribution ! ! !DESCRIPTION: ! ! This program shows examples of Field interfaces for ! sparse matrix multiplication of data. !----------------------------------------------------------------------------- #include "ESMF.h" #include "ESMF_Macros.inc" #undef ESMF_METHOD #define ESMF_METHOD "ESMF_FieldSMMEx" ! ESMF Framework module use ESMF use ESMF_TestMod implicit none !------------------------------------------------------------------------------ ! The following line turns the CVS identifier string into a printable variable. character(*), parameter :: version = & '$Id: ESMF_FieldSMMEx.F90,v 1.33 2011/07/11 21:52:34 svasquez Exp $' !------------------------------------------------------------------------------ ! Local variables integer :: rc, finalrc type(ESMF_Field) :: srcField, dstField type(ESMF_Field) :: srcFieldA, dstFieldA type(ESMF_Grid) :: grid, dstGrid type(ESMF_DistGrid) :: distgrid type(ESMF_VM) :: vm type(ESMF_RouteHandle) :: routehandle type(ESMF_Array) :: srcArray, dstArray type(ESMF_ArraySpec) :: arrayspec integer :: localrc, lpe, i integer, allocatable :: src_farray(:), dst_farray(:) integer :: fa_shape(1), tlb(1), tub(1) integer, pointer :: fptr(:) real(ESMF_KIND_R4), pointer :: fptr2d(:,:) real(ESMF_KIND_R4), pointer :: src_farray2(:) real(ESMF_KIND_R4), allocatable :: factorList(:) integer, allocatable :: factorIndexList(:,:) rc = ESMF_SUCCESS finalrc = ESMF_SUCCESS !------------------------------------------------------------------------------ call ESMF_Initialize(defaultlogfilename="FieldSMMEx.Log", & logkindflag=ESMF_LOGKIND_MULTI, rc=rc) if (rc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) if (.not. ESMF_TestMinPETs(4, ESMF_SRCLINE)) & call ESMF_Finalize(endflag=ESMF_END_ABORT) !------------------------------------------------------------------------------ !BOE ! \subsubsection{Sparse matrix multiplication from source Field to destination Field} ! \label{sec:field:usage:smm_1dptr} ! ! A user can use {\tt ESMF\_FieldSMM()} interface to perform sparse matrix multiplication ! from ! source Field to destination Field. This interface is overloaded by type and kind; ! ! In this example, we first create two 1D Fields, a source Field and a destination ! Field. Then we use {\tt ESMF\_FieldSMM} to ! perform sparse matrix multiplication from source Field to destination Field. ! ! The source and destination Field data are arranged such that each of the 4 PETs has 4 ! data elements. Moreover, the source Field has all its data elements initialized to a linear ! function based on local PET number. ! Then collectively on each PET, a SMM according to the following formula ! is preformed: \newline ! $dstField(i) = i * srcField(i), i = 1 ... 4$ \newline ! \newline ! ! Because source Field data are initialized to a linear function based on local PET number, ! the formula predicts that ! the result destination Field data on each PET is {1,2,3,4}. This is verified in the ! example. ! ! Section \ref{Array:SparseMatMul} provides a detailed discussion of the ! sparse matrix mulitiplication operation implemented in ESMF. ! !EOE !BOC ! Get current VM and pet number call ESMF_VMGetCurrent(vm, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_VMGet(vm, localPet=lpe, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! create distgrid and grid distgrid = ESMF_DistGridCreate(minIndex=(/1/), maxIndex=(/16/), & regDecomp=(/4/), & rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE grid = ESMF_GridCreate(distgrid=distgrid, & name="grid", rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_GridGetFieldBounds(grid, localDe=0, totalCount=fa_shape, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! create src\_farray, srcArray, and srcField ! +--------+--------+--------+--------+ ! 1 2 3 4 ! value ! 1 4 8 12 16 ! bounds allocate(src_farray(fa_shape(1)) ) src_farray = lpe+1 srcArray = ESMF_ArrayCreate(distgrid, src_farray, & indexflag=ESMF_INDEX_DELOCAL, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE srcField = ESMF_FieldCreate(grid, srcArray, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! create dst_farray, dstArray, and dstField ! +--------+--------+--------+--------+ ! 0 0 0 0 ! value ! 1 4 8 12 16 ! bounds allocate(dst_farray(fa_shape(1)) ) dst_farray = 0 dstArray = ESMF_ArrayCreate(distgrid, dst_farray, & indexflag=ESMF_INDEX_DELOCAL, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE dstField = ESMF_FieldCreate(grid, dstArray, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! perform sparse matrix multiplication ! 1. setup routehandle from source Field to destination Field ! initialize factorList and factorIndexList allocate(factorList(4)) allocate(factorIndexList(2,4)) factorList = (/1,2,3,4/) factorIndexList(1,:) = (/lpe*4+1,lpe*4+2,lpe*4+3,lpe*4+4/) factorIndexList(2,:) = (/lpe*4+1,lpe*4+2,lpe*4+3,lpe*4+4/) call ESMF_FieldSMMStore(srcField, dstField, routehandle, & factorList, factorIndexList, rc=localrc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! 2. use precomputed routehandle to perform SMM call ESMF_FieldSMM(srcfield, dstField, routehandle, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! verify sparse matrix multiplication call ESMF_FieldGet(dstField, localDe=0, farrayPtr=fptr, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! Verify that the result data in dstField is correct. ! Before the SMM op, the dst Field contains all 0. ! The SMM op reset the values to the index value, verify this is the case. ! +--------+--------+--------+--------+ ! 1 2 3 4 2 4 6 8 3 6 9 12 4 8 12 16 ! value ! 1 4 8 12 16 ! bounds do i = lbound(fptr, 1), ubound(fptr, 1) if(fptr(i) /= i*(lpe+1)) rc = ESMF_FAILURE enddo !EOC if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !BOE ! Field sparse matrix matmul can also be performed between weakly congruent Fields. ! In this case, source and destination Fields can have ungridded dimensions ! with size different from the Field pair used to compute the routehandle. !EOE !BOC call ESMF_ArraySpecSet(arrayspec, typekind=ESMF_TYPEKIND_I4, rank=2, rc=rc) !EOC if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !BOE ! Create two fields with ungridded dimensions using the Grid created previously. ! The new Field pair has matching number of elements. The ungridded dimension ! is mapped to the first dimension of either Field. !EOE !BOC srcFieldA = ESMF_FieldCreate(grid, arrayspec, gridToFieldMap=(/2/), & ungriddedLBound=(/1/), ungriddedUBound=(/10/), rc=rc) !EOC if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !BOC dstFieldA = ESMF_FieldCreate(grid, arrayspec, gridToFieldMap=(/2/), & ungriddedLBound=(/1/), ungriddedUBound=(/10/), rc=rc) !EOC if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !BOE ! Using the previously computed routehandle, weakly congruent Fields can perform ! sparse matrix matmul. !EOE !BOC call ESMF_FieldSMM(srcfieldA, dstFieldA, routehandle, rc=rc) !EOC if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !BOC ! release route handle call ESMF_FieldSMMRelease(routehandle, rc=rc) !EOC if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! destroy all objects created in this example to prevent memory leak call ESMF_FieldDestroy(srcField, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_FieldDestroy(dstField, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_FieldDestroy(srcFieldA, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_FieldDestroy(dstFieldA, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_ArrayDestroy(srcArray, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_ArrayDestroy(dstArray, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_GridDestroy(grid, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_DistGridDestroy(distgrid, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE deallocate(src_farray, dst_farray, factorList, factorIndexList) !BOE ! In the following discussion, we demonstrate how to set up a SMM routehandle ! between a pair of Fields that are different in number of gridded dimensions ! and the size of those gridded dimensions. The source Field has a 1D decomposition ! with 16 total elements; the destination Field has a 2D decomposition with ! 12 total elements. For ease of understanding of the actual matrix calculation, ! a global indexing scheme is used. !EOE !BOC distgrid = ESMF_DistGridCreate(minIndex=(/1/), maxIndex=(/16/), & indexflag=ESMF_INDEX_GLOBAL, & regDecomp=(/4/), & rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE grid = ESMF_GridCreate(distgrid=distgrid, & indexflag=ESMF_INDEX_GLOBAL, & name="grid", rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_GridGetFieldBounds(grid, localDe=0, totalLBound=tlb, & totalUBound=tub, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !EOC !BOE ! create 1D src\_farray, srcArray, and srcField !\begin{verbatim} ! + PET0 + PET1 + PET2 + PET3 + ! +--------+--------+--------+--------+ ! 1 2 3 4 ! value ! 1 4 8 12 16 ! bounds of seq indices !\end{verbatim} !EOE !BOC allocate(src_farray2(tlb(1):tub(1)) ) src_farray2 = lpe+1 srcArray = ESMF_ArrayCreate(distgrid, src_farray2, & indexflag=ESMF_INDEX_GLOBAL, & rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !print *, lpe, '+', tlb, tub, '+', src_farray2 srcField = ESMF_FieldCreate(grid, srcArray, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !EOC !BOE ! Create 2D dstField on the following distribution ! (numbers are the sequence indices): !\begin{verbatim} ! + PET0 + PET1 + PET2 + PET3 + ! +--------+--------+--------+--------+ ! | | | | | ! | 1 | 4 | 7 | 10 | ! | | | | | ! +--------+--------+--------+--------+ ! | | | | | ! | 2 | 5 | 8 | 11 | ! | | | | | ! +--------+--------+--------+--------+ ! | | | | | ! | 3 | 6 | 9 | 12 | ! | | | | | ! +--------+--------+--------+--------+ !\end{verbatim} !EOE !BOC ! Create the destination Grid dstGrid = ESMF_GridCreateNoPeriDim(minIndex=(/1,1/), maxIndex=(/3,4/), & indexflag = ESMF_INDEX_GLOBAL, & regDecomp = (/1,4/), & rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE dstField = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R4, & indexflag=ESMF_INDEX_GLOBAL, & rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !EOC !BOE ! Perform sparse matrix multiplication $dst_i$ = $M_{i,j}$ * $src_j$ ! First setup routehandle from source Field to destination Field ! with prescribed factorList and factorIndexList. ! ! The sparse matrix is of size 12x16, however only the following entries ! are filled: ! \begin{verbatim} ! M(3,1) = 0.1 ! M(3,10) = 0.4 ! M(8,2) = 0.25 ! M(8,16) = 0.5 ! M(12,1) = 0.3 ! M(12,16) = 0.7 ! \end{verbatim} ! ! By the definition of matrix calculation, the 8th element on PET2 in the ! dstField equals to 0.25*srcField(2) + 0.5*srcField(16) = 0.25*1+0.5*4=2.25 ! For simplicity, we will load the factorList and factorIndexList on ! PET 0 and 1, the SMMStore engine will load balance the parameters on all 4 ! PETs internally for optimal performance. !EOE !BOC if(lpe == 0) then allocate(factorList(3), factorIndexList(2,3)) factorList=(/0.1,0.4,0.25/) factorIndexList(1,:)=(/1,10,2/) factorIndexList(2,:)=(/3,3,8/) call ESMF_FieldSMMStore(srcField, dstField, routehandle=routehandle, & factorList=factorList, factorIndexList=factorIndexList, rc=localrc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE else if(lpe == 1) then allocate(factorList(3), factorIndexList(2,3)) factorList=(/0.5,0.3,0.7/) factorIndexList(1,:)=(/16,1,16/) factorIndexList(2,:)=(/8,12,12/) call ESMF_FieldSMMStore(srcField, dstField, routehandle=routehandle, & factorList=factorList, factorIndexList=factorIndexList, rc=localrc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE else call ESMF_FieldSMMStore(srcField, dstField, routehandle=routehandle, & rc=localrc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE endif ! 2. use precomputed routehandle to perform SMM call ESMF_FieldSMM(srcfield, dstField, routehandle=routehandle, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !EOC ! verify sparse matrix multiplication call ESMF_FieldGet(dstField, localDe=0, farrayPtr=fptr2d, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE ! Verify that the result data in dstField is correct. print *, lpe, '-', lbound(fptr2d), '-',ubound(fptr2d),'-', fptr2d ! destroy all objects created in this example to prevent memory leak call ESMF_FieldSMMRelease(routehandle, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_FieldDestroy(srcField, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_FieldDestroy(dstField, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !call ESMF_FieldDestroy(srcFieldA, rc=rc) !if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE !call ESMF_FieldDestroy(dstFieldA, rc=rc) !if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_ArrayDestroy(srcArray, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_GridDestroy(grid, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_GridDestroy(dstgrid, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE call ESMF_DistGridDestroy(distgrid, rc=rc) if(rc .ne. ESMF_SUCCESS) finalrc = ESMF_FAILURE deallocate(src_farray2) if(allocated(factorList)) deallocate(factorList, factorIndexList) call ESMF_Finalize(rc=rc) if (rc.NE.ESMF_SUCCESS) then finalrc = ESMF_FAILURE end if if (finalrc.EQ.ESMF_SUCCESS) then print *, "PASS: ESMF_FieldSMMEx.F90" else print *, "FAIL: ESMF_FieldSMMEx.F90" end if end program FieldSMMEx