Actual source code: mathematica.c


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

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

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

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

 19:   Level: developer

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

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

 35:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

265:   PetscFunctionBegin;
268:   PetscCall(PetscStrallocpy(name, &vmath->linkname));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

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

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

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

286:   PetscFunctionBegin;
289:   PetscCall(PetscStrallocpy(host, &vmath->linkhost));
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

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

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

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

306:   Collective

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

314:   Output Parameter:
315: . viewer  - The Mathematica viewer

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

325:   Level: intermediate

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

334:                 or

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

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

357: /*@C
358:   PetscViewerMathematicaGetLink - Returns the link to Mathematica from a `PETSCVIEWERMATHEMATICA`

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

364:   Level: intermediate

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

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

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

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

385:   Level: advanced

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

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

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

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

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

413:   Level: intermediate

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

421:   PetscFunctionBegin;
424:   *name = vmath->objName;
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }

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

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

435:   Level: intermediate

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

443:   PetscFunctionBegin;
446:   vmath->objName = name;
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

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

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

456:   Level: intermediate

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

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

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

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

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

479:   Level: intermediate

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

492:   PetscFunctionBegin;

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

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

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

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

522:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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