Actual source code: mathematica.c

petsc-3.3-p7 2013-05-11
  2: #include <petsc-private/viewerimpl.h>   /* "petscsys.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 = PETSC_NULL;
 12: static void *mathematicaEnv                   = PETSC_NULL;

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

 21:   Level: developer

 23: .keywords: Petsc, destroy, package, mathematica
 24: .seealso: PetscFinalize()
 25: @*/
 26: PetscErrorCode  PetscViewerMathematicaFinalizePackage(void)
 27: {
 29:   if (mathematicaEnv) MLDeinitialize((MLEnvironment) mathematicaEnv);
 30:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
 31:   return(0);
 32: }

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

 41:   Input Parameter:
 42:   path - The dynamic library path, or PETSC_NULL

 44:   Level: developer

 46: .keywords: Petsc, initialize, package, PLAPACK
 47: .seealso: PetscSysInitializePackage(), PetscInitialize()
 48: @*/
 49: PetscErrorCode  PetscViewerMathematicaInitializePackage(const char path[])
 50: {
 51:   PetscError ierr;

 54:   if (PetscViewerMathematicaPackageInitialized) return(0);
 55:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
 56:   mathematicaEnv = (void*) MLInitialize(0);
 57:   PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
 58:   return(0);
 59: }


 64: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
 65: {

 69:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return(0);
 70:   PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_NULL, PETSC_NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
 71:   return(0);
 72: }

 76: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
 77: {
 78:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
 79:   PetscErrorCode          ierr;

 82:   MLClose(vmath->link);
 83:   PetscFree(vmath->linkname);
 84:   PetscFree(vmath->linkhost);
 85:   PetscFree(vmath);
 86:   return(0);
 87: }

 91: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
 92: {

 96:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
 97:     PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
 98:   }
 99:   return(0);
100: }

104: PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
105: {
106:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
107: #ifdef MATHEMATICA_3_0
108:   int                     argc = 6;
109:   char                    *argv[6];
110: #else
111:   int                     argc = 5;
112:   char                    *argv[5];
113: #endif
114:   char                    hostname[256];
115:   long                    lerr;
116:   PetscErrorCode          ierr;

119:   /* Link name */
120:   argv[0] = "-linkname";
121:   if (!vmath->linkname) {
122:     argv[1] = "math -mathlink";
123:   } else {
124:     argv[1] = vmath->linkname;
125:   }

127:   /* Link host */
128:   argv[2] = "-linkhost";
129:   if (!vmath->linkhost) {
130:     PetscGetHostName(hostname, 255);
131:     argv[3] = hostname;
132:   } else {
133:     argv[3] = vmath->linkhost;
134:   }

136:   /* Link mode */
137: #ifdef MATHEMATICA_3_0
138:   argv[4] = "-linkmode";
139:   switch(vmath->linkmode) {
140:   case MATHEMATICA_LINK_CREATE:
141:     argv[5] = "Create";
142:     break;
143:   case MATHEMATICA_LINK_CONNECT:
144:     argv[5] = "Connect";
145:     break;
146:   case MATHEMATICA_LINK_LAUNCH:
147:     argv[5] = "Launch";
148:     break;
149:   }
150: #else
151:   switch(vmath->linkmode) {
152:   case MATHEMATICA_LINK_CREATE:
153:     argv[4] = "-linkcreate";
154:     break;
155:   case MATHEMATICA_LINK_CONNECT:
156:     argv[4] = "-linkconnect";
157:     break;
158:   case MATHEMATICA_LINK_LAUNCH:
159:     argv[4] = "-linklaunch";
160:     break;
161:   }
162: #endif
163:   vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
164: #endif
165:   return(0);
166: }

168: EXTERN_C_BEGIN
171: PetscErrorCode  PetscViewerCreate_Mathematica(PetscViewer v)
172: {
173:   PetscViewer_Mathematica *vmath;
174:   PetscErrorCode          ierr;

177: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
178:   PetscViewerMathematicaInitializePackage(PETSC_NULL);
179: #endif

181:   PetscNewLog(v,PetscViewer_Mathematica, &vmath);
182:   v->data         = (void*) vmath;
183:   v->ops->destroy = PetscViewerDestroy_Mathematica;
184:   v->ops->flush   = 0;
185:   PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);

187:   vmath->linkname         = PETSC_NULL;
188:   vmath->linkhost         = PETSC_NULL;
189:   vmath->linkmode         = MATHEMATICA_LINK_CONNECT;
190:   vmath->graphicsType     = GRAPHICS_MOTIF;
191:   vmath->plotType         = MATHEMATICA_TRIANGULATION_PLOT;
192:   vmath->objName          = PETSC_NULL;

194:   PetscViewerMathematicaSetFromOptions(v);
195:   PetscViewerMathematicaSetupConnection_Private(v);
196:   return(0);
197: }
198: EXTERN_C_END

202: PetscErrorCode PetscViewerMathematicaParseLinkMode_Private(char *modename, LinkMode *mode) {
203:   PetscBool      isCreate, isConnect, isLaunch;

207:   PetscStrcasecmp(modename, "Create",  &isCreate);
208:   PetscStrcasecmp(modename, "Connect", &isConnect);
209:   PetscStrcasecmp(modename, "Launch",  &isLaunch);
210:   if (isCreate) {
211:     *mode = MATHEMATICA_LINK_CREATE;
212:   } else if (isConnect) {
213:     *mode = MATHEMATICA_LINK_CONNECT;
214:   } else if (isLaunch) {
215:     *mode = MATHEMATICA_LINK_LAUNCH;
216:   } else {
217:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
218:   }
219:   return(0);
220: }

224: PetscErrorCode  PetscViewerMathematicaSetFromOptions(PetscViewer v)
225: {
226:   PetscViewer_Mathematica  *vmath = (PetscViewer_Mathematica *) v->data;
227:   char                     linkname[256];
228:   char                     modename[256];
229:   char                     hostname[256];
230:   char                     type[256];
231:   PetscInt                 numPorts;
232:   PetscInt                 *ports;
233:   PetscInt                 numHosts;
234:   int                      h;
235:   char                     **hosts;
236:   PetscMPIInt              size, rank;
237:   PetscBool                opt;
238:   PetscErrorCode           ierr;

241:   MPI_Comm_size(((PetscObject)v)->comm, &size);
242:   MPI_Comm_rank(((PetscObject)v)->comm, &rank);

244:   /* Get link name */
245:   PetscOptionsGetString("viewer_", "-math_linkname", linkname, 256, &opt);
246:   if (opt) {
247:     PetscViewerMathematicaSetLinkName(v, linkname);
248:   }
249:   /* Get link port */
250:   numPorts = size;
251:   PetscMalloc(size*sizeof(int), &ports);
252:   PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
253:   if (opt) {
254:     if (numPorts > rank) {
255:       snprintf(linkname, 255, "%6d", ports[rank]);
256:     } else {
257:       snprintf(linkname, 255, "%6d", ports[0]);
258:     }
259:     PetscViewerMathematicaSetLinkName(v, linkname);
260:   }
261:   PetscFree(ports);
262:   /* Get link host */
263:   numHosts = size;
264:   PetscMalloc(size*sizeof(char *), &hosts);
265:   PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
266:   if (opt) {
267:     if (numHosts > rank) {
268:       PetscStrncpy(hostname, hosts[rank], 255);
269:     } else {
270:       PetscStrncpy(hostname, hosts[0], 255);
271:     }
272:     PetscViewerMathematicaSetLinkHost(v, hostname);
273:   }
274:   for(h = 0; h < numHosts; h++) {
275:     PetscFree(hosts[h]);
276:   }
277:   PetscFree(hosts);
278:   /* Get link mode */
279:   PetscOptionsGetString("viewer_", "-math_linkmode", modename, 256, &opt);
280:   if (opt) {
281:     LinkMode mode;

283:     PetscViewerMathematicaParseLinkMode_Private(modename, &mode);
284:     PetscViewerMathematicaSetLinkMode(v, mode);
285:   }
286:   /* Get graphics type */
287:   PetscOptionsGetString("viewer_", "-math_graphics", type, 256, &opt);
288:   if (opt) {
289:     PetscBool  isMotif, isPS, isPSFile;

291:     PetscStrcasecmp(type, "Motif",  &isMotif);
292:     PetscStrcasecmp(type, "PS",     &isPS);
293:     PetscStrcasecmp(type, "PSFile", &isPSFile);
294:     if (isMotif) {
295:       vmath->graphicsType = GRAPHICS_MOTIF;
296:     } else if (isPS) {
297:       vmath->graphicsType = GRAPHICS_PS_STDOUT;
298:     } else if (isPSFile) {
299:       vmath->graphicsType = GRAPHICS_PS_FILE;
300:     }
301:   }
302:   /* Get plot type */
303:   PetscOptionsGetString("viewer_", "-math_type", type, 256, &opt);
304:   if (opt) {
305:     PetscBool  isTri, isVecTri, isVec, isSurface;

307:     PetscStrcasecmp(type, "Triangulation",       &isTri);
308:     PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
309:     PetscStrcasecmp(type, "Vector",              &isVec);
310:     PetscStrcasecmp(type, "Surface",             &isSurface);
311:     if (isTri) {
312:       vmath->plotType     = MATHEMATICA_TRIANGULATION_PLOT;
313:     } else if (isVecTri) {
314:       vmath->plotType     = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
315:     } else if (isVec) {
316:       vmath->plotType     = MATHEMATICA_VECTOR_PLOT;
317:     } else if (isSurface) {
318:       vmath->plotType     = MATHEMATICA_SURFACE_PLOT;
319:     }
320:   }
321:   return(0);
322: }

326: PetscErrorCode  PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name) {
327:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
328:   PetscErrorCode          ierr;

333:   PetscStrallocpy(name, &vmath->linkname);
334:   return(0);
335: }

339: PetscErrorCode  PetscViewerMathematicaSetLinkPort(PetscViewer v, int port) {
340:   char           name[16];

344:   snprintf(name, 16, "%6d", port);
345:   PetscViewerMathematicaSetLinkName(v, name);
346:   return(0);
347: }

351: PetscErrorCode  PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host) {
352:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
353:   PetscErrorCode          ierr;

358:   PetscStrallocpy(host, &vmath->linkhost);
359:   return(0);
360: }

364: PetscErrorCode  PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode) {
365:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;

368:   vmath->linkmode = mode;
369:   return(0);
370: }

372: /*----------------------------------------- Public Functions --------------------------------------------------------*/
375: /*@C
376:   PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.

378:   Collective on comm

380:   Input Parameters:
381: + comm    - The MPI communicator
382: . port    - [optional] The port to connect on, or PETSC_DECIDE
383: . machine - [optional] The machine to run Mathematica on, or PETSC_NULL
384: - mode    - [optional] The connection mode, or PETSC_NULL

386:   Output Parameter:
387: . viewer  - The Mathematica viewer

389:   Level: intermediate

391:   Notes:
392:   Most users should employ the following commands to access the 
393:   Mathematica viewers
394: $
395: $    PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
396: $    MatView(Mat matrix, PetscViewer viewer)
397: $
398: $                or
399: $
400: $    PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
401: $    VecView(Vec vector, PetscViewer viewer)

403:    Options Database Keys:
404: $    -viewer_math_linkhost <machine> - The host machine for the kernel
405: $    -viewer_math_linkname <name>    - The full link name for the connection
406: $    -viewer_math_linkport <port>    - The port for the connection
407: $    -viewer_math_mode <mode>        - The mode, e.g. Launch, Connect
408: $    -viewer_math_type <type>        - The plot type, e.g. Triangulation, Vector
409: $    -viewer_math_graphics <output>  - The output type, e.g. Motif, PS, PSFile

411: .keywords: PetscViewer, Mathematica, open

413: .seealso: MatView(), VecView()
414: @*/
415: PetscErrorCode  PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
416: {

420:   PetscViewerCreate(comm, v);
421: #if 0
422:   LinkMode linkmode;
423:   PetscViewerMathematicaSetLinkPort(*v, port);
424:   PetscViewerMathematicaSetLinkHost(*v, machine);
425:   PetscViewerMathematicaParseLinkMode_Private(mode, &linkmode);
426:   PetscViewerMathematicaSetLinkMode(*v, linkmode);
427: #endif
428:   PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
429:   return(0);
430: }

434: /*@C
435:   PetscViewerMathematicaGetLink - Returns the link to Mathematica

437:   Input Parameters:
438: . viewer - The Mathematica viewer
439: . link   - The link to Mathematica

441:   Level: intermediate

443: .keywords PetscViewer, Mathematica, link
444: .seealso PetscViewerMathematicaOpen()
445: @*/
446: PetscErrorCode  PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
447: {
448:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;

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

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

461:   Input Parameters:
462: . viewer - The Mathematica viewer
463: . type   - The packet type to search for, e.g RETURNPKT

465:   Level: advanced

467: .keywords PetscViewer, Mathematica, packets
468: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
469: @*/
470: PetscErrorCode  PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
471: {
472:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
473:   MLINK                   link  = vmath->link; /* The link to Mathematica */
474:   int                     pkt;                 /* The packet type */

477:   while((pkt = MLNextPacket(link)) && (pkt != type))
478:     MLNewPacket(link);
479:   if (!pkt) {
480:     MLClearError(link);
481:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char *) MLErrorMessage(link));
482:   }
483:   return(0);
484: }

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

491:   Input Parameter:
492: . viewer - The Mathematica viewer

494:   Output Parameter:
495: . name   - The name for new objects created in Mathematica

497:   Level: intermediate

499: .keywords PetscViewer, Mathematica, name
500: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
501: @*/
502: PetscErrorCode  PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
503: {
504:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;

509:   *name = vmath->objName;
510:   return(0);
511: }

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

518:   Input Parameters:
519: . viewer - The Mathematica viewer
520: . name   - The name for new objects created in Mathematica

522:   Level: intermediate

524: .keywords PetscViewer, Mathematica, name
525: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
526: @*/
527: PetscErrorCode  PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
528: {
529:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;

534:   vmath->objName = name;
535:   return(0);
536: }

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

543:   Input Parameter:
544: . viewer - The Mathematica viewer

546:   Level: intermediate

548: .keywords PetscViewer, Mathematica, name
549: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
550: @*/
551: PetscErrorCode  PetscViewerMathematicaClearName(PetscViewer viewer)
552: {
553:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;

557:   vmath->objName = PETSC_NULL;
558:   return(0);
559: }

563: /*@C
564:   PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica

566:   Input Parameter:
567: . viewer - The Mathematica viewer

569:   Output Parameter:
570: . v      - The vector

572:   Level: intermediate

574: .keywords PetscViewer, Mathematica, vector
575: .seealso VecView(), PetscViewerMathematicaPutVector()
576: @*/
577: PetscErrorCode  PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v) {
578:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
579:   MLINK                   link;   /* The link to Mathematica */
580:   char                    *name;
581:   PetscScalar             *mArray,*array;
582:   long                    mSize;
583:   int                     n;
584:   PetscErrorCode          ierr;


590:   /* Determine the object name */
591:   if (!vmath->objName) {
592:     name = "vec";
593:   } else {
594:     name = (char *) vmath->objName;
595:   }

597:   link = vmath->link;
598:   VecGetLocalSize(v, &n);
599:   VecGetArray(v, &array);
600:   MLPutFunction(link, "EvaluatePacket", 1);
601:     MLPutSymbol(link, name);
602:   MLEndPacket(link);
603:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
604:   MLGetRealList(link, &mArray, &mSize);
605:   if (n != mSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d",n,mSize);
606:   PetscMemcpy(array, mArray, mSize * sizeof(double));
607:   MLDisownRealList(link, mArray, mSize);
608:   VecRestoreArray(v, &array);

610:   return(0);
611: }

615: /*@C
616:   PetscViewerMathematicaPutVector - Send a vector to Mathematica

618:   Input Parameters:
619: + viewer - The Mathematica viewer
620: - v      - The vector

622:   Level: intermediate

624: .keywords PetscViewer, Mathematica, vector
625: .seealso VecView(), PetscViewerMathematicaGetVector()
626: @*/
627: PetscErrorCode  PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
628: {
629:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
630:   MLINK                   link  = vmath->link; /* The link to Mathematica */
631:   char                    *name;
632:   PetscScalar             *array;
633:   int                     n;
634:   PetscErrorCode          ierr;

637:   /* Determine the object name */
638:   if (!vmath->objName) {
639:     name = "vec";
640:   } else {
641:     name = (char *) vmath->objName;
642:   }
643:   VecGetLocalSize(v, &n);
644:   VecGetArray(v, &array);

646:   /* Send the Vector object */
647:   MLPutFunction(link, "EvaluatePacket", 1);
648:     MLPutFunction(link, "Set", 2);
649:       MLPutSymbol(link, name);
650:       MLPutRealList(link, array, n);
651:   MLEndPacket(link);
652:   /* Skip packets until ReturnPacket */
653:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
654:   /* Skip ReturnPacket */
655:   MLNewPacket(link);

657:   VecRestoreArray(v, &array);
658:   return(0);
659: }

661: PetscErrorCode  PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
662: {
663:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
664:   MLINK                   link  = vmath->link; /* The link to Mathematica */
665:   char                    *name;
666:   PetscErrorCode          ierr;

669:   /* Determine the object name */
670:   if (!vmath->objName) {
671:     name = "mat";
672:   } else {
673:     name = (char *) vmath->objName;
674:   }

676:   /* Send the dense matrix object */
677:   MLPutFunction(link, "EvaluatePacket", 1);
678:     MLPutFunction(link, "Set", 2);
679:       MLPutSymbol(link, name);
680:       MLPutFunction(link, "Transpose", 1);
681:         MLPutFunction(link, "Partition", 2);
682:           MLPutRealList(link, a, m*n);
683:           MLPutInteger(link, m);
684:   MLEndPacket(link);
685:   /* Skip packets until ReturnPacket */
686:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
687:   /* Skip ReturnPacket */
688:   MLNewPacket(link);

690:   return(0);
691: }

693: PetscErrorCode  PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
694: {
695:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
696:   MLINK                   link  = vmath->link; /* The link to Mathematica */
697:   const char              *symbol;
698:   char                    *name;
699:   PetscBool               match;
700:   PetscErrorCode          ierr;

703:   /* Determine the object name */
704:   if (!vmath->objName) {
705:     name = "mat";
706:   } else {
707:     name = (char *) vmath->objName;
708:   }

710:   /* Make sure Mathematica recognizes sparse matrices */
711:   MLPutFunction(link, "EvaluatePacket", 1);
712:     MLPutFunction(link, "Needs", 1);
713:       MLPutString(link, "LinearAlgebra`CSRMatrix`");
714:   MLEndPacket(link);
715:   /* Skip packets until ReturnPacket */
716:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
717:   /* Skip ReturnPacket */
718:   MLNewPacket(link);

720:   /* Send the CSRMatrix object */
721:   MLPutFunction(link, "EvaluatePacket", 1);
722:     MLPutFunction(link, "Set", 2);
723:       MLPutSymbol(link, name);
724:       MLPutFunction(link, "CSRMatrix", 5);
725:         MLPutInteger(link, m);
726:         MLPutInteger(link, n);
727:         MLPutFunction(link, "Plus", 2);
728:           MLPutIntegerList(link, i, m+1);
729:           MLPutInteger(link, 1);
730:         MLPutFunction(link, "Plus", 2);
731:           MLPutIntegerList(link, j, i[m]);
732:           MLPutInteger(link, 1);
733:         MLPutRealList(link, a, i[m]);
734:   MLEndPacket(link);
735:   /* Skip packets until ReturnPacket */
736:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
737:   /* Skip ReturnPacket */
738:   MLNewPacket(link);

740:   /* Check that matrix is valid */
741:   MLPutFunction(link, "EvaluatePacket", 1);
742:     MLPutFunction(link, "ValidQ", 1);
743:       MLPutSymbol(link, name);
744:   MLEndPacket(link);
745:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
746:   MLGetSymbol(link, &symbol);
747:   PetscStrcmp("True", (char *) symbol, &match);
748:   if (!match) {
749:     MLDisownSymbol(link, symbol);
750:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
751:   }
752:   MLDisownSymbol(link, symbol);
753:   /* Skip ReturnPacket */
754:   MLNewPacket(link);

756:   return(0);
757: }