#include <petsc/private/ftnimpl.h>
#include <petscdmda.h>

#if defined(PETSC_HAVE_FORTRAN_CAPS)
  #define dmdavecgetarray1_         DMDAVECGETARRAY1
  #define dmdavecrestorearray1_     DMDAVECRESTOREARRAY1
  #define dmdavecgetarray2_         DMDAVECGETARRAY2
  #define dmdavecrestorearray2_     DMDAVECRESTOREARRAY2
  #define dmdavecgetarray3_         DMDAVECGETARRAY3
  #define dmdavecrestorearray3_     DMDAVECRESTOREARRAY3
  #define dmdavecgetarray4_         DMDAVECGETARRAY4
  #define dmdavecrestorearray4_     DMDAVECRESTOREARRAY4
  #define dmdavecgetarrayread1_     DMDAVECGETARRAYREAD1
  #define dmdavecrestorearrayread1_ DMDAVECRESTOREARRAYREAD1
  #define dmdavecgetarrayread2_     DMDAVECGETARRAYREAD2
  #define dmdavecrestorearrayread2_ DMDAVECRESTOREARRAYREAD2
  #define dmdavecgetarrayread3_     DMDAVECGETARRAYREAD3
  #define dmdavecrestorearrayread3_ DMDAVECRESTOREARRAYREAD3
  #define dmdavecgetarrayread4_     DMDAVECGETARRAYREAD4
  #define dmdavecrestorearrayread4_ DMDAVECRESTOREARRAYREAD4
  #define dmdagetelements_          DMDAGETELEMENTS
  #define dmdarestoreelements_      DMDARESTOREELEMENTS
#elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
  #define dmdavecgetarray1_         dmdavecgetarray1
  #define dmdavecrestorearray1_     dmdavecrestorearray1
  #define dmdavecgetarray2_         dmdavecgetarray2
  #define dmdavecrestorearray2_     dmdavecrestorearray2
  #define dmdavecgetarray3_         dmdavecgetarray3
  #define dmdavecrestorearray3_     dmdavecrestorearray3
  #define dmdavecgetarray4_         dmdavecgetarray4
  #define dmdavecrestorearray4_     dmdavecrestorearray4
  #define dmdavecgetarrayread1_     dmdavecgetarrayread1
  #define dmdavecrestorearrayread1_ dmdavecrestorearrayread1
  #define dmdavecgetarrayread2_     dmdavecgetarrayread2
  #define dmdavecrestorearrayread2_ dmdavecrestorearrayread2
  #define dmdavecgetarrayread3_     dmdavecgetarrayread3
  #define dmdavecrestorearrayread3_ dmdavecrestorearrayread3
  #define dmdavecgetarrayread4_     dmdavecgetarrayread4
  #define dmdavecrestorearrayread4_ dmdavecrestorearrayread4
  #define dmdagetelements_          dmdagetelements
  #define dmdarestoreelements_      dmdarestoreelements
#endif

PETSC_EXTERN void dmdagetelements_(DM *dm, PetscInt *nel, PetscInt *nen, F90Array1d *e, int *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  const PetscInt *fa;

  if (!e) {
    *ierr = PetscError(((PetscObject)e)->comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__, PETSC_ERR_ARG_BADPTR, PETSC_ERROR_INITIAL, "e==NULL, maybe #include <petsc/finclude/petscvec.h> is missing?");
    return;
  }
  *ierr = DMDAGetElements(*dm, nel, nen, &fa);
  if (*ierr) return;
  *ierr = F90Array1dCreate((PetscInt *)fa, MPIU_INT, 1, (*nel) * (*nen), e PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdarestoreelements_(DM *dm, PetscInt *nel, PetscInt *nen, F90Array1d *e, int *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  if (!e) {
    *ierr = PetscError(((PetscObject)e)->comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__, PETSC_ERR_ARG_BADPTR, PETSC_ERROR_INITIAL, "e==NULL, maybe #include <petsc/finclude/petscvec.h> is missing?");
    return;
  }
  *ierr = F90Array1dDestroy(e, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarray1_(DM *da, Vec *v, F90Array1d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt     xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof;
  PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  *ierr = VecGetArray(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array1dCreate(aa, MPIU_SCALAR, gxs, gxm, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearray1_(DM *da, Vec *v, F90Array1d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscScalar *fa;
  *ierr = F90Array1dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArray(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array1dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarray2_(DM *da, Vec *v, F90Array2d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt     xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof;
  PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  if (dim == 1) {
    gys = gxs;
    gym = gxm;
    gxs = 0;
    gxm = dof;
  }
  *ierr = VecGetArray(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array2dCreate(aa, MPIU_SCALAR, gxs, gxm, gys, gym, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearray2_(DM *da, Vec *v, F90Array2d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscScalar *fa;
  *ierr = F90Array2dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArray(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array2dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarray3_(DM *da, Vec *v, F90Array3d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt     xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof;
  PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  if (dim == 2) {
    gzs = gys;
    gzm = gym;
    gys = gxs;
    gym = gxm;
    gxs = 0;
    gxm = dof;
  }
  *ierr = VecGetArray(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array3dCreate(aa, MPIU_SCALAR, gxs, gxm, gys, gym, gzs, gzm, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearray3_(DM *da, Vec *v, F90Array3d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscScalar *fa;
  *ierr = F90Array3dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArray(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array3dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarray4_(DM *da, Vec *v, F90Array4d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt     xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof, zero = 0;
  PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  *ierr = VecGetArray(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array4dCreate(aa, MPIU_SCALAR, zero, dof, gxs, gxm, gys, gym, gzs, gzm, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearray4_(DM *da, Vec *v, F90Array4d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscScalar *fa;
  /*
    F90Array4dAccess is not implemented, so the following call would fail
  */
  *ierr = F90Array4dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArray(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array4dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarrayread1_(DM *da, Vec *v, F90Array1d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt           xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof;
  const PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  *ierr = VecGetArrayRead(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array1dCreate((void *)aa, MPIU_SCALAR, gxs, gxm, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearrayread1_(DM *da, Vec *v, F90Array1d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  const PetscScalar *fa;
  *ierr = F90Array1dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArrayRead(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array1dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarrayread2_(DM *da, Vec *v, F90Array2d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt           xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof;
  const PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  if (dim == 1) {
    gys = gxs;
    gym = gxm;
    gxs = 0;
    gxm = dof;
  }
  *ierr = VecGetArrayRead(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array2dCreate((void *)aa, MPIU_SCALAR, gxs, gxm, gys, gym, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearrayread2_(DM *da, Vec *v, F90Array2d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  const PetscScalar *fa;
  *ierr = F90Array2dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArrayRead(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array2dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarrayread3_(DM *da, Vec *v, F90Array3d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt           xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof;
  const PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  if (dim == 2) {
    gzs = gys;
    gzm = gym;
    gys = gxs;
    gym = gxm;
    gxs = 0;
    gxm = dof;
  }
  *ierr = VecGetArrayRead(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array3dCreate((void *)aa, MPIU_SCALAR, gxs, gxm, gys, gym, gzs, gzm, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearrayread3_(DM *da, Vec *v, F90Array3d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  const PetscScalar *fa;
  *ierr = F90Array3dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArrayRead(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array3dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}

PETSC_EXTERN void dmdavecgetarrayread4_(DM *da, Vec *v, F90Array4d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  PetscInt           xs, ys, zs, xm, ym, zm, gxs, gys, gzs, gxm, gym, gzm, N, dim, dof, zero = 0;
  const PetscScalar *aa;

  *ierr = DMDAGetCorners(*da, &xs, &ys, &zs, &xm, &ym, &zm);
  if (*ierr) return;
  *ierr = DMDAGetGhostCorners(*da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  if (*ierr) return;
  *ierr = DMDAGetInfo(*da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
  if (*ierr) return;

  /* Handle case where user passes in global vector as opposed to local */
  *ierr = VecGetLocalSize(*v, &N);
  if (*ierr) return;
  if (N == xm * ym * zm * dof) {
    gxm = xm;
    gym = ym;
    gzm = zm;
    gxs = xs;
    gys = ys;
    gzs = zs;
  } else if (N != gxm * gym * gzm * dof) {
    *ierr = PETSC_ERR_ARG_INCOMP;
    return;
  }
  *ierr = VecGetArrayRead(*v, &aa);
  if (*ierr) return;
  *ierr = F90Array4dCreate((void *)aa, MPIU_SCALAR, zero, dof, gxs, gxm, gys, gym, gzs, gzm, a PETSC_F90_2PTR_PARAM(ptrd));
  if (*ierr) return;
}

PETSC_EXTERN void dmdavecrestorearrayread4_(DM *da, Vec *v, F90Array4d *a, PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptrd))
{
  const PetscScalar *fa;
  /*
    F90Array4dAccess is not implemented, so the following call would fail
  */
  *ierr = F90Array4dAccess(a, MPIU_SCALAR, (void **)&fa PETSC_F90_2PTR_PARAM(ptrd));
  *ierr = VecRestoreArrayRead(*v, &fa);
  if (*ierr) return;
  *ierr = F90Array4dDestroy(a, MPIU_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
}