Whether the metis_tests.f90 has bug?

[trimtrim]

at line 55 in metis_tests.f90 code: call write_graph("test1.graph",fxadj,fadjncy,1) fxadj and fadjncy are C style numbering. The program generates the following error message:

forrtl: severe (408): fort: (11): Subscript #1 of the array ADJNCY has value 0 which is less than the lower bound of 1

Image PC Routine Line Source libifcoremdd.dll 00007FFF8E4E5597 Unknown Unknown Unknown fmetis_test.exe 00007FF73BB920AE METIS_IO_mp_WRITE 45 metis_io.f90 fmetis_test.exe 00007FF73BBA225F METIS_TESTS_mp_TE 75 metis_tests.f90 fmetis_test.exe 00007FF73BBA433C MAIN__ 235 metis_tests.f90 fmetis_test.exe 00007FF73BBB8992 Unknown Unknown Unknown fmetis_test.exe 00007FF73BBE7994 Unknown Unknown Unknown fmetis_test.exe 00007FF73BBE78BE Unknown Unknown Unknown fmetis_test.exe 00007FF73BBE777E Unknown Unknown Unknown fmetis_test.exe 00007FF73BBE7A09 Unknown Unknown Unknown KERNEL32.DLL 00007FF80EB53DC4 Unknown Unknown Unknown ntdll.dll 00007FF80ECA3691 Unknown Unknown Unknown

[ivan-pi]

Hello trimtrim,

you are correct, the arrays fxadj and fadjncy assume C-style numbering in this example. The numflag dummy variable should be set to 0. Surprisingly, I got no error with gfortran. I'm assuming this is because the arrays are actually pointer to C memory, but would need to check first.

Please note that an advanced version of the write_graph subroutine can be found in the file metis_oo_interface.f90, that matches the entire METIS input file format and can export vertex weights, vertex sizes and edge weights. Most likely I will delete the file metis_io.f90 completely in the near future. Moreover IO functionality for graphs and meshes is already present in the METIS library (in the file io.c), and it might make more sense to just gain access to those functions (although they are not part of the "public" API) than rewrite everything on the Fortran side.

[trimtrim]

Hi Ivan, I am glad it helped. I have another question. For the METIS_MeshToNodal interface, the arguments "xadj", "adjncy" are declared as type(C_PTR) and it needs to use c_f_pointer to convert it to Fortran pointer. This might difficult for Fortran programmer. The "xadj", "adjncy" are 1 dimensional (1D) array in C. I am wondering whether the Fortran interface can declare the two arguments as Fortran 1D array? such as: INTEGER(C_INT),DIMENSION(*) ***** ::xadj,adjncy

Declare as 1D array can avoid use function c_f_pointer to convert the c pointer to Fortran pointer and simplify the coding?

[ivan-pi]

I agree that the need to explictly call the subroutine c_f_pointer to gain access to the values of the array is not very user friendly. Also the fact that the correct size of adjncy must be inferred from the last element of xadj is not that elegant and took me half an hour to figure it out.

What you suggest is using an assumed-size array and would most likely not work (give it a try?). To cite from this discussion: "An assumed-size array can only be used in whole array expressions when it's an actual argument in a procedure reference that does not require the array's shape."

The "quickest fix" at the moment is to write a wrapper which copies the array into Fortran allocated arrays. Something along the lines of:

subroutine ForMETIS_MeshToNodal(ne,nn,eptr,eind,numflag,xadj,adjncy,stat)
        use iso_c_binding, only : c_int, c_ptr, c_f_pointer
        integer, intent(in) :: ne
        integer, intent(in) :: nn
        integer, intent(in) :: eptr(ne+1)
        integer, intent(in) :: eind(:)
        integer, intent(in) :: numflag
        integer, intent(out), allocatable :: xadj(:)
        integer, intent(out), allocatable :: adjncy(:)
        integer, intent(out) :: stat

        type(c_ptr) :: c_xadj, c_adjncy

        integer(c_int), pointer :: f_xadj(:) => null(), f_adjncy(:) => null()

        stat = METIS_MeshToNodal(ne,nn,eptr,eind,numflag,c_xadj,c_adjncy)
        if (stat < 0) return

        call c_f_pointer(c_xadj,f_xadj,shape=[nn+1])

        select case(numflag)
        case(0) ! C-style numbering
            call c_f_pointer(c_adjncy,f_adjncy,shape=[f_xadj(nn+1)])
        case(1) ! Fortran-style numbering
            call c_f_pointer(c_adjncy,f_adjncy,shape=[f_xadj(nn+1)-1])
        case default
            write(*,*) "Wrong flag!"
            error stop 1
        end select

        allocate(xadj,source=f_xadj)
        allocate(adjncy,source=f_adjncy)

        stat = METIS_Free(c_xadj); f_xadj => null() 
        ! if (stat < 0) return
        stat = METIS_Free(c_adjncy); f_adjncy => null()
        ! if (stat < 0) return
end subroutine

By only working with the "public" METIS API, the price to pay then is the memory and time overhead needed to copy the values from C-allocated to Fortran-allocated entities.

The second (and better) way to fix this issue would be to modify the METIS source files and use the features defined in Technical Specification TS29113 "Further Interoperability C" which makes it possible to pass allocatable and pointer arrays to and from C, through a so-called "C descriptor". Specifically, this would mean going to the file "$METIS_FOLDER/libmetis/mesh.c" and modifying the function METIS_MeshToNodal and make it return a Fortran allocatable array. It would probably be an interesting project to learn how to use these new interoperability features, but at the moment is not on my priority list. (You can find a small example using C descriptors to allocate a Fortran array from C in this blog post by Steve Lionel from Intel.)

Note that on the main METIS home page a survey is open regarding the next release METIS 5.x. One of the questions in the survey is whether a Fortran interface should be included. Perhaps the original METIS developers are already planning to do something like this.

[ivan-pi]

Using the concept of C descriptors it is possible to copy the adjacency graph directly into Fortran allocatable arrays.

First I created a new file mesh_fortran.c in the metis-5.1.0/libmetis folder containing the following code:

#include "metislib.h"
#include "ISO_Fortran_binding.h"

inline void check_CFI_result(int res,int *stat) {
    if (res) { /* (de)allocation failed */
        if (stat) {
            *stat = res;
            return;
        } else { 
            exit(2);
        }
    }
}

int FMETIS_MeshToNodal(idx_t *ne, idx_t *nn, idx_t *eptr, idx_t *eind, 
                  idx_t *numflag, CFI_cdesc_t *fxadj, CFI_cdesc_t *fadjncy, 
                  int *stat /* NULL if not needed */)
{

    idx_t *xadj;
    idx_t *adjncy;

    /* Call METIS Core routine */
    int status = METIS_MeshToNodal(ne,nn,eptr,eind,numflag,&xadj,&adjncy);

    if(fxadj->base_addr) {
        int res = CFI_deallocate(fxadj);
        check_CFI_result(res,stat);
    }

    CFI_index_t lower[fxadj->rank], upper[fxadj->rank];
    lower[0] = 1;
    size_t dummy;

     /* Allocate memory for xadj, since you know its size. */
    upper[0] = ((CFI_index_t) *nn) + 1;
    int res = CFI_allocate(fxadj, lower, upper, dummy);
    check_CFI_result(res,stat);

    char *elt = (char *) (fxadj->base_addr); // first element
    for (idx_t i = 0; i <= *nn; i++) {
        * (int *) elt = xadj[i]; // copy value
        elt += fxadj->dim[0].sm; // increment by one memory stride
    }

    if(fadjncy->base_addr) {
        int res = CFI_deallocate(fadjncy);
        check_CFI_result(res,stat);
    }

    /* Allocate memory for adjncy, since you also know its size. */
    upper[0] = ((CFI_index_t) xadj[*nn]);
    res = CFI_allocate(fadjncy, lower, upper, dummy);
    check_CFI_result(res,stat);

    elt = (char *) (fadjncy->base_addr); // first element
    for (idx_t i = 0; i < xadj[*nn]; i++) {
        * (int *) elt = adjncy[i]; // copy value
        elt += fadjncy->dim[0].sm; // increment by one memory stride
    }

    free(xadj);
    free(adjncy);

    return status;
}

The following line are then added to the include header file metis-5.1.0/include/metis.h:

#include "ISO_Fortran_binding.h"

 /* ... Original code ... */

METIS_API(int) FMETIS_MeshToNodal(idx_t *ne, idx_t *nn, idx_t *eptr, idx_t *eind, 
                  idx_t *numflag, CFI_cdesc_t *fxadj, CFI_cdesc_t *fadjncy,
                  int *stat /* NULL if not needed */);

The METIS library should then be rebuilt. While gfortran 9 (install instructions here) added support for the "ISO_Fortran_binding.h" file, I get an internal compiler error. Building METIS with the Intel compiler works however.

Next I wrote the Fortran wrapper for the new MeshToNodal routine:

module fmetis

    use iso_c_binding
    implicit none

    integer(c_int), parameter :: METIS_OK = 1_c_int
    integer(c_int), parameter :: METIS_NOPTIONS = 40_c_int

    interface
    function FMETIS_MeshToNodal(ne,nn,eptr,eind,numflag,fxadj,fadjncy,stat) result(ierr) bind(C,name="FMETIS_MeshToNodal")
        import c_int, c_ptr

        ! Parameters
        integer(c_int), intent(in) :: ne
        integer(c_int), intent(in) :: nn
        integer(c_int), intent(in), dimension(*) :: eptr, eind
        integer(c_int), intent(in) :: numflag
        integer(c_int), intent(out), allocatable :: fxadj(:), fadjncy(:)
        integer(c_int), intent(out), optional :: stat ! Allocation status

        ! Returns
        integer(c_int) :: ierr
    end function
    end interface

end module

Finally, a usage example:

program test_fmetis_mesh2nodal

    use, intrinsic :: iso_c_binding
    use fmetis

    implicit none

    integer(c_int), parameter :: ne = 3 ! number of elements
    integer(c_int), parameter :: nn = 8 ! number of nodes
    integer(c_int), parameter :: npel = 4 ! nodes per element

    integer(c_int) :: eptr(ne+1)
    integer(c_int) :: eind(ne*npel)
    integer(c_int) :: epart(ne), npart(nn)
    integer(c_int) :: options(0:METIS_NOPTIONS-1)
    integer(c_int) :: ierr, objval, istat

    type(c_ptr) :: c_xadj, c_adjncy
    integer(c_int), pointer :: xadj(:) => null(), adjncy(:) => null()

    integer(c_int), allocatable :: fxadj(:), fadjncy(:)

    ! 0---1---4---6
    ! | 0 | 1 | 2 |
    ! 3---2---5---7

    eptr = [0,4,8,12]
    eind = [0,1,2,3,1,4,5,2,4,6,7,5]

    ierr = FMETIS_MeshToNodal(ne,nn,eptr,eind,numflag=0,fxadj=fxadj,fadjncy=fadjncy,stat=istat)
    if (ios /= METIS_OK) then
        write(*,*) "METIS_MeshToNodal failed with error: ", ios
        error stop 1
    end if
    if (istat /= 0_c_int) then
        write(*,*) "Fortran (de)allocation error: ", istat
        error stop 1
    end if

    write(*,'(A,*(I0,:,1X))') "fxadj = ", fxadj
    write(*,'(A,*(I1,:,1X))') "fadjncy = ", fadjncy

end program

The output of the program is:

fxadj = 0 3 8 13 16 21 26 29 32
fadjncy = 1 2 3 0 2 3 4 5 0 1 3 4 5 0 1 2 1 5 2 6 7 1 4 2 6 7 4 7 5 4 6 5

and can be checked against the drawing of the mesh for correctness. While the interface is now "clean" on the Fortran side, we have in fact achieved the same thing that could be done from the Fortran side. Bigger changes to METIS would be needed in order to directly return the adjacency graph as a Fortran array, without copying.