impls/matlab/vmatlab.c

#include <petsc/private/viewerimpl.h>
#include <mat.h>

typedef struct {
  MATFile      *ep;
  PetscMPIInt   rank;
  PetscFileMode btype;
} PetscViewer_Matlab;

/*@C
    PetscViewerMatlabPutArray - Puts an array into the `PETSCVIEWERMATLAB` viewer.

      Not collective: only processor zero saves the array

    Input Parameters:
+    mfile - the viewer
.    m,n - the dimensions of the array
.    array - the array (represented in one dimension)
-    name - the name of the array

   Level: advanced

    Note:
    Only writes array values on processor 0.

.seealso: `PETSCVIEWERMATLAB`, `PetscViewerMatlabGetArray()`
@*/
PetscErrorCode PetscViewerMatlabPutArray(PetscViewer mfile, int m, int n, const PetscScalar *array, const char *name) {
  PetscViewer_Matlab *ml;
  mxArray            *mat;

  PetscFunctionBegin;
  PetscCheck(mfile, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null argument: probably PETSC_VIEWER_MATLAB_() failed");
  ml = (PetscViewer_Matlab *)mfile->data;
  if (!ml->rank) {
    PetscCall(PetscInfo(mfile, "Putting MATLAB array %s\n", name));
#if !defined(PETSC_USE_COMPLEX)
    mat = mxCreateDoubleMatrix(m, n, mxREAL);
#else
    mat = mxCreateDoubleMatrix(m, n, mxCOMPLEX);
#endif
    PetscCall(PetscArraycpy(mxGetPr(mat), array, m * n));
    matPutVariable(ml->ep, name, mat);

    PetscCall(PetscInfo(mfile, "Put MATLAB array %s\n", name));
  }
  PetscFunctionReturn(0);
}

PetscErrorCode PetscViewerMatlabPutVariable(PetscViewer viewer, const char *name, void *mat) {
  PetscViewer_Matlab *ml = (PetscViewer_Matlab *)viewer->data;

  PetscFunctionBegin;
  matPutVariable(ml->ep, name, (mxArray *)mat);
  PetscFunctionReturn(0);
}

/*@C
    PetscViewerMatlabGetArray - Gets a variable from a `PETSCVIEWERMATLAB` viewer into an array

    Not Collective; only processor zero reads in the array

    Input Parameters:
+    mfile - the MATLAB file viewer
.    m,n - the dimensions of the array
.    array - the array (represented in one dimension)
-    name - the name of the array

   Level: advanced

    Note:
    Only reads in array values on processor 0.

.seealso: `PETSCVIEWERMATLAB`, `PetscViewerMatlabPutArray()`
@*/
PetscErrorCode PetscViewerMatlabGetArray(PetscViewer mfile, int m, int n, PetscScalar *array, const char *name) {
  PetscViewer_Matlab *ml;
  mxArray            *mat;

  PetscFunctionBegin;
  PetscCheck(mfile, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null argument: probably PETSC_VIEWER_MATLAB_() failed");
  ml = (PetscViewer_Matlab *)mfile->data;
  if (!ml->rank) {
    PetscCall(PetscInfo(mfile, "Getting MATLAB array %s\n", name));
    mat = matGetVariable(ml->ep, name);
    PetscCheck(mat, PETSC_COMM_SELF, PETSC_ERR_LIB, "Unable to get array %s from matlab", name);
    PetscCall(PetscArraycpy(array, mxGetPr(mat), m * n));
    PetscCall(PetscInfo(mfile, "Got MATLAB array %s\n", name));
  }
  PetscFunctionReturn(0);
}

PetscErrorCode PetscViewerFileSetMode_Matlab(PetscViewer viewer, PetscFileMode type) {
  PetscViewer_Matlab *vmatlab = (PetscViewer_Matlab *)viewer->data;

  PetscFunctionBegin;
  vmatlab->btype = type;
  PetscFunctionReturn(0);
}

/*
        Actually opens the file
*/
PetscErrorCode PetscViewerFileSetName_Matlab(PetscViewer viewer, const char name[]) {
  PetscViewer_Matlab *vmatlab = (PetscViewer_Matlab *)viewer->data;
  PetscFileMode       type    = vmatlab->btype;

  PetscFunctionBegin;
  PetscCheck(type != (PetscFileMode)-1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
  if (vmatlab->ep) matClose(vmatlab->ep);

  /* only first processor opens file */
  if (!vmatlab->rank) {
    if (type == FILE_MODE_READ) vmatlab->ep = matOpen(name, "r");
    else if (type == FILE_MODE_WRITE) vmatlab->ep = matOpen(name, "w");
    else {
      PetscCheck(type != FILE_MODE_UNDEFINED, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
      SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Unsupported file mode %s", PetscFileModes[type]);
    }
  }
  PetscFunctionReturn(0);
}

PetscErrorCode PetscViewerDestroy_Matlab(PetscViewer v) {
  PetscViewer_Matlab *vf = (PetscViewer_Matlab *)v->data;

  PetscFunctionBegin;
  if (vf->ep) matClose(vf->ep);
  PetscCall(PetscFree(vf));
  PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetName_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscViewerFileSetMode_C", NULL));
  PetscFunctionReturn(0);
}

/*MC
   PETSCVIEWERMATLAB - A viewer that saves the variables into a MATLAB .mat file that may be read into MATLAB
       with load('filename').

   Level: intermediate

       Notes:
             Currently can only save PETSc vectors to .mat files, not matrices (use the `PETSCVIEWERBINARY` and
             ${PETSC_DIR}/share/petsc/matlab/PetscBinaryRead.m to read matrices into MATLAB).

             For parallel vectors obtained with `DMCreateGlobalVector()` or `DMGetGlobalVector()` the vectors are saved to
             the .mat file in natural ordering. You can use DMView() to save the `DMDA` information to the .mat file
             the fields in the MATLAB loaded da variable give the array dimensions so you can reshape the MATLAB
             vector to the same multidimensional shape as it had in PETSc for plotting etc. For example,

             In your PETSc C/C++ code (assuming a two dimensional `DMDA` with one degree of freedom per node)
.vb
                PetscObjectSetName((PetscObject)x,"x");
                VecView(x,PETSC_VIEWER_MATLAB_WORLD);
                PetscObjectSetName((PetscObject)da,"da");
                DMView(x,PETSC_VIEWER_MATLAB_WORLD);
.ve
             Then from MATLAB
.vb
                load('matlaboutput.mat')   % matlaboutput.mat is the default filename
                xnew = zeros(da.n,da.m);
                xnew(:) = x;    % reshape one dimensional vector back to two dimensions
.ve

              If you wish to put the same variable into the .mat file several times you need to give it a new
              name before each call to view.

              Use `PetscViewerMatlabPutArray()` to just put an array of doubles into the .mat file

.seealso: `PETSC_VIEWER_MATLAB_()`, `PETSC_VIEWER_MATLAB_SELF`, `PETSC_VIEWER_MATLAB_WORLD`, `PetscViewerCreate()`,
          `PetscViewerMatlabOpen()`, `VecView()`, `DMView()`, `PetscViewerMatlabPutArray()`, `PETSCVIEWERBINARY`, `PETSCVIEWERASCII`, `PETSCVIEWERDRAW`,
          `PETSC_VIEWER_STDOUT_()`, `PetscViewerFileSetName()`, `PetscViewerFileSetMode()`, `PetscViewerFormat`, `PetscMatlabEngine`
M*/
PETSC_EXTERN PetscErrorCode PetscViewerCreate_Matlab(PetscViewer viewer) {
  PetscViewer_Matlab *e;

  PetscFunctionBegin;
  PetscCall(PetscNew(&e));
  PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &e->rank));
  e->btype     = FILE_MODE_UNDEFINED;
  viewer->data = (void *)e;

  PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetName_C", PetscViewerFileSetName_Matlab));
  PetscCall(PetscObjectComposeFunction((PetscObject)viewer, "PetscViewerFileSetMode_C", PetscViewerFileSetMode_Matlab));

  viewer->ops->destroy = PetscViewerDestroy_Matlab;
  PetscFunctionReturn(0);
}

/*@C
   PetscViewerMatlabOpen - Opens a Matlab .mat file for output

   Collective

   Input Parameters:
+  comm - MPI communicator
.  name - name of file
-  type - type of file
$    `FILE_MODE_WRITE` - create new file for MATLAB output
$    `FILE_MODE_READ` - open existing file for MATLAB input
$    `FILE_MODE_WRITE` - open existing file for MATLAB output

   Output Parameter:
.  binv - PetscViewer for MATLAB output to use with the specified file

   Level: beginner

   Notes:
   This `PetscViewer` should be destroyed with `PetscViewerDestroy()`.

   For writing files it only opens the file on processor 0 in the communicator.

   This only saves `Vec`s it cannot be used to save `Mat`s. We recommend using the `PETSCVIEWERBINARY` to save objects to be loaded into MATLAB
   instead of this routine.

   PETSc must be configured with the option -with-matlab for this functionality

.seealso: `PETSCVIEWERMATLAB`, `PetscViewerASCIIOpen()`, `PetscViewerPushFormat()`, `PetscViewerDestroy()`, `PETSCVIEWERBINARY`, `PetscViewerBinaryOpen()`
          `VecView()`, `MatView()`, `VecLoad()`, `MatLoad()`
@*/
PetscErrorCode PetscViewerMatlabOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *binv) {
  PetscFunctionBegin;
  PetscCall(PetscViewerCreate(comm, binv));
  PetscCall(PetscViewerSetType(*binv, PETSCVIEWERMATLAB));
  PetscCall(PetscViewerFileSetMode(*binv, type));
  PetscCall(PetscViewerFileSetName(*binv, name));
  PetscFunctionReturn(0);
}

static PetscMPIInt Petsc_Viewer_Matlab_keyval = MPI_KEYVAL_INVALID;

/*@C
     PETSC_VIEWER_MATLAB_ - Creates a `PETSCVIEWERMATLAB` `PetscViewer` shared by all processors
                     in a communicator.

     Collective

     Input Parameter:
.    comm - the MPI communicator to share the Matlab `PetscViewer`

   Options Database Key:
.    -viewer_matlab_filename <name> - name of the Matlab file

   Environmental variable:
.   `PETSC_VIEWER_MATLAB_FILENAME` - name of the Matlab file

     Level: intermediate

     Note:
     Unlike almost all other PETSc routines, `PETSC_VIEWER_MATLAB_()` does not return
     an error code.  The matlab PetscViewer is usually used in the form
$       XXXView(XXX object,PETSC_VIEWER_MATLAB_(comm));

     Use `PETSC_VIEWER_SOCKET_()` or `PetscViewerSocketOpen()` to communicator with an interactive MATLAB session.

.seealso: `PETSC_VIEWER_MATLAB_WORLD`, `PETSC_VIEWER_MATLAB_SELF`, `PetscViewerMatlabOpen()`, `PetscViewerCreate()`,
          `PetscViewerDestroy()`
@*/
PetscViewer PETSC_VIEWER_MATLAB_(MPI_Comm comm) {
  PetscErrorCode ierr;
  PetscBool      flg;
  PetscViewer    viewer;
  char           fname[PETSC_MAX_PATH_LEN];
  MPI_Comm       ncomm;

  PetscFunctionBegin;
  ierr = PetscCommDuplicate(comm, &ncomm, NULL);
  if (ierr) {
    PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
    PetscFunctionReturn(0);
  }
  if (Petsc_Viewer_Matlab_keyval == MPI_KEYVAL_INVALID) {
    ierr = MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, MPI_COMM_NULL_DELETE_FN, &Petsc_Viewer_Matlab_keyval, 0);
    if (ierr) {
      PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
      PetscFunctionReturn(NULL);
    }
  }
  ierr = MPI_Comm_get_attr(ncomm, Petsc_Viewer_Matlab_keyval, (void **)&viewer, (int *)&flg);
  if (ierr) {
    PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
    PetscFunctionReturn(NULL);
  }
  if (!flg) { /* PetscViewer not yet created */
    ierr = PetscOptionsGetenv(ncomm, "PETSC_VIEWER_MATLAB_FILENAME", fname, PETSC_MAX_PATH_LEN, &flg);
    if (ierr) {
      PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_REPEAT, " ");
      PetscFunctionReturn(NULL);
    }
    if (!flg) {
      ierr = PetscStrcpy(fname, "matlaboutput.mat");
      if (ierr) {
        PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_REPEAT, " ");
        PetscFunctionReturn(NULL);
      }
    }
    ierr = PetscViewerMatlabOpen(ncomm, fname, FILE_MODE_WRITE, &viewer);
    if (ierr) {
      PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_REPEAT, " ");
      PetscFunctionReturn(NULL);
    }
    ierr = PetscObjectRegisterDestroy((PetscObject)viewer);
    if (ierr) {
      PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_REPEAT, " ");
      PetscFunctionReturn(NULL);
    }
    ierr = MPI_Comm_set_attr(ncomm, Petsc_Viewer_Matlab_keyval, (void *)viewer);
    if (ierr) {
      PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_INITIAL, " ");
      PetscFunctionReturn(NULL);
    }
  }
  ierr = PetscCommDestroy(&ncomm);
  if (ierr) {
    PetscError(PETSC_COMM_SELF, __LINE__, "PETSC_VIEWER_MATLAB_", __FILE__, PETSC_ERR_PLIB, PETSC_ERROR_REPEAT, " ");
    PetscFunctionReturn(NULL);
  }
  PetscFunctionReturn(viewer);
}