Actual source code: mathematica.c


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

 31: /*@C
 32:   PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
 33:   called from PetscViewerInitializePackage().

 35:   Level: developer

 37: .seealso: PetscSysInitializePackage(), PetscInitialize()
 38: @*/
 39: PetscErrorCode  PetscViewerMathematicaInitializePackage(void)
 40: {
 41:   PetscError ierr;

 44:   if (PetscViewerMathematicaPackageInitialized) return(0);
 45:   PetscViewerMathematicaPackageInitialized = PETSC_TRUE;

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

 49:   PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
 50:   return(0);
 51: }

 53: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
 54: {

 58:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return(0);
 59:   PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
 60:   return(0);
 61: }

 63: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
 64: {
 65:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
 66:   PetscErrorCode          ierr;

 69:   MLClose(vmath->link);
 70:   PetscFree(vmath->linkname);
 71:   PetscFree(vmath->linkhost);
 72:   PetscFree(vmath);
 73:   return(0);
 74: }

 76: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
 77: {

 81:   if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
 82:     PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
 83:   }
 84:   return(0);
 85: }

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

102:   /* Link name */
103:   argv[0] = "-linkname";
104:   if (!vmath->linkname) argv[1] = "math -mathlink";
105:   else                  argv[1] = vmath->linkname;

107:   /* Link host */
108:   argv[2] = "-linkhost";
109:   if (!vmath->linkhost) {
110:     PetscGetHostName(hostname, sizeof(hostname));
111:     argv[3] = hostname;
112:   } else argv[3] = vmath->linkhost;

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

146: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
147: {
148:   PetscViewer_Mathematica *vmath;
149:   PetscErrorCode          ierr;

152:   PetscViewerMathematicaInitializePackage();

154:   PetscNewLog(v,&vmath);
155:   v->data         = (void*) vmath;
156:   v->ops->destroy = PetscViewerDestroy_Mathematica;
157:   v->ops->flush   = 0;
158:   PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);

160:   vmath->linkname     = NULL;
161:   vmath->linkhost     = NULL;
162:   vmath->linkmode     = MATHEMATICA_LINK_CONNECT;
163:   vmath->graphicsType = GRAPHICS_MOTIF;
164:   vmath->plotType     = MATHEMATICA_TRIANGULATION_PLOT;
165:   vmath->objName      = NULL;

167:   PetscViewerMathematicaSetFromOptions(v);
168:   PetscViewerMathematicaSetupConnection_Private(v);
169:   return(0);
170: }

172: static PetscErrorCode PetscViewerMathematicaParseLinkMode(char *modename, LinkMode *mode)
173: {
174:   PetscBool      isCreate, isConnect, isLaunch;

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

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

205:   MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
206:   MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);

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

244:     PetscViewerMathematicaParseLinkMode(modename, &mode);
245:     PetscViewerMathematicaSetLinkMode(v, mode);
246:   }
247:   /* Get graphics type */
248:   PetscOptionsGetString("viewer_", "-math_graphics", type, sizeof(type), &opt);
249:   if (opt) {
250:     PetscBool isMotif, isPS, isPSFile;

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

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

276: PetscErrorCode  PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
277: {
278:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
279:   PetscErrorCode          ierr;

284:   PetscStrallocpy(name, &vmath->linkname);
285:   return(0);
286: }

288: PetscErrorCode  PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
289: {
290:   char           name[16];

294:   snprintf(name, 16, "%6d", port);
295:   PetscViewerMathematicaSetLinkName(v, name);
296:   return(0);
297: }

299: PetscErrorCode  PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
300: {
301:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
302:   PetscErrorCode          ierr;

307:   PetscStrallocpy(host, &vmath->linkhost);
308:   return(0);
309: }

311: PetscErrorCode  PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
312: {
313:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;

316:   vmath->linkmode = mode;
317:   return(0);
318: }

320: /*----------------------------------------- Public Functions --------------------------------------------------------*/
321: /*@C
322:   PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.

324:   Collective

326:   Input Parameters:
327: + comm    - The MPI communicator
328: . port    - [optional] The port to connect on, or PETSC_DECIDE
329: . machine - [optional] The machine to run Mathematica on, or NULL
330: - mode    - [optional] The connection mode, or NULL

332:   Output Parameter:
333: . viewer  - The Mathematica viewer

335:   Level: intermediate

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

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

357: .seealso: MatView(), VecView()
358: @*/
359: PetscErrorCode  PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
360: {

364:   PetscViewerCreate(comm, v);
365: #if 0
366:   LinkMode linkmode;
367:   PetscViewerMathematicaSetLinkPort(*v, port);
368:   PetscViewerMathematicaSetLinkHost(*v, machine);
369:   PetscViewerMathematicaParseLinkMode(mode, &linkmode);
370:   PetscViewerMathematicaSetLinkMode(*v, linkmode);
371: #endif
372:   PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
373:   return(0);
374: }

376: /*@C
377:   PetscViewerMathematicaGetLink - Returns the link to Mathematica

379:   Input Parameters:
380: + viewer - The Mathematica viewer
381: - link   - The link to Mathematica

383:   Level: intermediate

385: .keywords PetscViewer, Mathematica, link
386: .seealso PetscViewerMathematicaOpen()
387: @*/
388: PetscErrorCode  PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
389: {
390:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

394:   *link = vmath->link;
395:   return(0);
396: }

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

401:   Input Parameters:
402: + viewer - The Mathematica viewer
403: - type   - The packet type to search for, e.g RETURNPKT

405:   Level: advanced

407: .keywords PetscViewer, Mathematica, packets
408: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
409: @*/
410: PetscErrorCode  PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
411: {
412:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
413:   MLINK                   link   = vmath->link; /* The link to Mathematica */
414:   int                     pkt;                 /* The packet type */

417:   while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
418:   if (!pkt) {
419:     MLClearError(link);
420:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char*) MLErrorMessage(link));
421:   }
422:   return(0);
423: }

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

428:   Input Parameter:
429: . viewer - The Mathematica viewer

431:   Output Parameter:
432: . name   - The name for new objects created in Mathematica

434:   Level: intermediate

436: .keywords PetscViewer, Mathematica, name
437: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
438: @*/
439: PetscErrorCode  PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
440: {
441:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

446:   *name = vmath->objName;
447:   return(0);
448: }

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

453:   Input Parameters:
454: + viewer - The Mathematica viewer
455: - name   - The name for new objects created in Mathematica

457:   Level: intermediate

459: .keywords PetscViewer, Mathematica, name
460: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
461: @*/
462: PetscErrorCode  PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
463: {
464:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

469:   vmath->objName = name;
470:   return(0);
471: }

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

476:   Input Parameter:
477: . viewer - The Mathematica viewer

479:   Level: intermediate

481: .keywords PetscViewer, Mathematica, name
482: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
483: @*/
484: PetscErrorCode  PetscViewerMathematicaClearName(PetscViewer viewer)
485: {
486:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;

490:   vmath->objName = NULL;
491:   return(0);
492: }

494: /*@C
495:   PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica

497:   Input Parameter:
498: . viewer - The Mathematica viewer

500:   Output Parameter:
501: . v      - The vector

503:   Level: intermediate

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


522:   /* Determine the object name */
523:   if (!vmath->objName) name = "vec";
524:   else                 name = (char*) vmath->objName;

526:   link = vmath->link;
527:   VecGetLocalSize(v, &n);
528:   VecGetArray(v, &array);
529:   MLPutFunction(link, "EvaluatePacket", 1);
530:   MLPutSymbol(link, name);
531:   MLEndPacket(link);
532:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
533:   MLGetRealList(link, &mArray, &mSize);
534:   if (n != mSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d",n,mSize);
535:   PetscArraycpy(array, mArray, mSize);
536:   MLDisownRealList(link, mArray, mSize);
537:   VecRestoreArray(v, &array);
538:   return(0);
539: }

541: /*@C
542:   PetscViewerMathematicaPutVector - Send a vector to Mathematica

544:   Input Parameters:
545: + viewer - The Mathematica viewer
546: - v      - The vector

548:   Level: intermediate

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

563:   /* Determine the object name */
564:   if (!vmath->objName) name = "vec";
565:   else                 name = (char*) vmath->objName;

567:   VecGetLocalSize(v, &n);
568:   VecGetArray(v, &array);

570:   /* Send the Vector object */
571:   MLPutFunction(link, "EvaluatePacket", 1);
572:   MLPutFunction(link, "Set", 2);
573:   MLPutSymbol(link, name);
574:   MLPutRealList(link, array, n);
575:   MLEndPacket(link);
576:   /* Skip packets until ReturnPacket */
577:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
578:   /* Skip ReturnPacket */
579:   MLNewPacket(link);

581:   VecRestoreArray(v, &array);
582:   return(0);
583: }

585: PetscErrorCode  PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
586: {
587:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
588:   MLINK                   link   = vmath->link; /* The link to Mathematica */
589:   char                    *name;
590:   PetscErrorCode          ierr;

593:   /* Determine the object name */
594:   if (!vmath->objName) name = "mat";
595:   else                 name = (char*) vmath->objName;

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

613: PetscErrorCode  PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
614: {
615:   PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
616:   MLINK                   link   = vmath->link; /* The link to Mathematica */
617:   const char              *symbol;
618:   char                    *name;
619:   PetscBool               match;
620:   PetscErrorCode          ierr;

623:   /* Determine the object name */
624:   if (!vmath->objName) name = "mat";
625:   else                 name = (char*) vmath->objName;

627:   /* Make sure Mathematica recognizes sparse matrices */
628:   MLPutFunction(link, "EvaluatePacket", 1);
629:   MLPutFunction(link, "Needs", 1);
630:   MLPutString(link, "LinearAlgebra`CSRMatrix`");
631:   MLEndPacket(link);
632:   /* Skip packets until ReturnPacket */
633:   PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
634:   /* Skip ReturnPacket */
635:   MLNewPacket(link);

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

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