xref: /petsc/src/sys/classes/viewer/impls/mathematica/mathematica.c (revision 03047865b8d8757cf1cf9cda45785c1537b01dc1) !
1 #include <petsc/private/viewerimpl.h> /* "petscsys.h" */
2 #include <mathematica.h>
3 
4 #if defined(PETSC_HAVE__SNPRINTF) && !defined(PETSC_HAVE_SNPRINTF)
5   #define snprintf _snprintf
6 #endif
7 
8 PetscViewer  PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE = NULL;
9 static void *mathematicaEnv                         = NULL;
10 
11 static PetscBool PetscViewerMathematicaPackageInitialized = PETSC_FALSE;
12 /*@C
13   PetscViewerMathematicaFinalizePackage - This function destroys everything in the PETSc interface to Mathematica. It is
14   called from PetscFinalize().
15 
16   Level: developer
17 
18 .seealso: `PetscFinalize()`
19 @*/
PetscViewerMathematicaFinalizePackage(void)20 PetscErrorCode PetscViewerMathematicaFinalizePackage(void)
21 {
22   PetscFunctionBegin;
23   if (mathematicaEnv) MLDeinitialize((MLEnvironment)mathematicaEnv);
24   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
25   PetscFunctionReturn(PETSC_SUCCESS);
26 }
27 
28 /*@C
29   PetscViewerMathematicaInitializePackage - This function initializes everything in the PETSc interface to Mathematica. It is
30   called from `PetscViewerInitializePackage()`.
31 
32   Level: developer
33 
34 .seealso: `PetscSysInitializePackage()`, `PetscInitialize()`
35 @*/
PetscViewerMathematicaInitializePackage(void)36 PetscErrorCode PetscViewerMathematicaInitializePackage(void)
37 {
38   PetscError ierr;
39 
40   PetscFunctionBegin;
41   if (PetscViewerMathematicaPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
42   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
43 
44   mathematicaEnv = (void *)MLInitialize(0);
45 
46   PetscCall(PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage));
47   PetscFunctionReturn(PETSC_SUCCESS);
48 }
49 
PetscViewerInitializeMathematicaWorld_Private()50 static PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
51 {
52   PetscFunctionBegin;
53   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) PetscFunctionReturn(PETSC_SUCCESS);
54   PetscCall(PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE));
55   PetscFunctionReturn(PETSC_SUCCESS);
56 }
57 
PetscViewerDestroy_Mathematica(PetscViewer viewer)58 static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
59 {
60   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
61 
62   PetscFunctionBegin;
63   MLClose(vmath->link);
64   PetscCall(PetscFree(vmath->linkname));
65   PetscCall(PetscFree(vmath->linkhost));
66   PetscCall(PetscFree(vmath));
67   PetscFunctionReturn(PETSC_SUCCESS);
68 }
69 
PetscViewerDestroyMathematica_Private(void)70 static PetscErrorCode PetscViewerDestroyMathematica_Private(void)
71 {
72   PetscFunctionBegin;
73   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) PetscCall(PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE));
74   PetscFunctionReturn(PETSC_SUCCESS);
75 }
76 
PetscViewerMathematicaSetupConnection_Private(PetscViewer v)77 static PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
78 {
79   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
80 #if defined(MATHEMATICA_3_0)
81   int   argc = 6;
82   char *argv[6];
83 #else
84   int   argc = 5;
85   char *argv[5];
86 #endif
87   char hostname[256];
88   long lerr;
89 
90   PetscFunctionBegin;
91   /* Link name */
92   argv[0] = "-linkname";
93   if (!vmath->linkname) argv[1] = "math -mathlink";
94   else argv[1] = vmath->linkname;
95 
96   /* Link host */
97   argv[2] = "-linkhost";
98   if (!vmath->linkhost) {
99     PetscCall(PetscGetHostName(hostname, sizeof(hostname)));
100     argv[3] = hostname;
101   } else argv[3] = vmath->linkhost;
102 
103   /* Link mode */
104 #if defined(MATHEMATICA_3_0)
105   argv[4] = "-linkmode";
106   switch (vmath->linkmode) {
107   case MATHEMATICA_LINK_CREATE:
108     argv[5] = "Create";
109     break;
110   case MATHEMATICA_LINK_CONNECT:
111     argv[5] = "Connect";
112     break;
113   case MATHEMATICA_LINK_LAUNCH:
114     argv[5] = "Launch";
115     break;
116   }
117 #else
118   switch (vmath->linkmode) {
119   case MATHEMATICA_LINK_CREATE:
120     argv[4] = "-linkcreate";
121     break;
122   case MATHEMATICA_LINK_CONNECT:
123     argv[4] = "-linkconnect";
124     break;
125   case MATHEMATICA_LINK_LAUNCH:
126     argv[4] = "-linklaunch";
127     break;
128   }
129 #endif
130   vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
131 #endif
132   PetscFunctionReturn(PETSC_SUCCESS);
133 }
134 
PetscViewerCreate_Mathematica(PetscViewer v)135 PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
136 {
137   PetscViewer_Mathematica *vmath;
138 
139   PetscFunctionBegin;
140   PetscCall(PetscViewerMathematicaInitializePackage());
141 
142   PetscCall(PetscNew(&vmath));
143   v->data         = (void *)vmath;
144   v->ops->destroy = PetscViewerDestroy_Mathematica;
145   v->ops->flush   = 0;
146   PetscCall(PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name));
147 
148   vmath->linkname     = NULL;
149   vmath->linkhost     = NULL;
150   vmath->linkmode     = MATHEMATICA_LINK_CONNECT;
151   vmath->graphicsType = GRAPHICS_MOTIF;
152   vmath->plotType     = MATHEMATICA_TRIANGULATION_PLOT;
153   vmath->objName      = NULL;
154 
155   PetscCall(PetscViewerMathematicaSetFromOptions(v));
156   PetscCall(PetscViewerMathematicaSetupConnection_Private(v));
157   PetscFunctionReturn(PETSC_SUCCESS);
158 }
159 
PetscViewerMathematicaParseLinkMode(char * modename,LinkMode * mode)160 static PetscErrorCode PetscViewerMathematicaParseLinkMode(char *modename, LinkMode *mode)
161 {
162   PetscBool isCreate, isConnect, isLaunch;
163 
164   PetscFunctionBegin;
165   PetscCall(PetscStrcasecmp(modename, "Create", &isCreate));
166   PetscCall(PetscStrcasecmp(modename, "Connect", &isConnect));
167   PetscCall(PetscStrcasecmp(modename, "Launch", &isLaunch));
168   if (isCreate) *mode = MATHEMATICA_LINK_CREATE;
169   else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
170   else if (isLaunch) *mode = MATHEMATICA_LINK_LAUNCH;
171   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
172   PetscFunctionReturn(PETSC_SUCCESS);
173 }
174 
PetscViewerMathematicaSetFromOptions(PetscViewer v)175 PetscErrorCode PetscViewerMathematicaSetFromOptions(PetscViewer v)
176 {
177   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
178   char                     linkname[256];
179   char                     modename[256];
180   char                     hostname[256];
181   char                     type[256];
182   PetscInt                 numPorts;
183   PetscInt                *ports;
184   PetscInt                 numHosts;
185   int                      h;
186   char                   **hosts;
187   PetscMPIInt              size, rank;
188   PetscBool                opt;
189 
190   PetscFunctionBegin;
191   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
192   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
193 
194   /* Get link name */
195   PetscCall(PetscOptionsGetString("viewer_", "-math_linkname", linkname, sizeof(linkname), &opt));
196   if (opt) PetscCall(PetscViewerMathematicaSetLinkName(v, linkname));
197   /* Get link port */
198   numPorts = size;
199   PetscCall(PetscMalloc1(size, &ports));
200   PetscCall(PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt));
201   if (opt) {
202     if (numPorts > rank) snprintf(linkname, sizeof(linkname), "%6d", ports[rank]);
203     else snprintf(linkname, sizeof(linkname), "%6d", ports[0]);
204     PetscCall(PetscViewerMathematicaSetLinkName(v, linkname));
205   }
206   PetscCall(PetscFree(ports));
207   /* Get link host */
208   numHosts = size;
209   PetscCall(PetscMalloc1(size, &hosts));
210   PetscCall(PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt));
211   if (opt) {
212     if (numHosts > rank) {
213       PetscCall(PetscStrncpy(hostname, hosts[rank], sizeof(hostname)));
214     } else {
215       PetscCall(PetscStrncpy(hostname, hosts[0], sizeof(hostname)));
216     }
217     PetscCall(PetscViewerMathematicaSetLinkHost(v, hostname));
218   }
219   for (h = 0; h < numHosts; h++) PetscCall(PetscFree(hosts[h]));
220   PetscCall(PetscFree(hosts));
221   /* Get link mode */
222   PetscCall(PetscOptionsGetString("viewer_", "-math_linkmode", modename, sizeof(modename), &opt));
223   if (opt) {
224     LinkMode mode;
225 
226     PetscCall(PetscViewerMathematicaParseLinkMode(modename, &mode));
227     PetscCall(PetscViewerMathematicaSetLinkMode(v, mode));
228   }
229   /* Get graphics type */
230   PetscCall(PetscOptionsGetString("viewer_", "-math_graphics", type, sizeof(type), &opt));
231   if (opt) {
232     PetscBool isMotif, isPS, isPSFile;
233 
234     PetscCall(PetscStrcasecmp(type, "Motif", &isMotif));
235     PetscCall(PetscStrcasecmp(type, "PS", &isPS));
236     PetscCall(PetscStrcasecmp(type, "PSFile", &isPSFile));
237     if (isMotif) vmath->graphicsType = GRAPHICS_MOTIF;
238     else if (isPS) vmath->graphicsType = GRAPHICS_PS_STDOUT;
239     else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
240   }
241   /* Get plot type */
242   PetscCall(PetscOptionsGetString("viewer_", "-math_type", type, sizeof(type), &opt));
243   if (opt) {
244     PetscBool isTri, isVecTri, isVec, isSurface;
245 
246     PetscCall(PetscStrcasecmp(type, "Triangulation", &isTri));
247     PetscCall(PetscStrcasecmp(type, "VectorTriangulation", &isVecTri));
248     PetscCall(PetscStrcasecmp(type, "Vector", &isVec));
249     PetscCall(PetscStrcasecmp(type, "Surface", &isSurface));
250     if (isTri) vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
251     else if (isVecTri) vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
252     else if (isVec) vmath->plotType = MATHEMATICA_VECTOR_PLOT;
253     else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
254   }
255   PetscFunctionReturn(PETSC_SUCCESS);
256 }
257 
PetscViewerMathematicaSetLinkName(PetscViewer v,const char * name)258 PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
259 {
260   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
261 
262   PetscFunctionBegin;
263   PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 1);
264   PetscAssertPointer(name, 2);
265   PetscCall(PetscStrallocpy(name, &vmath->linkname));
266   PetscFunctionReturn(PETSC_SUCCESS);
267 }
268 
PetscViewerMathematicaSetLinkPort(PetscViewer v,int port)269 PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
270 {
271   char name[16];
272 
273   PetscFunctionBegin;
274   snprintf(name, 16, "%6d", port);
275   PetscCall(PetscViewerMathematicaSetLinkName(v, name));
276   PetscFunctionReturn(PETSC_SUCCESS);
277 }
278 
PetscViewerMathematicaSetLinkHost(PetscViewer v,const char * host)279 PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
280 {
281   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
282 
283   PetscFunctionBegin;
284   PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 1);
285   PetscAssertPointer(host, 2);
286   PetscCall(PetscStrallocpy(host, &vmath->linkhost));
287   PetscFunctionReturn(PETSC_SUCCESS);
288 }
289 
PetscViewerMathematicaSetLinkMode(PetscViewer v,LinkMode mode)290 PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
291 {
292   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
293 
294   PetscFunctionBegin;
295   vmath->linkmode = mode;
296   PetscFunctionReturn(PETSC_SUCCESS);
297 }
298 
299 /*@
300   PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
301 
302   Collective
303 
304   Input Parameters:
305 + comm    - The MPI communicator
306 . port    - [optional] The port to connect on, or PETSC_DECIDE
307 . machine - [optional] The machine to run Mathematica on, or NULL
308 - mode    - [optional] The connection mode, or NULL
309 
310   Output Parameter:
311 . v - The Mathematica viewer
312 
313   Options Database Keys:
314 + -viewer_math_linkhost <machine> - The host machine for the kernel
315 . -viewer_math_linkname <name>    - The full link name for the connection
316 . -viewer_math_linkport <port>    - The port for the connection
317 . -viewer_math_mode <mode>        - The mode, e.g. Launch, Connect
318 . -viewer_math_type <type>        - The plot type, e.g. Triangulation, Vector
319 - -viewer_math_graphics <output>  - The output type, e.g. Motif, PS, PSFile
320 
321   Level: intermediate
322 
323   Note:
324   Most users should employ the following commands to access the
325   Mathematica viewers
326 .vb
327     PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
328     MatView(Mat matrix, PetscViewer viewer)
329 
330                 or
331 
332     PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
333     VecView(Vec vector, PetscViewer viewer)
334 .ve
335 
336 .seealso: `PETSCVIEWERMATHEMATICA`, `MatView()`, `VecView()`
337 @*/
PetscViewerMathematicaOpen(MPI_Comm comm,int port,const char machine[],const char mode[],PetscViewer * v)338 PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
339 {
340   PetscFunctionBegin;
341   PetscCall(PetscViewerCreate(comm, v));
342 #if 0
343   LinkMode linkmode;
344   PetscCall(PetscViewerMathematicaSetLinkPort(*v, port));
345   PetscCall(PetscViewerMathematicaSetLinkHost(*v, machine));
346   PetscCall(PetscViewerMathematicaParseLinkMode(mode, &linkmode));
347   PetscCall(PetscViewerMathematicaSetLinkMode(*v, linkmode));
348 #endif
349   PetscCall(PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA));
350   PetscFunctionReturn(PETSC_SUCCESS);
351 }
352 
353 /*
354   PetscViewerMathematicaGetLink - Returns the link to Mathematica from a `PETSCVIEWERMATHEMATICA`
355 
356   Input Parameters:
357 + viewer - The Mathematica viewer
358 - link   - The link to Mathematica
359 
360   Level: intermediate
361 
362 .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaOpen()`
363 */
PetscViewerMathematicaGetLink(PetscViewer viewer,MLINK * link)364 static PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
365 {
366   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
367 
368   PetscFunctionBegin;
369   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1);
370   *link = vmath->link;
371   PetscFunctionReturn(PETSC_SUCCESS);
372 }
373 
374 /*@C
375   PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
376 
377   Input Parameters:
378 + viewer - The Mathematica viewer
379 - type   - The packet type to search for, e.g RETURNPKT
380 
381   Level: advanced
382 
383 .seealso: `PetscViewerMathematicaSetName()`, `PetscViewerMathematicaGetVector()`
384 @*/
PetscViewerMathematicaSkipPackets(PetscViewer viewer,int type)385 PetscErrorCode PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
386 {
387   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
388   MLINK                    link  = vmath->link; /* The link to Mathematica */
389   int                      pkt;                 /* The packet type */
390 
391   PetscFunctionBegin;
392   while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
393   if (!pkt) {
394     MLClearError(link);
395     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, (char *)MLErrorMessage(link));
396   }
397   PetscFunctionReturn(PETSC_SUCCESS);
398 }
399 
400 /*@C
401   PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica via `PETSCVIEWERMATHEMATICA`
402 
403   Input Parameter:
404 . viewer - The Mathematica viewer
405 
406   Output Parameter:
407 . name - The name for new objects created in Mathematica
408 
409   Level: intermediate
410 
411 .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaSetName()`, `PetscViewerMathematicaClearName()`
412 @*/
PetscViewerMathematicaGetName(PetscViewer viewer,const char ** name)413 PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
414 {
415   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
416 
417   PetscFunctionBegin;
418   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1);
419   PetscAssertPointer(name, 2);
420   *name = vmath->objName;
421   PetscFunctionReturn(PETSC_SUCCESS);
422 }
423 
424 /*@C
425   PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica via `PETSCVIEWERMATHEMATICA`
426 
427   Input Parameters:
428 + viewer - The Mathematica viewer
429 - name   - The name for new objects created in Mathematica
430 
431   Level: intermediate
432 
433 .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaClearName()`
434 @*/
PetscViewerMathematicaSetName(PetscViewer viewer,const char name[])435 PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
436 {
437   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
438 
439   PetscFunctionBegin;
440   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1);
441   PetscAssertPointer(name, 2);
442   vmath->objName = name;
443   PetscFunctionReturn(PETSC_SUCCESS);
444 }
445 
446 /*@
447   PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
448 
449   Input Parameter:
450 . viewer - The Mathematica viewer
451 
452   Level: intermediate
453 
454 .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaGetName()`, `PetscViewerMathematicaSetName()`
455 @*/
PetscViewerMathematicaClearName(PetscViewer viewer)456 PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
457 {
458   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
459 
460   PetscFunctionBegin;
461   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1);
462   vmath->objName = NULL;
463   PetscFunctionReturn(PETSC_SUCCESS);
464 }
465 
PetscViewerMathematicaPutMatrix(PetscViewer viewer,int m,int n,PetscReal * a)466 PetscErrorCode PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
467 {
468   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
469   MLINK                    link  = vmath->link; /* The link to Mathematica */
470   char                    *name;
471 
472   PetscFunctionBegin;
473   /* Determine the object name */
474   if (!vmath->objName) name = "mat";
475   else name = (char *)vmath->objName;
476 
477   /* Send the dense matrix object */
478   MLPutFunction(link, "EvaluatePacket", 1);
479   MLPutFunction(link, "Set", 2);
480   MLPutSymbol(link, name);
481   MLPutFunction(link, "Transpose", 1);
482   MLPutFunction(link, "Partition", 2);
483   MLPutRealList(link, a, m * n);
484   MLPutInteger(link, m);
485   MLEndPacket(link);
486   /* Skip packets until ReturnPacket */
487   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
488   /* Skip ReturnPacket */
489   MLNewPacket(link);
490   PetscFunctionReturn(PETSC_SUCCESS);
491 }
492 
PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer,int m,int n,int * i,int * j,PetscReal * a)493 PetscErrorCode PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
494 {
495   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
496   MLINK                    link  = vmath->link; /* The link to Mathematica */
497   const char              *symbol;
498   char                    *name;
499   PetscBool                match;
500 
501   PetscFunctionBegin;
502   /* Determine the object name */
503   if (!vmath->objName) name = "mat";
504   else name = (char *)vmath->objName;
505 
506   /* Make sure Mathematica recognizes sparse matrices */
507   MLPutFunction(link, "EvaluatePacket", 1);
508   MLPutFunction(link, "Needs", 1);
509   MLPutString(link, "LinearAlgebra`CSRMatrix`");
510   MLEndPacket(link);
511   /* Skip packets until ReturnPacket */
512   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
513   /* Skip ReturnPacket */
514   MLNewPacket(link);
515 
516   /* Send the CSRMatrix object */
517   MLPutFunction(link, "EvaluatePacket", 1);
518   MLPutFunction(link, "Set", 2);
519   MLPutSymbol(link, name);
520   MLPutFunction(link, "CSRMatrix", 5);
521   MLPutInteger(link, m);
522   MLPutInteger(link, n);
523   MLPutFunction(link, "Plus", 2);
524   MLPutIntegerList(link, i, m + 1);
525   MLPutInteger(link, 1);
526   MLPutFunction(link, "Plus", 2);
527   MLPutIntegerList(link, j, i[m]);
528   MLPutInteger(link, 1);
529   MLPutRealList(link, a, i[m]);
530   MLEndPacket(link);
531   /* Skip packets until ReturnPacket */
532   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
533   /* Skip ReturnPacket */
534   MLNewPacket(link);
535 
536   /* Check that matrix is valid */
537   MLPutFunction(link, "EvaluatePacket", 1);
538   MLPutFunction(link, "ValidQ", 1);
539   MLPutSymbol(link, name);
540   MLEndPacket(link);
541   PetscCall(PetscViewerMathematicaSkipPackets(viewer, RETURNPKT));
542   MLGetSymbol(link, &symbol);
543   PetscCall(PetscStrcmp("True", (char *)symbol, &match));
544   if (!match) {
545     MLDisownSymbol(link, symbol);
546     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
547   }
548   MLDisownSymbol(link, symbol);
549   /* Skip ReturnPacket */
550   MLNewPacket(link);
551   PetscFunctionReturn(PETSC_SUCCESS);
552 }
553