Actual source code: mathematica.c

  1: #define PETSC_DLL
  2: /* 
  3:         Written by Matt Knepley, knepley@cs.purdue.edu 7/23/97
  4:         Major overhall for interactivity               11/14/97
  5:         Reorganized                                    11/8/98
  6: */
 7:  #include src/sys/viewer/viewerimpl.h
 8:  #include private/pcimpl.h
 9:  #include src/mat/impls/aij/seq/aij.h
 10:  #include mathematica.h
 11: #include "petscfix.h"

 13: #if defined (PETSC_HAVE__SNPRINTF) && !defined(PETSC_HAVE_SNPRINTF)
 14: #define snprintf _snprintf
 15: #endif

 17: PetscViewer  PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE = PETSC_NULL;
 18: static void *mathematicaEnv                   = PETSC_NULL;

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

 27:   Input Parameter:
 28:   path - The dynamic library path, or PETSC_NULL

 30:   Level: developer

 32: .keywords: Petsc, initialize, package, PLAPACK
 33: .seealso: PetscInitializePackage(), PetscInitialize()
 34: @*/
 35: PetscErrorCode  PetscViewerMathematicaInitializePackage(char *path)
 36: {
 37:   static PetscTruth initialized = PETSC_FALSE;

 40:   if (initialized) return(0);
 41:   initialized = PETSC_TRUE;
 42:   mathematicaEnv = (void*) MLInitialize(0);
 43:   return(0);
 44: }

 48: /*@C
 49:   PetscViewerMathematicaDestroyPackage - This function destroys everything in the Petsc interface to Mathematica. It is
 50:   called from PetscFinalize().

 52:   Level: developer

 54: .keywords: Petsc, destroy, package, mathematica
 55: .seealso: PetscFinalize()
 56: @*/
 57: PetscErrorCode  PetscViewerMathematicaFinalizePackage(void)
 58: {
 60:   if (mathematicaEnv) MLDeinitialize((MLEnvironment) mathematicaEnv);
 61:   return(0);
 62: }

 66: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
 67: {

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

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

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

 93: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
 94: {

 98:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
 99:     PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
100:   }
101:   return(0);
102: }

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

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

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

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

173: PetscErrorCode  PetscViewerCreate_Mathematica(PetscViewer v)
174: {
175:   PetscViewer_Mathematica *vmath;
176:   PetscErrorCode          ierr;


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

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

193:   PetscViewerMathematicaSetFromOptions(v);
194:   PetscViewerMathematicaSetupConnection_Private(v);
195:   return(0);
196: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

377:   Collective on comm

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

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

388:   Level: intermediate

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

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

410: .keywords: PetscViewer, Mathematica, open

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

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

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

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

440:   Level: intermediate

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

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

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

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

464:   Level: advanced

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

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

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

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

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

496:   Level: intermediate

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

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

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

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

521:   Level: intermediate

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

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

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

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

545:   Level: intermediate

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

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

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

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

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

571:   Level: intermediate

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


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

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

609:   return(0);
610: }

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

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

621:   Level: intermediate

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

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

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

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

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

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

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

689:   return(0);
690: }

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

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

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

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

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

755:   return(0);
756: }