Actual source code: mathematica.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  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 = NULL;
 12: static void *mathematicaEnv                        = 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:   Level: developer

 43: .keywords: Petsc, initialize, package
 44: .seealso: PetscSysInitializePackage(), PetscInitialize()
 45: @*/
 46: PetscErrorCode  PetscViewerMathematicaInitializePackage(void)
 47: {
 48:   PetscError ierr;

 51:   if (PetscViewerMathematicaPackageInitialized) return(0);
 52:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;

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

 56:   PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
 57:   return(0);
 58: }


 63: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
 64: {

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

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

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

 90: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
 91: {

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

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

118:   /* Link name */
119:   argv[0] = "-linkname";
120:   if (!vmath->linkname) argv[1] = "math -mathlink";
121:   else                  argv[1] = vmath->linkname;

123:   /* Link host */
124:   argv[2] = "-linkhost";
125:   if (!vmath->linkhost) {
126:     PetscGetHostName(hostname, 255);
127:     argv[3] = hostname;
128:   } else argv[3] = vmath->linkhost;

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

164: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
165: {
166:   PetscViewer_Mathematica *vmath;
167:   PetscErrorCode          ierr;

170:   PetscViewerMathematicaInitializePackage();

172:   PetscNewLog(v,&vmath);
173:   v->data         = (void*) vmath;
174:   v->ops->destroy = PetscViewerDestroy_Mathematica;
175:   v->ops->flush   = 0;
176:   PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);

178:   vmath->linkname     = NULL;
179:   vmath->linkhost     = NULL;
180:   vmath->linkmode     = MATHEMATICA_LINK_CONNECT;
181:   vmath->graphicsType = GRAPHICS_MOTIF;
182:   vmath->plotType     = MATHEMATICA_TRIANGULATION_PLOT;
183:   vmath->objName      = NULL;

185:   PetscViewerMathematicaSetFromOptions(v);
186:   PetscViewerMathematicaSetupConnection_Private(v);
187:   return(0);
188: }

192: PetscErrorCode PetscViewerMathematicaParseLinkMode_Private(char *modename, LinkMode *mode)
193: {
194:   PetscBool      isCreate, isConnect, isLaunch;

198:   PetscStrcasecmp(modename, "Create",  &isCreate);
199:   PetscStrcasecmp(modename, "Connect", &isConnect);
200:   PetscStrcasecmp(modename, "Launch",  &isLaunch);
201:   if (isCreate)       *mode = MATHEMATICA_LINK_CREATE;
202:   else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
203:   else if (isLaunch)  *mode = MATHEMATICA_LINK_LAUNCH;
204:   else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
205:   return(0);
206: }

210: PetscErrorCode  PetscViewerMathematicaSetFromOptions(PetscViewer v)
211: {
212:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
213:   char                    linkname[256];
214:   char                    modename[256];
215:   char                    hostname[256];
216:   char                    type[256];
217:   PetscInt                numPorts;
218:   PetscInt                *ports;
219:   PetscInt                numHosts;
220:   int                     h;
221:   char                    **hosts;
222:   PetscMPIInt             size, rank;
223:   PetscBool               opt;
224:   PetscErrorCode          ierr;

227:   MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
228:   MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);

230:   /* Get link name */
231:   PetscOptionsGetString("viewer_", "-math_linkname", linkname, 256, &opt);
232:   if (opt) {
233:     PetscViewerMathematicaSetLinkName(v, linkname);
234:   }
235:   /* Get link port */
236:   numPorts = size;
237:   PetscMalloc1(size, &ports);
238:   PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
239:   if (opt) {
240:     if (numPorts > rank) snprintf(linkname, 255, "%6d", ports[rank]);
241:     else                 snprintf(linkname, 255, "%6d", ports[0]);
242:     PetscViewerMathematicaSetLinkName(v, linkname);
243:   }
244:   PetscFree(ports);
245:   /* Get link host */
246:   numHosts = size;
247:   PetscMalloc1(size, &hosts);
248:   PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
249:   if (opt) {
250:     if (numHosts > rank) {
251:       PetscStrncpy(hostname, hosts[rank], 255);
252:     } else {
253:       PetscStrncpy(hostname, hosts[0], 255);
254:     }
255:     PetscViewerMathematicaSetLinkHost(v, hostname);
256:   }
257:   for (h = 0; h < numHosts; h++) {
258:     PetscFree(hosts[h]);
259:   }
260:   PetscFree(hosts);
261:   /* Get link mode */
262:   PetscOptionsGetString("viewer_", "-math_linkmode", modename, 256, &opt);
263:   if (opt) {
264:     LinkMode mode;

266:     PetscViewerMathematicaParseLinkMode_Private(modename, &mode);
267:     PetscViewerMathematicaSetLinkMode(v, mode);
268:   }
269:   /* Get graphics type */
270:   PetscOptionsGetString("viewer_", "-math_graphics", type, 256, &opt);
271:   if (opt) {
272:     PetscBool isMotif, isPS, isPSFile;

274:     PetscStrcasecmp(type, "Motif",  &isMotif);
275:     PetscStrcasecmp(type, "PS",     &isPS);
276:     PetscStrcasecmp(type, "PSFile", &isPSFile);
277:     if (isMotif)       vmath->graphicsType = GRAPHICS_MOTIF;
278:     else if (isPS)     vmath->graphicsType = GRAPHICS_PS_STDOUT;
279:     else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
280:   }
281:   /* Get plot type */
282:   PetscOptionsGetString("viewer_", "-math_type", type, 256, &opt);
283:   if (opt) {
284:     PetscBool isTri, isVecTri, isVec, isSurface;

286:     PetscStrcasecmp(type, "Triangulation",       &isTri);
287:     PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
288:     PetscStrcasecmp(type, "Vector",              &isVec);
289:     PetscStrcasecmp(type, "Surface",             &isSurface);
290:     if (isTri)          vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
291:     else if (isVecTri)  vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
292:     else if (isVec)     vmath->plotType = MATHEMATICA_VECTOR_PLOT;
293:     else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
294:   }
295:   return(0);
296: }

300: PetscErrorCode  PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
301: {
302:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
303:   PetscErrorCode          ierr;

308:   PetscStrallocpy(name, &vmath->linkname);
309:   return(0);
310: }

314: PetscErrorCode  PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
315: {
316:   char           name[16];

320:   snprintf(name, 16, "%6d", port);
321:   PetscViewerMathematicaSetLinkName(v, name);
322:   return(0);
323: }

327: PetscErrorCode  PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
328: {
329:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
330:   PetscErrorCode          ierr;

335:   PetscStrallocpy(host, &vmath->linkhost);
336:   return(0);
337: }

341: PetscErrorCode  PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
342: {
343:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;

346:   vmath->linkmode = mode;
347:   return(0);
348: }

350: /*----------------------------------------- Public Functions --------------------------------------------------------*/
353: /*@C
354:   PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.

356:   Collective on comm

358:   Input Parameters:
359: + comm    - The MPI communicator
360: . port    - [optional] The port to connect on, or PETSC_DECIDE
361: . machine - [optional] The machine to run Mathematica on, or NULL
362: - mode    - [optional] The connection mode, or NULL

364:   Output Parameter:
365: . viewer  - The Mathematica viewer

367:   Level: intermediate

369:   Notes:
370:   Most users should employ the following commands to access the
371:   Mathematica viewers
372: $
373: $    PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
374: $    MatView(Mat matrix, PetscViewer viewer)
375: $
376: $                or
377: $
378: $    PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
379: $    VecView(Vec vector, PetscViewer viewer)

381:    Options Database Keys:
382: $    -viewer_math_linkhost <machine> - The host machine for the kernel
383: $    -viewer_math_linkname <name>    - The full link name for the connection
384: $    -viewer_math_linkport <port>    - The port for the connection
385: $    -viewer_math_mode <mode>        - The mode, e.g. Launch, Connect
386: $    -viewer_math_type <type>        - The plot type, e.g. Triangulation, Vector
387: $    -viewer_math_graphics <output>  - The output type, e.g. Motif, PS, PSFile

389: .keywords: PetscViewer, Mathematica, open

391: .seealso: MatView(), VecView()
392: @*/
393: PetscErrorCode  PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
394: {

398:   PetscViewerCreate(comm, v);
399: #if 0
400:   LinkMode linkmode;
401:   PetscViewerMathematicaSetLinkPort(*v, port);
402:   PetscViewerMathematicaSetLinkHost(*v, machine);
403:   PetscViewerMathematicaParseLinkMode_Private(mode, &linkmode);
404:   PetscViewerMathematicaSetLinkMode(*v, linkmode);
405: #endif
406:   PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
407:   return(0);
408: }

412: /*@C
413:   PetscViewerMathematicaGetLink - Returns the link to Mathematica

415:   Input Parameters:
416: . viewer - The Mathematica viewer
417: . link   - The link to Mathematica

419:   Level: intermediate

421: .keywords PetscViewer, Mathematica, link
422: .seealso PetscViewerMathematicaOpen()
423: @*/
424: PetscErrorCode  PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
425: {
426:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

430:   *link = vmath->link;
431:   return(0);
432: }

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

439:   Input Parameters:
440: . viewer - The Mathematica viewer
441: . type   - The packet type to search for, e.g RETURNPKT

443:   Level: advanced

445: .keywords PetscViewer, Mathematica, packets
446: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
447: @*/
448: PetscErrorCode  PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
449: {
450:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
451:   MLINK                   link   = vmath->link; /* The link to Mathematica */
452:   int                     pkt;                 /* The packet type */

455:   while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
456:   if (!pkt) {
457:     MLClearError(link);
458:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char*) MLErrorMessage(link));
459:   }
460:   return(0);
461: }

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

468:   Input Parameter:
469: . viewer - The Mathematica viewer

471:   Output Parameter:
472: . name   - The name for new objects created in Mathematica

474:   Level: intermediate

476: .keywords PetscViewer, Mathematica, name
477: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
478: @*/
479: PetscErrorCode  PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
480: {
481:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

486:   *name = vmath->objName;
487:   return(0);
488: }

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

495:   Input Parameters:
496: . viewer - The Mathematica viewer
497: . name   - The name for new objects created in Mathematica

499:   Level: intermediate

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

511:   vmath->objName = name;
512:   return(0);
513: }

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

520:   Input Parameter:
521: . viewer - The Mathematica viewer

523:   Level: intermediate

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

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

540: /*@C
541:   PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica

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

546:   Output Parameter:
547: . v      - The vector

549:   Level: intermediate

551: .keywords PetscViewer, Mathematica, vector
552: .seealso VecView(), PetscViewerMathematicaPutVector()
553: @*/
554: PetscErrorCode  PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
555: {
556:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
557:   MLINK                   link;   /* The link to Mathematica */
558:   char                    *name;
559:   PetscScalar             *mArray,*array;
560:   long                    mSize;
561:   int                     n;
562:   PetscErrorCode          ierr;


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

572:   link = vmath->link;
573:   VecGetLocalSize(v, &n);
574:   VecGetArray(v, &array);
575:   MLPutFunction(link, "EvaluatePacket", 1);
576:   MLPutSymbol(link, name);
577:   MLEndPacket(link);
578:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
579:   MLGetRealList(link, &mArray, &mSize);
580:   if (n != mSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d",n,mSize);
581:   PetscMemcpy(array, mArray, mSize * sizeof(double));
582:   MLDisownRealList(link, mArray, mSize);
583:   VecRestoreArray(v, &array);
584:   return(0);
585: }

589: /*@C
590:   PetscViewerMathematicaPutVector - Send a vector to Mathematica

592:   Input Parameters:
593: + viewer - The Mathematica viewer
594: - v      - The vector

596:   Level: intermediate

598: .keywords PetscViewer, Mathematica, vector
599: .seealso VecView(), PetscViewerMathematicaGetVector()
600: @*/
601: PetscErrorCode  PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
602: {
603:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
604:   MLINK                   link   = vmath->link; /* The link to Mathematica */
605:   char                    *name;
606:   PetscScalar             *array;
607:   int                     n;
608:   PetscErrorCode          ierr;

611:   /* Determine the object name */
612:   if (!vmath->objName) name = "vec";
613:   else                 name = (char*) vmath->objName;

615:   VecGetLocalSize(v, &n);
616:   VecGetArray(v, &array);

618:   /* Send the Vector object */
619:   MLPutFunction(link, "EvaluatePacket", 1);
620:   MLPutFunction(link, "Set", 2);
621:   MLPutSymbol(link, name);
622:   MLPutRealList(link, array, n);
623:   MLEndPacket(link);
624:   /* Skip packets until ReturnPacket */
625:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
626:   /* Skip ReturnPacket */
627:   MLNewPacket(link);

629:   VecRestoreArray(v, &array);
630:   return(0);
631: }

633: PetscErrorCode  PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
634: {
635:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
636:   MLINK                   link   = vmath->link; /* The link to Mathematica */
637:   char                    *name;
638:   PetscErrorCode          ierr;

641:   /* Determine the object name */
642:   if (!vmath->objName) name = "mat";
643:   else                 name = (char*) vmath->objName;

645:   /* Send the dense matrix object */
646:   MLPutFunction(link, "EvaluatePacket", 1);
647:   MLPutFunction(link, "Set", 2);
648:   MLPutSymbol(link, name);
649:   MLPutFunction(link, "Transpose", 1);
650:   MLPutFunction(link, "Partition", 2);
651:   MLPutRealList(link, a, m*n);
652:   MLPutInteger(link, m);
653:   MLEndPacket(link);
654:   /* Skip packets until ReturnPacket */
655:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
656:   /* Skip ReturnPacket */
657:   MLNewPacket(link);
658:   return(0);
659: }

661: PetscErrorCode  PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
662: {
663:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
664:   MLINK                   link   = vmath->link; /* The link to Mathematica */
665:   const char              *symbol;
666:   char                    *name;
667:   PetscBool               match;
668:   PetscErrorCode          ierr;

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

675:   /* Make sure Mathematica recognizes sparse matrices */
676:   MLPutFunction(link, "EvaluatePacket", 1);
677:   MLPutFunction(link, "Needs", 1);
678:   MLPutString(link, "LinearAlgebra`CSRMatrix`");
679:   MLEndPacket(link);
680:   /* Skip packets until ReturnPacket */
681:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
682:   /* Skip ReturnPacket */
683:   MLNewPacket(link);

685:   /* Send the CSRMatrix object */
686:   MLPutFunction(link, "EvaluatePacket", 1);
687:   MLPutFunction(link, "Set", 2);
688:   MLPutSymbol(link, name);
689:   MLPutFunction(link, "CSRMatrix", 5);
690:   MLPutInteger(link, m);
691:   MLPutInteger(link, n);
692:   MLPutFunction(link, "Plus", 2);
693:   MLPutIntegerList(link, i, m+1);
694:   MLPutInteger(link, 1);
695:   MLPutFunction(link, "Plus", 2);
696:   MLPutIntegerList(link, j, i[m]);
697:   MLPutInteger(link, 1);
698:   MLPutRealList(link, a, i[m]);
699:   MLEndPacket(link);
700:   /* Skip packets until ReturnPacket */
701:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
702:   /* Skip ReturnPacket */
703:   MLNewPacket(link);

705:   /* Check that matrix is valid */
706:   MLPutFunction(link, "EvaluatePacket", 1);
707:   MLPutFunction(link, "ValidQ", 1);
708:   MLPutSymbol(link, name);
709:   MLEndPacket(link);
710:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
711:   MLGetSymbol(link, &symbol);
712:   PetscStrcmp("True", (char*) symbol, &match);
713:   if (!match) {
714:     MLDisownSymbol(link, symbol);
715:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
716:   }
717:   MLDisownSymbol(link, symbol);
718:   /* Skip ReturnPacket */
719:   MLNewPacket(link);
720:   return(0);
721: }