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: }