1
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>
6
7 #if defined(PETSC_HAVE__SNPRINTF) && !defined(PETSC_HAVE_SNPRINTF)
8 #define snprintf _snprintf
9 #endif
10
11 PetscViewer PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE = NULL;
12 static void *mathematicaEnv = NULL;
13
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().
18
19 Level: developer
20
21 .seealso: PetscFinalize()
22 @*/
PetscViewerMathematicaFinalizePackage(void)23 PetscErrorCode PetscViewerMathematicaFinalizePackage(void)
24 {
25 PetscFunctionBegin;
26 if (mathematicaEnv) MLDeinitialize((MLEnvironment) mathematicaEnv);
27 PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
28 PetscFunctionReturn(0);
29 }
30
31 /*@C
32 PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
33 called from PetscViewerInitializePackage().
34
35 Level: developer
36
37 .seealso: PetscSysInitializePackage(), PetscInitialize()
38 @*/
PetscViewerMathematicaInitializePackage(void)39 PetscErrorCode PetscViewerMathematicaInitializePackage(void)
40 {
41 PetscError ierr;
42
43 PetscFunctionBegin;
44 if (PetscViewerMathematicaPackageInitialized) PetscFunctionReturn(0);
45 PetscViewerMathematicaPackageInitialized = PETSC_TRUE;
46
47 mathematicaEnv = (void*) MLInitialize(0);
48
49 ierr = PetscRegisterFinalize(PetscViewerMathematicaFinalizePackage);CHKERRQ(ierr);
50 PetscFunctionReturn(0);
51 }
52
53
PetscViewerInitializeMathematicaWorld_Private()54 PetscErrorCode PetscViewerInitializeMathematicaWorld_Private()
55 {
56 PetscErrorCode ierr;
57
58 PetscFunctionBegin;
59 if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) PetscFunctionReturn(0);
60 ierr = PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, NULL, NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);CHKERRQ(ierr);
61 PetscFunctionReturn(0);
62 }
63
PetscViewerDestroy_Mathematica(PetscViewer viewer)64 static PetscErrorCode PetscViewerDestroy_Mathematica(PetscViewer viewer)
65 {
66 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
67 PetscErrorCode ierr;
68
69 PetscFunctionBegin;
70 MLClose(vmath->link);
71 ierr = PetscFree(vmath->linkname);CHKERRQ(ierr);
72 ierr = PetscFree(vmath->linkhost);CHKERRQ(ierr);
73 ierr = PetscFree(vmath);CHKERRQ(ierr);
74 PetscFunctionReturn(0);
75 }
76
PetscViewerDestroyMathematica_Private(void)77 PetscErrorCode PetscViewerDestroyMathematica_Private(void)
78 {
79 PetscErrorCode ierr;
80
81 PetscFunctionBegin;
82 if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
83 ierr = PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);CHKERRQ(ierr);
84 }
85 PetscFunctionReturn(0);
86 }
87
PetscViewerMathematicaSetupConnection_Private(PetscViewer v)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;
101
102 PetscFunctionBegin;
103 /* Link name */
104 argv[0] = "-linkname";
105 if (!vmath->linkname) argv[1] = "math -mathlink";
106 else argv[1] = vmath->linkname;
107
108 /* Link host */
109 argv[2] = "-linkhost";
110 if (!vmath->linkhost) {
111 ierr = PetscGetHostName(hostname, sizeof(hostname));CHKERRQ(ierr);
112 argv[3] = hostname;
113 } else argv[3] = vmath->linkhost;
114
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 PetscFunctionReturn(0);
145 }
146
PetscViewerCreate_Mathematica(PetscViewer v)147 PETSC_EXTERN PetscErrorCode PetscViewerCreate_Mathematica(PetscViewer v)
148 {
149 PetscViewer_Mathematica *vmath;
150 PetscErrorCode ierr;
151
152 PetscFunctionBegin;
153 ierr = PetscViewerMathematicaInitializePackage();CHKERRQ(ierr);
154
155 ierr = PetscNewLog(v,&vmath);CHKERRQ(ierr);
156 v->data = (void*) vmath;
157 v->ops->destroy = PetscViewerDestroy_Mathematica;
158 v->ops->flush = 0;
159 ierr = PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &((PetscObject)v)->type_name);CHKERRQ(ierr);
160
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;
167
168 ierr = PetscViewerMathematicaSetFromOptions(v);CHKERRQ(ierr);
169 ierr = PetscViewerMathematicaSetupConnection_Private(v);CHKERRQ(ierr);
170 PetscFunctionReturn(0);
171 }
172
PetscViewerMathematicaParseLinkMode(char * modename,LinkMode * mode)173 static PetscErrorCode PetscViewerMathematicaParseLinkMode(char *modename, LinkMode *mode)
174 {
175 PetscBool isCreate, isConnect, isLaunch;
176 PetscErrorCode ierr;
177
178 PetscFunctionBegin;
179 ierr = PetscStrcasecmp(modename, "Create", &isCreate);CHKERRQ(ierr);
180 ierr = PetscStrcasecmp(modename, "Connect", &isConnect);CHKERRQ(ierr);
181 ierr = PetscStrcasecmp(modename, "Launch", &isLaunch);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
187 }
188
PetscViewerMathematicaSetFromOptions(PetscViewer v)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;
204
205 PetscFunctionBegin;
206 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);CHKERRQ(ierr);
207 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank);CHKERRQ(ierr);
208
209 /* Get link name */
210 ierr = PetscOptionsGetString("viewer_", "-math_linkname", linkname, sizeof(linkname), &opt);CHKERRQ(ierr);
211 if (opt) {
212 ierr = PetscViewerMathematicaSetLinkName(v, linkname);CHKERRQ(ierr);
213 }
214 /* Get link port */
215 numPorts = size;
216 ierr = PetscMalloc1(size, &ports);CHKERRQ(ierr);
217 ierr = PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);CHKERRQ(ierr);
218 if (opt) {
219 if (numPorts > rank) snprintf(linkname, sizeof(linkname), "%6d", ports[rank]);
220 else snprintf(linkname, sizeof(linkname), "%6d", ports[0]);
221 ierr = PetscViewerMathematicaSetLinkName(v, linkname);CHKERRQ(ierr);
222 }
223 ierr = PetscFree(ports);CHKERRQ(ierr);
224 /* Get link host */
225 numHosts = size;
226 ierr = PetscMalloc1(size, &hosts);CHKERRQ(ierr);
227 ierr = PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);CHKERRQ(ierr);
228 if (opt) {
229 if (numHosts > rank) {
230 ierr = PetscStrncpy(hostname, hosts[rank], sizeof(hostname));CHKERRQ(ierr);
231 } else {
232 ierr = PetscStrncpy(hostname, hosts[0], sizeof(hostname));CHKERRQ(ierr);
233 }
234 ierr = PetscViewerMathematicaSetLinkHost(v, hostname);CHKERRQ(ierr);
235 }
236 for (h = 0; h < numHosts; h++) {
237 ierr = PetscFree(hosts[h]);CHKERRQ(ierr);
238 }
239 ierr = PetscFree(hosts);CHKERRQ(ierr);
240 /* Get link mode */
241 ierr = PetscOptionsGetString("viewer_", "-math_linkmode", modename, sizeof(modename), &opt);CHKERRQ(ierr);
242 if (opt) {
243 LinkMode mode;
244
245 ierr = PetscViewerMathematicaParseLinkMode(modename, &mode);CHKERRQ(ierr);
246 ierr = PetscViewerMathematicaSetLinkMode(v, mode);CHKERRQ(ierr);
247 }
248 /* Get graphics type */
249 ierr = PetscOptionsGetString("viewer_", "-math_graphics", type, sizeof(type), &opt);CHKERRQ(ierr);
250 if (opt) {
251 PetscBool isMotif, isPS, isPSFile;
252
253 ierr = PetscStrcasecmp(type, "Motif", &isMotif);CHKERRQ(ierr);
254 ierr = PetscStrcasecmp(type, "PS", &isPS);CHKERRQ(ierr);
255 ierr = PetscStrcasecmp(type, "PSFile", &isPSFile);CHKERRQ(ierr);
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 ierr = PetscOptionsGetString("viewer_", "-math_type", type, sizeof(type), &opt);CHKERRQ(ierr);
262 if (opt) {
263 PetscBool isTri, isVecTri, isVec, isSurface;
264
265 ierr = PetscStrcasecmp(type, "Triangulation", &isTri);CHKERRQ(ierr);
266 ierr = PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);CHKERRQ(ierr);
267 ierr = PetscStrcasecmp(type, "Vector", &isVec);CHKERRQ(ierr);
268 ierr = PetscStrcasecmp(type, "Surface", &isSurface);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
275 }
276
PetscViewerMathematicaSetLinkName(PetscViewer v,const char * name)277 PetscErrorCode PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name)
278 {
279 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
280 PetscErrorCode ierr;
281
282 PetscFunctionBegin;
283 PetscValidHeaderSpecific(v,PETSC_VIEWER_CLASSID,1);
284 PetscValidCharPointer(name,2);
285 ierr = PetscStrallocpy(name, &vmath->linkname);CHKERRQ(ierr);
286 PetscFunctionReturn(0);
287 }
288
PetscViewerMathematicaSetLinkPort(PetscViewer v,int port)289 PetscErrorCode PetscViewerMathematicaSetLinkPort(PetscViewer v, int port)
290 {
291 char name[16];
292 PetscErrorCode ierr;
293
294 PetscFunctionBegin;
295 snprintf(name, 16, "%6d", port);
296 ierr = PetscViewerMathematicaSetLinkName(v, name);CHKERRQ(ierr);
297 PetscFunctionReturn(0);
298 }
299
PetscViewerMathematicaSetLinkHost(PetscViewer v,const char * host)300 PetscErrorCode PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host)
301 {
302 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
303 PetscErrorCode ierr;
304
305 PetscFunctionBegin;
306 PetscValidHeaderSpecific(v,PETSC_VIEWER_CLASSID,1);
307 PetscValidCharPointer(host,2);
308 ierr = PetscStrallocpy(host, &vmath->linkhost);CHKERRQ(ierr);
309 PetscFunctionReturn(0);
310 }
311
PetscViewerMathematicaSetLinkMode(PetscViewer v,LinkMode mode)312 PetscErrorCode PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode)
313 {
314 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) v->data;
315
316 PetscFunctionBegin;
317 vmath->linkmode = mode;
318 PetscFunctionReturn(0);
319 }
320
321 /*----------------------------------------- Public Functions --------------------------------------------------------*/
322 /*@C
323 PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
324
325 Collective
326
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
332
333 Output Parameter:
334 . viewer - The Mathematica viewer
335
336 Level: intermediate
337
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)
349
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
357
358 .seealso: MatView(), VecView()
359 @*/
PetscViewerMathematicaOpen(MPI_Comm comm,int port,const char machine[],const char mode[],PetscViewer * v)360 PetscErrorCode PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v)
361 {
362 PetscErrorCode ierr;
363
364 PetscFunctionBegin;
365 ierr = PetscViewerCreate(comm, v);CHKERRQ(ierr);
366 #if 0
367 LinkMode linkmode;
368 ierr = PetscViewerMathematicaSetLinkPort(*v, port);CHKERRQ(ierr);
369 ierr = PetscViewerMathematicaSetLinkHost(*v, machine);CHKERRQ(ierr);
370 ierr = PetscViewerMathematicaParseLinkMode(mode, &linkmode);CHKERRQ(ierr);
371 ierr = PetscViewerMathematicaSetLinkMode(*v, linkmode);CHKERRQ(ierr);
372 #endif
373 ierr = PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);CHKERRQ(ierr);
374 PetscFunctionReturn(0);
375 }
376
377 /*@C
378 PetscViewerMathematicaGetLink - Returns the link to Mathematica
379
380 Input Parameters:
381 + viewer - The Mathematica viewer
382 - link - The link to Mathematica
383
384 Level: intermediate
385
386 .keywords PetscViewer, Mathematica, link
387 .seealso PetscViewerMathematicaOpen()
388 @*/
PetscViewerMathematicaGetLink(PetscViewer viewer,MLINK * link)389 PetscErrorCode PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
390 {
391 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
392
393 PetscFunctionBegin;
394 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID,1);
395 *link = vmath->link;
396 PetscFunctionReturn(0);
397 }
398
399 /*@C
400 PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
401
402 Input Parameters:
403 + viewer - The Mathematica viewer
404 - type - The packet type to search for, e.g RETURNPKT
405
406 Level: advanced
407
408 .keywords PetscViewer, Mathematica, packets
409 .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
410 @*/
PetscViewerMathematicaSkipPackets(PetscViewer viewer,int type)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 */
416
417 PetscFunctionBegin;
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 PetscFunctionReturn(0);
424 }
425
426 /*@C
427 PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica
428
429 Input Parameter:
430 . viewer - The Mathematica viewer
431
432 Output Parameter:
433 . name - The name for new objects created in Mathematica
434
435 Level: intermediate
436
437 .keywords PetscViewer, Mathematica, name
438 .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
439 @*/
PetscViewerMathematicaGetName(PetscViewer viewer,const char ** name)440 PetscErrorCode PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
441 {
442 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
443
444 PetscFunctionBegin;
445 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID,1);
446 PetscValidPointer(name,2);
447 *name = vmath->objName;
448 PetscFunctionReturn(0);
449 }
450
451 /*@C
452 PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica
453
454 Input Parameters:
455 + viewer - The Mathematica viewer
456 - name - The name for new objects created in Mathematica
457
458 Level: intermediate
459
460 .keywords PetscViewer, Mathematica, name
461 .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
462 @*/
PetscViewerMathematicaSetName(PetscViewer viewer,const char name[])463 PetscErrorCode PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
464 {
465 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
466
467 PetscFunctionBegin;
468 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID,1);
469 PetscValidPointer(name,2);
470 vmath->objName = name;
471 PetscFunctionReturn(0);
472 }
473
474 /*@C
475 PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
476
477 Input Parameter:
478 . viewer - The Mathematica viewer
479
480 Level: intermediate
481
482 .keywords PetscViewer, Mathematica, name
483 .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
484 @*/
PetscViewerMathematicaClearName(PetscViewer viewer)485 PetscErrorCode PetscViewerMathematicaClearName(PetscViewer viewer)
486 {
487 PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica*) viewer->data;
488
489 PetscFunctionBegin;
490 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID,1);
491 vmath->objName = NULL;
492 PetscFunctionReturn(0);
493 }
494
495 /*@C
496 PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica
497
498 Input Parameter:
499 . viewer - The Mathematica viewer
500
501 Output Parameter:
502 . v - The vector
503
504 Level: intermediate
505
506 .keywords PetscViewer, Mathematica, vector
507 .seealso VecView(), PetscViewerMathematicaPutVector()
508 @*/
PetscViewerMathematicaGetVector(PetscViewer viewer,Vec v)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;
518
519 PetscFunctionBegin;
520 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID,1);
521 PetscValidHeaderSpecific(v, VEC_CLASSID,2);
522
523 /* Determine the object name */
524 if (!vmath->objName) name = "vec";
525 else name = (char*) vmath->objName;
526
527 link = vmath->link;
528 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
529 ierr = VecGetArray(v, &array);CHKERRQ(ierr);
530 MLPutFunction(link, "EvaluatePacket", 1);
531 MLPutSymbol(link, name);
532 MLEndPacket(link);
533 ierr = PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);CHKERRQ(ierr);
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 ierr = PetscArraycpy(array, mArray, mSize);CHKERRQ(ierr);
537 MLDisownRealList(link, mArray, mSize);
538 ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
539 PetscFunctionReturn(0);
540 }
541
542 /*@C
543 PetscViewerMathematicaPutVector - Send a vector to Mathematica
544
545 Input Parameters:
546 + viewer - The Mathematica viewer
547 - v - The vector
548
549 Level: intermediate
550
551 .keywords PetscViewer, Mathematica, vector
552 .seealso VecView(), PetscViewerMathematicaGetVector()
553 @*/
PetscViewerMathematicaPutVector(PetscViewer viewer,Vec v)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;
562
563 PetscFunctionBegin;
564 /* Determine the object name */
565 if (!vmath->objName) name = "vec";
566 else name = (char*) vmath->objName;
567
568 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
569 ierr = VecGetArray(v, &array);CHKERRQ(ierr);
570
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 ierr = PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);CHKERRQ(ierr);
579 /* Skip ReturnPacket */
580 MLNewPacket(link);
581
582 ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
583 PetscFunctionReturn(0);
584 }
585
PetscViewerMathematicaPutMatrix(PetscViewer viewer,int m,int n,PetscReal * a)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;
592
593 PetscFunctionBegin;
594 /* Determine the object name */
595 if (!vmath->objName) name = "mat";
596 else name = (char*) vmath->objName;
597
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 ierr = PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);CHKERRQ(ierr);
609 /* Skip ReturnPacket */
610 MLNewPacket(link);
611 PetscFunctionReturn(0);
612 }
613
PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer,int m,int n,int * i,int * j,PetscReal * a)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;
622
623 PetscFunctionBegin;
624 /* Determine the object name */
625 if (!vmath->objName) name = "mat";
626 else name = (char*) vmath->objName;
627
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 ierr = PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);CHKERRQ(ierr);
635 /* Skip ReturnPacket */
636 MLNewPacket(link);
637
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 ierr = PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);CHKERRQ(ierr);
655 /* Skip ReturnPacket */
656 MLNewPacket(link);
657
658 /* Check that matrix is valid */
659 MLPutFunction(link, "EvaluatePacket", 1);
660 MLPutFunction(link, "ValidQ", 1);
661 MLPutSymbol(link, name);
662 MLEndPacket(link);
663 ierr = PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);CHKERRQ(ierr);
664 MLGetSymbol(link, &symbol);
665 ierr = PetscStrcmp("True", (char*) symbol, &match);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
674 }
675
676