Actual source code: mathematica.c
petsc-3.11.4 2019-09-28
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: .keywords: Petsc, destroy, package, mathematica
22: .seealso: PetscFinalize()
23: @*/
24: PetscErrorCode PetscViewerMathematicaFinalizePackage(void)
25: {
27: if (mathematicaEnv) MLDeinitialize((MLEnvironment) mathematicaEnv);
28: PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
29: return(0);
30: }
32: /*@C
33: PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
34: called from PetscViewerInitializePackage().
36: Level: developer
38: .keywords: Petsc, initialize, package
39: .seealso: PetscSysInitializePackage(), PetscInitialize()
40: @*/
41: PetscErrorCode PetscViewerMathematicaInitializePackage(void)
42: {
43: PetscError ierr;
46: if (PetscViewerMathematicaPackageInitialized) return(0);
47: PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
49: mathematicaEnv = (void*) MLInitialize(0);
51: PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
52: return(0);
53: }
56: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
57: {
61: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return(0);
62: PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
63: return(0);
64: }
66: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
67: {
68: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
69: PetscErrorCode ierr;
72: MLClose(vmath->link);
73: PetscFree(vmath->linkname);
74: PetscFree(vmath->linkhost);
75: PetscFree(vmath);
76: return(0);
77: }
79: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
80: {
84: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
85: PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
86: }
87: return(0);
88: }
90: PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
91: {
92: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
93: #if defined(MATHEMATICA_3_0)
94: int argc = 6;
95: char *argv[6];
96: #else
97: int argc = 5;
98: char *argv[5];
99: #endif
100: char hostname[256];
101: long lerr;
102: PetscErrorCode ierr;
105: /* Link name */
106: argv[0] = "-linkname";
107: if (!vmath->linkname) argv[1] = "math -mathlink";
108: else argv[1] = vmath->linkname;
110: /* Link host */
111: argv[2] = "-linkhost";
112: if (!vmath->linkhost) {
113: PetscGetHostName(hostname, 255);
114: argv[3] = hostname;
115: } else argv[3] = vmath->linkhost;
117: /* Link mode */
118: #if defined(MATHEMATICA_3_0)
119: argv[4] = "-linkmode";
120: switch (vmath->linkmode) {
121: case MATHEMATICA_LINK_CREATE:
122: argv[5] = "Create";
123: break;
124: case MATHEMATICA_LINK_CONNECT:
125: argv[5] = "Connect";
126: break;
127: case MATHEMATICA_LINK_LAUNCH:
128: argv[5] = "Launch";
129: break;
130: }
131: #else
132: switch (vmath->linkmode) {
133: case MATHEMATICA_LINK_CREATE:
134: argv[4] = "-linkcreate";
135: break;
136: case MATHEMATICA_LINK_CONNECT:
137: argv[4] = "-linkconnect";
138: break;
139: case MATHEMATICA_LINK_LAUNCH:
140: argv[4] = "-linklaunch";
141: break;
142: }
143: #endif
144: vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
145: #endif
146: return(0);
147: }
149: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
150: {
151: PetscViewer_Mathematica *vmath;
152: PetscErrorCode ierr;
155: PetscViewerMathematicaInitializePackage();
157: PetscNewLog(v,&vmath);
158: v->data = (void*) vmath;
159: v->ops->destroy = PetscViewerDestroy_Mathematica;
160: v->ops->flush = 0;
161: PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);
163: vmath->linkname = NULL;
164: vmath->linkhost = NULL;
165: vmath->linkmode = MATHEMATICA_LINK_CONNECT;
166: vmath->graphicsType = GRAPHICS_MOTIF;
167: vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
168: vmath->objName = NULL;
170: PetscViewerMathematicaSetFromOptions(v);
171: PetscViewerMathematicaSetupConnection_Private(v);
172: return(0);
173: }
175: static PetscErrorCode PetscViewerMathematicaParseLinkMode(char *modename, LinkMode *mode)
176: {
177: PetscBool isCreate, isConnect, isLaunch;
181: PetscStrcasecmp(modename, "Create", &isCreate);
182: PetscStrcasecmp(modename, "Connect", &isConnect);
183: PetscStrcasecmp(modename, "Launch", &isLaunch);
184: if (isCreate) *mode = MATHEMATICA_LINK_CREATE;
185: else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
186: else if (isLaunch) *mode = MATHEMATICA_LINK_LAUNCH;
187: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
188: return(0);
189: }
191: PetscErrorCode PetscViewerMathematicaSetFromOptions(PetscViewer v)
192: {
193: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
194: char linkname[256];
195: char modename[256];
196: char hostname[256];
197: char type[256];
198: PetscInt numPorts;
199: PetscInt *ports;
200: PetscInt numHosts;
201: int h;
202: char **hosts;
203: PetscMPIInt size, rank;
204: PetscBool opt;
205: PetscErrorCode ierr;
208: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
209: MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);
211: /* Get link name */
212: PetscOptionsGetString("viewer_", "-math_linkname", linkname, 256, &opt);
213: if (opt) {
214: PetscViewerMathematicaSetLinkName(v, linkname);
215: }
216: /* Get link port */
217: numPorts = size;
218: PetscMalloc1(size, &ports);
219: PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
220: if (opt) {
221: if (numPorts > rank) snprintf(linkname, 255, "%6d", ports[rank]);
222: else snprintf(linkname, 255, "%6d", ports[0]);
223: PetscViewerMathematicaSetLinkName(v, linkname);
224: }
225: PetscFree(ports);
226: /* Get link host */
227: numHosts = size;
228: PetscMalloc1(size, &hosts);
229: PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
230: if (opt) {
231: if (numHosts > rank) {
232: PetscStrncpy(hostname, hosts[rank], 255);
233: } else {
234: PetscStrncpy(hostname, hosts[0], 255);
235: }
236: PetscViewerMathematicaSetLinkHost(v, hostname);
237: }
238: for (h = 0; h < numHosts; h++) {
239: PetscFree(hosts[h]);
240: }
241: PetscFree(hosts);
242: /* Get link mode */
243: PetscOptionsGetString("viewer_", "-math_linkmode", modename, 256, &opt);
244: if (opt) {
245: LinkMode mode;
247: PetscViewerMathematicaParseLinkMode(modename, &mode);
248: PetscViewerMathematicaSetLinkMode(v, mode);
249: }
250: /* Get graphics type */
251: PetscOptionsGetString("viewer_", "-math_graphics", type, 256, &opt);
252: if (opt) {
253: PetscBool isMotif, isPS, isPSFile;
255: PetscStrcasecmp(type, "Motif", &isMotif);
256: PetscStrcasecmp(type, "PS", &isPS);
257: PetscStrcasecmp(type, "PSFile", &isPSFile);
258: if (isMotif) vmath->graphicsType = GRAPHICS_MOTIF;
259: else if (isPS) vmath->graphicsType = GRAPHICS_PS_STDOUT;
260: else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
261: }
262: /* Get plot type */
263: PetscOptionsGetString("viewer_", "-math_type", type, 256, &opt);
264: if (opt) {
265: PetscBool isTri, isVecTri, isVec, isSurface;
267: PetscStrcasecmp(type, "Triangulation", &isTri);
268: PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
269: PetscStrcasecmp(type, "Vector", &isVec);
270: PetscStrcasecmp(type, "Surface", &isSurface);
271: if (isTri) vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
272: else if (isVecTri) vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
273: else if (isVec) vmath->plotType = MATHEMATICA_VECTOR_PLOT;
274: else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
275: }
276: return(0);
277: }
279: PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
280: {
281: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
282: PetscErrorCode ierr;
287: PetscStrallocpy(name, &vmath->linkname);
288: return(0);
289: }
291: PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
292: {
293: char name[16];
297: snprintf(name, 16, "%6d", port);
298: PetscViewerMathematicaSetLinkName(v, name);
299: return(0);
300: }
302: PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
303: {
304: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
305: PetscErrorCode ierr;
310: PetscStrallocpy(host, &vmath->linkhost);
311: return(0);
312: }
314: PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
315: {
316: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
319: vmath->linkmode = mode;
320: return(0);
321: }
323: /*----------------------------------------- Public Functions --------------------------------------------------------*/
324: /*@C
325: PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
327: Collective on comm
329: Input Parameters:
330: + comm - The MPI communicator
331: . port - [optional] The port to connect on, or PETSC_DECIDE
332: . machine - [optional] The machine to run Mathematica on, or NULL
333: - mode - [optional] The connection mode, or NULL
335: Output Parameter:
336: . viewer - The Mathematica viewer
338: Level: intermediate
340: Notes:
341: Most users should employ the following commands to access the
342: Mathematica viewers
343: $
344: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
345: $ MatView(Mat matrix, PetscViewer viewer)
346: $
347: $ or
348: $
349: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
350: $ VecView(Vec vector, PetscViewer viewer)
352: Options Database Keys:
353: + -viewer_math_linkhost <machine> - The host machine for the kernel
354: . -viewer_math_linkname <name> - The full link name for the connection
355: . -viewer_math_linkport <port> - The port for the connection
356: . -viewer_math_mode <mode> - The mode, e.g. Launch, Connect
357: . -viewer_math_type <type> - The plot type, e.g. Triangulation, Vector
358: - -viewer_math_graphics <output> - The output type, e.g. Motif, PS, PSFile
360: .keywords: PetscViewer, Mathematica, open
362: .seealso: MatView(), VecView()
363: @*/
364: PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
365: {
369: PetscViewerCreate(comm, v);
370: #if 0
371: LinkMode linkmode;
372: PetscViewerMathematicaSetLinkPort(*v, port);
373: PetscViewerMathematicaSetLinkHost(*v, machine);
374: PetscViewerMathematicaParseLinkMode(mode, &linkmode);
375: PetscViewerMathematicaSetLinkMode(*v, linkmode);
376: #endif
377: PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
378: return(0);
379: }
381: /*@C
382: PetscViewerMathematicaGetLink - Returns the link to Mathematica
384: Input Parameters:
385: . viewer - The Mathematica viewer
386: . link - The link to Mathematica
388: Level: intermediate
390: .keywords PetscViewer, Mathematica, link
391: .seealso PetscViewerMathematicaOpen()
392: @*/
393: PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
394: {
395: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
399: *link = vmath->link;
400: return(0);
401: }
403: /*@C
404: PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
406: Input Parameters:
407: . viewer - The Mathematica viewer
408: . type - The packet type to search for, e.g RETURNPKT
410: Level: advanced
412: .keywords PetscViewer, Mathematica, packets
413: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
414: @*/
415: PetscErrorCode PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
416: {
417: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
418: MLINK link = vmath->link; /* The link to Mathematica */
419: int pkt; /* The packet type */
422: while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
423: if (!pkt) {
424: MLClearError(link);
425: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char*) MLErrorMessage(link));
426: }
427: return(0);
428: }
430: /*@C
431: PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica
433: Input Parameter:
434: . viewer - The Mathematica viewer
436: Output Parameter:
437: . name - The name for new objects created in Mathematica
439: Level: intermediate
441: .keywords PetscViewer, Mathematica, name
442: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
443: @*/
444: PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
445: {
446: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
451: *name = vmath->objName;
452: return(0);
453: }
455: /*@C
456: PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica
458: Input Parameters:
459: . viewer - The Mathematica viewer
460: . name - The name for new objects created in Mathematica
462: Level: intermediate
464: .keywords PetscViewer, Mathematica, name
465: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
466: @*/
467: PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
468: {
469: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
474: vmath->objName = name;
475: return(0);
476: }
478: /*@C
479: PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
481: Input Parameter:
482: . viewer - The Mathematica viewer
484: Level: intermediate
486: .keywords PetscViewer, Mathematica, name
487: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
488: @*/
489: PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
490: {
491: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
495: vmath->objName = NULL;
496: return(0);
497: }
499: /*@C
500: PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica
502: Input Parameter:
503: . viewer - The Mathematica viewer
505: Output Parameter:
506: . v - The vector
508: Level: intermediate
510: .keywords PetscViewer, Mathematica, vector
511: .seealso VecView(), PetscViewerMathematicaPutVector()
512: @*/
513: PetscErrorCode PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
514: {
515: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
516: MLINK link; /* The link to Mathematica */
517: char *name;
518: PetscScalar *mArray,*array;
519: long mSize;
520: int n;
521: PetscErrorCode ierr;
527: /* Determine the object name */
528: if (!vmath->objName) name = "vec";
529: else name = (char*) vmath->objName;
531: link = vmath->link;
532: VecGetLocalSize(v, &n);
533: VecGetArray(v, &array);
534: MLPutFunction(link, "EvaluatePacket", 1);
535: MLPutSymbol(link, name);
536: MLEndPacket(link);
537: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
538: MLGetRealList(link, &mArray, &mSize);
539: if (n != mSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d",n,mSize);
540: PetscMemcpy(array, mArray, mSize * sizeof(double));
541: MLDisownRealList(link, mArray, mSize);
542: VecRestoreArray(v, &array);
543: return(0);
544: }
546: /*@C
547: PetscViewerMathematicaPutVector - Send a vector to Mathematica
549: Input Parameters:
550: + viewer - The Mathematica viewer
551: - v - The vector
553: Level: intermediate
555: .keywords PetscViewer, Mathematica, vector
556: .seealso VecView(), PetscViewerMathematicaGetVector()
557: @*/
558: PetscErrorCode PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
559: {
560: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
561: MLINK link = vmath->link; /* The link to Mathematica */
562: char *name;
563: PetscScalar *array;
564: int n;
565: PetscErrorCode ierr;
568: /* Determine the object name */
569: if (!vmath->objName) name = "vec";
570: else name = (char*) vmath->objName;
572: VecGetLocalSize(v, &n);
573: VecGetArray(v, &array);
575: /* Send the Vector object */
576: MLPutFunction(link, "EvaluatePacket", 1);
577: MLPutFunction(link, "Set", 2);
578: MLPutSymbol(link, name);
579: MLPutRealList(link, array, n);
580: MLEndPacket(link);
581: /* Skip packets until ReturnPacket */
582: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
583: /* Skip ReturnPacket */
584: MLNewPacket(link);
586: VecRestoreArray(v, &array);
587: return(0);
588: }
590: PetscErrorCode PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
591: {
592: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
593: MLINK link = vmath->link; /* The link to Mathematica */
594: char *name;
595: PetscErrorCode ierr;
598: /* Determine the object name */
599: if (!vmath->objName) name = "mat";
600: else name = (char*) vmath->objName;
602: /* Send the dense matrix object */
603: MLPutFunction(link, "EvaluatePacket", 1);
604: MLPutFunction(link, "Set", 2);
605: MLPutSymbol(link, name);
606: MLPutFunction(link, "Transpose", 1);
607: MLPutFunction(link, "Partition", 2);
608: MLPutRealList(link, a, m*n);
609: MLPutInteger(link, m);
610: MLEndPacket(link);
611: /* Skip packets until ReturnPacket */
612: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
613: /* Skip ReturnPacket */
614: MLNewPacket(link);
615: return(0);
616: }
618: PetscErrorCode PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
619: {
620: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
621: MLINK link = vmath->link; /* The link to Mathematica */
622: const char *symbol;
623: char *name;
624: PetscBool match;
625: PetscErrorCode ierr;
628: /* Determine the object name */
629: if (!vmath->objName) name = "mat";
630: else name = (char*) vmath->objName;
632: /* Make sure Mathematica recognizes sparse matrices */
633: MLPutFunction(link, "EvaluatePacket", 1);
634: MLPutFunction(link, "Needs", 1);
635: MLPutString(link, "LinearAlgebra`CSRMatrix`");
636: MLEndPacket(link);
637: /* Skip packets until ReturnPacket */
638: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
639: /* Skip ReturnPacket */
640: MLNewPacket(link);
642: /* Send the CSRMatrix object */
643: MLPutFunction(link, "EvaluatePacket", 1);
644: MLPutFunction(link, "Set", 2);
645: MLPutSymbol(link, name);
646: MLPutFunction(link, "CSRMatrix", 5);
647: MLPutInteger(link, m);
648: MLPutInteger(link, n);
649: MLPutFunction(link, "Plus", 2);
650: MLPutIntegerList(link, i, m+1);
651: MLPutInteger(link, 1);
652: MLPutFunction(link, "Plus", 2);
653: MLPutIntegerList(link, j, i[m]);
654: MLPutInteger(link, 1);
655: MLPutRealList(link, a, i[m]);
656: MLEndPacket(link);
657: /* Skip packets until ReturnPacket */
658: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
659: /* Skip ReturnPacket */
660: MLNewPacket(link);
662: /* Check that matrix is valid */
663: MLPutFunction(link, "EvaluatePacket", 1);
664: MLPutFunction(link, "ValidQ", 1);
665: MLPutSymbol(link, name);
666: MLEndPacket(link);
667: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
668: MLGetSymbol(link, &symbol);
669: PetscStrcmp("True", (char*) symbol, &match);
670: if (!match) {
671: MLDisownSymbol(link, symbol);
672: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
673: }
674: MLDisownSymbol(link, symbol);
675: /* Skip ReturnPacket */
676: MLNewPacket(link);
677: return(0);
678: }