xref: /petsc/src/mat/interface/ftn-custom/zmatrixf.c (revision 5d83a8b16d06840f96948f1a43aa9c83c769a60a)

#include <petsc/private/fortranimpl.h>
#include <petsc/private/f90impl.h>
#include <petscmat.h>
#include <petscviewer.h>

#if defined(PETSC_HAVE_FORTRAN_CAPS)
  #define matdestroymatrices_       MATDESTROYMATRICES
  #define matdestroysubmatrices_    MATDESTROYSUBMATRICES
  #define matgetrowij_              MATGETROWIJ
  #define matrestorerowij_          MATRESTOREROWIJ
  #define matgetrow_                MATGETROW
  #define matrestorerow_            MATRESTOREROW
  #define matseqaijgetarray_        MATSEQAIJGETARRAY
  #define matseqaijrestorearray_    MATSEQAIJRESTOREARRAY
  #define matdensegetarray_         MATDENSEGETARRAY
  #define matdensegetarrayread_     MATDENSEGETARRAYREAD
  #define matdenserestorearray_     MATDENSERESTOREARRAY
  #define matdenserestorearrayread_ MATDENSERESTOREARRAYREAD
  #define matcreatesubmatrices_     MATCREATESUBMATRICES
  #define matcreatesubmatricesmpi_  MATCREATESUBMATRICESMPI
  #define matnullspacesetfunction_  MATNULLSPACESETFUNCTION
  #define matfindnonzerorows_       MATFINDNONZEROROWS
#elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
  #define matdestroymatrices_       matdestroymatrices
  #define matdestroysubmatrices_    matdestroysubmatrices
  #define matgetrowij_              matgetrowij
  #define matrestorerowij_          matrestorerowij
  #define matgetrow_                matgetrow
  #define matrestorerow_            matrestorerow
  #define matseqaijgetarray_        matseqaijgetarray
  #define matseqaijrestorearray_    matseqaijrestorearray
  #define matdensegetarray_         matdensegetarray
  #define matdensegetarrayread_     matdensegetarrayread
  #define matdenserestorearray_     matdenserestorearray
  #define matdenserestorearrayread_ matdenserestorearrayread
  #define matcreatesubmatrices_     matcreatesubmatrices
  #define matcreatesubmatricesmpi_  matcreatesubmatricesmpi
  #define matnullspacesetfunction_  matnullspacesetfunction
  #define matfindnonzerorows_       matfindnonzerorows
#endif

static PetscErrorCode ournullfunction(MatNullSpace sp, Vec x, void *ctx)
{
  PetscCallFortranVoidFunction((*(void (*)(MatNullSpace *, Vec *, void *, PetscErrorCode *))(((PetscObject)sp)->fortran_func_pointers[0]))(&sp, &x, ctx, &ierr));
  return PETSC_SUCCESS;
}

PETSC_EXTERN void matnullspacesetfunction_(MatNullSpace *sp, PetscErrorCode (*rem)(MatNullSpace, Vec, void *), void *ctx, PetscErrorCode *ierr)
{
  PetscObjectAllocateFortranPointers(*sp, 1);
  ((PetscObject)*sp)->fortran_func_pointers[0] = (PetscVoidFn *)rem;

  *ierr = MatNullSpaceSetFunction(*sp, ournullfunction, ctx);
}

PETSC_EXTERN void matgetrowij_(Mat *B, PetscInt *shift, PetscBool *sym, PetscBool *blockcompressed, PetscInt *n, PetscInt *ia, size_t *iia, PetscInt *ja, size_t *jja, PetscBool *done, PetscErrorCode *ierr)
{
  const PetscInt *IA, *JA;
  *ierr = MatGetRowIJ(*B, *shift, *sym, *blockcompressed, n, &IA, &JA, done);
  if (*ierr) return;
  *iia = PetscIntAddressToFortran(ia, (PetscInt *)IA);
  *jja = PetscIntAddressToFortran(ja, (PetscInt *)JA);
}

PETSC_EXTERN void matrestorerowij_(Mat *B, PetscInt *shift, PetscBool *sym, PetscBool *blockcompressed, PetscInt *n, PetscInt *ia, size_t *iia, PetscInt *ja, size_t *jja, PetscBool *done, PetscErrorCode *ierr)
{
  const PetscInt *IA = PetscIntAddressFromFortran(ia, *iia), *JA = PetscIntAddressFromFortran(ja, *jja);
  *ierr = MatRestoreRowIJ(*B, *shift, *sym, *blockcompressed, n, &IA, &JA, done);
}

/*
   This is a poor way of storing the column and value pointers
  generated by MatGetRow() to be returned with MatRestoreRow()
  but there is not natural,good place else to store them. Hence
  Fortran programmers can only have one outstanding MatGetRows()
  at a time.
*/
static int                matgetrowactive = 0;
static const PetscInt    *my_ocols        = NULL;
static const PetscScalar *my_ovals        = NULL;

PETSC_EXTERN void matgetrow_(Mat *mat, PetscInt *row, PetscInt *ncols, PetscInt *cols, PetscScalar *vals, PetscErrorCode *ierr)
{
  const PetscInt    **oocols = &my_ocols;
  const PetscScalar **oovals = &my_ovals;

  if (matgetrowactive) {
    *ierr = PetscError(PETSC_COMM_SELF, __LINE__, "MatGetRow_Fortran", __FILE__, PETSC_ERR_ARG_WRONGSTATE, PETSC_ERROR_INITIAL, "Cannot have two MatGetRow() active simultaneously\n\
               call MatRestoreRow() before calling MatGetRow() a second time");
    *ierr = PETSC_ERR_ARG_WRONGSTATE;
    return;
  }

  CHKFORTRANNULLINTEGER(cols);
  if (!cols) oocols = NULL;
  CHKFORTRANNULLSCALAR(vals);
  if (!vals) oovals = NULL;

  *ierr = MatGetRow(*mat, *row, ncols, oocols, oovals);
  if (*ierr) return;

  if (oocols) {
    *ierr = PetscArraycpy(cols, my_ocols, *ncols);
    if (*ierr) return;
  }
  if (oovals) {
    *ierr = PetscArraycpy(vals, my_ovals, *ncols);
    if (*ierr) return;
  }
  matgetrowactive = 1;
}

PETSC_EXTERN void matrestorerow_(Mat *mat, PetscInt *row, PetscInt *ncols, PetscInt *cols, PetscScalar *vals, PetscErrorCode *ierr)
{
  const PetscInt    **oocols = &my_ocols;
  const PetscScalar **oovals = &my_ovals;

  if (!matgetrowactive) {
    *ierr = PetscError(PETSC_COMM_SELF, __LINE__, "MatRestoreRow_Fortran", __FILE__, PETSC_ERR_ARG_WRONGSTATE, PETSC_ERROR_INITIAL, "Must call MatGetRow() first");
    *ierr = PETSC_ERR_ARG_WRONGSTATE;
    return;
  }
  CHKFORTRANNULLINTEGER(cols);
  if (!cols) oocols = NULL;
  CHKFORTRANNULLSCALAR(vals);
  if (!vals) oovals = NULL;

  *ierr           = MatRestoreRow(*mat, *row, ncols, oocols, oovals);
  matgetrowactive = 0;
}

PETSC_EXTERN void matseqaijgetarray_(Mat *mat, PetscScalar *fa, size_t *ia, PetscErrorCode *ierr)
{
  PetscScalar *mm;
  PetscInt     m, n;

  *ierr = MatSeqAIJGetArray(*mat, &mm);
  if (*ierr) return;
  *ierr = MatGetSize(*mat, &m, &n);
  if (*ierr) return;
  *ierr = PetscScalarAddressToFortran((PetscObject)*mat, 1, fa, mm, m * n, ia);
  if (*ierr) return;
}

PETSC_EXTERN void matseqaijrestorearray_(Mat *mat, PetscScalar *fa, size_t *ia, PetscErrorCode *ierr)
{
  PetscScalar *lx;
  PetscInt     m, n;

  *ierr = MatGetSize(*mat, &m, &n);
  if (*ierr) return;
  *ierr = PetscScalarAddressFromFortran((PetscObject)*mat, fa, *ia, m * n, &lx);
  if (*ierr) return;
  *ierr = MatSeqAIJRestoreArray(*mat, &lx);
  if (*ierr) return;
}

PETSC_EXTERN void matdensegetarray_(Mat *mat, PetscScalar *fa, size_t *ia, PetscErrorCode *ierr)
{
  PetscScalar *mm;
  PetscInt     m, n;

  *ierr = MatDenseGetArray(*mat, &mm);
  if (*ierr) return;
  *ierr = MatGetSize(*mat, &m, &n);
  if (*ierr) return;
  *ierr = PetscScalarAddressToFortran((PetscObject)*mat, 1, fa, mm, m * n, ia);
  if (*ierr) return;
}

PETSC_EXTERN void matdenserestorearray_(Mat *mat, PetscScalar *fa, size_t *ia, PetscErrorCode *ierr)
{
  PetscScalar *lx;
  PetscInt     m, n;

  *ierr = MatGetSize(*mat, &m, &n);
  if (*ierr) return;
  *ierr = PetscScalarAddressFromFortran((PetscObject)*mat, fa, *ia, m * n, &lx);
  if (*ierr) return;
  *ierr = MatDenseRestoreArray(*mat, &lx);
  if (*ierr) return;
}

PETSC_EXTERN void matdensegetarrayread_(Mat *mat, PetscScalar *fa, size_t *ia, PetscErrorCode *ierr)
{
  const PetscScalar *mm;
  PetscInt           m, n;

  *ierr = MatDenseGetArrayRead(*mat, &mm);
  if (*ierr) return;
  *ierr = MatGetSize(*mat, &m, &n);
  if (*ierr) return;
  *ierr = PetscScalarAddressToFortran((PetscObject)*mat, 1, fa, (PetscScalar *)mm, m * n, ia);
  if (*ierr) return;
}

PETSC_EXTERN void matdenserestorearrayread_(Mat *mat, PetscScalar *fa, size_t *ia, PetscErrorCode *ierr)
{
  const PetscScalar *lx;
  PetscInt           m, n;

  *ierr = MatGetSize(*mat, &m, &n);
  if (*ierr) return;
  *ierr = PetscScalarAddressFromFortran((PetscObject)*mat, fa, *ia, m * n, (PetscScalar **)&lx);
  if (*ierr) return;
  *ierr = MatDenseRestoreArrayRead(*mat, &lx);
  if (*ierr) return;
}

/*
    MatCreateSubmatrices() is slightly different from C since the
    Fortran provides the array to hold the submatrix objects,while in C that
    array is allocated by the MatCreateSubmatrices()
*/
PETSC_EXTERN void matcreatesubmatrices_(Mat *mat, PetscInt *n, IS *isrow, IS *iscol, MatReuse *scall, Mat *smat, PetscErrorCode *ierr)
{
  Mat     *lsmat;
  PetscInt i;

  if (*scall == MAT_INITIAL_MATRIX) {
    *ierr = MatCreateSubMatrices(*mat, *n, isrow, iscol, *scall, &lsmat);
    for (i = 0; i <= *n; i++) { /* lsmat[*n] might be a dummy matrix for saving data structure */
      smat[i] = lsmat[i];
    }
    *ierr = PetscFree(lsmat);
  } else {
    *ierr = MatCreateSubMatrices(*mat, *n, isrow, iscol, *scall, &smat);
  }
}

/*
    MatCreateSubmatrices() is slightly different from C since the
    Fortran provides the array to hold the submatrix objects,while in C that
    array is allocated by the MatCreateSubmatrices()
*/
PETSC_EXTERN void matcreatesubmatricesmpi_(Mat *mat, PetscInt *n, IS *isrow, IS *iscol, MatReuse *scall, Mat *smat, PetscErrorCode *ierr)
{
  Mat     *lsmat;
  PetscInt i;

  if (*scall == MAT_INITIAL_MATRIX) {
    *ierr = MatCreateSubMatricesMPI(*mat, *n, isrow, iscol, *scall, &lsmat);
    for (i = 0; i <= *n; i++) { /* lsmat[*n] might be a dummy matrix for saving data structure */
      smat[i] = lsmat[i];
    }
    *ierr = PetscFree(lsmat);
  } else {
    *ierr = MatCreateSubMatricesMPI(*mat, *n, isrow, iscol, *scall, &smat);
  }
}

/*
    MatDestroyMatrices() is slightly different from C since the
    Fortran does not free the array of matrix objects, while in C that
    the array is freed
*/
PETSC_EXTERN void matdestroymatrices_(PetscInt *n, Mat *smat, PetscErrorCode *ierr)
{
  PetscInt i;

  for (i = 0; i < *n; i++) {
    PETSC_FORTRAN_OBJECT_F_DESTROYED_TO_C_NULL(&smat[i]);
    *ierr = MatDestroy(&smat[i]);
    if (*ierr) return;
    PETSC_FORTRAN_OBJECT_C_NULL_TO_F_DESTROYED(&smat[i]);
  }
}

/*
    MatDestroySubMatrices() is slightly different from C since the
    Fortran provides the array to hold the submatrix objects, while in C that
    array is allocated by the MatCreateSubmatrices()

    An extra matrix may be stored at the end of the array, hence the check see
    MatDestroySubMatrices_Dummy()
*/
PETSC_EXTERN void matdestroysubmatrices_(PetscInt *n, Mat *smat, PetscErrorCode *ierr)
{
  Mat     *lsmat;
  PetscInt i;

  if (*n == 0) return;
  *ierr = PetscMalloc1(*n + 1, &lsmat);
  if (*ierr) return;
  for (i = 0; i <= *n; i++) { lsmat[i] = smat[i]; }
  *ierr = MatDestroySubMatrices(*n, &lsmat);
  if (*ierr) return;
  for (i = 0; i <= *n; i++) { PETSC_FORTRAN_OBJECT_C_NULL_TO_F_DESTROYED(&smat[i]); }
}