Actual source code: mathematica.c

petsc-master 2019-10-23
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: .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: }


 54: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
 55: {

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

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

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

 77: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
 78: {

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

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

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

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

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

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

153:   PetscViewerMathematicaInitializePackage();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

325:   Collective

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

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

336:   Level: intermediate

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

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

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

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

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

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

384:   Level: intermediate

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

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

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

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

406:   Level: advanced

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

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

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

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

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

435:   Level: intermediate

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

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

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

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

458:   Level: intermediate

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

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

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

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

480:   Level: intermediate

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

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

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

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

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

504:   Level: intermediate

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


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

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

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

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

549:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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