Actual source code: mathematica.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/viewerimpl.h> /* "petscsys.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;
17: /*@C
18: PetscViewerMathematicaFinalizePackage - This function destroys everything in the Petsc interface to Mathematica. It is
19: called from PetscFinalize().
21: Level: developer
23: .keywords: Petsc, destroy, package, mathematica
24: .seealso: PetscFinalize()
25: @*/
26: PetscErrorCode PetscViewerMathematicaFinalizePackage(void)
27: {
29: if (mathematicaEnv) MLDeinitialize((MLEnvironment) mathematicaEnv);
30: PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
31: return(0);
32: }
36: /*@C
37: PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
38: called from PetscDLLibraryRegister() when using dynamic libraries, and on the call to PetscInitialize()
39: when using static libraries.
41: Level: developer
43: .keywords: Petsc, initialize, package
44: .seealso: PetscSysInitializePackage(), PetscInitialize()
45: @*/
46: PetscErrorCode PetscViewerMathematicaInitializePackage(void)
47: {
48: PetscError ierr;
51: if (PetscViewerMathematicaPackageInitialized) return(0);
52: PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
54: mathematicaEnv = (void*) MLInitialize(0);
56: PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);
57: return(0);
58: }
63: PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
64: {
68: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return(0);
69: PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
70: return(0);
71: }
75: static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
76: {
77: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
78: PetscErrorCode ierr;
81: MLClose(vmath->link);
82: PetscFree(vmath->linkname);
83: PetscFree(vmath->linkhost);
84: PetscFree(vmath);
85: return(0);
86: }
90: PetscErrorCode PetscViewerDestroyMathematica_Private(void)
91: {
95: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
96: PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
97: }
98: return(0);
99: }
103: PetscErrorCode PetscViewerMathematicaSetupConnection_Private(PetscViewer v)
104: {
105: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
106: #if defined(MATHEMATICA_3_0)
107: int argc = 6;
108: char *argv[6];
109: #else
110: int argc = 5;
111: char *argv[5];
112: #endif
113: char hostname[256];
114: long lerr;
115: PetscErrorCode ierr;
118: /* Link name */
119: argv[0] = "-linkname";
120: if (!vmath->linkname) argv[1] = "math -mathlink";
121: else argv[1] = vmath->linkname;
123: /* Link host */
124: argv[2] = "-linkhost";
125: if (!vmath->linkhost) {
126: PetscGetHostName(hostname, 255);
127: argv[3] = hostname;
128: } else argv[3] = vmath->linkhost;
130: /* Link mode */
131: #if defined(MATHEMATICA_3_0)
132: argv[4] = "-linkmode";
133: switch (vmath->linkmode) {
134: case MATHEMATICA_LINK_CREATE:
135: argv[5] = "Create";
136: break;
137: case MATHEMATICA_LINK_CONNECT:
138: argv[5] = "Connect";
139: break;
140: case MATHEMATICA_LINK_LAUNCH:
141: argv[5] = "Launch";
142: break;
143: }
144: #else
145: switch (vmath->linkmode) {
146: case MATHEMATICA_LINK_CREATE:
147: argv[4] = "-linkcreate";
148: break;
149: case MATHEMATICA_LINK_CONNECT:
150: argv[4] = "-linkconnect";
151: break;
152: case MATHEMATICA_LINK_LAUNCH:
153: argv[4] = "-linklaunch";
154: break;
155: }
156: #endif
157: vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
158: #endif
159: return(0);
160: }
164: PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
165: {
166: PetscViewer_Mathematica *vmath;
167: PetscErrorCode ierr;
170: PetscViewerMathematicaInitializePackage();
172: PetscNewLog(v,&vmath);
173: v->data = (void*) vmath;
174: v->ops->destroy = PetscViewerDestroy_Mathematica;
175: v->ops->flush = 0;
176: PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);
178: vmath->linkname = NULL;
179: vmath->linkhost = NULL;
180: vmath->linkmode = MATHEMATICA_LINK_CONNECT;
181: vmath->graphicsType = GRAPHICS_MOTIF;
182: vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
183: vmath->objName = NULL;
185: PetscViewerMathematicaSetFromOptions(v);
186: PetscViewerMathematicaSetupConnection_Private(v);
187: return(0);
188: }
192: PetscErrorCode PetscViewerMathematicaParseLinkMode_Private(char *modename, LinkMode *mode)
193: {
194: PetscBool isCreate, isConnect, isLaunch;
198: PetscStrcasecmp(modename, "Create", &isCreate);
199: PetscStrcasecmp(modename, "Connect", &isConnect);
200: PetscStrcasecmp(modename, "Launch", &isLaunch);
201: if (isCreate) *mode = MATHEMATICA_LINK_CREATE;
202: else if (isConnect) *mode = MATHEMATICA_LINK_CONNECT;
203: else if (isLaunch) *mode = MATHEMATICA_LINK_LAUNCH;
204: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
205: return(0);
206: }
210: PetscErrorCode PetscViewerMathematicaSetFromOptions(PetscViewer v)
211: {
212: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
213: char linkname[256];
214: char modename[256];
215: char hostname[256];
216: char type[256];
217: PetscInt numPorts;
218: PetscInt *ports;
219: PetscInt numHosts;
220: int h;
221: char **hosts;
222: PetscMPIInt size, rank;
223: PetscBool opt;
224: PetscErrorCode ierr;
227: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
228: MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);
230: /* Get link name */
231: PetscOptionsGetString("viewer_", "-math_linkname", linkname, 256, &opt);
232: if (opt) {
233: PetscViewerMathematicaSetLinkName(v, linkname);
234: }
235: /* Get link port */
236: numPorts = size;
237: PetscMalloc1(size, &ports);
238: PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
239: if (opt) {
240: if (numPorts > rank) snprintf(linkname, 255, "%6d", ports[rank]);
241: else snprintf(linkname, 255, "%6d", ports[0]);
242: PetscViewerMathematicaSetLinkName(v, linkname);
243: }
244: PetscFree(ports);
245: /* Get link host */
246: numHosts = size;
247: PetscMalloc1(size, &hosts);
248: PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
249: if (opt) {
250: if (numHosts > rank) {
251: PetscStrncpy(hostname, hosts[rank], 255);
252: } else {
253: PetscStrncpy(hostname, hosts[0], 255);
254: }
255: PetscViewerMathematicaSetLinkHost(v, hostname);
256: }
257: for (h = 0; h < numHosts; h++) {
258: PetscFree(hosts[h]);
259: }
260: PetscFree(hosts);
261: /* Get link mode */
262: PetscOptionsGetString("viewer_", "-math_linkmode", modename, 256, &opt);
263: if (opt) {
264: LinkMode mode;
266: PetscViewerMathematicaParseLinkMode_Private(modename, &mode);
267: PetscViewerMathematicaSetLinkMode(v, mode);
268: }
269: /* Get graphics type */
270: PetscOptionsGetString("viewer_", "-math_graphics", type, 256, &opt);
271: if (opt) {
272: PetscBool isMotif, isPS, isPSFile;
274: PetscStrcasecmp(type, "Motif", &isMotif);
275: PetscStrcasecmp(type, "PS", &isPS);
276: PetscStrcasecmp(type, "PSFile", &isPSFile);
277: if (isMotif) vmath->graphicsType = GRAPHICS_MOTIF;
278: else if (isPS) vmath->graphicsType = GRAPHICS_PS_STDOUT;
279: else if (isPSFile) vmath->graphicsType = GRAPHICS_PS_FILE;
280: }
281: /* Get plot type */
282: PetscOptionsGetString("viewer_", "-math_type", type, 256, &opt);
283: if (opt) {
284: PetscBool isTri, isVecTri, isVec, isSurface;
286: PetscStrcasecmp(type, "Triangulation", &isTri);
287: PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
288: PetscStrcasecmp(type, "Vector", &isVec);
289: PetscStrcasecmp(type, "Surface", &isSurface);
290: if (isTri) vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
291: else if (isVecTri) vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
292: else if (isVec) vmath->plotType = MATHEMATICA_VECTOR_PLOT;
293: else if (isSurface) vmath->plotType = MATHEMATICA_SURFACE_PLOT;
294: }
295: return(0);
296: }
300: PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
301: {
302: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
303: PetscErrorCode ierr;
308: PetscStrallocpy(name, &vmath->linkname);
309: return(0);
310: }
314: PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
315: {
316: char name[16];
320: snprintf(name, 16, "%6d", port);
321: PetscViewerMathematicaSetLinkName(v, name);
322: return(0);
323: }
327: PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
328: {
329: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
330: PetscErrorCode ierr;
335: PetscStrallocpy(host, &vmath->linkhost);
336: return(0);
337: }
341: PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
342: {
343: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
346: vmath->linkmode = mode;
347: return(0);
348: }
350: /*----------------------------------------- Public Functions --------------------------------------------------------*/
353: /*@C
354: PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
356: Collective on comm
358: Input Parameters:
359: + comm - The MPI communicator
360: . port - [optional] The port to connect on, or PETSC_DECIDE
361: . machine - [optional] The machine to run Mathematica on, or NULL
362: - mode - [optional] The connection mode, or NULL
364: Output Parameter:
365: . viewer - The Mathematica viewer
367: Level: intermediate
369: Notes:
370: Most users should employ the following commands to access the
371: Mathematica viewers
372: $
373: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
374: $ MatView(Mat matrix, PetscViewer viewer)
375: $
376: $ or
377: $
378: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
379: $ VecView(Vec vector, PetscViewer viewer)
381: Options Database Keys:
382: + -viewer_math_linkhost <machine> - The host machine for the kernel
383: . -viewer_math_linkname <name> - The full link name for the connection
384: . -viewer_math_linkport <port> - The port for the connection
385: . -viewer_math_mode <mode> - The mode, e.g. Launch, Connect
386: . -viewer_math_type <type> - The plot type, e.g. Triangulation, Vector
387: - -viewer_math_graphics <output> - The output type, e.g. Motif, PS, PSFile
389: .keywords: PetscViewer, Mathematica, open
391: .seealso: MatView(), VecView()
392: @*/
393: PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
394: {
398: PetscViewerCreate(comm, v);
399: #if 0
400: LinkMode linkmode;
401: PetscViewerMathematicaSetLinkPort(*v, port);
402: PetscViewerMathematicaSetLinkHost(*v, machine);
403: PetscViewerMathematicaParseLinkMode_Private(mode, &linkmode);
404: PetscViewerMathematicaSetLinkMode(*v, linkmode);
405: #endif
406: PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
407: return(0);
408: }
412: /*@C
413: PetscViewerMathematicaGetLink - Returns the link to Mathematica
415: Input Parameters:
416: . viewer - The Mathematica viewer
417: . link - The link to Mathematica
419: Level: intermediate
421: .keywords PetscViewer, Mathematica, link
422: .seealso PetscViewerMathematicaOpen()
423: @*/
424: PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
425: {
426: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
430: *link = vmath->link;
431: return(0);
432: }
436: /*@C
437: PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
439: Input Parameters:
440: . viewer - The Mathematica viewer
441: . type - The packet type to search for, e.g RETURNPKT
443: Level: advanced
445: .keywords PetscViewer, Mathematica, packets
446: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
447: @*/
448: PetscErrorCode PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
449: {
450: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
451: MLINK link = vmath->link; /* The link to Mathematica */
452: int pkt; /* The packet type */
455: while ((pkt = MLNextPacket(link)) && (pkt != type)) MLNewPacket(link);
456: if (!pkt) {
457: MLClearError(link);
458: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, (char*) MLErrorMessage(link));
459: }
460: return(0);
461: }
465: /*@C
466: PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica
468: Input Parameter:
469: . viewer - The Mathematica viewer
471: Output Parameter:
472: . name - The name for new objects created in Mathematica
474: Level: intermediate
476: .keywords PetscViewer, Mathematica, name
477: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
478: @*/
479: PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
480: {
481: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
486: *name = vmath->objName;
487: return(0);
488: }
492: /*@C
493: PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica
495: Input Parameters:
496: . viewer - The Mathematica viewer
497: . name - The name for new objects created in Mathematica
499: Level: intermediate
501: .keywords PetscViewer, Mathematica, name
502: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
503: @*/
504: PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
505: {
506: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
511: vmath->objName = name;
512: return(0);
513: }
517: /*@C
518: PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
520: Input Parameter:
521: . viewer - The Mathematica viewer
523: Level: intermediate
525: .keywords PetscViewer, Mathematica, name
526: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
527: @*/
528: PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
529: {
530: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
534: vmath->objName = NULL;
535: return(0);
536: }
540: /*@C
541: PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica
543: Input Parameter:
544: . viewer - The Mathematica viewer
546: Output Parameter:
547: . v - The vector
549: Level: intermediate
551: .keywords PetscViewer, Mathematica, vector
552: .seealso VecView(), PetscViewerMathematicaPutVector()
553: @*/
554: PetscErrorCode PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v)
555: {
556: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
557: MLINK link; /* The link to Mathematica */
558: char *name;
559: PetscScalar *mArray,*array;
560: long mSize;
561: int n;
562: PetscErrorCode ierr;
568: /* Determine the object name */
569: if (!vmath->objName) name = "vec";
570: else name = (char*) vmath->objName;
572: link = vmath->link;
573: VecGetLocalSize(v, &n);
574: VecGetArray(v, &array);
575: MLPutFunction(link, "EvaluatePacket", 1);
576: MLPutSymbol(link, name);
577: MLEndPacket(link);
578: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
579: MLGetRealList(link, &mArray, &mSize);
580: if (n != mSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Incompatible vector sizes %d %d",n,mSize);
581: PetscMemcpy(array, mArray, mSize * sizeof(double));
582: MLDisownRealList(link, mArray, mSize);
583: VecRestoreArray(v, &array);
584: return(0);
585: }
589: /*@C
590: PetscViewerMathematicaPutVector - Send a vector to Mathematica
592: Input Parameters:
593: + viewer - The Mathematica viewer
594: - v - The vector
596: Level: intermediate
598: .keywords PetscViewer, Mathematica, vector
599: .seealso VecView(), PetscViewerMathematicaGetVector()
600: @*/
601: PetscErrorCode PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v)
602: {
603: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
604: MLINK link = vmath->link; /* The link to Mathematica */
605: char *name;
606: PetscScalar *array;
607: int n;
608: PetscErrorCode ierr;
611: /* Determine the object name */
612: if (!vmath->objName) name = "vec";
613: else name = (char*) vmath->objName;
615: VecGetLocalSize(v, &n);
616: VecGetArray(v, &array);
618: /* Send the Vector object */
619: MLPutFunction(link, "EvaluatePacket", 1);
620: MLPutFunction(link, "Set", 2);
621: MLPutSymbol(link, name);
622: MLPutRealList(link, array, n);
623: MLEndPacket(link);
624: /* Skip packets until ReturnPacket */
625: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
626: /* Skip ReturnPacket */
627: MLNewPacket(link);
629: VecRestoreArray(v, &array);
630: return(0);
631: }
633: PetscErrorCode PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a)
634: {
635: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
636: MLINK link = vmath->link; /* The link to Mathematica */
637: char *name;
638: PetscErrorCode ierr;
641: /* Determine the object name */
642: if (!vmath->objName) name = "mat";
643: else name = (char*) vmath->objName;
645: /* Send the dense matrix object */
646: MLPutFunction(link, "EvaluatePacket", 1);
647: MLPutFunction(link, "Set", 2);
648: MLPutSymbol(link, name);
649: MLPutFunction(link, "Transpose", 1);
650: MLPutFunction(link, "Partition", 2);
651: MLPutRealList(link, a, m*n);
652: MLPutInteger(link, m);
653: MLEndPacket(link);
654: /* Skip packets until ReturnPacket */
655: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
656: /* Skip ReturnPacket */
657: MLNewPacket(link);
658: return(0);
659: }
661: PetscErrorCode PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a)
662: {
663: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
664: MLINK link = vmath->link; /* The link to Mathematica */
665: const char *symbol;
666: char *name;
667: PetscBool match;
668: PetscErrorCode ierr;
671: /* Determine the object name */
672: if (!vmath->objName) name = "mat";
673: else name = (char*) vmath->objName;
675: /* Make sure Mathematica recognizes sparse matrices */
676: MLPutFunction(link, "EvaluatePacket", 1);
677: MLPutFunction(link, "Needs", 1);
678: MLPutString(link, "LinearAlgebra`CSRMatrix`");
679: MLEndPacket(link);
680: /* Skip packets until ReturnPacket */
681: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
682: /* Skip ReturnPacket */
683: MLNewPacket(link);
685: /* Send the CSRMatrix object */
686: MLPutFunction(link, "EvaluatePacket", 1);
687: MLPutFunction(link, "Set", 2);
688: MLPutSymbol(link, name);
689: MLPutFunction(link, "CSRMatrix", 5);
690: MLPutInteger(link, m);
691: MLPutInteger(link, n);
692: MLPutFunction(link, "Plus", 2);
693: MLPutIntegerList(link, i, m+1);
694: MLPutInteger(link, 1);
695: MLPutFunction(link, "Plus", 2);
696: MLPutIntegerList(link, j, i[m]);
697: MLPutInteger(link, 1);
698: MLPutRealList(link, a, i[m]);
699: MLEndPacket(link);
700: /* Skip packets until ReturnPacket */
701: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
702: /* Skip ReturnPacket */
703: MLNewPacket(link);
705: /* Check that matrix is valid */
706: MLPutFunction(link, "EvaluatePacket", 1);
707: MLPutFunction(link, "ValidQ", 1);
708: MLPutSymbol(link, name);
709: MLEndPacket(link);
710: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
711: MLGetSymbol(link, &symbol);
712: PetscStrcmp("True", (char*) symbol, &match);
713: if (!match) {
714: MLDisownSymbol(link, symbol);
715: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
716: }
717: MLDisownSymbol(link, symbol);
718: /* Skip ReturnPacket */
719: MLNewPacket(link);
720: return(0);
721: }