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: {
25: if (mathematicaEnv) MLDeinitialize((MLEnvironment)mathematicaEnv);
26: PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
27: return 0;
28: }
30: /*@C
31: PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
32: called from `PetscViewerInitializePackage()`.
34: Level: developer
36: .seealso: `PetscSysInitializePackage()`, `PetscInitialize()`
37: @*/
38: PetscErrorCode PetscViewerMathematicaInitializePackage(void)
39: {
40: PetscError ierr;
42: if (PetscViewerMathematicaPackageInitialized) return 0;
43: PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
45: mathematicaEnv = (void *)MLInitialize(0);
47: PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
48: return 0;
49: }
51: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
52: {
53: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return 0;
54: PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
55: return 0;
56: }
58: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
59: {
60: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
62: MLClose(vmath->link);
63: PetscFree(vmath->linkname);
64: PetscFree(vmath->linkhost);
65: PetscFree(vmath);
66: return 0;
67: }
69: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
70: {
71: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
72: return 0;
73: }
75: PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
76: {
77: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
78: #if defined(MATHEMATICA_3_0)
79: int argc = 6;
80: char *argv[6];
81: #else
82: int argc = 5;
83: char *argv[5];
84: #endif
85: char hostname[256];
86: long lerr;
88: /* Link name */
89: argv[0] = "-linkname";
90: if (!vmath->linkname) argv[1] = "math -mathlink";
91: else argv[1] = vmath->linkname;
93: /* Link host */
94: argv[2] = "-linkhost";
95: if (!vmath->linkhost) {
96: PetscGetHostName(hostname, sizeof(hostname));
97: argv[3] = hostname;
98: } else argv[3] = vmath->linkhost;
100: /* Link mode */
101: #if defined(MATHEMATICA_3_0)
102: argv[4] = "-linkmode";
103: switch (vmath->linkmode) {
104: case MATHEMATICA_LINK_CREATE:
105: argv[5] = "Create";
106: break;
107: case MATHEMATICA_LINK_CONNECT:
108: argv[5] = "Connect";
109: break;
110: case MATHEMATICA_LINK_LAUNCH:
111: argv[5] = "Launch";
112: break;
113: }
114: #else
115: switch (vmath->linkmode) {
116: case MATHEMATICA_LINK_CREATE:
117: argv[4] = "-linkcreate";
118: break;
119: case MATHEMATICA_LINK_CONNECT:
120: argv[4] = "-linkconnect";
121: break;
122: case MATHEMATICA_LINK_LAUNCH:
123: argv[4] = "-linklaunch";
124: break;
125: }
126: #endif
127: vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
128: #endif
129: return 0;
130: }
132: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
133: {
134: PetscViewer_Mathematica *vmath;
136: PetscViewerMathematicaInitializePackage();
138: PetscNew(&vmath);
139: v->data = (void *)vmath;
140: v->ops->destroy = PetscViewerDestroy_Mathematica;
141: v->ops->flush = 0;
142: PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);
144: vmath->linkname = NULL;
145: vmath->linkhost = NULL;
146: vmath->linkmode = MATHEMATICA_LINK_CONNECT;
147: vmath->graphicsType = GRAPHICS_MOTIF;
148: vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
149: vmath->objName = NULL;
151: PetscViewerMathematicaSetFromOptions(v);
152: PetscViewerMathematicaSetupConnection_Private(v);
153: return 0;
154: }
156: static PetscErrorCode PetscViewerMathematicaParseLinkMode(char *modename, LinkMode *mode)
157: {
158: PetscBool isCreate, isConnect, isLaunch;
160: PetscStrcasecmp(modename, "Create", &isCreate);
161: PetscStrcasecmp(modename, "Connect", &isConnect);
162: PetscStrcasecmp(modename, "Launch", &isLaunch);
163: if (isCreate) *mode = MATHEMATICA_LINK_CREATE;
164: else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
165: else if (isLaunch) *mode = MATHEMATICA_LINK_LAUNCH;
166: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
167: return 0;
168: }
170: PetscErrorCode PetscViewerMathematicaSetFromOptions(PetscViewer v)
171: {
172: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
173: char linkname[256];
174: char modename[256];
175: char hostname[256];
176: char type[256];
177: PetscInt numPorts;
178: PetscInt *ports;
179: PetscInt numHosts;
180: int h;
181: char **hosts;
182: PetscMPIInt size, rank;
183: PetscBool opt;
185: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
186: MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);
188: /* Get link name */
189: PetscOptionsGetString("viewer_", "-math_linkname", linkname, sizeof(linkname), &opt);
190: if (opt) PetscViewerMathematicaSetLinkName(v, linkname);
191: /* Get link port */
192: numPorts = size;
193: PetscMalloc1(size, &ports);
194: PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
195: if (opt) {
196: if (numPorts > rank) snprintf(linkname, sizeof(linkname), "%6d", ports[rank]);
197: else snprintf(linkname, sizeof(linkname), "%6d", ports[0]);
198: PetscViewerMathematicaSetLinkName(v, linkname);
199: }
200: PetscFree(ports);
201: /* Get link host */
202: numHosts = size;
203: PetscMalloc1(size, &hosts);
204: PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
205: if (opt) {
206: if (numHosts > rank) {
207: PetscStrncpy(hostname, hosts[rank], sizeof(hostname));
208: } else {
209: PetscStrncpy(hostname, hosts[0], sizeof(hostname));
210: }
211: PetscViewerMathematicaSetLinkHost(v, hostname);
212: }
213: for (h = 0; h < numHosts; h++) PetscFree(hosts[h]);
214: PetscFree(hosts);
215: /* Get link mode */
216: PetscOptionsGetString("viewer_", "-math_linkmode", modename, sizeof(modename), &opt);
217: if (opt) {
218: LinkMode mode;
220: PetscViewerMathematicaParseLinkMode(modename, &mode);
221: PetscViewerMathematicaSetLinkMode(v, mode);
222: }
223: /* Get graphics type */
224: PetscOptionsGetString("viewer_", "-math_graphics", type, sizeof(type), &opt);
225: if (opt) {
226: PetscBool isMotif, isPS, isPSFile;
228: PetscStrcasecmp(type, "Motif", &isMotif);
229: PetscStrcasecmp(type, "PS", &isPS);
230: PetscStrcasecmp(type, "PSFile", &isPSFile);
231: if (isMotif) vmath->graphicsType = GRAPHICS_MOTIF;
232: else if (isPS) vmath->graphicsType = GRAPHICS_PS_STDOUT;
233: else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
234: }
235: /* Get plot type */
236: PetscOptionsGetString("viewer_", "-math_type", type, sizeof(type), &opt);
237: if (opt) {
238: PetscBool isTri, isVecTri, isVec, isSurface;
240: PetscStrcasecmp(type, "Triangulation", &isTri);
241: PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
242: PetscStrcasecmp(type, "Vector", &isVec);
243: PetscStrcasecmp(type, "Surface", &isSurface);
244: if (isTri) vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
245: else if (isVecTri) vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
246: else if (isVec) vmath->plotType = MATHEMATICA_VECTOR_PLOT;
247: else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
248: }
249: return 0;
250: }
252: PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
253: {
254: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
258: PetscStrallocpy(name, &vmath->linkname);
259: return 0;
260: }
262: PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
263: {
264: char name[16];
266: snprintf(name, 16, "%6d", port);
267: PetscViewerMathematicaSetLinkName(v, name);
268: return 0;
269: }
271: PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
272: {
273: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
277: PetscStrallocpy(host, &vmath->linkhost);
278: return 0;
279: }
281: PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
282: {
283: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)v->data;
285: vmath->linkmode = mode;
286: return 0;
287: }
289: /*----------------------------------------- Public Functions --------------------------------------------------------*/
290: /*@C
291: PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
293: Collective
295: Input Parameters:
296: + comm - The MPI communicator
297: . port - [optional] The port to connect on, or PETSC_DECIDE
298: . machine - [optional] The machine to run Mathematica on, or NULL
299: - mode - [optional] The connection mode, or NULL
301: Output Parameter:
302: . viewer - The Mathematica viewer
304: Options Database Keys:
305: + -viewer_math_linkhost <machine> - The host machine for the kernel
306: . -viewer_math_linkname <name> - The full link name for the connection
307: . -viewer_math_linkport <port> - The port for the connection
308: . -viewer_math_mode <mode> - The mode, e.g. Launch, Connect
309: . -viewer_math_type <type> - The plot type, e.g. Triangulation, Vector
310: - -viewer_math_graphics <output> - The output type, e.g. Motif, PS, PSFile
312: Level: intermediate
314: Note:
315: Most users should employ the following commands to access the
316: Mathematica viewers
317: .vb
318: PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
319: MatView(Mat matrix, PetscViewer viewer)
321: or
323: PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
324: VecView(Vec vector, PetscViewer viewer)
325: .ve
327: .seealso: `PETSCVIEWERMATHEMATICA`, `MatView()`, `VecView()`
328: @*/
329: PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
330: {
331: PetscViewerCreate(comm, v);
332: #if 0
333: LinkMode linkmode;
334: PetscViewerMathematicaSetLinkPort(*v, port);
335: PetscViewerMathematicaSetLinkHost(*v, machine);
336: PetscViewerMathematicaParseLinkMode(mode, &linkmode);
337: PetscViewerMathematicaSetLinkMode(*v, linkmode);
338: #endif
339: PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
340: return 0;
341: }
343: /*@C
344: PetscViewerMathematicaGetLink - Returns the link to Mathematica from a `PETSCVIEWERMATHEMATICA`
346: Input Parameters:
347: + viewer - The Mathematica viewer
348: - link - The link to Mathematica
350: Level: intermediate
352: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaOpen()`
353: @*/
354: PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
355: {
356: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
359: *link = vmath->link;
360: return 0;
361: }
363: /*@C
364: PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
366: Input Parameters:
367: + viewer - The Mathematica viewer
368: - type - The packet type to search for, e.g RETURNPKT
370: Level: advanced
372: .seealso: `PetscViewerMathematicaSetName()`, `PetscViewerMathematicaGetVector()`
373: @*/
374: PetscErrorCode PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
375: {
376: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
377: MLINK link = vmath->link; /* The link to Mathematica */
378: int pkt; /* The packet type */
380: while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
381: if (!pkt) {
382: MLClearError(link);
383: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, (char *)MLErrorMessage(link));
384: }
385: return 0;
386: }
388: /*@C
389: PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica via `PETSCVIEWERMATHEMATICA`
391: Input Parameter:
392: . viewer - The Mathematica viewer
394: Output Parameter:
395: . name - The name for new objects created in Mathematica
397: Level: intermediate
399: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaSetName()`, `PetscViewerMathematicaClearName()`
400: @*/
401: PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
402: {
403: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
407: *name = vmath->objName;
408: return 0;
409: }
411: /*@C
412: PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica via `PETSCVIEWERMATHEMATICA`
414: Input Parameters:
415: + viewer - The Mathematica viewer
416: - name - The name for new objects created in Mathematica
418: Level: intermediate
420: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaSetName()`, `PetscViewerMathematicaClearName()`
421: @*/
422: PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
423: {
424: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
428: vmath->objName = name;
429: return 0;
430: }
432: /*@
433: PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
435: Input Parameter:
436: . viewer - The Mathematica viewer
438: Level: intermediate
440: .seealso: `PETSCVIEWERMATHEMATICA`, `PetscViewerMathematicaGetName()`, `PetscViewerMathematicaSetName()`
441: @*/
442: PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
443: {
444: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
447: vmath->objName = NULL;
448: return 0;
449: }
451: /*@
452: PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica via a `PETSCVIEWERMATHEMATICA`
454: Input Parameter:
455: . viewer - The Mathematica viewer
457: Output Parameter:
458: . v - The vector
460: Level: intermediate
462: .seealso: `PETSCVIEWERMATHEMATICA`, `VecView()`, `PetscViewerMathematicaPutVector()`
463: @*/
464: PetscErrorCode PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
465: {
466: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
467: MLINK link; /* The link to Mathematica */
468: char *name;
469: PetscScalar *mArray, *array;
470: long mSize;
471: int n;
476: /* Determine the object name */
477: if (!vmath->objName) name = "vec";
478: else name = (char *)vmath->objName;
480: link = vmath->link;
481: VecGetLocalSize(v, &n);
482: VecGetArray(v, &array);
483: MLPutFunction(link, "EvaluatePacket", 1);
484: MLPutSymbol(link, name);
485: MLEndPacket(link);
486: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
487: MLGetRealList(link, &mArray, &mSize);
489: PetscArraycpy(array, mArray, mSize);
490: MLDisownRealList(link, mArray, mSize);
491: VecRestoreArray(v, &array);
492: return 0;
493: }
495: /*@
496: PetscViewerMathematicaPutVector - Send a vector to Mathematica via a `PETSCVIEWERMATHEMATICA` `PetscViewer`
498: Input Parameters:
499: + viewer - The Mathematica viewer
500: - v - The vector
502: Level: intermediate
504: .seealso: `PETSCVIEWERMATHEMATICA`, `VecView()`, `PetscViewerMathematicaGetVector()`
505: @*/
506: PetscErrorCode PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
507: {
508: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
509: MLINK link = vmath->link; /* The link to Mathematica */
510: char *name;
511: PetscScalar *array;
512: int n;
514: /* Determine the object name */
515: if (!vmath->objName) name = "vec";
516: else name = (char *)vmath->objName;
518: VecGetLocalSize(v, &n);
519: VecGetArray(v, &array);
521: /* Send the Vector object */
522: MLPutFunction(link, "EvaluatePacket", 1);
523: MLPutFunction(link, "Set", 2);
524: MLPutSymbol(link, name);
525: MLPutRealList(link, array, n);
526: MLEndPacket(link);
527: /* Skip packets until ReturnPacket */
528: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
529: /* Skip ReturnPacket */
530: MLNewPacket(link);
532: VecRestoreArray(v, &array);
533: return 0;
534: }
536: PetscErrorCode PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
537: {
538: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
539: MLINK link = vmath->link; /* The link to Mathematica */
540: char *name;
542: /* Determine the object name */
543: if (!vmath->objName) name = "mat";
544: else name = (char *)vmath->objName;
546: /* Send the dense matrix object */
547: MLPutFunction(link, "EvaluatePacket", 1);
548: MLPutFunction(link, "Set", 2);
549: MLPutSymbol(link, name);
550: MLPutFunction(link, "Transpose", 1);
551: MLPutFunction(link, "Partition", 2);
552: MLPutRealList(link, a, m * n);
553: MLPutInteger(link, m);
554: MLEndPacket(link);
555: /* Skip packets until ReturnPacket */
556: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
557: /* Skip ReturnPacket */
558: MLNewPacket(link);
559: return 0;
560: }
562: PetscErrorCode PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
563: {
564: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *)viewer->data;
565: MLINK link = vmath->link; /* The link to Mathematica */
566: const char *symbol;
567: char *name;
568: PetscBool match;
570: /* Determine the object name */
571: if (!vmath->objName) name = "mat";
572: else name = (char *)vmath->objName;
574: /* Make sure Mathematica recognizes sparse matrices */
575: MLPutFunction(link, "EvaluatePacket", 1);
576: MLPutFunction(link, "Needs", 1);
577: MLPutString(link, "LinearAlgebra`CSRMatrix`");
578: MLEndPacket(link);
579: /* Skip packets until ReturnPacket */
580: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
581: /* Skip ReturnPacket */
582: MLNewPacket(link);
584: /* Send the CSRMatrix object */
585: MLPutFunction(link, "EvaluatePacket", 1);
586: MLPutFunction(link, "Set", 2);
587: MLPutSymbol(link, name);
588: MLPutFunction(link, "CSRMatrix", 5);
589: MLPutInteger(link, m);
590: MLPutInteger(link, n);
591: MLPutFunction(link, "Plus", 2);
592: MLPutIntegerList(link, i, m + 1);
593: MLPutInteger(link, 1);
594: MLPutFunction(link, "Plus", 2);
595: MLPutIntegerList(link, j, i[m]);
596: MLPutInteger(link, 1);
597: MLPutRealList(link, a, i[m]);
598: MLEndPacket(link);
599: /* Skip packets until ReturnPacket */
600: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
601: /* Skip ReturnPacket */
602: MLNewPacket(link);
604: /* Check that matrix is valid */
605: MLPutFunction(link, "EvaluatePacket", 1);
606: MLPutFunction(link, "ValidQ", 1);
607: MLPutSymbol(link, name);
608: MLEndPacket(link);
609: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
610: MLGetSymbol(link, &symbol);
611: PetscStrcmp("True", (char *)symbol, &match);
612: if (!match) {
613: MLDisownSymbol(link, symbol);
614: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
615: }
616: MLDisownSymbol(link, symbol);
617: /* Skip ReturnPacket */
618: MLNewPacket(link);
619: return 0;
620: }