Actual source code: mathematica.c
petsc-3.14.6 2021-03-30
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, sizeof(hostname));
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, sizeof(linkname), &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, sizeof(linkname), "%6d", ports[rank]);
220: else snprintf(linkname, sizeof(linkname), "%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], sizeof(hostname));
231: } else {
232: PetscStrncpy(hostname, hosts[0], sizeof(hostname));
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, sizeof(modename), &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, sizeof(type), &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, sizeof(type), &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: }