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) {
72: PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
73: }
74: return 0;
75: }
77: PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
78: {
79: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
80: #if defined(MATHEMATICA_3_0)
81: int argc = 6;
82: char *argv[6];
83: #else
84: int argc = 5;
85: char *argv[5];
86: #endif
87: char hostname[256];
88: long lerr;
90: /* Link name */
91: argv[0] = "-linkname";
92: if (!vmath->linkname) argv[1] = "math -mathlink";
93: else argv[1] = vmath->linkname;
95: /* Link host */
96: argv[2] = "-linkhost";
97: if (!vmath->linkhost) {
98: PetscGetHostName(hostname, sizeof(hostname));
99: argv[3] = hostname;
100: } else argv[3] = vmath->linkhost;
102: /* Link mode */
103: #if defined(MATHEMATICA_3_0)
104: argv[4] = "-linkmode";
105: switch (vmath->linkmode) {
106: case MATHEMATICA_LINK_CREATE:
107: argv[5] = "Create";
108: break;
109: case MATHEMATICA_LINK_CONNECT:
110: argv[5] = "Connect";
111: break;
112: case MATHEMATICA_LINK_LAUNCH:
113: argv[5] = "Launch";
114: break;
115: }
116: #else
117: switch (vmath->linkmode) {
118: case MATHEMATICA_LINK_CREATE:
119: argv[4] = "-linkcreate";
120: break;
121: case MATHEMATICA_LINK_CONNECT:
122: argv[4] = "-linkconnect";
123: break;
124: case MATHEMATICA_LINK_LAUNCH:
125: argv[4] = "-linklaunch";
126: break;
127: }
128: #endif
129: vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
130: #endif
131: return 0;
132: }
134: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
135: {
136: PetscViewer_Mathematica *vmath;
138: PetscViewerMathematicaInitializePackage();
140: PetscNewLog(v,&vmath);
141: v->data = (void*) vmath;
142: v->ops->destroy = PetscViewerDestroy_Mathematica;
143: v->ops->flush = 0;
144: PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);
146: vmath->linkname = NULL;
147: vmath->linkhost = NULL;
148: vmath->linkmode = MATHEMATICA_LINK_CONNECT;
149: vmath->graphicsType = GRAPHICS_MOTIF;
150: vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
151: vmath->objName = NULL;
153: PetscViewerMathematicaSetFromOptions(v);
154: PetscViewerMathematicaSetupConnection_Private(v);
155: return 0;
156: }
158: static PetscErrorCode PetscViewerMathematicaParseLinkMode(char *modename, LinkMode *mode)
159: {
160: PetscBool isCreate, isConnect, isLaunch;
162: PetscStrcasecmp(modename, "Create", &isCreate);
163: PetscStrcasecmp(modename, "Connect", &isConnect);
164: PetscStrcasecmp(modename, "Launch", &isLaunch);
165: if (isCreate) *mode = MATHEMATICA_LINK_CREATE;
166: else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
167: else if (isLaunch) *mode = MATHEMATICA_LINK_LAUNCH;
168: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
169: return 0;
170: }
172: PetscErrorCode PetscViewerMathematicaSetFromOptions(PetscViewer v)
173: {
174: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
175: char linkname[256];
176: char modename[256];
177: char hostname[256];
178: char type[256];
179: PetscInt numPorts;
180: PetscInt *ports;
181: PetscInt numHosts;
182: int h;
183: char **hosts;
184: PetscMPIInt size, rank;
185: PetscBool opt;
187: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
188: MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);
190: /* Get link name */
191: PetscOptionsGetString("viewer_", "-math_linkname", linkname, sizeof(linkname), &opt);
192: if (opt) {
193: PetscViewerMathematicaSetLinkName(v, linkname);
194: }
195: /* Get link port */
196: numPorts = size;
197: PetscMalloc1(size, &ports);
198: PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
199: if (opt) {
200: if (numPorts > rank) snprintf(linkname, sizeof(linkname), "%6d", ports[rank]);
201: else snprintf(linkname, sizeof(linkname), "%6d", ports[0]);
202: PetscViewerMathematicaSetLinkName(v, linkname);
203: }
204: PetscFree(ports);
205: /* Get link host */
206: numHosts = size;
207: PetscMalloc1(size, &hosts);
208: PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
209: if (opt) {
210: if (numHosts > rank) {
211: PetscStrncpy(hostname, hosts[rank], sizeof(hostname));
212: } else {
213: PetscStrncpy(hostname, hosts[0], sizeof(hostname));
214: }
215: PetscViewerMathematicaSetLinkHost(v, hostname);
216: }
217: for (h = 0; h < numHosts; h++) {
218: PetscFree(hosts[h]);
219: }
220: PetscFree(hosts);
221: /* Get link mode */
222: PetscOptionsGetString("viewer_", "-math_linkmode", modename, sizeof(modename), &opt);
223: if (opt) {
224: LinkMode mode;
226: PetscViewerMathematicaParseLinkMode(modename, &mode);
227: PetscViewerMathematicaSetLinkMode(v, mode);
228: }
229: /* Get graphics type */
230: PetscOptionsGetString("viewer_", "-math_graphics", type, sizeof(type), &opt);
231: if (opt) {
232: PetscBool isMotif, isPS, isPSFile;
234: PetscStrcasecmp(type, "Motif", &isMotif);
235: PetscStrcasecmp(type, "PS", &isPS);
236: PetscStrcasecmp(type, "PSFile", &isPSFile);
237: if (isMotif) vmath->graphicsType = GRAPHICS_MOTIF;
238: else if (isPS) vmath->graphicsType = GRAPHICS_PS_STDOUT;
239: else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
240: }
241: /* Get plot type */
242: PetscOptionsGetString("viewer_", "-math_type", type, sizeof(type), &opt);
243: if (opt) {
244: PetscBool isTri, isVecTri, isVec, isSurface;
246: PetscStrcasecmp(type, "Triangulation", &isTri);
247: PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
248: PetscStrcasecmp(type, "Vector", &isVec);
249: PetscStrcasecmp(type, "Surface", &isSurface);
250: if (isTri) vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
251: else if (isVecTri) vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
252: else if (isVec) vmath->plotType = MATHEMATICA_VECTOR_PLOT;
253: else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
254: }
255: return 0;
256: }
258: PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
259: {
260: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
264: PetscStrallocpy(name, &vmath->linkname);
265: return 0;
266: }
268: PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
269: {
270: char name[16];
272: snprintf(name, 16, "%6d", port);
273: PetscViewerMathematicaSetLinkName(v, name);
274: return 0;
275: }
277: PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
278: {
279: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
283: PetscStrallocpy(host, &vmath->linkhost);
284: return 0;
285: }
287: PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
288: {
289: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
291: vmath->linkmode = mode;
292: return 0;
293: }
295: /*----------------------------------------- Public Functions --------------------------------------------------------*/
296: /*@C
297: PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
299: Collective
301: Input Parameters:
302: + comm - The MPI communicator
303: . port - [optional] The port to connect on, or PETSC_DECIDE
304: . machine - [optional] The machine to run Mathematica on, or NULL
305: - mode - [optional] The connection mode, or NULL
307: Output Parameter:
308: . viewer - The Mathematica viewer
310: Level: intermediate
312: Notes:
313: Most users should employ the following commands to access the
314: Mathematica viewers
315: $
316: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
317: $ MatView(Mat matrix, PetscViewer viewer)
318: $
319: $ or
320: $
321: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
322: $ VecView(Vec vector, PetscViewer viewer)
324: Options Database Keys:
325: + -viewer_math_linkhost <machine> - The host machine for the kernel
326: . -viewer_math_linkname <name> - The full link name for the connection
327: . -viewer_math_linkport <port> - The port for the connection
328: . -viewer_math_mode <mode> - The mode, e.g. Launch, Connect
329: . -viewer_math_type <type> - The plot type, e.g. Triangulation, Vector
330: - -viewer_math_graphics <output> - The output type, e.g. Motif, PS, PSFile
332: .seealso: MatView(), VecView()
333: @*/
334: PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
335: {
336: PetscViewerCreate(comm, v);
337: #if 0
338: LinkMode linkmode;
339: PetscViewerMathematicaSetLinkPort(*v, port);
340: PetscViewerMathematicaSetLinkHost(*v, machine);
341: PetscViewerMathematicaParseLinkMode(mode, &linkmode);
342: PetscViewerMathematicaSetLinkMode(*v, linkmode);
343: #endif
344: PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
345: return 0;
346: }
348: /*@C
349: PetscViewerMathematicaGetLink - Returns the link to Mathematica
351: Input Parameters:
352: + viewer - The Mathematica viewer
353: - link - The link to Mathematica
355: Level: intermediate
357: .keywords PetscViewer, Mathematica, link
358: .seealso PetscViewerMathematicaOpen()
359: @*/
360: PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
361: {
362: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
365: *link = vmath->link;
366: return 0;
367: }
369: /*@C
370: PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
372: Input Parameters:
373: + viewer - The Mathematica viewer
374: - type - The packet type to search for, e.g RETURNPKT
376: Level: advanced
378: .keywords PetscViewer, Mathematica, packets
379: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
380: @*/
381: PetscErrorCode PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
382: {
383: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
384: MLINK link = vmath->link; /* The link to Mathematica */
385: int pkt; /* The packet type */
387: while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
388: if (!pkt) {
389: MLClearError(link);
390: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char*) MLErrorMessage(link));
391: }
392: return 0;
393: }
395: /*@C
396: PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica
398: Input Parameter:
399: . viewer - The Mathematica viewer
401: Output Parameter:
402: . name - The name for new objects created in Mathematica
404: Level: intermediate
406: .keywords PetscViewer, Mathematica, name
407: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
408: @*/
409: PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
410: {
411: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
415: *name = vmath->objName;
416: return 0;
417: }
419: /*@C
420: PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica
422: Input Parameters:
423: + viewer - The Mathematica viewer
424: - name - The name for new objects created in Mathematica
426: Level: intermediate
428: .keywords PetscViewer, Mathematica, name
429: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
430: @*/
431: PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
432: {
433: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
437: vmath->objName = name;
438: return 0;
439: }
441: /*@C
442: PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
444: Input Parameter:
445: . viewer - The Mathematica viewer
447: Level: intermediate
449: .keywords PetscViewer, Mathematica, name
450: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
451: @*/
452: PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
453: {
454: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
457: vmath->objName = NULL;
458: return 0;
459: }
461: /*@C
462: PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica
464: Input Parameter:
465: . viewer - The Mathematica viewer
467: Output Parameter:
468: . v - The vector
470: Level: intermediate
472: .keywords PetscViewer, Mathematica, vector
473: .seealso VecView(), PetscViewerMathematicaPutVector()
474: @*/
475: PetscErrorCode PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
476: {
477: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
478: MLINK link; /* The link to Mathematica */
479: char *name;
480: PetscScalar *mArray,*array;
481: long mSize;
482: int n;
487: /* Determine the object name */
488: if (!vmath->objName) name = "vec";
489: else name = (char*) vmath->objName;
491: link = vmath->link;
492: VecGetLocalSize(v, &n);
493: VecGetArray(v, &array);
494: MLPutFunction(link, "EvaluatePacket", 1);
495: MLPutSymbol(link, name);
496: MLEndPacket(link);
497: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
498: MLGetRealList(link, &mArray, &mSize);
500: PetscArraycpy(array, mArray, mSize);
501: MLDisownRealList(link, mArray, mSize);
502: VecRestoreArray(v, &array);
503: return 0;
504: }
506: /*@C
507: PetscViewerMathematicaPutVector - Send a vector to Mathematica
509: Input Parameters:
510: + viewer - The Mathematica viewer
511: - v - The vector
513: Level: intermediate
515: .keywords PetscViewer, Mathematica, vector
516: .seealso VecView(), PetscViewerMathematicaGetVector()
517: @*/
518: PetscErrorCode PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
519: {
520: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
521: MLINK link = vmath->link; /* The link to Mathematica */
522: char *name;
523: PetscScalar *array;
524: int n;
526: /* Determine the object name */
527: if (!vmath->objName) name = "vec";
528: else name = (char*) vmath->objName;
530: VecGetLocalSize(v, &n);
531: VecGetArray(v, &array);
533: /* Send the Vector object */
534: MLPutFunction(link, "EvaluatePacket", 1);
535: MLPutFunction(link, "Set", 2);
536: MLPutSymbol(link, name);
537: MLPutRealList(link, array, n);
538: MLEndPacket(link);
539: /* Skip packets until ReturnPacket */
540: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
541: /* Skip ReturnPacket */
542: MLNewPacket(link);
544: VecRestoreArray(v, &array);
545: return 0;
546: }
548: PetscErrorCode PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
549: {
550: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
551: MLINK link = vmath->link; /* The link to Mathematica */
552: char *name;
554: /* Determine the object name */
555: if (!vmath->objName) name = "mat";
556: else name = (char*) vmath->objName;
558: /* Send the dense matrix object */
559: MLPutFunction(link, "EvaluatePacket", 1);
560: MLPutFunction(link, "Set", 2);
561: MLPutSymbol(link, name);
562: MLPutFunction(link, "Transpose", 1);
563: MLPutFunction(link, "Partition", 2);
564: MLPutRealList(link, a, m*n);
565: MLPutInteger(link, m);
566: MLEndPacket(link);
567: /* Skip packets until ReturnPacket */
568: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
569: /* Skip ReturnPacket */
570: MLNewPacket(link);
571: return 0;
572: }
574: PetscErrorCode PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
575: {
576: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
577: MLINK link = vmath->link; /* The link to Mathematica */
578: const char *symbol;
579: char *name;
580: PetscBool match;
582: /* Determine the object name */
583: if (!vmath->objName) name = "mat";
584: else name = (char*) vmath->objName;
586: /* Make sure Mathematica recognizes sparse matrices */
587: MLPutFunction(link, "EvaluatePacket", 1);
588: MLPutFunction(link, "Needs", 1);
589: MLPutString(link, "LinearAlgebra`CSRMatrix`");
590: MLEndPacket(link);
591: /* Skip packets until ReturnPacket */
592: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
593: /* Skip ReturnPacket */
594: MLNewPacket(link);
596: /* Send the CSRMatrix object */
597: MLPutFunction(link, "EvaluatePacket", 1);
598: MLPutFunction(link, "Set", 2);
599: MLPutSymbol(link, name);
600: MLPutFunction(link, "CSRMatrix", 5);
601: MLPutInteger(link, m);
602: MLPutInteger(link, n);
603: MLPutFunction(link, "Plus", 2);
604: MLPutIntegerList(link, i, m+1);
605: MLPutInteger(link, 1);
606: MLPutFunction(link, "Plus", 2);
607: MLPutIntegerList(link, j, i[m]);
608: MLPutInteger(link, 1);
609: MLPutRealList(link, a, i[m]);
610: MLEndPacket(link);
611: /* Skip packets until ReturnPacket */
612: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
613: /* Skip ReturnPacket */
614: MLNewPacket(link);
616: /* Check that matrix is valid */
617: MLPutFunction(link, "EvaluatePacket", 1);
618: MLPutFunction(link, "ValidQ", 1);
619: MLPutSymbol(link, name);
620: MLEndPacket(link);
621: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
622: MLGetSymbol(link, &symbol);
623: PetscStrcmp("True", (char*) symbol, &match);
624: if (!match) {
625: MLDisownSymbol(link, symbol);
626: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
627: }
628: MLDisownSymbol(link, symbol);
629: /* Skip ReturnPacket */
630: MLNewPacket(link);
631: return 0;
632: }