Actual source code: mathematica.c

petsc-3.8.3 2017-12-09
Report Typos and Errors

  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: .keywords: Petsc, destroy, package, mathematica
 22: .seealso: PetscFinalize()
 23: @*/
 24: PetscErrorCode  PetscViewerMathematicaFinalizePackage(void)
 25: {
 27:   if (mathematicaEnv) MLDeinitialize((MLEnvironment) mathematicaEnv);
 28:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
 29:   return(0);
 30: }

 32: /*@C
 33:   PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
 34:   called from PetscDLLibraryRegister() when using dynamic libraries, and on the call to PetscInitialize()
 35:   when using static libraries.

 37:   Level: developer

 39: .keywords: Petsc, initialize, package
 40: .seealso: PetscSysInitializePackage(), PetscInitialize()
 41: @*/
 42: PetscErrorCode  PetscViewerMathematicaInitializePackage(void)
 43: {
 44:   PetscError ierr;

 47:   if (PetscViewerMathematicaPackageInitialized) return(0);
 48:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;

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

 52:   PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
 53:   return(0);
 54: }


 57: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
 58: {

 62:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return(0);
 63:   PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
 64:   return(0);
 65: }

 67: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
 68: {
 69:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
 70:   PetscErrorCode          ierr;

 73:   MLClose(vmath->link);
 74:   PetscFree(vmath->linkname);
 75:   PetscFree(vmath->linkhost);
 76:   PetscFree(vmath);
 77:   return(0);
 78: }

 80: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
 81: {

 85:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
 86:     PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
 87:   }
 88:   return(0);
 89: }

 91: PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
 92: {
 93:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
 94: #if defined(MATHEMATICA_3_0)
 95:   int                     argc = 6;
 96:   char                    *argv[6];
 97: #else
 98:   int                     argc = 5;
 99:   char                    *argv[5];
100: #endif
101:   char                    hostname[256];
102:   long                    lerr;
103:   PetscErrorCode          ierr;

106:   /* Link name */
107:   argv[0] = "-linkname";
108:   if (!vmath->linkname) argv[1] = "math -mathlink";
109:   else                  argv[1] = vmath->linkname;

111:   /* Link host */
112:   argv[2] = "-linkhost";
113:   if (!vmath->linkhost) {
114:     PetscGetHostName(hostname, 255);
115:     argv[3] = hostname;
116:   } else argv[3] = vmath->linkhost;

118:   /* Link mode */
119: #if defined(MATHEMATICA_3_0)
120:   argv[4] = "-linkmode";
121:   switch (vmath->linkmode) {
122:   case MATHEMATICA_LINK_CREATE:
123:     argv[5] = "Create";
124:     break;
125:   case MATHEMATICA_LINK_CONNECT:
126:     argv[5] = "Connect";
127:     break;
128:   case MATHEMATICA_LINK_LAUNCH:
129:     argv[5] = "Launch";
130:     break;
131:   }
132: #else
133:   switch (vmath->linkmode) {
134:   case MATHEMATICA_LINK_CREATE:
135:     argv[4] = "-linkcreate";
136:     break;
137:   case MATHEMATICA_LINK_CONNECT:
138:     argv[4] = "-linkconnect";
139:     break;
140:   case MATHEMATICA_LINK_LAUNCH:
141:     argv[4] = "-linklaunch";
142:     break;
143:   }
144: #endif
145:   vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
146: #endif
147:   return(0);
148: }

150: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
151: {
152:   PetscViewer_Mathematica *vmath;
153:   PetscErrorCode          ierr;

156:   PetscViewerMathematicaInitializePackage();

158:   PetscNewLog(v,&vmath);
159:   v->data         = (void*) vmath;
160:   v->ops->destroy = PetscViewerDestroy_Mathematica;
161:   v->ops->flush   = 0;
162:   PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);

164:   vmath->linkname     = NULL;
165:   vmath->linkhost     = NULL;
166:   vmath->linkmode     = MATHEMATICA_LINK_CONNECT;
167:   vmath->graphicsType = GRAPHICS_MOTIF;
168:   vmath->plotType     = MATHEMATICA_TRIANGULATION_PLOT;
169:   vmath->objName      = NULL;

171:   PetscViewerMathematicaSetFromOptions(v);
172:   PetscViewerMathematicaSetupConnection_Private(v);
173:   return(0);
174: }

176: PetscErrorCode PetscViewerMathematicaParseLinkMode_Private(char *modename, LinkMode *mode)
177: {
178:   PetscBool      isCreate, isConnect, isLaunch;

182:   PetscStrcasecmp(modename, "Create",  &isCreate);
183:   PetscStrcasecmp(modename, "Connect", &isConnect);
184:   PetscStrcasecmp(modename, "Launch",  &isLaunch);
185:   if (isCreate)       *mode = MATHEMATICA_LINK_CREATE;
186:   else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
187:   else if (isLaunch)  *mode = MATHEMATICA_LINK_LAUNCH;
188:   else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
189:   return(0);
190: }

192: PetscErrorCode  PetscViewerMathematicaSetFromOptions(PetscViewer v)
193: {
194:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
195:   char                    linkname[256];
196:   char                    modename[256];
197:   char                    hostname[256];
198:   char                    type[256];
199:   PetscInt                numPorts;
200:   PetscInt                *ports;
201:   PetscInt                numHosts;
202:   int                     h;
203:   char                    **hosts;
204:   PetscMPIInt             size, rank;
205:   PetscBool               opt;
206:   PetscErrorCode          ierr;

209:   MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
210:   MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);

212:   /* Get link name */
213:   PetscOptionsGetString("viewer_", "-math_linkname", linkname, 256, &opt);
214:   if (opt) {
215:     PetscViewerMathematicaSetLinkName(v, linkname);
216:   }
217:   /* Get link port */
218:   numPorts = size;
219:   PetscMalloc1(size, &ports);
220:   PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
221:   if (opt) {
222:     if (numPorts > rank) snprintf(linkname, 255, "%6d", ports[rank]);
223:     else                 snprintf(linkname, 255, "%6d", ports[0]);
224:     PetscViewerMathematicaSetLinkName(v, linkname);
225:   }
226:   PetscFree(ports);
227:   /* Get link host */
228:   numHosts = size;
229:   PetscMalloc1(size, &hosts);
230:   PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
231:   if (opt) {
232:     if (numHosts > rank) {
233:       PetscStrncpy(hostname, hosts[rank], 255);
234:     } else {
235:       PetscStrncpy(hostname, hosts[0], 255);
236:     }
237:     PetscViewerMathematicaSetLinkHost(v, hostname);
238:   }
239:   for (h = 0; h < numHosts; h++) {
240:     PetscFree(hosts[h]);
241:   }
242:   PetscFree(hosts);
243:   /* Get link mode */
244:   PetscOptionsGetString("viewer_", "-math_linkmode", modename, 256, &opt);
245:   if (opt) {
246:     LinkMode mode;

248:     PetscViewerMathematicaParseLinkMode_Private(modename, &mode);
249:     PetscViewerMathematicaSetLinkMode(v, mode);
250:   }
251:   /* Get graphics type */
252:   PetscOptionsGetString("viewer_", "-math_graphics", type, 256, &opt);
253:   if (opt) {
254:     PetscBool isMotif, isPS, isPSFile;

256:     PetscStrcasecmp(type, "Motif",  &isMotif);
257:     PetscStrcasecmp(type, "PS",     &isPS);
258:     PetscStrcasecmp(type, "PSFile", &isPSFile);
259:     if (isMotif)       vmath->graphicsType = GRAPHICS_MOTIF;
260:     else if (isPS)     vmath->graphicsType = GRAPHICS_PS_STDOUT;
261:     else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
262:   }
263:   /* Get plot type */
264:   PetscOptionsGetString("viewer_", "-math_type", type, 256, &opt);
265:   if (opt) {
266:     PetscBool isTri, isVecTri, isVec, isSurface;

268:     PetscStrcasecmp(type, "Triangulation",       &isTri);
269:     PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
270:     PetscStrcasecmp(type, "Vector",              &isVec);
271:     PetscStrcasecmp(type, "Surface",             &isSurface);
272:     if (isTri)          vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
273:     else if (isVecTri)  vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
274:     else if (isVec)     vmath->plotType = MATHEMATICA_VECTOR_PLOT;
275:     else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
276:   }
277:   return(0);
278: }

280: PetscErrorCode  PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
281: {
282:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
283:   PetscErrorCode          ierr;

288:   PetscStrallocpy(name, &vmath->linkname);
289:   return(0);
290: }

292: PetscErrorCode  PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
293: {
294:   char           name[16];

298:   snprintf(name, 16, "%6d", port);
299:   PetscViewerMathematicaSetLinkName(v, name);
300:   return(0);
301: }

303: PetscErrorCode  PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
304: {
305:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
306:   PetscErrorCode          ierr;

311:   PetscStrallocpy(host, &vmath->linkhost);
312:   return(0);
313: }

315: PetscErrorCode  PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
316: {
317:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;

320:   vmath->linkmode = mode;
321:   return(0);
322: }

324: /*----------------------------------------- Public Functions --------------------------------------------------------*/
325: /*@C
326:   PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.

328:   Collective on comm

330:   Input Parameters:
331: + comm    - The MPI communicator
332: . port    - [optional] The port to connect on, or PETSC_DECIDE
333: . machine - [optional] The machine to run Mathematica on, or NULL
334: - mode    - [optional] The connection mode, or NULL

336:   Output Parameter:
337: . viewer  - The Mathematica viewer

339:   Level: intermediate

341:   Notes:
342:   Most users should employ the following commands to access the
343:   Mathematica viewers
344: $
345: $    PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
346: $    MatView(Mat matrix, PetscViewer viewer)
347: $
348: $                or
349: $
350: $    PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
351: $    VecView(Vec vector, PetscViewer viewer)

353:    Options Database Keys:
354: +    -viewer_math_linkhost <machine> - The host machine for the kernel
355: .    -viewer_math_linkname <name>    - The full link name for the connection
356: .    -viewer_math_linkport <port>    - The port for the connection
357: .    -viewer_math_mode <mode>        - The mode, e.g. Launch, Connect
358: .    -viewer_math_type <type>        - The plot type, e.g. Triangulation, Vector
359: -    -viewer_math_graphics <output>  - The output type, e.g. Motif, PS, PSFile

361: .keywords: PetscViewer, Mathematica, open

363: .seealso: MatView(), VecView()
364: @*/
365: PetscErrorCode  PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
366: {

370:   PetscViewerCreate(comm, v);
371: #if 0
372:   LinkMode linkmode;
373:   PetscViewerMathematicaSetLinkPort(*v, port);
374:   PetscViewerMathematicaSetLinkHost(*v, machine);
375:   PetscViewerMathematicaParseLinkMode_Private(mode, &linkmode);
376:   PetscViewerMathematicaSetLinkMode(*v, linkmode);
377: #endif
378:   PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
379:   return(0);
380: }

382: /*@C
383:   PetscViewerMathematicaGetLink - Returns the link to Mathematica

385:   Input Parameters:
386: . viewer - The Mathematica viewer
387: . link   - The link to Mathematica

389:   Level: intermediate

391: .keywords PetscViewer, Mathematica, link
392: .seealso PetscViewerMathematicaOpen()
393: @*/
394: PetscErrorCode  PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
395: {
396:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

400:   *link = vmath->link;
401:   return(0);
402: }

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

407:   Input Parameters:
408: . viewer - The Mathematica viewer
409: . type   - The packet type to search for, e.g RETURNPKT

411:   Level: advanced

413: .keywords PetscViewer, Mathematica, packets
414: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
415: @*/
416: PetscErrorCode  PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
417: {
418:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
419:   MLINK                   link   = vmath->link; /* The link to Mathematica */
420:   int                     pkt;                 /* The packet type */

423:   while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
424:   if (!pkt) {
425:     MLClearError(link);
426:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char*) MLErrorMessage(link));
427:   }
428:   return(0);
429: }

431: /*@C
432:   PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica

434:   Input Parameter:
435: . viewer - The Mathematica viewer

437:   Output Parameter:
438: . name   - The name for new objects created in Mathematica

440:   Level: intermediate

442: .keywords PetscViewer, Mathematica, name
443: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
444: @*/
445: PetscErrorCode  PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
446: {
447:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

452:   *name = vmath->objName;
453:   return(0);
454: }

456: /*@C
457:   PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica

459:   Input Parameters:
460: . viewer - The Mathematica viewer
461: . name   - The name for new objects created in Mathematica

463:   Level: intermediate

465: .keywords PetscViewer, Mathematica, name
466: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
467: @*/
468: PetscErrorCode  PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
469: {
470:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

475:   vmath->objName = name;
476:   return(0);
477: }

479: /*@C
480:   PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica

482:   Input Parameter:
483: . viewer - The Mathematica viewer

485:   Level: intermediate

487: .keywords PetscViewer, Mathematica, name
488: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
489: @*/
490: PetscErrorCode  PetscViewerMathematicaClearName(PetscViewer viewer)
491: {
492:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

496:   vmath->objName = NULL;
497:   return(0);
498: }

500: /*@C
501:   PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica

503:   Input Parameter:
504: . viewer - The Mathematica viewer

506:   Output Parameter:
507: . v      - The vector

509:   Level: intermediate

511: .keywords PetscViewer, Mathematica, vector
512: .seealso VecView(), PetscViewerMathematicaPutVector()
513: @*/
514: PetscErrorCode  PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
515: {
516:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
517:   MLINK                   link;   /* The link to Mathematica */
518:   char                    *name;
519:   PetscScalar             *mArray,*array;
520:   long                    mSize;
521:   int                     n;
522:   PetscErrorCode          ierr;


528:   /* Determine the object name */
529:   if (!vmath->objName) name = "vec";
530:   else                 name = (char*) vmath->objName;

532:   link = vmath->link;
533:   VecGetLocalSize(v, &n);
534:   VecGetArray(v, &array);
535:   MLPutFunction(link, "EvaluatePacket", 1);
536:   MLPutSymbol(link, name);
537:   MLEndPacket(link);
538:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
539:   MLGetRealList(link, &mArray, &mSize);
540:   if (n != mSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d",n,mSize);
541:   PetscMemcpy(array, mArray, mSize * sizeof(double));
542:   MLDisownRealList(link, mArray, mSize);
543:   VecRestoreArray(v, &array);
544:   return(0);
545: }

547: /*@C
548:   PetscViewerMathematicaPutVector - Send a vector to Mathematica

550:   Input Parameters:
551: + viewer - The Mathematica viewer
552: - v      - The vector

554:   Level: intermediate

556: .keywords PetscViewer, Mathematica, vector
557: .seealso VecView(), PetscViewerMathematicaGetVector()
558: @*/
559: PetscErrorCode  PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
560: {
561:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
562:   MLINK                   link   = vmath->link; /* The link to Mathematica */
563:   char                    *name;
564:   PetscScalar             *array;
565:   int                     n;
566:   PetscErrorCode          ierr;

569:   /* Determine the object name */
570:   if (!vmath->objName) name = "vec";
571:   else                 name = (char*) vmath->objName;

573:   VecGetLocalSize(v, &n);
574:   VecGetArray(v, &array);

576:   /* Send the Vector object */
577:   MLPutFunction(link, "EvaluatePacket", 1);
578:   MLPutFunction(link, "Set", 2);
579:   MLPutSymbol(link, name);
580:   MLPutRealList(link, array, n);
581:   MLEndPacket(link);
582:   /* Skip packets until ReturnPacket */
583:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
584:   /* Skip ReturnPacket */
585:   MLNewPacket(link);

587:   VecRestoreArray(v, &array);
588:   return(0);
589: }

591: PetscErrorCode  PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
592: {
593:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
594:   MLINK                   link   = vmath->link; /* The link to Mathematica */
595:   char                    *name;
596:   PetscErrorCode          ierr;

599:   /* Determine the object name */
600:   if (!vmath->objName) name = "mat";
601:   else                 name = (char*) vmath->objName;

603:   /* Send the dense matrix object */
604:   MLPutFunction(link, "EvaluatePacket", 1);
605:   MLPutFunction(link, "Set", 2);
606:   MLPutSymbol(link, name);
607:   MLPutFunction(link, "Transpose", 1);
608:   MLPutFunction(link, "Partition", 2);
609:   MLPutRealList(link, a, m*n);
610:   MLPutInteger(link, m);
611:   MLEndPacket(link);
612:   /* Skip packets until ReturnPacket */
613:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
614:   /* Skip ReturnPacket */
615:   MLNewPacket(link);
616:   return(0);
617: }

619: PetscErrorCode  PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
620: {
621:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
622:   MLINK                   link   = vmath->link; /* The link to Mathematica */
623:   const char              *symbol;
624:   char                    *name;
625:   PetscBool               match;
626:   PetscErrorCode          ierr;

629:   /* Determine the object name */
630:   if (!vmath->objName) name = "mat";
631:   else                 name = (char*) vmath->objName;

633:   /* Make sure Mathematica recognizes sparse matrices */
634:   MLPutFunction(link, "EvaluatePacket", 1);
635:   MLPutFunction(link, "Needs", 1);
636:   MLPutString(link, "LinearAlgebra`CSRMatrix`");
637:   MLEndPacket(link);
638:   /* Skip packets until ReturnPacket */
639:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
640:   /* Skip ReturnPacket */
641:   MLNewPacket(link);

643:   /* Send the CSRMatrix object */
644:   MLPutFunction(link, "EvaluatePacket", 1);
645:   MLPutFunction(link, "Set", 2);
646:   MLPutSymbol(link, name);
647:   MLPutFunction(link, "CSRMatrix", 5);
648:   MLPutInteger(link, m);
649:   MLPutInteger(link, n);
650:   MLPutFunction(link, "Plus", 2);
651:   MLPutIntegerList(link, i, m+1);
652:   MLPutInteger(link, 1);
653:   MLPutFunction(link, "Plus", 2);
654:   MLPutIntegerList(link, j, i[m]);
655:   MLPutInteger(link, 1);
656:   MLPutRealList(link, a, i[m]);
657:   MLEndPacket(link);
658:   /* Skip packets until ReturnPacket */
659:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
660:   /* Skip ReturnPacket */
661:   MLNewPacket(link);

663:   /* Check that matrix is valid */
664:   MLPutFunction(link, "EvaluatePacket", 1);
665:   MLPutFunction(link, "ValidQ", 1);
666:   MLPutSymbol(link, name);
667:   MLEndPacket(link);
668:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
669:   MLGetSymbol(link, &symbol);
670:   PetscStrcmp("True", (char*) symbol, &match);
671:   if (!match) {
672:     MLDisownSymbol(link, symbol);
673:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
674:   }
675:   MLDisownSymbol(link, symbol);
676:   /* Skip ReturnPacket */
677:   MLNewPacket(link);
678:   return(0);
679: }