Bugfix/fortran nvgetarraypointer (#214)

Temporary fix for issue #160 
* add GetLocalLength operation
* add GetSubvectorLocalLength to manyvector

Co-authored-by: David Gardner <[email protected]>
Co-authored-by: Daniel R. Reynolds <[email protected]>

examples/arkode/F2003_parallel/CMakeLists.txt

@@ -16,12 +16,14 @@
 
 # Example lists are tuples "name\;nodes\;tasks\;type" where the
 # type is develop for examples excluded from 'make test' in releases
-
-# Examples using SUNDIALS linear solvers
-set(FARKODE_examples
-  "ark_brusselator1D_task_local_nls_f2003\;--monitor\;1\;4\;develop\;2"
-  "ark_brusselator1D_task_local_nls_f2003\;--monitor --global-nls\;1\;4\;develop\;2"
-  "ark_brusselator1D_task_local_nls_f2003\;--monitor --explicit --tf 3\;1\;4\;develop\;2")
+if(SUNDIALS_LOGGING_ENABLE_MPI)
+  set(FARKODE_examples
+    "ark_brusselator1D_task_local_nls_f2003\;--monitor\;1\;4\;develop\;2"
+    "ark_brusselator1D_task_local_nls_f2003\;--monitor --global-nls\;1\;4\;develop\;2"
+    "ark_brusselator1D_task_local_nls_f2003\;--monitor --explicit --tf 3\;1\;4\;develop\;2")
+else()
+  set(FARKODE_examples )
+endif()
 
 if(MPI_Fortran_COMPILER)
   # use MPI wrapper as the compiler

include/nvector/nvector_manyvector.h

@@ -49,10 +49,10 @@ extern "C" {
    ----------------------------------------------------------------- */
 
 struct _N_VectorContent_ManyVector {
-  sunindextype  num_subvectors;  /* number of vectors attached      */
-  sunindextype  global_length;   /* overall manyvector length       */
-  N_Vector*     subvec_array;    /* pointer to N_Vector array       */
-  booleantype   own_data;        /* flag indicating data ownership  */
+  sunindextype  num_subvectors;  /* number of vectors attached       */
+  sunindextype  global_length;   /* overall global manyvector length */
+  N_Vector*     subvec_array;    /* pointer to N_Vector array        */
+  booleantype   own_data;        /* flag indicating data ownership   */
 };
 
 typedef struct _N_VectorContent_ManyVector *N_VectorContent_ManyVector;
@@ -86,6 +86,7 @@ SUNDIALS_EXPORT void N_VDestroy_ManyVector(N_Vector v);
 SUNDIALS_EXPORT void N_VSpace_ManyVector(N_Vector v, sunindextype *lrw,
                                          sunindextype *liw);
 SUNDIALS_EXPORT sunindextype N_VGetLength_ManyVector(N_Vector v);
+SUNDIALS_EXPORT sunindextype N_VGetSubvectorLocalLength_ManyVector(N_Vector v, sunindextype vec_num);
 SUNDIALS_EXPORT void N_VLinearSum_ManyVector(realtype a, N_Vector x,
                                              realtype b, N_Vector y,
                                              N_Vector z);

include/nvector/nvector_mpimanyvector.h

@@ -52,10 +52,10 @@ extern "C" {
 
 struct _N_VectorContent_MPIManyVector {
   MPI_Comm      comm;            /* overall MPI communicator        */
-  sunindextype  num_subvectors;  /* number of vectors attached      */
-  sunindextype  global_length;   /* overall manyvector length       */
-  N_Vector*     subvec_array;    /* pointer to N_Vector array       */
-  booleantype   own_data;        /* flag indicating data ownership  */
+  sunindextype  num_subvectors;  /* number of vectors attached       */
+  sunindextype  global_length;   /* overall global manyvector length */
+  N_Vector*     subvec_array;    /* pointer to N_Vector array        */
+  booleantype   own_data;        /* flag indicating data ownership   */
 };
 
 typedef struct _N_VectorContent_MPIManyVector *N_VectorContent_MPIManyVector;
@@ -95,6 +95,7 @@ SUNDIALS_EXPORT void N_VSpace_MPIManyVector(N_Vector v, sunindextype *lrw,
                                             sunindextype *liw);
 SUNDIALS_EXPORT void *N_VGetCommunicator_MPIManyVector(N_Vector v);
 SUNDIALS_EXPORT sunindextype N_VGetLength_MPIManyVector(N_Vector v);
+SUNDIALS_EXPORT sunindextype N_VGetSubvectorLocalLength_MPIManyVector(N_Vector v, sunindextype vec_num);
 SUNDIALS_EXPORT void N_VLinearSum_MPIManyVector(realtype a, N_Vector x,
                                                 realtype b, N_Vector y,
                                                 N_Vector z);

include/sundials/sundials_nvector.h

@@ -182,6 +182,10 @@ struct _generic_N_Vector_Ops
   void (*nvprint)(N_Vector);
   void (*nvprintfile)(N_Vector, FILE*);
 
+  /* WARNING: this function should not be used, it is here as a temporary
+      fix for https://github.com/LLNL/sundials/issues/160 (fortran getarraypointer bug). */
+  sunindextype (*nvgetlocallength)(N_Vector);
+
 #ifdef __cplusplus
   _generic_N_Vector_Ops() = default;
 #endif
@@ -315,6 +319,11 @@ SUNDIALS_EXPORT void N_VSetVecAtIndexVectorArray(N_Vector* vs, int index, N_Vect
 SUNDIALS_EXPORT void N_VPrint(N_Vector v);
 SUNDIALS_EXPORT void N_VPrintFile(N_Vector v, FILE* outfile);
 
+
+  /* WARNING: this function should not be used, it is here as a temporary
+      fix for https://github.com/LLNL/sundials/issues/160 (fortran getarraypointer bug). */
+  SUNDIALS_DEPRECATED_EXPORT sunindextype N_VGetLocalLength(N_Vector v);
+
 #ifdef __cplusplus
 }
 #endif

src/arkode/fmod/farkode_erkstep_mod.f90

@@ -1097,18 +1097,18 @@ function FERKStepSetTable(arkode_mem, b) &
 swig_result = fresult
 end function
 
-function FERKStepSetTableNum(arkode_mem, itable) &
+function FERKStepSetTableNum(arkode_mem, etable) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
 integer(C_INT) :: swig_result
 type(C_PTR) :: arkode_mem
-integer(ARKODE_ERKTableID), intent(in) :: itable
+integer(ARKODE_ERKTableID), intent(in) :: etable
 integer(C_INT) :: fresult 
 type(C_PTR) :: farg1 
 integer(C_INT) :: farg2 
 
 farg1 = arkode_mem
-farg2 = itable
+farg2 = etable
 fresult = swigc_FERKStepSetTableNum(farg1, farg2)
 swig_result = fresult
 end function
@@ -1131,19 +1131,19 @@ subroutine SWIG_string_to_chararray(string, chars, wrap)
   wrap%size = len(string)
 end subroutine
 
-function FERKStepSetTableName(arkode_mem, table) &
+function FERKStepSetTableName(arkode_mem, etable) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
 integer(C_INT) :: swig_result
 type(C_PTR) :: arkode_mem
-character(kind=C_CHAR, len=*), target :: table
+character(kind=C_CHAR, len=*), target :: etable
 character(kind=C_CHAR), dimension(:), allocatable, target :: farg2_chars
 integer(C_INT) :: fresult 
 type(C_PTR) :: farg1 
 type(SwigArrayWrapper) :: farg2 
 
 farg1 = arkode_mem
-call SWIG_string_to_chararray(table, farg2_chars, farg2)
+call SWIG_string_to_chararray(etable, farg2_chars, farg2)
 fresult = swigc_FERKStepSetTableName(farg1, farg2)
 swig_result = fresult
 end function

src/nvector/manyvector/fmod/fnvector_manyvector_mod.c

@@ -366,6 +366,20 @@ SWIGEXPORT int64_t _wrap_FN_VGetLength_ManyVector(N_Vector farg1) {
 }
 
 
+SWIGEXPORT int64_t _wrap_FN_VGetSubvectorLocalLength_ManyVector(N_Vector farg1, int64_t const *farg2) {
+  int64_t fresult ;
+  N_Vector arg1 = (N_Vector) 0 ;
+  sunindextype arg2 ;
+  sunindextype result;
+
+  arg1 = (N_Vector)(farg1);
+  arg2 = (sunindextype)(*farg2);
+  result = N_VGetSubvectorLocalLength_ManyVector(arg1,arg2);
+  fresult = (sunindextype)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT void _wrap_FN_VLinearSum_ManyVector(double const *farg1, N_Vector farg2, double const *farg3, N_Vector farg4, N_Vector farg5) {
   realtype arg1 ;
   N_Vector arg2 = (N_Vector) 0 ;

src/nvector/manyvector/fmod/fnvector_manyvector_mod.f90

@@ -40,6 +40,7 @@ module fnvector_manyvector_mod
  public :: FN_VDestroy_ManyVector
  public :: FN_VSpace_ManyVector
  public :: FN_VGetLength_ManyVector
+ public :: FN_VGetSubvectorLocalLength_ManyVector
  public :: FN_VLinearSum_ManyVector
  public :: FN_VConst_ManyVector
  public :: FN_VProd_ManyVector
@@ -191,6 +192,15 @@ function swigc_FN_VGetLength_ManyVector(farg1) &
 integer(C_INT64_T) :: fresult
 end function
 
+function swigc_FN_VGetSubvectorLocalLength_ManyVector(farg1, farg2) &
+bind(C, name="_wrap_FN_VGetSubvectorLocalLength_ManyVector") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+integer(C_INT64_T), intent(in) :: farg2
+integer(C_INT64_T) :: fresult
+end function
+
 subroutine swigc_FN_VLinearSum_ManyVector(farg1, farg2, farg3, farg4, farg5) &
 bind(C, name="_wrap_FN_VLinearSum_ManyVector")
 use, intrinsic :: ISO_C_BINDING
@@ -642,7 +652,7 @@ function FN_VGetSubvectorArrayPointer_ManyVector(v, vec_num) &
 farg1 = c_loc(v)
 farg2 = vec_num
 fresult = swigc_FN_VGetSubvectorArrayPointer_ManyVector(farg1, farg2)
-call c_f_pointer(fresult, swig_result, [1])
+call c_f_pointer(fresult, swig_result, [FN_VGetSubvectorLocalLength_ManyVector(v, vec_num)])
 end function
 
 function FN_VSetSubvectorArrayPointer_ManyVector(v_data, v, vec_num) &
@@ -774,6 +784,22 @@ function FN_VGetLength_ManyVector(v) &
 swig_result = fresult
 end function
 
+function FN_VGetSubvectorLocalLength_ManyVector(v, vec_num) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT64_T) :: swig_result
+type(N_Vector), target, intent(inout) :: v
+integer(C_INT64_T), intent(in) :: vec_num
+integer(C_INT64_T) :: fresult 
+type(C_PTR) :: farg1 
+integer(C_INT64_T) :: farg2 
+
+farg1 = c_loc(v)
+farg2 = vec_num
+fresult = swigc_FN_VGetSubvectorLocalLength_ManyVector(farg1, farg2)
+swig_result = fresult
+end function
+
 subroutine FN_VLinearSum_ManyVector(a, x, b, y, z)
 use, intrinsic :: ISO_C_BINDING
 real(C_DOUBLE), intent(in) :: a

src/nvector/manyvector/fmod/fnvector_mpimanyvector_mod.c

@@ -400,6 +400,20 @@ SWIGEXPORT int64_t _wrap_FN_VGetLength_MPIManyVector(N_Vector farg1) {
 }
 
 
+SWIGEXPORT int64_t _wrap_FN_VGetSubvectorLocalLength_MPIManyVector(N_Vector farg1, int64_t const *farg2) {
+  int64_t fresult ;
+  N_Vector arg1 = (N_Vector) 0 ;
+  sunindextype arg2 ;
+  sunindextype result;
+
+  arg1 = (N_Vector)(farg1);
+  arg2 = (sunindextype)(*farg2);
+  result = N_VGetSubvectorLocalLength_MPIManyVector(arg1,arg2);
+  fresult = (sunindextype)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT void _wrap_FN_VLinearSum_MPIManyVector(double const *farg1, N_Vector farg2, double const *farg3, N_Vector farg4, N_Vector farg5) {
   realtype arg1 ;
   N_Vector arg2 = (N_Vector) 0 ;

src/nvector/manyvector/fmod/fnvector_mpimanyvector_mod.f90

@@ -42,6 +42,7 @@ module fnvector_mpimanyvector_mod
  public :: FN_VSpace_MPIManyVector
  public :: FN_VGetCommunicator_MPIManyVector
  public :: FN_VGetLength_MPIManyVector
+ public :: FN_VGetSubvectorLocalLength_MPIManyVector
  public :: FN_VLinearSum_MPIManyVector
  public :: FN_VConst_MPIManyVector
  public :: FN_VProd_MPIManyVector
@@ -220,6 +221,15 @@ function swigc_FN_VGetLength_MPIManyVector(farg1) &
 integer(C_INT64_T) :: fresult
 end function
 
+function swigc_FN_VGetSubvectorLocalLength_MPIManyVector(farg1, farg2) &
+bind(C, name="_wrap_FN_VGetSubvectorLocalLength_MPIManyVector") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+integer(C_INT64_T), intent(in) :: farg2
+integer(C_INT64_T) :: fresult
+end function
+
 subroutine swigc_FN_VLinearSum_MPIManyVector(farg1, farg2, farg3, farg4, farg5) &
 bind(C, name="_wrap_FN_VLinearSum_MPIManyVector")
 use, intrinsic :: ISO_C_BINDING
@@ -764,7 +774,7 @@ function FN_VGetSubvectorArrayPointer_MPIManyVector(v, vec_num) &
 farg1 = c_loc(v)
 farg2 = vec_num
 fresult = swigc_FN_VGetSubvectorArrayPointer_MPIManyVector(farg1, farg2)
-call c_f_pointer(fresult, swig_result, [1])
+call c_f_pointer(fresult, swig_result, [FN_VGetSubvectorLocalLength_MPIManyVector(v, vec_num)])
 end function
 
 function FN_VSetSubvectorArrayPointer_MPIManyVector(v_data, v, vec_num) &
@@ -909,6 +919,22 @@ function FN_VGetLength_MPIManyVector(v) &
 swig_result = fresult
 end function
 
+function FN_VGetSubvectorLocalLength_MPIManyVector(v, vec_num) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT64_T) :: swig_result
+type(N_Vector), target, intent(inout) :: v
+integer(C_INT64_T), intent(in) :: vec_num
+integer(C_INT64_T) :: fresult 
+type(C_PTR) :: farg1 
+integer(C_INT64_T) :: farg2 
+
+farg1 = c_loc(v)
+farg2 = vec_num
+fresult = swigc_FN_VGetSubvectorLocalLength_MPIManyVector(farg1, farg2)
+swig_result = fresult
+end function
+
 subroutine FN_VLinearSum_MPIManyVector(a, x, b, y, z)
 use, intrinsic :: ISO_C_BINDING
 real(C_DOUBLE), intent(in) :: a

src/nvector/manyvector/nvector_manyvector.c

@@ -18,6 +18,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <sundials/sundials_math.h>
+#include "sundials/sundials_nvector.h"
 #ifdef MANYVECTOR_BUILD_WITH_MPI
 #include <nvector/nvector_mpimanyvector.h>
 #else
@@ -271,7 +272,7 @@ N_Vector N_VNew_ManyVector(sunindextype num_subvectors,
 {
   N_Vector v;
   N_VectorContent_ManyVector content;
-  sunindextype i;
+  sunindextype i, local_length;
 
   /* Check that input N_Vectors are non-NULL */
   if (vec_array == NULL)  return(NULL);
@@ -368,15 +369,16 @@ N_Vector N_VNew_ManyVector(sunindextype num_subvectors,
     content->subvec_array[i] = vec_array[i];
 
   /* Determine overall ManyVector length: sum contributions from all subvectors */
-  content->global_length = 0;
+  local_length = 0;
   for (i=0; i<num_subvectors; i++) {
     if (vec_array[i]->ops->nvgetlength) {
-      content->global_length += N_VGetLength(vec_array[i]);
+      local_length += N_VGetLength(vec_array[i]);
     } else {
       N_VDestroy(v);
       return(NULL);
     }
   }
+  content->global_length = local_length;
 
   return(v);
 }
@@ -560,6 +562,10 @@ sunindextype MVAPPEND(N_VGetLength)(N_Vector v)
   return(MANYVECTOR_GLOBLENGTH(v));
 }
 
+sunindextype MVAPPEND(N_VGetSubvectorLocalLength)(N_Vector v, sunindextype vec_num)
+{
+  return(N_VGetLocalLength(MVAPPEND(N_VGetSubvector)(v, vec_num)));
+}
 
 /* Performs the linear sum z = a*x + b*y by calling N_VLinearSum on all subvectors;
    this routine does not check that x, y and z are ManyVectors, if they have the

src/nvector/mpiplusx/fmod/fnvector_mpiplusx_mod.f90

@@ -157,7 +157,7 @@ function FN_VGetArrayPointer_MPIPlusX(v) &
 
 farg1 = c_loc(v)
 fresult = swigc_FN_VGetArrayPointer_MPIPlusX(farg1)
-call c_f_pointer(fresult, swig_result, [1])
+call c_f_pointer(fresult, swig_result, [FN_VGetLocalLength_MPIPlusX(v)])
 end function
 
 subroutine FN_VSetArrayPointer_MPIPlusX(vdata, v)

src/nvector/mpiplusx/nvector_mpiplusx.c

@@ -36,6 +36,7 @@ N_Vector N_VMake_MPIPlusX(MPI_Comm comm, N_Vector X, SUNContext sunctx)
   v->ops->nvgetvectorid = N_VGetVectorID_MPIPlusX;
   v->ops->nvgetarraypointer = N_VGetArrayPointer_MPIPlusX;
   v->ops->nvsetarraypointer = N_VSetArrayPointer_MPIPlusX;
+  v->ops->nvgetlocallength = N_VGetLocalLength_MPIPlusX;
 
   /* debugging functions */
   if (X->ops->nvprint)

src/nvector/openmp/fmod/fnvector_openmp_mod.f90

@@ -767,7 +767,7 @@ function FN_VGetArrayPointer_OpenMP(v) &
 
 farg1 = c_loc(v)
 fresult = swigc_FN_VGetArrayPointer_OpenMP(farg1)
-call c_f_pointer(fresult, swig_result, [1])
+call c_f_pointer(fresult, swig_result, [FN_VGetLength_OpenMP(v)])
 end function
 
 subroutine FN_VSetArrayPointer_OpenMP(v_data, v)

src/nvector/openmp/nvector_openmp.c

@@ -94,6 +94,7 @@ N_Vector N_VNewEmpty_OpenMP(sunindextype length, int num_threads, SUNContext sun
   v->ops->nvgetarraypointer = N_VGetArrayPointer_OpenMP;
   v->ops->nvsetarraypointer = N_VSetArrayPointer_OpenMP;
   v->ops->nvgetlength       = N_VGetLength_OpenMP;
+  v->ops->nvgetlocallength  = N_VGetLength_OpenMP;
 
   /* standard vector operations */
   v->ops->nvlinearsum    = N_VLinearSum_OpenMP;

src/nvector/parallel/fmod/fnvector_parallel_mod.f90

@@ -911,7 +911,7 @@ function FN_VGetArrayPointer_Parallel(v) &
 
 farg1 = c_loc(v)
 fresult = swigc_FN_VGetArrayPointer_Parallel(farg1)
-call c_f_pointer(fresult, swig_result, [1])
+call c_f_pointer(fresult, swig_result, [FN_VGetLocalLength_Parallel(v)])
 end function
 
 subroutine FN_VSetArrayPointer_Parallel(v_data, v)

src/nvector/parallel/nvector_parallel.c

@@ -107,6 +107,7 @@ N_Vector N_VNewEmpty_Parallel(MPI_Comm comm,
   v->ops->nvsetarraypointer = N_VSetArrayPointer_Parallel;
   v->ops->nvgetcommunicator = N_VGetCommunicator_Parallel;
   v->ops->nvgetlength       = N_VGetLength_Parallel;
+  v->ops->nvgetlocallength  = N_VGetLocalLength_Parallel;
 
   /* standard vector operations */
   v->ops->nvlinearsum    = N_VLinearSum_Parallel;

src/nvector/parhyp/nvector_parhyp.c

@@ -23,6 +23,7 @@
 
 #include <nvector/nvector_parhyp.h>
 #include <sundials/sundials_math.h>
+#include "sundials/sundials_nvector.h"
 
 #define ZERO   RCONST(0.0)
 #define HALF   RCONST(0.5)