Actual source code: mathematica.c
petsc-3.8.4 2018-03-24
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 PetscDLLibraryRegister() when using dynamic libraries, and on the call to PetscInitialize()
35: when using static libraries.
37: Level: developer
39: .keywords: Petsc, initialize, package
40: .seealso: PetscSysInitializePackage(), PetscInitialize()
41: @*/
42: PetscErrorCode PetscViewerMathematicaInitializePackage(void)
43: {
44: PetscError ierr;
47: if (PetscViewerMathematicaPackageInitialized) return(0);
48: PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
50: mathematicaEnv = (void*) MLInitialize(0);
52: PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
53: return(0);
54: }
57: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
58: {
62: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return(0);
63: PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
64: return(0);
65: }
67: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
68: {
69: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
70: PetscErrorCode ierr;
73: MLClose(vmath->link);
74: PetscFree(vmath->linkname);
75: PetscFree(vmath->linkhost);
76: PetscFree(vmath);
77: return(0);
78: }
80: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
81: {
85: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
86: PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
87: }
88: return(0);
89: }
91: PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
92: {
93: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
94: #if defined(MATHEMATICA_3_0)
95: int argc = 6;
96: char *argv[6];
97: #else
98: int argc = 5;
99: char *argv[5];
100: #endif
101: char hostname[256];
102: long lerr;
103: PetscErrorCode ierr;
106: /* Link name */
107: argv[0] = "-linkname";
108: if (!vmath->linkname) argv[1] = "math -mathlink";
109: else argv[1] = vmath->linkname;
111: /* Link host */
112: argv[2] = "-linkhost";
113: if (!vmath->linkhost) {
114: PetscGetHostName(hostname, 255);
115: argv[3] = hostname;
116: } else argv[3] = vmath->linkhost;
118: /* Link mode */
119: #if defined(MATHEMATICA_3_0)
120: argv[4] = "-linkmode";
121: switch (vmath->linkmode) {
122: case MATHEMATICA_LINK_CREATE:
123: argv[5] = "Create";
124: break;
125: case MATHEMATICA_LINK_CONNECT:
126: argv[5] = "Connect";
127: break;
128: case MATHEMATICA_LINK_LAUNCH:
129: argv[5] = "Launch";
130: break;
131: }
132: #else
133: switch (vmath->linkmode) {
134: case MATHEMATICA_LINK_CREATE:
135: argv[4] = "-linkcreate";
136: break;
137: case MATHEMATICA_LINK_CONNECT:
138: argv[4] = "-linkconnect";
139: break;
140: case MATHEMATICA_LINK_LAUNCH:
141: argv[4] = "-linklaunch";
142: break;
143: }
144: #endif
145: vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
146: #endif
147: return(0);
148: }
150: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
151: {
152: PetscViewer_Mathematica *vmath;
153: PetscErrorCode ierr;
156: PetscViewerMathematicaInitializePackage();
158: PetscNewLog(v,&vmath);
159: v->data = (void*) vmath;
160: v->ops->destroy = PetscViewerDestroy_Mathematica;
161: v->ops->flush = 0;
162: PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);
164: vmath->linkname = NULL;
165: vmath->linkhost = NULL;
166: vmath->linkmode = MATHEMATICA_LINK_CONNECT;
167: vmath->graphicsType = GRAPHICS_MOTIF;
168: vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
169: vmath->objName = NULL;
171: PetscViewerMathematicaSetFromOptions(v);
172: PetscViewerMathematicaSetupConnection_Private(v);
173: return(0);
174: }
176: PetscErrorCode PetscViewerMathematicaParseLinkMode_Private(char *modename, LinkMode *mode)
177: {
178: PetscBool isCreate, isConnect, isLaunch;
182: PetscStrcasecmp(modename, "Create", &isCreate);
183: PetscStrcasecmp(modename, "Connect", &isConnect);
184: PetscStrcasecmp(modename, "Launch", &isLaunch);
185: if (isCreate) *mode = MATHEMATICA_LINK_CREATE;
186: else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
187: else if (isLaunch) *mode = MATHEMATICA_LINK_LAUNCH;
188: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
189: return(0);
190: }
192: PetscErrorCode PetscViewerMathematicaSetFromOptions(PetscViewer v)
193: {
194: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
195: char linkname[256];
196: char modename[256];
197: char hostname[256];
198: char type[256];
199: PetscInt numPorts;
200: PetscInt *ports;
201: PetscInt numHosts;
202: int h;
203: char **hosts;
204: PetscMPIInt size, rank;
205: PetscBool opt;
206: PetscErrorCode ierr;
209: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
210: MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);
212: /* Get link name */
213: PetscOptionsGetString("viewer_", "-math_linkname", linkname, 256, &opt);
214: if (opt) {
215: PetscViewerMathematicaSetLinkName(v, linkname);
216: }
217: /* Get link port */
218: numPorts = size;
219: PetscMalloc1(size, &ports);
220: PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
221: if (opt) {
222: if (numPorts > rank) snprintf(linkname, 255, "%6d", ports[rank]);
223: else snprintf(linkname, 255, "%6d", ports[0]);
224: PetscViewerMathematicaSetLinkName(v, linkname);
225: }
226: PetscFree(ports);
227: /* Get link host */
228: numHosts = size;
229: PetscMalloc1(size, &hosts);
230: PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
231: if (opt) {
232: if (numHosts > rank) {
233: PetscStrncpy(hostname, hosts[rank], 255);
234: } else {
235: PetscStrncpy(hostname, hosts[0], 255);
236: }
237: PetscViewerMathematicaSetLinkHost(v, hostname);
238: }
239: for (h = 0; h < numHosts; h++) {
240: PetscFree(hosts[h]);
241: }
242: PetscFree(hosts);
243: /* Get link mode */
244: PetscOptionsGetString("viewer_", "-math_linkmode", modename, 256, &opt);
245: if (opt) {
246: LinkMode mode;
248: PetscViewerMathematicaParseLinkMode_Private(modename, &mode);
249: PetscViewerMathematicaSetLinkMode(v, mode);
250: }
251: /* Get graphics type */
252: PetscOptionsGetString("viewer_", "-math_graphics", type, 256, &opt);
253: if (opt) {
254: PetscBool isMotif, isPS, isPSFile;
256: PetscStrcasecmp(type, "Motif", &isMotif);
257: PetscStrcasecmp(type, "PS", &isPS);
258: PetscStrcasecmp(type, "PSFile", &isPSFile);
259: if (isMotif) vmath->graphicsType = GRAPHICS_MOTIF;
260: else if (isPS) vmath->graphicsType = GRAPHICS_PS_STDOUT;
261: else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
262: }
263: /* Get plot type */
264: PetscOptionsGetString("viewer_", "-math_type", type, 256, &opt);
265: if (opt) {
266: PetscBool isTri, isVecTri, isVec, isSurface;
268: PetscStrcasecmp(type, "Triangulation", &isTri);
269: PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
270: PetscStrcasecmp(type, "Vector", &isVec);
271: PetscStrcasecmp(type, "Surface", &isSurface);
272: if (isTri) vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
273: else if (isVecTri) vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
274: else if (isVec) vmath->plotType = MATHEMATICA_VECTOR_PLOT;
275: else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
276: }
277: return(0);
278: }
280: PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
281: {
282: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
283: PetscErrorCode ierr;
288: PetscStrallocpy(name, &vmath->linkname);
289: return(0);
290: }
292: PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
293: {
294: char name[16];
298: snprintf(name, 16, "%6d", port);
299: PetscViewerMathematicaSetLinkName(v, name);
300: return(0);
301: }
303: PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
304: {
305: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
306: PetscErrorCode ierr;
311: PetscStrallocpy(host, &vmath->linkhost);
312: return(0);
313: }
315: PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
316: {
317: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
320: vmath->linkmode = mode;
321: return(0);
322: }
324: /*----------------------------------------- Public Functions --------------------------------------------------------*/
325: /*@C
326: PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
328: Collective on comm
330: Input Parameters:
331: + comm - The MPI communicator
332: . port - [optional] The port to connect on, or PETSC_DECIDE
333: . machine - [optional] The machine to run Mathematica on, or NULL
334: - mode - [optional] The connection mode, or NULL
336: Output Parameter:
337: . viewer - The Mathematica viewer
339: Level: intermediate
341: Notes:
342: Most users should employ the following commands to access the
343: Mathematica viewers
344: $
345: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
346: $ MatView(Mat matrix, PetscViewer viewer)
347: $
348: $ or
349: $
350: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
351: $ VecView(Vec vector, PetscViewer viewer)
353: Options Database Keys:
354: + -viewer_math_linkhost <machine> - The host machine for the kernel
355: . -viewer_math_linkname <name> - The full link name for the connection
356: . -viewer_math_linkport <port> - The port for the connection
357: . -viewer_math_mode <mode> - The mode, e.g. Launch, Connect
358: . -viewer_math_type <type> - The plot type, e.g. Triangulation, Vector
359: - -viewer_math_graphics <output> - The output type, e.g. Motif, PS, PSFile
361: .keywords: PetscViewer, Mathematica, open
363: .seealso: MatView(), VecView()
364: @*/
365: PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
366: {
370: PetscViewerCreate(comm, v);
371: #if 0
372: LinkMode linkmode;
373: PetscViewerMathematicaSetLinkPort(*v, port);
374: PetscViewerMathematicaSetLinkHost(*v, machine);
375: PetscViewerMathematicaParseLinkMode_Private(mode, &linkmode);
376: PetscViewerMathematicaSetLinkMode(*v, linkmode);
377: #endif
378: PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
379: return(0);
380: }
382: /*@C
383: PetscViewerMathematicaGetLink - Returns the link to Mathematica
385: Input Parameters:
386: . viewer - The Mathematica viewer
387: . link - The link to Mathematica
389: Level: intermediate
391: .keywords PetscViewer, Mathematica, link
392: .seealso PetscViewerMathematicaOpen()
393: @*/
394: PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
395: {
396: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
400: *link = vmath->link;
401: return(0);
402: }
404: /*@C
405: PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
407: Input Parameters:
408: . viewer - The Mathematica viewer
409: . type - The packet type to search for, e.g RETURNPKT
411: Level: advanced
413: .keywords PetscViewer, Mathematica, packets
414: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
415: @*/
416: PetscErrorCode PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
417: {
418: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
419: MLINK link = vmath->link; /* The link to Mathematica */
420: int pkt; /* The packet type */
423: while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
424: if (!pkt) {
425: MLClearError(link);
426: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char*) MLErrorMessage(link));
427: }
428: return(0);
429: }
431: /*@C
432: PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica
434: Input Parameter:
435: . viewer - The Mathematica viewer
437: Output Parameter:
438: . name - The name for new objects created in Mathematica
440: Level: intermediate
442: .keywords PetscViewer, Mathematica, name
443: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
444: @*/
445: PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
446: {
447: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
452: *name = vmath->objName;
453: return(0);
454: }
456: /*@C
457: PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica
459: Input Parameters:
460: . viewer - The Mathematica viewer
461: . name - The name for new objects created in Mathematica
463: Level: intermediate
465: .keywords PetscViewer, Mathematica, name
466: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
467: @*/
468: PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
469: {
470: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
475: vmath->objName = name;
476: return(0);
477: }
479: /*@C
480: PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
482: Input Parameter:
483: . viewer - The Mathematica viewer
485: Level: intermediate
487: .keywords PetscViewer, Mathematica, name
488: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
489: @*/
490: PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
491: {
492: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
496: vmath->objName = NULL;
497: return(0);
498: }
500: /*@C
501: PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica
503: Input Parameter:
504: . viewer - The Mathematica viewer
506: Output Parameter:
507: . v - The vector
509: Level: intermediate
511: .keywords PetscViewer, Mathematica, vector
512: .seealso VecView(), PetscViewerMathematicaPutVector()
513: @*/
514: PetscErrorCode PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
515: {
516: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
517: MLINK link; /* The link to Mathematica */
518: char *name;
519: PetscScalar *mArray,*array;
520: long mSize;
521: int n;
522: PetscErrorCode ierr;
528: /* Determine the object name */
529: if (!vmath->objName) name = "vec";
530: else name = (char*) vmath->objName;
532: link = vmath->link;
533: VecGetLocalSize(v, &n);
534: VecGetArray(v, &array);
535: MLPutFunction(link, "EvaluatePacket", 1);
536: MLPutSymbol(link, name);
537: MLEndPacket(link);
538: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
539: MLGetRealList(link, &mArray, &mSize);
540: if (n != mSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d",n,mSize);
541: PetscMemcpy(array, mArray, mSize * sizeof(double));
542: MLDisownRealList(link, mArray, mSize);
543: VecRestoreArray(v, &array);
544: return(0);
545: }
547: /*@C
548: PetscViewerMathematicaPutVector - Send a vector to Mathematica
550: Input Parameters:
551: + viewer - The Mathematica viewer
552: - v - The vector
554: Level: intermediate
556: .keywords PetscViewer, Mathematica, vector
557: .seealso VecView(), PetscViewerMathematicaGetVector()
558: @*/
559: PetscErrorCode PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
560: {
561: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
562: MLINK link = vmath->link; /* The link to Mathematica */
563: char *name;
564: PetscScalar *array;
565: int n;
566: PetscErrorCode ierr;
569: /* Determine the object name */
570: if (!vmath->objName) name = "vec";
571: else name = (char*) vmath->objName;
573: VecGetLocalSize(v, &n);
574: VecGetArray(v, &array);
576: /* Send the Vector object */
577: MLPutFunction(link, "EvaluatePacket", 1);
578: MLPutFunction(link, "Set", 2);
579: MLPutSymbol(link, name);
580: MLPutRealList(link, array, n);
581: MLEndPacket(link);
582: /* Skip packets until ReturnPacket */
583: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
584: /* Skip ReturnPacket */
585: MLNewPacket(link);
587: VecRestoreArray(v, &array);
588: return(0);
589: }
591: PetscErrorCode PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
592: {
593: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
594: MLINK link = vmath->link; /* The link to Mathematica */
595: char *name;
596: PetscErrorCode ierr;
599: /* Determine the object name */
600: if (!vmath->objName) name = "mat";
601: else name = (char*) vmath->objName;
603: /* Send the dense matrix object */
604: MLPutFunction(link, "EvaluatePacket", 1);
605: MLPutFunction(link, "Set", 2);
606: MLPutSymbol(link, name);
607: MLPutFunction(link, "Transpose", 1);
608: MLPutFunction(link, "Partition", 2);
609: MLPutRealList(link, a, m*n);
610: MLPutInteger(link, m);
611: MLEndPacket(link);
612: /* Skip packets until ReturnPacket */
613: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
614: /* Skip ReturnPacket */
615: MLNewPacket(link);
616: return(0);
617: }
619: PetscErrorCode PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
620: {
621: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
622: MLINK link = vmath->link; /* The link to Mathematica */
623: const char *symbol;
624: char *name;
625: PetscBool match;
626: PetscErrorCode ierr;
629: /* Determine the object name */
630: if (!vmath->objName) name = "mat";
631: else name = (char*) vmath->objName;
633: /* Make sure Mathematica recognizes sparse matrices */
634: MLPutFunction(link, "EvaluatePacket", 1);
635: MLPutFunction(link, "Needs", 1);
636: MLPutString(link, "LinearAlgebra`CSRMatrix`");
637: MLEndPacket(link);
638: /* Skip packets until ReturnPacket */
639: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
640: /* Skip ReturnPacket */
641: MLNewPacket(link);
643: /* Send the CSRMatrix object */
644: MLPutFunction(link, "EvaluatePacket", 1);
645: MLPutFunction(link, "Set", 2);
646: MLPutSymbol(link, name);
647: MLPutFunction(link, "CSRMatrix", 5);
648: MLPutInteger(link, m);
649: MLPutInteger(link, n);
650: MLPutFunction(link, "Plus", 2);
651: MLPutIntegerList(link, i, m+1);
652: MLPutInteger(link, 1);
653: MLPutFunction(link, "Plus", 2);
654: MLPutIntegerList(link, j, i[m]);
655: MLPutInteger(link, 1);
656: MLPutRealList(link, a, i[m]);
657: MLEndPacket(link);
658: /* Skip packets until ReturnPacket */
659: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
660: /* Skip ReturnPacket */
661: MLNewPacket(link);
663: /* Check that matrix is valid */
664: MLPutFunction(link, "EvaluatePacket", 1);
665: MLPutFunction(link, "ValidQ", 1);
666: MLPutSymbol(link, name);
667: MLEndPacket(link);
668: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
669: MLGetSymbol(link, &symbol);
670: PetscStrcmp("True", (char*) symbol, &match);
671: if (!match) {
672: MLDisownSymbol(link, symbol);
673: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
674: }
675: MLDisownSymbol(link, symbol);
676: /* Skip ReturnPacket */
677: MLNewPacket(link);
678: return(0);
679: }