Actual source code: mathematica.c

      1: #include <petsc/private/viewerimpl.h>
      2: #include <petsc/private/pcimpl.h>
      3: #include <../src/mat/impls/aij/seq/aij.h>
      4: #include <mathematica.h>

      6: #if defined(PETSC_HAVE__SNPRINTF) && !defined(PETSC_HAVE_SNPRINTF)
      7:   #define snprintf _snprintf
      8: #endif

     10: PetscViewer  PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE = NULL;
     11: static void *mathematicaEnv                         = NULL;

     13: static PetscBool PetscViewerMathematicaPackageInitialized = PETSC_FALSE;
     14: /*@C
     15:   PetscViewerMathematicaFinalizePackage - This function destroys everything in the Petsc interface to Mathematica. It is
     16:   called from PetscFinalize().

     18:   Level: developer

     20: .seealso: `PetscFinalize()`
     21: @*/
     22: PetscErrorCode PetscViewerMathematicaFinalizePackage(void)
     23: {
     24:   PetscFunctionBegin;
     25:   if (mathematicaEnv) MLDeinitialize((MLEnvironment)mathematicaEnv);
     26:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
     27:   PetscFunctionReturn(PETSC_SUCCESS);
     28: }

     30: /*@C
     31:   PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
     32:   called from `PetscViewerInitializePackage()`.

     34:   Level: developer

     36: .seealso: `PetscSysInitializePackage()`, `PetscInitialize()`
     37: @*/
     38: PetscErrorCode PetscViewerMathematicaInitializePackage(void)
     39: {
     40:   PetscError ierr;

     42:   PetscFunctionBegin;
     43:   if (PetscViewerMathematicaPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
     44:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;

     46:   mathematicaEnv = (void *)MLInitialize(0);

     48:   PetscCall(PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage));
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: static PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
     53: {
     54:   PetscFunctionBegin;
     55:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) PetscFunctionReturn(PETSC_SUCCESS);
     56:   PetscCall(PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE));
     57:   PetscFunctionReturn(PETSC_SUCCESS);
     58: }

     60: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
     61: {
     62:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;

     64:   PetscFunctionBegin;
     65:   MLClose(vmath->link);
     66:   PetscCall(PetscFree(vmath->linkname));
     67:   PetscCall(PetscFree(vmath->linkhost));
     68:   PetscCall(PetscFree(vmath));
     69:   PetscFunctionReturn(PETSC_SUCCESS);
     70: }

     72: static PetscErrorCode PetscViewerDestroyMathematica_Private(void)
     73: {
     74:   PetscFunctionBegin;
     75:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) PetscCall(PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE));
     76:   PetscFunctionReturn(PETSC_SUCCESS);
     77: }

     79: static PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
     80: {
     81:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
     82: #if defined(MATHEMATICA_3_0)
     83:   int   argc = 6;
     84:   char *argv[6];
     85: #else
     86:   int   argc = 5;
     87:   char *argv[5];
     88: #endif
     89:   char hostname[256];
     90:   long lerr;

     92:   PetscFunctionBegin;
     93:   /* Link name */
     94:   argv[0] = "-linkname";
     95:   if (!vmath->linkname) argv[1] = "math -mathlink";
     96:   else argv[1] = vmath->linkname;

     98:   /* Link host */
     99:   argv[2] = "-linkhost";
    100:   if (!vmath->linkhost) {
    101:     PetscCall(PetscGetHostName(hostname, sizeof(hostname)));
    102:     argv[3] = hostname;
    103:   } else argv[3] = vmath->linkhost;

    105:   /* Link mode */
    106: #if defined(MATHEMATICA_3_0)
    107:   argv[4] = "-linkmode";
    108:   switch (vmath->linkmode) {
    109:   case MATHEMATICA_LINK_CREATE:
    110:     argv[5] = "Create";
    111:     break;
    112:   case MATHEMATICA_LINK_CONNECT:
    113:     argv[5] = "Connect";
    114:     break;
    115:   case MATHEMATICA_LINK_LAUNCH:
    116:     argv[5] = "Launch";
    117:     break;
    118:   }
    119: #else
    120:   switch (vmath->linkmode) {
    121:   case MATHEMATICA_LINK_CREATE:
    122:     argv[4] = "-linkcreate";
    123:     break;
    124:   case MATHEMATICA_LINK_CONNECT:
    125:     argv[4] = "-linkconnect";
    126:     break;
    127:   case MATHEMATICA_LINK_LAUNCH:
    128:     argv[4] = "-linklaunch";
    129:     break;
    130:   }
    131: #endif
    132:   vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
    133: #endif
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
    138: {
    139:   PetscViewer_Mathematica *vmath;

    141:   PetscFunctionBegin;
    142:   PetscCall(PetscViewerMathematicaInitializePackage());

    144:   PetscCall(PetscNew(&vmath));
    145:   v->data         = (void *)vmath;
    146:   v->ops->destroy = PetscViewerDestroy_Mathematica;
    147:   v->ops->flush   = 0;
    148:   PetscCall(PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name));

    150:   vmath->linkname     = NULL;
    151:   vmath->linkhost     = NULL;
    152:   vmath->linkmode     = MATHEMATICA_LINK_CONNECT;
    153:   vmath->graphicsType = GRAPHICS_MOTIF;
    154:   vmath->plotType     = MATHEMATICA_TRIANGULATION_PLOT;
    155:   vmath->objName      = NULL;

    157:   PetscCall(PetscViewerMathematicaSetFromOptions(v));
    158:   PetscCall(PetscViewerMathematicaSetupConnection_Private(v));
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: static PetscErrorCode PetscViewerMathematicaParseLinkMode(char *modename, LinkMode *mode)
    163: {
    164:   PetscBool isCreate, isConnect, isLaunch;

    166:   PetscFunctionBegin;
    167:   PetscCall(PetscStrcasecmp(modename, "Create", &isCreate));
    168:   PetscCall(PetscStrcasecmp(modename, "Connect", &isConnect));
    169:   PetscCall(PetscStrcasecmp(modename, "Launch", &isLaunch));
    170:   if (isCreate) *mode = MATHEMATICA_LINK_CREATE;
    171:   else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
    172:   else if (isLaunch) *mode = MATHEMATICA_LINK_LAUNCH;
    173:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
    174:   PetscFunctionReturn(PETSC_SUCCESS);
    175: }

    177: PetscErrorCode PetscViewerMathematicaSetFromOptions(PetscViewer v)
    178: {
    179:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
    180:   char                     linkname[256];
    181:   char                     modename[256];
    182:   char                     hostname[256];
    183:   char                     type[256];
    184:   PetscInt                 numPorts;
    185:   PetscInt                *ports;
    186:   PetscInt                 numHosts;
    187:   int                      h;
    188:   char                   **hosts;
    189:   PetscMPIInt              size, rank;
    190:   PetscBool                opt;

    192:   PetscFunctionBegin;
    193:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
    194:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));

    196:   /* Get link name */
    197:   PetscCall(PetscOptionsGetString("viewer_", "-math_linkname", linkname, sizeof(linkname), &opt));
    198:   if (opt) PetscCall(PetscViewerMathematicaSetLinkName(v, linkname));
    199:   /* Get link port */
    200:   numPorts = size;
    201:   PetscCall(PetscMalloc1(size, &ports));
    202:   PetscCall(PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt));
    203:   if (opt) {
    204:     if (numPorts > rank) snprintf(linkname, sizeof(linkname), "%6d", ports[rank]);
    205:     else snprintf(linkname, sizeof(linkname), "%6d", ports[0]);
    206:     PetscCall(PetscViewerMathematicaSetLinkName(v, linkname));
    207:   }
    208:   PetscCall(PetscFree(ports));
    209:   /* Get link host */
    210:   numHosts = size;
    211:   PetscCall(PetscMalloc1(size, &hosts));
    212:   PetscCall(PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt));
    213:   if (opt) {
    214:     if (numHosts > rank) {
    215:       PetscCall(PetscStrncpy(hostname, hosts[rank], sizeof(hostname)));
    216:     } else {
    217:       PetscCall(PetscStrncpy(hostname, hosts[0], sizeof(hostname)));
    218:     }
    219:     PetscCall(PetscViewerMathematicaSetLinkHost(v, hostname));
    220:   }
    221:   for (h = 0; h < numHosts; h++) PetscCall(PetscFree(hosts[h]));
    222:   PetscCall(PetscFree(hosts));
    223:   /* Get link mode */
    224:   PetscCall(PetscOptionsGetString("viewer_", "-math_linkmode", modename, sizeof(modename), &opt));
    225:   if (opt) {
    226:     LinkMode mode;

    228:     PetscCall(PetscViewerMathematicaParseLinkMode(modename, &mode));
    229:     PetscCall(PetscViewerMathematicaSetLinkMode(v, mode));
    230:   }
    231:   /* Get graphics type */
    232:   PetscCall(PetscOptionsGetString("viewer_", "-math_graphics", type, sizeof(type), &opt));
    233:   if (opt) {
    234:     PetscBool isMotif, isPS, isPSFile;

    236:     PetscCall(PetscStrcasecmp(type, "Motif", &isMotif));
    237:     PetscCall(PetscStrcasecmp(type, "PS", &isPS));
    238:     PetscCall(PetscStrcasecmp(type, "PSFile", &isPSFile));
    239:     if (isMotif) vmath->graphicsType = GRAPHICS_MOTIF;
    240:     else if (isPS) vmath->graphicsType = GRAPHICS_PS_STDOUT;
    241:     else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
    242:   }
    243:   /* Get plot type */
    244:   PetscCall(PetscOptionsGetString("viewer_", "-math_type", type, sizeof(type), &opt));
    245:   if (opt) {
    246:     PetscBool isTri, isVecTri, isVec, isSurface;

    248:     PetscCall(PetscStrcasecmp(type, "Triangulation", &isTri));
    249:     PetscCall(PetscStrcasecmp(type, "VectorTriangulation", &isVecTri));
    250:     PetscCall(PetscStrcasecmp(type, "Vector", &isVec));
    251:     PetscCall(PetscStrcasecmp(type, "Surface", &isSurface));
    252:     if (isTri) vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
    253:     else if (isVecTri) vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
    254:     else if (isVec) vmath->plotType = MATHEMATICA_VECTOR_PLOT;
    255:     else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
    256:   }
    257:   PetscFunctionReturn(PETSC_SUCCESS);
    258: }

    260: PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
    261: {
    262:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;

    264:   PetscFunctionBegin;
    266:   PetscAssertPointer(name, 2);
    267:   PetscCall(PetscStrallocpy(name, &vmath->linkname));
    268:   PetscFunctionReturn(PETSC_SUCCESS);
    269: }

    271: PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
    272: {
    273:   char name[16];

    275:   PetscFunctionBegin;
    276:   snprintf(name, 16, "%6d", port);
    277:   PetscCall(PetscViewerMathematicaSetLinkName(v, name));
    278:   PetscFunctionReturn(PETSC_SUCCESS);
    279: }

    281: PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
    282: {
    283:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;

    285:   PetscFunctionBegin;
    287:   PetscAssertPointer(host, 2);
    288:   PetscCall(PetscStrallocpy(host, &vmath->linkhost));
    289:   PetscFunctionReturn(PETSC_SUCCESS);
    290: }

    292: PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
    293: {
    294:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;

    296:   PetscFunctionBegin;
    297:   vmath->linkmode = mode;
    298:   PetscFunctionReturn(PETSC_SUCCESS);
    299: }

    301: /*----------------------------------------- Public Functions --------------------------------------------------------*/
    302: /*@
    303:   PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.

    305:   Collective

    307:   Input Parameters:
    308: + comm    - The MPI communicator
    309: . port    - [optional] The port to connect on, or PETSC_DECIDE
    310: . machine - [optional] The machine to run Mathematica on, or NULL
    311: - mode    - [optional] The connection mode, or NULL

    313:   Output Parameter:
    314: . v - The Mathematica viewer

    316:   Options Database Keys:
    317: + -viewer_math_linkhost <machine> - The host machine for the kernel
    318: . -viewer_math_linkname <name>    - The full link name for the connection
    319: . -viewer_math_linkport <port>    - The port for the connection
    320: . -viewer_math_mode <mode>        - The mode, e.g. Launch, Connect
    321: . -viewer_math_type <type>        - The plot type, e.g. Triangulation, Vector
    322: - -viewer_math_graphics <output>  - The output type, e.g. Motif, PS, PSFile

    324:   Level: intermediate

    326:   Note:
    327:   Most users should employ the following commands to access the
    328:   Mathematica viewers
    329: .vb
    330:     PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
    331:     MatView(Mat matrix, PetscViewer viewer)

    333:                 or

    335:     PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
    336:     VecView(Vec vector, PetscViewer viewer)
    337: .ve

    339: .seealso: `PETSCVIEWERMATHEMATICA`, `MatView()`, `VecView()`
    340: @*/
    341: PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
    342: {
    343:   PetscFunctionBegin;
    344:   PetscCall(PetscViewerCreate(comm, v));
    345: #if 0
    346:   LinkMode linkmode;
    347:   PetscCall(PetscViewerMathematicaSetLinkPort(*v, port));
    348:   PetscCall(PetscViewerMathematicaSetLinkHost(*v, machine));
    349:   PetscCall(PetscViewerMathematicaParseLinkMode(mode, &linkmode));
    350:   PetscCall(PetscViewerMathematicaSetLinkMode(*v, linkmode));
    351: #endif
    352:   PetscCall(PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA));
    353:   PetscFunctionReturn(PETSC_SUCCESS);
    354: }

    356: /*
    357:   PetscViewerMathematicaGetLink - Returns the link to Mathematica from a `PETSCVIEWERMATHEMATICA`

    359:   Input Parameters:
    360: + viewer - The Mathematica viewer
    361: - link   - The link to Mathematica

    363:   Level: intermediate

    365: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaOpen()`
    366: */
    367: static PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
    368: {
    369:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;

    371:   PetscFunctionBegin;
    373:   *link = vmath->link;
    374:   PetscFunctionReturn(PETSC_SUCCESS);
    375: }

    377: /*@C
    378:   PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received

    380:   Input Parameters:
    381: + viewer - The Mathematica viewer
    382: - type   - The packet type to search for, e.g RETURNPKT

    384:   Level: advanced

    386: .seealso: `PetscViewerMathematicaSetName()`, `PetscViewerMathematicaGetVector()`
    387: @*/
    388: PetscErrorCode PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
    389: {
    390:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
    391:   MLINK                    link  = vmath->link; /* The link to Mathematica */
    392:   int                      pkt;                 /* The packet type */

    394:   PetscFunctionBegin;
    395:   while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
    396:   if (!pkt) {
    397:     MLClearError(link);
    398:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, (char *)MLErrorMessage(link));
    399:   }
    400:   PetscFunctionReturn(PETSC_SUCCESS);
    401: }

    403: /*@C
    404:   PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica via `PETSCVIEWERMATHEMATICA`

    406:   Input Parameter:
    407: . viewer - The Mathematica viewer

    409:   Output Parameter:
    410: . name - The name for new objects created in Mathematica

    412:   Level: intermediate

    414: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaSetName()`, `PetscViewerMathematicaClearName()`
    415: @*/
    416: PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
    417: {
    418:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;

    420:   PetscFunctionBegin;
    422:   PetscAssertPointer(name, 2);
    423:   *name = vmath->objName;
    424:   PetscFunctionReturn(PETSC_SUCCESS);
    425: }

    427: /*@C
    428:   PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica via `PETSCVIEWERMATHEMATICA`

    430:   Input Parameters:
    431: + viewer - The Mathematica viewer
    432: - name   - The name for new objects created in Mathematica

    434:   Level: intermediate

    436: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaClearName()`
    437: @*/
    438: PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
    439: {
    440:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;

    442:   PetscFunctionBegin;
    444:   PetscAssertPointer(name, 2);
    445:   vmath->objName = name;
    446:   PetscFunctionReturn(PETSC_SUCCESS);
    447: }

    449: /*@
    450:   PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica

    452:   Input Parameter:
    453: . viewer - The Mathematica viewer

    455:   Level: intermediate

    457: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaGetName()`, `PetscViewerMathematicaSetName()`
    458: @*/
    459: PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
    460: {
    461:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;

    463:   PetscFunctionBegin;
    465:   vmath->objName = NULL;
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: /*@
    470:   PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica via a `PETSCVIEWERMATHEMATICA`

    472:   Input Parameter:
    473: . viewer - The Mathematica viewer

    475:   Output Parameter:
    476: . v - The vector

    478:   Level: intermediate

    480: .seealso: `PETSCVIEWERMATHEMATICA`, `VecView()`, `PetscViewerMathematicaPutVector()`
    481: @*/
    482: PetscErrorCode PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
    483: {
    484:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
    485:   MLINK                    link; /* The link to Mathematica */
    486:   char                    *name;
    487:   PetscScalar             *mArray, *array;
    488:   long                     mSize;
    489:   int                      n;

    491:   PetscFunctionBegin;

    495:   /* Determine the object name */
    496:   if (!vmath->objName) name = "vec";
    497:   else name = (char *)vmath->objName;

    499:   link = vmath->link;
    500:   PetscCall(VecGetLocalSize(v, &n));
    501:   PetscCall(VecGetArray(v, &array));
    502:   MLPutFunction(link, "EvaluatePacket", 1);
    503:   MLPutSymbol(link, name);
    504:   MLEndPacket(link);
    505:   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
    506:   MLGetRealList(link, &mArray, &mSize);
    507:   PetscCheck(n == mSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d", n, mSize);
    508:   PetscCall(PetscArraycpy(array, mArray, mSize));
    509:   MLDisownRealList(link, mArray, mSize);
    510:   PetscCall(VecRestoreArray(v, &array));
    511:   PetscFunctionReturn(PETSC_SUCCESS);
    512: }

    514: /*@
    515:   PetscViewerMathematicaPutVector - Send a vector to Mathematica via a `PETSCVIEWERMATHEMATICA` `PetscViewer`

    517:   Input Parameters:
    518: + viewer - The Mathematica viewer
    519: - v      - The vector

    521:   Level: intermediate

    523: .seealso: `PETSCVIEWERMATHEMATICA`, `VecView()`, `PetscViewerMathematicaGetVector()`
    524: @*/
    525: PetscErrorCode PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
    526: {
    527:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
    528:   MLINK                    link  = vmath->link; /* The link to Mathematica */
    529:   char                    *name;
    530:   PetscScalar             *array;
    531:   int                      n;

    533:   PetscFunctionBegin;
    534:   /* Determine the object name */
    535:   if (!vmath->objName) name = "vec";
    536:   else name = (char *)vmath->objName;

    538:   PetscCall(VecGetLocalSize(v, &n));
    539:   PetscCall(VecGetArray(v, &array));

    541:   /* Send the Vector object */
    542:   MLPutFunction(link, "EvaluatePacket", 1);
    543:   MLPutFunction(link, "Set", 2);
    544:   MLPutSymbol(link, name);
    545:   MLPutRealList(link, array, n);
    546:   MLEndPacket(link);
    547:   /* Skip packets until ReturnPacket */
    548:   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
    549:   /* Skip ReturnPacket */
    550:   MLNewPacket(link);

    552:   PetscCall(VecRestoreArray(v, &array));
    553:   PetscFunctionReturn(PETSC_SUCCESS);
    554: }

    556: PetscErrorCode PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
    557: {
    558:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
    559:   MLINK                    link  = vmath->link; /* The link to Mathematica */
    560:   char                    *name;

    562:   PetscFunctionBegin;
    563:   /* Determine the object name */
    564:   if (!vmath->objName) name = "mat";
    565:   else name = (char *)vmath->objName;

    567:   /* Send the dense matrix object */
    568:   MLPutFunction(link, "EvaluatePacket", 1);
    569:   MLPutFunction(link, "Set", 2);
    570:   MLPutSymbol(link, name);
    571:   MLPutFunction(link, "Transpose", 1);
    572:   MLPutFunction(link, "Partition", 2);
    573:   MLPutRealList(link, a, m * n);
    574:   MLPutInteger(link, m);
    575:   MLEndPacket(link);
    576:   /* Skip packets until ReturnPacket */
    577:   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
    578:   /* Skip ReturnPacket */
    579:   MLNewPacket(link);
    580:   PetscFunctionReturn(PETSC_SUCCESS);
    581: }

    583: PetscErrorCode PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
    584: {
    585:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
    586:   MLINK                    link  = vmath->link; /* The link to Mathematica */
    587:   const char              *symbol;
    588:   char                    *name;
    589:   PetscBool                match;

    591:   PetscFunctionBegin;
    592:   /* Determine the object name */
    593:   if (!vmath->objName) name = "mat";
    594:   else name = (char *)vmath->objName;

    596:   /* Make sure Mathematica recognizes sparse matrices */
    597:   MLPutFunction(link, "EvaluatePacket", 1);
    598:   MLPutFunction(link, "Needs", 1);
    599:   MLPutString(link, "LinearAlgebra`CSRMatrix`");
    600:   MLEndPacket(link);
    601:   /* Skip packets until ReturnPacket */
    602:   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
    603:   /* Skip ReturnPacket */
    604:   MLNewPacket(link);

    606:   /* Send the CSRMatrix object */
    607:   MLPutFunction(link, "EvaluatePacket", 1);
    608:   MLPutFunction(link, "Set", 2);
    609:   MLPutSymbol(link, name);
    610:   MLPutFunction(link, "CSRMatrix", 5);
    611:   MLPutInteger(link, m);
    612:   MLPutInteger(link, n);
    613:   MLPutFunction(link, "Plus", 2);
    614:   MLPutIntegerList(link, i, m + 1);
    615:   MLPutInteger(link, 1);
    616:   MLPutFunction(link, "Plus", 2);
    617:   MLPutIntegerList(link, j, i[m]);
    618:   MLPutInteger(link, 1);
    619:   MLPutRealList(link, a, i[m]);
    620:   MLEndPacket(link);
    621:   /* Skip packets until ReturnPacket */
    622:   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
    623:   /* Skip ReturnPacket */
    624:   MLNewPacket(link);

    626:   /* Check that matrix is valid */
    627:   MLPutFunction(link, "EvaluatePacket", 1);
    628:   MLPutFunction(link, "ValidQ", 1);
    629:   MLPutSymbol(link, name);
    630:   MLEndPacket(link);
    631:   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
    632:   MLGetSymbol(link, &symbol);
    633:   PetscCall(PetscStrcmp("True", (char *)symbol, &match));
    634:   if (!match) {
    635:     MLDisownSymbol(link, symbol);
    636:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
    637:   }
    638:   MLDisownSymbol(link, symbol);
    639:   /* Skip ReturnPacket */
    640:   MLNewPacket(link);
    641:   PetscFunctionReturn(PETSC_SUCCESS);
    642: }