1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
|
<html><head><title></title></head>
<body>
<a name="934888">
<h1>1.0 </a>Introduction</h1>
</a>
<a name="934254">
There are many mathematics libraries currently enjoying widespread use. At the time that </a>PML was begun, however, there were very few written in C and fewer still written with </a>portability in mind. The temptation has all too often been to program to a specific hardware platform and forget portability. More recently some very good work has been published by Press, Flannery, Teulkovsky, and Vetterling, </a>Numerical Recipes in C, which is highly portable. The reason for writing PML hasn’t changed in spite this. No one library is apt to have everything that a given application needs and there are other considerations such as interfaces and related functionality to drive the development of a math library.<p>
</a>
<a name="934220">
PML (</a>Portable Mathematics Library) was developed initially as a repository for any routines of a mathematical nature in what has become the </a>PACT system. More recently, I have tried to systematize the library and organize it for future growth.<p>
</a>
<a name="934237">
The library has two fundamental parts. First are the many pure math routines such as equation solvers. Second are the structure definitions and some related functions. The structures which PML supports are provided to aid C programmers in organizing their mathematical computations in a mathematical way. That is, they are encouraged to give primary consideration to the mathematics not the implementation of the mathematics. This also is intended to add to the mathematical rigor of </a>simulations and promote the use of some of the more abstract but critical ideas of </a>mathematics by giving application developers concrete representations of </a>abstract structures.<p>
</a>
<a name="934250">
The data structures for PML are presented first since they support most of the library. A summary of the mathematical functions comes next. A set of examples completes this manual.<p>
</a>
<a name="934889">
I would like to acknowledge the work of Charles F. McMillan in PML. He devised the original </a>matrix structure and wrote many of the routines to manipulate them. The final form that they take in PML is very close to his original conception.<p>
</a>
<a name="934239">
<h1>2.0 </a>Data Structures</h1>
</a>
<a name="934890">
</a>PML currently supports the following structures: </a>complex, </a>matrix; </a>set; </a>mapping; and </a>field. The </a>complex structure and its associated functions provide a rudimentary facility to do arithmetic with complex numbers. The </a>matrix structure contains “all” of the information needed to deal with first and second rank </a>tensors (i.e. </a>vectors and </a>matrices). The </a>set structure is intended to describe a discrete set of </a>elements. It includes information about the type of elements, the dimensionality, topology, and so on. Once the notion of a set has been built up, it is natural to proceed to discuss mappings of sets. This is obviously useful for such tasks as </a>interpolations between </a>computational </a>meshes and </a>visualization. The </a>mapping structure includes a domain and range set and information describing the mapping between them. To support the mapping ideas, it is necessary to consider whether the sets involved are </a></a>algebras, </a>rings, </a>fields, or whatever. To facilitate this thinking, a </a>field structure is defined to provide the set of operations which can be carried out on the elements of sets. This also makes the set and mapping structures very general (which it should be) and still reasonably efficient.<p>
</a>
<a name="934248">
<h2>2.1 </a>Complex</h2>
</a>
<a name="934276">
The </a>complex structure is the natural representation of a complex double precision number.<p>
</a>
<A NAME="934277"><PRE>
</PRE><A NAME="934309"><PRE> struct s_complex
</PRE><A NAME="934278"><PRE> {double real;
</PRE><A NAME="934279"><PRE> double imag;};
</PRE><A NAME="934266"><PRE> typedef struct s_complex complex;
</PRE><a name="934384">
The </a>real part of the complex number is kept in the member </a>real and the </a>imaginary part of the complex number is kept in the member </a>imag.<p>
</a>
<a name="934252">
<h2>2.2 </a>PM_matrix</h2>
</a>
<a name="934227">
The </a>PM_matrix structure is a useful way to group a contiguous block of memory together with two integers which describe the </a>logical </a>shape of the space. <p>
</a>
<A NAME="934264"><PRE>
</PRE><A NAME="934244"><PRE> struct s_PM_matrix
</PRE><A NAME="934246"><PRE> {int nrow;
</PRE><A NAME="934253"><PRE> int ncol;
</PRE><A NAME="934241"><PRE> REAL *array;};
</PRE><A NAME="934263"><PRE> typedef struct s_PM_matrix PM_matrix;
</PRE><a name="934256">
The </a>nrow member is the </a>number of </a>rows of the matrix and the </a>ncol member is the </a>number of </a>columns of the matrix. The </a>array member is a pointer to the actual space. The type declaration </a>REAL is set in the PACT system to be either </a>float or </a>double. Typically it is set to double, and float is used mainly on the smaller machines to conserve memory.<p>
</a>
<a name="934243">
In this implementation, matrices are assumed to have a </a>one based </a>indexing scheme. That is, all matrix indices go from 1 to ncol or from 1 to nrow. <p>
</a>
<a name="934329">
When operating on matrices it is necessary only to pass the matrix as an argument since the number of rows and columns is contained in the structure. A </a>vector is taken to be simply a matrix with either one row or one column.<p>
</a>
<a name="934229">
<h2>2.3 </a>Set</h2>
</a>
<a name="934330">
The </a>PM_set structure is an attempt to describe a usefully general mathematical set for such purposes as </a>interpolation and </a>visualization. It encapsulates a collection of information about a discrete set of </a>elements of arbitrary type.<p>
</a>
<A NAME="934273"><PRE>
</PRE><A NAME="934300"><PRE> struct s_PM_set
</PRE><A NAME="934332"><PRE> {char *name;
</PRE><A NAME="934333"><PRE> char *element_type;
</PRE><A NAME="934334"><PRE> int dimension;
</PRE><A NAME="934335"><PRE> int *max_index;
</PRE><A NAME="934336"><PRE> int dimension_elem;
</PRE><A NAME="934337"><PRE> long n_elements;
</PRE><A NAME="934232"><PRE> byte *elements;
</PRE><A NAME="934236"><PRE> char *es_type;
</PRE><A NAME="934257"><PRE> byte *extrema;
</PRE><A NAME="934260"><PRE> byte *scales;
</PRE><A NAME="934261"><PRE> PM_field *opers;
</PRE><A NAME="934339"><PRE> REAL *metric;
</PRE><A NAME="934342"><PRE> char *symmetry_type;
</PRE><A NAME="934343"><PRE> byte *symmetry;
</PRE><A NAME="934344"><PRE> char *topology_type;
</PRE><A NAME="934345"><PRE> byte *topology;
</PRE><A NAME="934349"><PRE> char *info_type;
</PRE><A NAME="934350"><PRE> byte *info;};
</PRE><A NAME="934352"><PRE> typedef struct s_PM_set PM_set;
</PRE><a name="934265">
The </a>PM_set members are: </a>name, every set has a unique identifier; </a>element_type, the type of element in the set; </a>dimension, the dimension of the set; </a>max_index, an array of the lengths along each dimension (it is dimension long); </a>dimension_elem, the dimensionality of the elements of the set</a>; n_elements, the </a>number of elements in the set; </a>elements, a pointer to the actual elements of the set; </a>es_type, the type of the extrema and the scales arrays (typically, element_type is double ** while es_type is double *); </a>extrema, an array of minimum and maximum values for each component (2*dimension_elem long); </a>scales, an array of scales for each dimension; </a>opers, specifies the </a>field for the set; </a>metric, an appropriate array of metric coefficients for the set if it is a metric set</a>; symmetry_type, an ASCII string naming the type of the data in the symmetry member; </a>symmetry, a pointer to data describing symmetries of the set, if any exist; </a>topology_type, an ASCII string naming the type of the data in the </a>topology member; </a>topology, a pointer to data describing the </a>connectivity of the elements of the set; </a>opers, a collection of functions describing defined </a>operations on elements of the set; </a>info_type, an ASCII string naming the type of the data in the info member; and </a>info, a pointer to data describing any additional information about the set.<p>
</a>
<a name="934395">
In order to be general and useful, many types of information are described by two members of the </a>PM_set. For example, an application may have a mesh (handled as a set) with an array of integers to define the </a>connectivity. The </a>topology member should point to that array and the </a>topology_type member should be set to integer or int. The application would have easy access to the connectivity and another application which might work with sets of other kinds could be programmed to distinguish between the </a>topological information associated with these mesh sets and that of other </a>sets it knows about.<p>
</a>
<a name="934394">
<h3>2.3.1 </a>Example of </a>building a </a>set</h3>
</a>
<a name="934385">
This example is taken from the source for PML and illustrates how a set can be built. The routine uses the </a>SCORE functions </a>MAKE, </a>MAKE_N, </a>SC_arrlen, and </a>SC_strsave as well as the SCORE variable </a>SC_DOUBLE_S.<p>
</a>
<A NAME="934393"><PRE>
</PRE><A NAME="934242"><PRE> #include “pml.h”
</PRE><A NAME="934292"><PRE>
</PRE><A NAME="934267"><PRE> main()
</PRE><A NAME="934270"><PRE> {int i, k, l, kmax, lmax, kxl;
</PRE><A NAME="934285"><PRE> REAL *x, *y, *u, *v, r, t;
</PRE><A NAME="934286"><PRE> REAL xmin, xmax, ymin, ymax, km, lm;
</PRE><A NAME="934288"><PRE> PM_set *domain, *range;
</PRE><A NAME="934341"><PRE>
</PRE><A NAME="934346"><PRE> /* set up data */
</PRE><A NAME="934347"><PRE> kmax = 20;
</PRE><A NAME="934408"><PRE> lmax = 20;
</PRE><A NAME="934409"><PRE> xmin = -5.0;
</PRE><A NAME="934410"><PRE> xmax = 5.0;
</PRE><A NAME="934411"><PRE> ymin = -5.0;
</PRE><A NAME="934413"><PRE> ymax = 5.0;
</PRE><A NAME="934420"><PRE>
</PRE><A NAME="934421"><PRE> kxl = kmax*lmax;
</PRE><A NAME="934422"><PRE> x = MAKE_N(REAL, kxl);
</PRE><A NAME="934423"><PRE> y = MAKE_N(REAL, kxl);
</PRE><A NAME="934424"><PRE> u = MAKE_N(REAL, kxl);
</PRE><A NAME="934425"><PRE> v = MAKE_N(REAL, kxl);
</PRE><A NAME="934428"><PRE>
</PRE><A NAME="934429"><PRE> km = 2.0*PI/((double) (kmax - 1));
</PRE><A NAME="934430"><PRE> lm = 2.0*PI/((double) (lmax - 1));
</PRE><A NAME="934431"><PRE> for (k = 0; k < kmax; k++)
</PRE><A NAME="934435"><PRE> for (l = 0; l < lmax; l++)
</PRE><A NAME="934436"><PRE> {i = l*kmax + k;
</PRE><A NAME="934437"><PRE> x[i] = k*km;
</PRE><A NAME="934439"><PRE> y[i] = l*lm;
</PRE><A NAME="934443"><PRE> u[i] = sin(k*km);
</PRE><A NAME="934459"><PRE> v[i] = sin(l*lm);};
</PRE><A NAME="934480"><PRE>
</PRE><A NAME="934483"><PRE> /* build the domain set */
</PRE><A NAME="934487"><PRE> domain = PM_make_set(“{x, y}”, SC_REAL_S, FALSE,
</PRE><A NAME="934491"><PRE> 2, kmax, lmax, 2, x, y);
</PRE><A NAME="934498"><PRE>
</PRE><A NAME="934500"><PRE> /* build the range set */
</PRE><A NAME="934502"><PRE> range = PM_make_set(“{u, v}”, SC_REAL_S, FALSE,
</PRE><A NAME="934504"><PRE> 2, kmax, lmax, 2, u, v);
</PRE><A NAME="934511"><PRE>
</PRE><A NAME="934597"><PRE> .
</PRE><A NAME="934601"><PRE> .
</PRE><A NAME="934271"><PRE> .
</PRE><A NAME="934287"><PRE>
</PRE><A NAME="934607"><PRE> return(0);}
</PRE><A NAME="934615"><PRE>
</PRE><a name="934233">
<h2>2.4 </a>Mapping</h2>
</a>
<a name="934584">
Once the </a>PM_set structure is there to organize data into some semblance of a mathematical </a>set, it is quite natural to express objects such as </a>computational meshes as sets. In fact, considering the needs of </a>numerical </a>simulation codes and </a>visualization codes, the concepts of </a>domain sets and </a>range sets arise and lead to the need for a structure to characterize a mathematical mapping.<p>
</a>
<A NAME="934371"><PRE>
</PRE><A NAME="934353"><PRE> struct s_PM_mapping
</PRE><A NAME="934255"><PRE> {char *name;
</PRE><A NAME="934268"><PRE> char *category;
</PRE><A NAME="934272"><PRE> PM_set *domain;
</PRE><A NAME="934374"><PRE> PM_set *range;
</PRE><A NAME="934375"><PRE> char *map_type;
</PRE><A NAME="934274"><PRE> byte *map;
</PRE><A NAME="934280"><PRE> in file_type;
</PRE><A NAME="934282"><PRE> byte *file_info;
</PRE><A NAME="934283"><PRE> char *file;
</PRE><A NAME="934380"><PRE> struct s_PM_mapping *next;};
</PRE><A NAME="934382"><PRE> typedef struct s_PM_mapping PM_mapping;
</PRE><a name="934574">
The members of the </a>PM_mapping are: </a>name, a name for the mapping; </a>category, the category to which the mapping belongs (e.g., Logical-Rectangular, Arbitrarily-Connected); </a>domain, the mapping’s domain; </a>range, (a subset of) the image of domain under map; </a>map_type, an ASCII string naming the type of the data pointed to by map; </a>map, a pointer to data describing the mapping between domain and range; </a>file_type, file type </a>ASCII, </a>BINARY, or </a>PDB; </a>file_info, </a>file information (cast to some struct with information); </a>file, the file name for the mapping; and </a>next; the next PM_mapping in the </a>chain.<p>
</a>
<a name="934377">
At this point, many of these ideas are evolving. To insulate the developer somewhat against these changes, the map member was used. It allows different structures to be used to capture the description of the mapping from domain to range. For example, to handle the case when one has a domain and a function, the range could be NULL. Another potential issue to handle is whether the mapping is a </a>bijection, </a>surjection, or </a>injection.<p>
</a>
<a name="934379">
A very immediate issue that arises instantly is the </a>centering of the range set relative to the domain set. For example, the domain might be a </a>mesh with </a>coordinates defining </a>node centers while the range set might be a </a>zone centered quantity. One small structure is provided to address this issue:<p>
</a>
<A NAME="934569"><PRE>
</PRE><A NAME="934284"><PRE> struct s_PM_map_info
</PRE><A NAME="934366"><PRE> {char *name;
</PRE><A NAME="934368"><PRE> int centering;};
</PRE><A NAME="934370"><PRE> typedef struct s_PM_map_info PM_map_info;
</PRE><a name="934580">
The </a>PM_map_info structure only has a </a>name member and an integer member, </a>centering, to specify the relative centering as being one of: </a>Z_CENT for </a>zone centered; N_CENT for </a>node centered; </a>F_CENT for </a>face centered; </a>E_CENT for </a>edge centered; or </a>U_CENT for </a>uncentered.<p>
</a>
<a name="934400">
<h2>2.5 Connectivity</h2>
</a>
<a name="934402">
To support general computational meshes PML has a model for arbitrarily connected meshes which is independent of dimensionality and a priori assumptions about shapes of zones or elements. We give a complete discussion below including definitions of terms. Hopefully this will suffice to give a useful account of this topic.<p>
</a>
<a name="934427">
<h3>2.5.1 Concepts</h3>
</a>
<a name="934434">
The goal is to capture non-trivial topological information describing computational meshes which is efficient in space usage and amenable to fast algorithms for visualization rendering and analytical operations.<p>
</a>
<a name="934447">
The representation of the connectivity scheme must accomodate a heirarchy of information. At the bottom of the heirarchy is the minimal information needed by a minimal connectivity scheme. At the top is a maximal set of information which is computable from the minimal information but is kept in memory rather than recomputing it.<p>
</a>
<a name="934471">
<h3>2.5.2 Definitions</h3>
</a>
<A NAME="934479"><PRE> </PRE> Computational mesh or mesh - a collection of discrete points in some vector space
along with the specification of relationships between points especially the neighbor
relationship.
<BR><A NAME="934489"><PRE> </PRE>Connectivity - the pattern of specifications of node neighbors.
<BR><A NAME="934497"><PRE> </PRE>N-cell - an N dimensional volume unit. Example: 3 dimensional mesh. A zone is a 3-
cell, a face is a 2-cell, an edge is a 1-cell, and a node is a 0-cell.
<BR><A NAME="934503"><PRE> </PRE>N-Center - the center of an N-cell is defined by some weighted average of the real
nodes which imply the N-cell. A 0-Center corresponds to a real node.
<BR><A NAME="934514"><PRE> </PRE>Boundary - the boundary of an N-cell is a set of (N-1)-cells such that the boundary of
the boundary is 0. Although each (N-1)-cell may have a boundary, the set of them may
not. This is a closure property. Each boundary (N-1)-cell has an orientation. The outward normal is the canonical orientation.
<BR><A NAME="934528"><PRE> </PRE>Zone - for a mesh in an N dimensional vector space a zone is the smallest N dimensional volume unit which is delineated by real nodes. In an N dimensional mesh, a zone
is an N-cell.
<BR><A NAME="934532"><PRE> </PRE>Node - for a mesh in an N dimensional vector space a node is a single point in the N
dimensional vector space. A node is always a 0-cell. The collection of nodes is one of
the defining elements of the mesh.
<BR><A NAME="934541"><PRE> </PRE>Real Node - a node in the mesh which is explicitly specified and is a part of the connectivity of the mesh.
<BR><A NAME="934544"><PRE> </PRE>Virtual Node - a node which is implicitly defined by the mesh connectivity scheme. For
example, the points defining the zone centers may be virtual nodes.
<BR><A NAME="934561"><PRE> </PRE>Decomposition - an N-cell which is implied by a set of real nodes can be expressed as a
sum of N-cells which are implied by the set of real nodes and at least one virtual node.
<BR><A NAME="934568"><PRE> </PRE>Irreducible Decomposition - an N-cell has a special, unique decomposition which is the
sum of N-cells each of which has: 1 0-Center, 1 1-Center, ..., 1 N-Center in their boundary and which collectively cover the N-cell.
<BR><a name="934269">
<p>
</a>
<a name="934290">
<p>
</a>
<a name="934291">
<p>
</a>
<a name="934293">
<p>
</a>
<a name="934294">
<p>
</a>
<a name="934295">
<p>
</a>
<a name="934296">
<p>
</a>
<a name="934297">
<p>
</a>
<a name="934298">
<p>
</a>
<a name="934299">
<p>
</a>
<A NAME="934314"><B>Figure: The cell A is one cell in the irreducible decomposition of the rectangular 2-Cell shown
</B><HR><a name="934575">
<h3>2.5.3 Axioms</h3>
</a>
<ul><a name="934577">
<li>The mimimum information necessary to describe a computational mesh is: 1) a list of nodes; and 2) for each node a list of its neighboring nodes.
</a>
<a name="934583">
<li>The information that will typically be available in a mesh description is: 1) a list of nodes; 2) a list of zones; and 3) for each node a list of its neighboring nodes.
</a>
</ul><a name="934599">
<h3>2.5.4 Representation Requirements</h3>
</a>
<a name="934602">
To represent an N-cell obviously depends on how much information you want to maintain.<p>
</a>
<a name="934405">
<p>
</a>
<ul><a name="934608">
<li>Minimal:
</a>
<A NAME="934404"><PRE> Boundary specification (necessary and sufficient)
</PRE><a name="934609">
<li>Maximal:
</a>
<A NAME="934403"><PRE> Center node
</PRE><A NAME="934610"><PRE> (N+1)-cell for which this N-cell is a boundary segment
</PRE><A NAME="934611"><PRE> Neighboring (N+1)-cell
</PRE><A NAME="934612"><PRE> Decomposition specification
</PRE><A NAME="934613"><PRE> Irreducible Decomposition specification
</PRE><a name="934614">
A flexible approach is taken by representing a cell as a block of contiguous integers in an array of cells. This means that an array of cells is in effect a two dimensional array of integers. The fastest changing dimension is bounded by the number of parameters used in representing the cells and the slowest changing dimension is the cell number.<p>
</a>
<a name="934623">
PML defines the following constants to name the cell description parameters<p>
</a>
<A NAME="934462">#define BND_CELL_MIN 0
<P><A NAME="934418">#define BND_CELL_MAX 1
<P><A NAME="934433">#define OPPOSITE_CELL 2
<P><A NAME="934441">#define PARENT_CELL 3
<P><A NAME="934444">#define NGB_CELL 4
<P><A NAME="934446">#define CENTER_CELL 5
<P><A NAME="934451">#define DEC_CELL_MIN 6
<P><A NAME="934456">#define DEC_CELL_MAX 7
<P></ul><a name="934626">
<h3>2.5.5 Representation of Connectivity</h3>
</a>
<a name="934627">
The PM_mesh_topology structure attempts to describe in complete generality a mesh of points and their neighbors as a heirarchy in which the highest dimensional volumes sit at the top and the line segment descriptions sit at the bottom.<p>
</a>
<a name="934463">
At each level in the heirarchy an N-cell is described by a contiguous block of (N-1)-cell boundary sections (next lower down in the heirarchy).<p>
</a>
<a name="934507">
Assuming the dimensionality to be n the heirarchy goes like:<p>
</a>
<A NAME="934515">boundaries[n] the description for the n-dim surfaces
<P><A NAME="934516">boundaries[n-1] the description for the (n-1)-dim surfaces
<P><A NAME="934518"> .
<P><A NAME="934468"> .
<P><A NAME="934476"> .
<P><A NAME="934530">boundaries[1] the description for the 1-dim surfaces, i.e. line segments whose endpoints are the node indices
<P><A NAME="934538">boundaries[0] a mimimal description for the nodes
<P><a name="934542">
<p>
</a>
<A NAME="934543"><PRE> struct s_PM_mesh_topology
</PRE><A NAME="934550"><PRE> {int n_dimensions; /* number of dimensions */
</PRE><A NAME="934555"><PRE> int *n_bound_params; /* number of bound params at each level */
</PRE><A NAME="934557"><PRE> int *n_cells; /* number of cells at each level */
</PRE><A NAME="934562"><PRE> long **boundaries;}; /* boundary info array at each level */
</PRE><a name="934563">
<p>
</a>
<a name="934565">
typedef struct s_PM_mesh_topology PM_mesh_topology;<p>
</a>
<a name="934570">
For example, in a 3 dimensional problem, the “zones” at the top of the hierarchy are 3-cells. Each “zone” has some number of “faces” or 2-cells bounding it. Each “face” has some number of “edges” or 1-cells bounding it. Implicitly, each “edge” has two points or 0-cells bounding it. Then boundaries[3] would contain the indices into boundaries[2] of boundary faces and, if that level of information is kept, indices of other zones in boundaries[3] which are neighbors.<p>
</a>
<a name="934638">
<p>
</a>
<a name="934315">
<p>
</a>
<a name="934316">
<p>
</a>
<a name="934317">
<p>
</a>
<a name="934318">
<p>
</a>
<a name="934319">
<p>
</a>
<a name="934320">
<p>
</a>
<a name="934414">
<p>
</a>
<a name="934415">
<p>
</a>
<a name="934416">
<p>
</a>
<a name="934417">
<p>
</a>
<a name="934419">
<p>
</a>
<A NAME="934440"><B>Figure: Diagram of hierarchy for 3 dimensional mesh connectivity
</B><HR><a name="934519">
<p>
</a>
<a name="934240">
<h2>2.6 </a>Field</h2>
</a>
<a name="934258">
Depending on which </a>operations are defined on a </a>set of </a>elements (in addition to other properties), one has a </a>vector space, an </a>algebra, a </a>ring, a </a>group, a </a>field, or variations on that same theme. As a practical matter, given a set of elements, applications need to know how to combine them (add, subtract, etc.). The </a>PM_field structure provides a way to make those connections in a generic way.<p>
</a>
<A NAME="934355"><PRE>
</PRE><A NAME="934259"><PRE> struct s_PM_field
</PRE><A NAME="934262"><PRE> {byte (*add)();
</PRE><A NAME="934356"><PRE> byte (*sub)();
</PRE><A NAME="934357"><PRE> byte (*scalar_mult)();
</PRE><A NAME="934358"><PRE> byte (*mult)();
</PRE><A NAME="934359"><PRE> byte (*div)();};
</PRE><A NAME="934579"><PRE> typedef struct s_PM_field PM_field;
</PRE><a name="934249">
<h1>3.0 </a>Summary of the PML C API</h1>
</a>
<a name="934624">
Here is a brief summary of the routines in </a>PML. They are broken down by category as: matrix manipulation routines; other equation solvers; geometry routines; mathematical functions; constructors and destructors; field operators; and sorting routines.<p>
</a>
<a name="934281">
<h2>3.1 </a>MATRIX Manipulation Routines</h2>
</a>
<A NAME="934310"><BR><B>
</B><BR><A NAME="934289"><BR><B>PM_matrix *</a>PM_create(int nrow, int ncol)
</B><BR><a name="934587">
</a>Create and return a new PM_matrix with nrow </a>rows and ncol </a>columns. If space cannot be allocated, NULL is returned.<p>
</a>
<A NAME="934661"><BR><B>PM_matrix *</a>PM_decompose(PM_matrix *a, int *ip, int flag)
</B><BR><a name="934481">
Perform an </a>LU </a>decomposition on the PM_matrix a. The </a>pivoting information is returned in ip. The space for ip must be passed into this function. If flag is TRUE a new matrix is allocated to hold the decomposed matrix. Otherwise the original matrix is overwritten by the decomposed matrix. NULL is returned if the routine fails.<p>
</a>
<A NAME="934666"><BR><B>PM_matrix *</a>PM_decb(int n, int ml, int mu, PM_matrix *b, int *ip)
</B><BR><a name="934677">
Perform an </a>LU decomposition on the </a>banded matrix b. The arguments n, ml, and mu are the </a>order of the matrix, the </a>lower </a>band width, and the </a>upper band width respectively. The </a>pivoting information is returned in ip. The space for ip must be passed into the routine. This routine replaces the contents of b with the decomposed matrix and returns a pointer to b if successful. If the routine fails it returns NULL.<p>
</a>
<A NAME="934674"><BR><B>int </a>PM_destroy(PM_matrix *m)
</B><BR><a name="934678">
</a>Release the storage associated with the matrix m. Both the array and the matrix structure are freed. TRUE is returned if successful otherwise FALSE is returned.<p>
</a>
<A NAME="934486"><BR><B></a>PM_element(PM_matrix *m, int row, int column)
</B><BR><a name="934449">
This </a>macro accesses the </a>matrix </a>element, m(row, column). Its value can either be read or set.<p>
</a>
<A NAME="934663"><BR><B>PM_matrix *</a>PM_inverse(PM_matrix *a)
</B><BR><a name="934478">
The matrix a is inverted via an </a>LU decomposition. A new space is allocated to hold the </a>inverse. If the inversion is successful the inverse matrix is returned otherwise NULL is returned.<p>
</a>
<A NAME="934668"><BR><B>int </a>PM_inv_times_ds(PM_matrix *b, PM_matrix *c, PM_matrix *t, int *ip)
</B><BR><a name="934485">
Using an </a>LU decomposed matrix b compute (b-1).c and return it in t. The </a>pivoting information from the decomposition is supplied in ip.<p>
</a>
<A NAME="934665"><BR><B>PM_matrix *</a>PM_lower(PM_matrix *a)
</B><BR><a name="934470">
</a>Create and fill a matrix consisting of the </a>sub-diagonal part of the input matrix a. The new matrix is returned.<p>
</a>
<A NAME="934659"><BR><B>PM_matrix *</a>PM_minus(PM_matrix *a, PM_matrix *b)
</B><BR><a name="934474">
</a>Create and fill a matrix consisting of the difference, a-b. The new matrix is returned.<p>
</a>
<A NAME="934499"><BR><B>PM_matrix *</a>PM_minus_s(PM_matrix *a, PM_matrix *b, PM_matrix *c)
</B><BR><a name="934473">
Compute the </a>difference, b-c. The result is returned in the matrix a.<p>
</a>
<A NAME="934658"><BR><B>PM_matrix *</a>PM_plus(PM_matrix *a, PM_matrix *b)
</B><BR><a name="934582">
</a>Create and fill a matrix consisting of the </a>sum, a+b. The new matrix is returned.<p>
</a>
<A NAME="934670"><BR><B>PM_matrix *</a>PM_plus_s(PM_matrix *a, PM_matrix *b, PM_matrix *c)
</B><BR><a name="934472">
Compute the </a>sum, b+c. The result is returned in the matrix a.<p>
</a>
<A NAME="934667"><BR><B>PM_matrix *</a>PM_print(PM_matrix *a)
</B><BR><a name="934495">
</a>Print the matrix a to stdout. The matrix a is returned<p>
</a>
<A NAME="934662"><BR><B>PM_matrix *</a>PM_sol(PM_matrix *ul, PM_matrix *b, int *ip, int flag)
</B><BR><a name="934669">
Using the </a>LU decomposed matrix ul and the </a>pivoting information in ip as computed by PM_decompose, </a>solve for the </a>solution to the equation A.x = b. If flag is TRUE a new matrix is allocated to contain the solution and is returned. If flag is FALSE the solution is returned in b and a pointer to b is the return value of the function.<p>
</a>
<A NAME="934482"><BR><B>PM_matrix *</a>PM_solb(int n, int ml, int mu, PM_matrix *b, PM_matrix *y, int *ip)
</B><BR><a name="934671">
Using the </a>LU decomposed </a>matrix b and the </a>pivoting information in ip as computed by PM_decb, </a>solve for the </a>solution to the equation A.x = y where A is a </a>banded matrix. The arguments n, ml, and mu are the </a>order of the matrix, the </a>lower </a>band width, and the </a>upper band width respectively. The solution is returned in y and a pointer to y is the return value of the function.<p>
</a>
<A NAME="934660"><BR><B>PM_matrix *</a>PM_solve(PM_matrix *a, PM_matrix *b)
</B><BR><a name="934475">
The equation a.x = b is </a>solved and the result returned in b.<p>
</a>
<A NAME="934657"><BR><B>PM_matrix *</a>PM_times(PM_matrix *a, PM_matrix *b)
</B><BR><a name="934588">
</a>Create and fill a matrix consisting of the </a>product, a.b. The new matrix is returned.<p>
</a>
<A NAME="934496"><BR><B>PM_matrix *</a>PM_times_s(PM_matrix *a, PM_matrix *b, PM_matrix *c)
</B><BR><a name="934675">
Compute the </a>product, b.c. The result is returned in the matrix a.<p>
</a>
<A NAME="934676"><BR><B>PM_matrix *</a>PM_transpose(PM_matrix *a)
</B><BR><a name="934477">
Compute the </a>transpose of the matrix a. A new matrix is created to contain the transpose and is returned.<p>
</a>
<A NAME="934664"><BR><B>PM_matrix *</a>PM_upper(PM_matrix *a)
</B><BR><a name="934648">
</a>Create and fill a </a>matrix consisting of the </a>super-diagonal part of the input matrix a. The new matrix is returned.<p>
</a>
<a name="934656">
<h2>3.2 Other </a>Equation </a>Solvers</h2>
</a>
<a name="934354">
The </a>matrix solvers presented here are for </a>sparse matrices. Sparse matrices do not fit very naturally into the </a>matrix structure discussed above. Hence the solvers here do not use matrix structures.<p>
</a>
<A NAME="934450"><PRE>
</PRE><A NAME="934452"><BR><B>int </a>PM_block_tridi(REAL *a, REAL *b, REAL *c, REAL *u, int ns, int nb)
</B><BR><a name="934457">
</a>Solve the equation M.x = u where M is </a>block </a>tridiagonal. There are nb </a>blocks each ns squared in size. The matrices a, b, and c are the </a>sub-diagonal, </a>diagonal, and </a>super-diagonal parts respectively. They are (ns*ns)*nb long and both the a[0] and c[nb-1] blocks are unused since the </a>off-diagonal parts are one block shorter than the diagonal. The result is returned in u. All </a>vectors are passed into the routine which returns TRUE if successful and FALSE otherwise.<p>
</a>
<A NAME="934454"><PRE>
</PRE><A NAME="934585"><BR><B>REAL </a>PM_iccg(int km, int lm, double eps, int ks, int maxit, REAL *a0,
REAL *a1, REAL *b0, REAL *b1, REAL *bm1,
REAL *x, REAL *y)
</B><BR><a name="934455">
Use the </a>ICCG method to </a>solve the </a>system of equations, M.x = y. The array M is </a>symmetric and </a>positive definite. Its components are passed in as the arrays: a0, a1, b0, b1, and bm1. The known values are passed in via the y array. The input arrays are all twice as long as needed to specify the problem; the second half is workspace. That is, the input arrays should all have dimension at least 2*(km*lm).The matrix M to be solved is symmetric, and has the structure:<p>
</a>
<A NAME="934458"><PRE> a01
</PRE><A NAME="934460"><PRE> a11 a02
</PRE><A NAME="934461"><PRE> a12 a03
</PRE><A NAME="934488"><PRE> .
</PRE><A NAME="934490"><PRE> .
</PRE><A NAME="934492"><PRE> .
</PRE><A NAME="934493"><PRE> a1km-2 a0km-1
</PRE><A NAME="934494"><PRE> a1km-1 a0km
</PRE><A NAME="934586"><PRE>
</PRE><A NAME="934589"><PRE> b01 bm11 a0km+1
</PRE><A NAME="934590"><PRE>
</PRE><A NAME="934655"><PRE> b11 b02 bm12 a1km+1 a0km+2
</PRE><A NAME="934672"><PRE>
</PRE><A NAME="934673"><PRE> b12 b03 bm13 a1km+2 a0km+3
</PRE><A NAME="934301"><PRE> . .
</PRE><a name="934302">
Note that elements a1km etc. are zero due to the </a>block structure. The parameters km and lm are the dimensions of the underlying </a>computational mesh. The </a>index 0 <= k < km is the rapidly varying index of the </a>mesh array. The </a>convergence criterion is controlled by the </a>dimensionless parameter eps. When the </a>solution changes by less than eps, </a>convergence is assumed. When </a>PML runs on a </a>vector machine, the </a>incomplete </a>Cholesky </a>decomposition is done by </a>cyclic reduction. To control the </a>number of levels of cyclic reduction the ks parameter is supplied. One more level than specified will be performed. The minimum possible value of ks is 0, and the maximum possible value of ks is kp-1, where kp is the highest power of 2 with 2kp <= lm. A good choice for ks is 4 for “most” problems. On a </a>scalar machine this parameter is ignored. The </a>number of iterations in the </a>conjugate gradient step is controlled by the parameter maxit. The return value, ret, is the actual number of </a>conjugate gradient </a>iterations plus the </a>norm of the </a>residual. The value is multiplied by -1 if certain exceptional conditions arise. Failure of the algorithm then is indicated if any of the following is true: <p>
</a>
<a name="934484">
ret < 0<p>
</a>
<a name="934467">
maxit <= int(ret)<p>
</a>
<a name="934469">
eps < frac(ret)<p>
</a>
<A NAME="934891"><PRE> The x array contains the </a>solution on return.
</PRE><A NAME="934308"><PRE>
</PRE><A NAME="934307"><BR><B>int </a>PM_tridi(REAL *a, REAL *b, REAL *c, REAL *r, REAL *u, int n)
</B><BR><a name="934649">
</a>Solve the equation M.u = r where M is </a>tridiagonal. The matrices a, b, and c are the </a>sub-diagonal, </a>diagonal, and </a>super-diagonal respectively and both a[0] and c[n-1] are unused since the </a>off-diagonal v</a>ectors are one element shorter than the diagonal. The number of equations is n. The result is returned in u. All vectors are passed into the routine which returns TRUE if successful and FALSE otherwise.<p>
</a>
<a name="934621">
<h2>3.3 </a>Geometry Routines</h2>
</a>
<A NAME="934501"><PRE>
</PRE><A NAME="934465"><BR><B>REAL </a>PM_dot(REAL *x, REAL *y, int n)
</B><BR><a name="934464">
Take the </a>inner product of two </a>vectors of length n and return the result.<p>
</a>
<A NAME="934466"><PRE>
</PRE><A NAME="934505"><BR><B>int </a>PM_cross(double x1, double y1, double x2, double y2, double x3, double y3,
double x4, double y4, double *px0, double *py0)
</B><BR><a name="934652">
Return the </a>intersection point of the </a>line segment defined by (X1 - X2) and the </a>ray defined by (X3 - X4). The point X3 is the end point of the ray. Each </a>vector Xi has the components Xi = (xi, yi). The intersection point is returned in (px0, py0). If the ray does not cross the </a>segment (not the line, the segment!) FALSE is returned and the vector (</a>HUGE, HUGE) is passed back in px0 and py0. If ray does cross TRUE is returned. HUGE is a </a>#</a>define’d constant signifying a large floating point number<p>
</a>
<a name="934388">
<p>
</a>
<A NAME="934397"><BR><B>int </a>PM_cross_seg(double x1, double y1, double x2, double y2, double x3, double
y3, double x4, double y4, double *px0, double *py0)
</B><BR><a name="934398">
Lik PM_cross this function will return the intersection point of the lines implied by the two segments given. However, it returns TRUE if and only if the segments cross.<p>
</a>
<A NAME="934401"><BR><B>int </a>PM_cross_line_plane(double x1, double y1, double z1, double x2, double y2,
double z2, REAL *px, REAL *py, REAL *pz, REAL *px0, REAL
*py0, REAL *pz0, int line)
</B><BR><a name="934389">
This function tests for the intersection of lines and a plane. The value of line may be 0, 1, 2 for segment, ray, or line. The line is defined by the vectors: X1 and X2 the plane is defined by the vectors: X3, X4, X5. For rays X1 is the terminations of the ray (i.e. tail of vector) and X2 is X1 - dX (dX defines the direction of the ray). Vectors are defined as X = (x, y). TRUE is returned if and only if the segment, ray, or line crosses the plan.<p>
</a>
<A NAME="934406"><BR><B>int </a>PM_colinear_2d(REAL *px, REAL *py, int n)
</B><BR><a name="934390">
Return TRUE if and only if the n points from the px and py arrays are colinear in 2 dimensions.<p>
</a>
<A NAME="934407"><BR><B>int </a>PM_colinear_3d(REAL *px, REAL *py, REAL *pz, int n)
</B><BR><a name="934391">
Return TRUE if and only if the n points from the px, py, and pz arrays are colinear in 3 dimensions.<p>
</a>
<A NAME="934412"><BR><B>int </a>PM_contains_2d(double x, double y, REAL *px, REAL *py, int n)
</B><BR><a name="934392">
Return TRUE if and only if the point (x, y) is contained in the plane bounded by the polygon defined by the n points in the px and py arrays.<p>
</a>
<A NAME="934426"><BR><B>int </a>PM_contains_3d(double x, double y, double z, REAL *px, REAL *py, REAL
*pz, int bnd)
</B><BR><a name="934399">
Return TRUE if and only if the point (x, y, z) is contained in the plane bounded by the polygon defined by the n points in the px, py, and pz arrays. If bnd it TRUE points on the boundary are included.<p>
</a>
<A NAME="934432"><BR><B>int </a>PM_intersect_line_polygon(REAL *pxmn, REAL *pymn, REAL *pxmx,
REAL *pymx, REAL *ppx, REAL *ppy, int np, int *pic)
</B><BR><a name="934438">
This routine computes the intersection points of the line given by (pxmn, pymn) and (pxmx, pymx) and the given polygon as defined by the np points in ppx and ppy. The number of intersection points is returned in pic and the intersection points are returned in (pxmn, pymn) and (pxmx, pymx). FALSE is returned if the line segment is completely outside the polygon.<p>
</a>
<A NAME="934442"><BR><B>void </a>PM_convex_hull(REAL *px, REAL *py, int nh, REAL **ppx, REAL **ppy,
int *pnp)
</B><BR><a name="934396">
This routine allocates and returns arrays (ppx, ppy) containing pnp points defining the convex polygon which contains the nh points specified in the px and py arrays.<p>
</a>
<A NAME="934445"><BR><B>void </a>PM_nearest_point(REAL *px, REAL *py, int n, double xs, double ys, REAL
*pxt, REAL *pyt, int *pi)
</B><BR><a name="934453">
This routine finds the point in the curve defined by the n points in the px and py arrays which is nearest to the point (xs, ys). The point is returned in (pxt, pyt) and the index of that point in pi.<p>
</a>
<a name="934506">
<h2>3.4 </a>Mathematical Functions</h2>
</a>
<A NAME="934508"><PRE>
</PRE><A NAME="934510"><BR><B>double </a>PM_</a>fix(double value)
</B><BR><a name="934509">
Return the </a>integer part of value as a REAL number. Unlike the C library function floor, PM_fix(-2.3) = -2.0;<p>
</a>
<A NAME="934512"><BR><B>double </a>PM_frac(double x)
</B><BR><a name="934303">
Return the </a>fractional part of x as a REAL number.<p>
</a>
<A NAME="934304"><BR><B>double </a>PM_</a>power(double x, int p)
</B><BR><a name="934305">
Return xp.<p>
</a>
<A NAME="934517"><BR><B>double </a>PM_random(double x)
</B><BR><a name="934596">
Return a </a>random number in the range -1.0 to 1.0.<p>
</a>
<A NAME="934592"><BR><B>double </a>PM_round(double x)
</B><BR><a name="934598">
</a>Round the number x to the nearest integer returned as a REAL number.<p>
</a>
<A NAME="934594"><BR><B>double </a>PM_sgn(double value, double sign)
</B><BR><a name="934600">
Return </a>sign(sign)*abs(value).<p>
</a>
<A NAME="934513"><BR><B>double </a>PM_sign(double x)
</B><BR><a name="934533">
Return the </a>sign of x (-1.0 or 1.0) or 0 if x = 0.<p>
</a>
<A NAME="934603"><BR><B>double </a>PM_sqr(double x)
</B><BR><A NAME="934653"><PRE> Return x2.
</PRE><a name="934650">
<h2>3.5 </a>Constructors and </a>Destructors</h2>
</a>
<A NAME="934534"><PRE>
</PRE><A NAME="934525"><BR><B>PM_set *</a>PM_make_set(char *name, char *type, int cp, int nd, ...)
</B><BR><a name="934529">
</a>Allocate and return a </a>PM_set with name name and data of type type. The cp flag specifies the data given will be copied if the value is TRUE otherwise the data given will be used directly. This can be very important because some applications require that the data in the set have been dynamically allocated by SCORE and the copy would guarantee that. Also this feature gives the application more control of data flow. The parameter nd is the </a>dimension of the set. It is followed by nd integers representing the maximum value of an index for that dimension. The indices start at 0. Next comes the dimensionality of the </a>elements of the </a>set, nde. This is followed by nde arrays whose values are double precision floating point numbers. Each of these arrays must be as long as the product of the nd maximum indices.<p>
</a>
<a name="934531">
This is not remotely the most general PM_set that can be represented, but it is one of the most frequently encountered types of set.<p>
</a>
<a name="934537">
See the examples in the discussion of the </a>PM_set structure.<p>
</a>
<A NAME="934540"><PRE>
</PRE><A NAME="934571"><BR><B>void </a>PM_rel_set(PM_set *set, int rel)
</B><BR><a name="934572">
</a>Release the space associated with the specified </a>set. If rel is TRUE the data arrays will be released as well as the set structure.<p>
</a>
<a name="934306">
<p>
</a>
<A NAME="934536"><BR><B>PM_mapping *</a>PM_make_mapping(char *name, char *cat, PM_set *domain,
PM_set *range, int centering, PM_mapping *next)
</B><BR><a name="934521">
</a>Allocate and return a </a>PM_mapping with name name, type cat, domain, range, and centering. The mapping type is one of “Logical-Rectangular” or “Arbitrarily-Connected”. If a linked list of mappings is desired the mappings to which the new one points is passed in next.<p>
</a>
<a name="934522">
See the examples in the discussion of the PM_mapping structure.<p>
</a>
<A NAME="934520"><PRE>
</PRE><A NAME="934539"><BR><B>void </a>PM_rel_mapping(PM_mapping *f, int rld, int rlr)
</B><BR><a name="934679">
</a>Release the space associated with the specified mapping. If rld is TRUE the data arrays of the domain set will be released and if rlr is TRUE the data arrays of the range set will be released.<p>
</a>
<a name="934523">
<h2>3.6 </a>Field </a>Operators</h2>
</a>
<a name="934526">
The following functions come in groups. They are provided to have at hand the functions for the </a>PM_field structure which is used in connection with </a>PM_set’s and </a>PM_mapping’s. They are also useful in applications where a </a>pointer to a function is needed and the function in question is </a>addition, </a>subtraction, </a>multiplication, or </a>division. <p>
</a>
<A NAME="934524"><BR><B>int </a>PM_iplus(int x, int y)
</B><BR><A NAME="934545"><BR><B>int </a>PM_iminus(int x, int y)
</B><BR><A NAME="934546"><BR><B>int </a>PM_itimes(int x, int y)
</B><BR><A NAME="934547"><BR><B>int </a>PM_idivide(int x, int y)
</B><BR><a name="934548">
Add, subtract, multiply, or divide integer x and y values and return an integer result.<p>
</a>
<A NAME="934535"><BR><B>int </a>PM_imodulo(int x, int y)
</B><BR><a name="934549">
Return </a>mod(x, y) as an integer value.<p>
</a>
<A NAME="934551"><BR><B>long </a>PM_lplus(long x, long y)
</B><BR><A NAME="934552"><BR><B>long </a>PM_lminus(long x, long y)
</B><BR><A NAME="934553"><BR><B>long </a>PM_ltimes(long x, long y)
</B><BR><A NAME="934554"><BR><B>long </a>PM_ldivide(long x, long y)
</B><BR><a name="934604">
Add, subtract, multiply, or divide long integer x and y values and return a long integer result.<p>
</a>
<A NAME="934605"><BR><B>long </a>PM_lmodulo(long x, long y)
</B><BR><a name="934556">
Return </a>mod(x, y) as a long integer value.<p>
</a>
<A NAME="934558"><BR><B>double </a>PM_fplus(double x, double y)
</B><BR><A NAME="934560"><BR><B>double </a>PM_fminus(double x, double y)
</B><BR><A NAME="934559"><BR><B>double </a>PM_ftimes(double x, double y)
</B><BR><A NAME="934606"><BR><B>double </a>PM_fdivide(double x, double y)
</B><BR><A NAME="934680"><PRE> Add, subtract, multiply, or divide double precision floating point x and y values and
return a double result.
</PRE><a name="934527">
<h2>3.7 </a>Sorting Routines</h2>
</a>
<A NAME="934564"><PRE>
</PRE><A NAME="934566"><BR><B>int *</a>PM_t_sort(int *in, int n_dep, int n_pts, int *ord)
</B><BR><a name="934654">
This function implements the </a>topological sort routine described in Knuth’s The Art of Computer Programming, Vol 1. Input is an array of indices in, which is 2*n_dep long, containing n_dep pairs of integers specifying the topology of the set to be sorted in terms of relations, j < k (“j precedes k”). The number of points in the set is n_pts. If ord is NULL this routine will allocate space to return the ordering, otherwise it builds the ordering into ord. The array ord must be n_pts long if it is passed in. In any case, if successful, a pointer to the ordering array is returned. If there are loops in the topology specified by in, or there are other errors NULL is returned.<p>
</a>
<a name="934313">
<h1>4.0 </a>Summary of the PML FORTRAN API</h1>
</a>
<a name="934245">
Here is a brief summary of the routines in </a>PML. They are broken down by category as: matrix manipulation routines; other equation solvers; geometry routines; mathematical functions; constructors and destructors; field operators; and sorting routines.<p>
</a>
<A NAME="934324"><BR><B></a>pmszft(integer n)
</B><BR><a name="934322">
Returns the nearest integer power of 2 in the argument n. This is used to check the size of the output arrays passed into </a>pmrfft or </a>pmcfft. The actual array sizes must be one larger than this value (this is for the zero frequence component).<p>
</a>
<A NAME="934327"><BR><B></a>pmrfft(real outyr, real outyi, real outx, integer n, real inx, real iny, real xn, real xx,
integer o)
</B><BR><a name="934323">
Performs an FFT on n non-evenly spaced real data points.<p>
</a>
<a name="934326">
The input arguments are: n an integer number of data points; inx, an array of n real x values; iny, an array of n real y values; xn, the minimum x value; xx, the maximum real value; and o, an integer flag signalling the ordering of the transformed data. If o is 1 then the FFT data is returned in increasing order of “frequency”. If o is 0 then the FFT data is returned with frequency starting at 0 and increasing followed by the most negative frequency increasing to 0.<p>
</a>
<a name="934325">
The output arguments are: outyr, an array of real values constituting the real part of the transformed range; outyi, an array of real values constituting the imaginary part of the transformed range; and outx, an array of real values constituting the transformed domain (e.g. frequency). The calling routine must supply the output arrays. The size of the arrays is one greater than the value return by a call to pmszft with the size of the input arrays.<p>
</a>
<A NAME="934361"><BR><B></a>pmcfft(real outyr, real outyi, real outx, integer n, real inx, real inyr, real inyi, real
xn, real xx, integer f, integer o)
</B><BR><a name="934362">
Performs an FFT on n non-evenly spaced complex data points.<p>
</a>
<a name="934331">
The input arguments are: n an integer number of data points; inx, an array of n real x values; inyr, an array of n real y values constituting the real part of the input range; inyi, an array of n real y values constituting the imaginary part of the input range; xn, the minimum x value; xx, the maximum real value; f, an integer 1 for an FFT and -1 for in inverse FFT; and o, an integer flag signalling the ordering of the transformed data. If o is 1 then the FFT data is returned in increasing order of “frequency”. If o is 0 then the FFT data is returned with frequency starting at 0 and increasing followed by the most negative frequency increasing to 0.<p>
</a>
<a name="934367">
The output arguments are: outyr, an array of real values constituting the real part of the transformed range; outyi, an array of real values constituting the imaginary part of the transformed range; and outx, an array of real values constituting the transformed domain (e.g. frequency). The calling routine must supply the output arrays. The size of the arrays is one greater than the value return by a call to pmszft with the size of the input arrays.<p>
</a>
<A NAME="934573"><BR><B>integer </a>pmbset(integer nn, char *name, integer nt, char *type, integer cp, integer
nd, integer nde, integer maxes, integer topid, integer nextid)
</B><BR><a name="934576">
Start building a </a>PM_set. This function builds a set upto the elements of the set which are added with the </a>pmaset routine. After all the components have been added a call to </a>pmeset must be made to complete the set. An integer identifier for the set is returned. If its value is -1 then an error has occured and the set does not exist. The arguments are:<p>
</a>
<A NAME="934578">nn number of characters in the set name
<P><A NAME="934581">name the name of the set
<P><A NAME="934591">nt number of characters in the set type
<P><A NAME="934593">type the type of the elements of the set
<P><A NAME="934595">cp if 1 the elements added with pmaset will be copied
<P><A NAME="934616">nd the number of dimensions of the set
<P><A NAME="934617">nde the number of dimensions of the set elements
<P><A NAME="934618">maxes an array of dimensions for logical rectangular sets
<P><A NAME="934619">topid an identifier for a PM_mesh_topology structure
<P><A NAME="934620">nextid an identifier for the next set in a linke list
<P><a name="934328">
If the set is </a>logically rectangular then maxes should be an array containing the lengths of each dimension and topid should be 0. If the set is </a>arbitrarily connected maxes should be 0 and topid should contain a valid mesh topology structure.<p>
</a>
<a name="934645">
The return value is an integer which is used as an identifier in functions which require a reference to a set. If there is an error -1 will be returned.<p>
</a>
<A NAME="934625"><BR><B>integer </a>pmaset(integer setid, integer component, void element)
</B><BR><a name="934622">
Add </a>component number component to the set. The data in element must be an array of numbers of the type defined in the set and with length matching the number of elements of the set. Typically this will be called nde times where nde is the number of dimensions of the set elements as supplied to pmbset. If the set was made with the copy flag on (1), then the values in element will be copied into a newly allocated space, otherwise the array element will be used directly. This point should always be very carefully considered since dynamic memory heaps can be corrupted if this is used improperly. The argument setid is the value returned from a call to </a>pmbset.<p>
</a>
<A NAME="934629"><BR><B>integer </a>pmeset(integer setid)
</B><BR><a name="934628">
This completes the construction of the set begun with </a>pmbset and </a>pmaset. The argument setid is the value return from a call to </a>pmbset.<p>
</a>
<A NAME="934631"><BR><B>integer </a>pmmtop(integer nd, integer nc, integer bp, integer bnd)
</B><BR><a name="934632">
Make a </a>PM_mesh_topology structure from the given </a>connectivity informations. You should be thoroughly familiar with the discussion of connectivity given earlier. The arguments are:<p>
</a>
<A NAME="934630">nd the number of dimensions
<P><A NAME="934633">nc an array of numbers of cells of each dimension
<P><A NAME="934634">bp an array of numbers of boundary parameters at each dimension
<P><A NAME="934635">bnd an array containing the cell descriptions
<P><a name="934636">
The layout of nc has nc(1) being the number of 0-cells (nodes), nc(2) being the number of 1-cells (edges), nc(3) being the number of 2-cells (faces), and so on. The layout of bp has bp(1) is 0, bp(2) the number of parameters describing all 1-cells, bp(3) the number of parameters describing all 2-cells, and so on. The layout of bnd is as follows:<p>
</a>
<A NAME="934637">bnd(n1) thru bnd(n2-1) contains all of the 1-cell information
<P><A NAME="934639">bnd(n2) thru bnd(n3-1) contains all of the 2-cell information
<P><A NAME="934640">bnd(n3) thru bnd(n4-1) contains all of the 3-cell information
<P><A NAME="934643"> .
<P><A NAME="934641"> .
<P><A NAME="934644"> .
<P><a name="934646">
where<p>
</a>
<A NAME="934642">ni = 1 + nc(i)*bp(i)
<P><a name="934647">
It is crucial to remember that the values in bnd are indeces which are 0 based! Be sure to review the section on connectivity.<p>
</a>
<a name="934321">
<h1>5.0 </a>EXAMPLES</h1>
</a>
<a name="934729">
To provide some examples of the use of these routines, the source code for some of the </a>PML test programs is provided here. The first piece of code demonstrates the use of most of the matrix routines. This is followed by the ICCG test program and the test program for the topological sort routine. The routine to generate a </a>PM_set was included in the discussion of the </a>PM_set structure. Users may contact the author for further examples.<p>
</a>
<a name="934730">
<h2>5.1 </a>Examples of MATRIX Functions</h2>
</a>
<A NAME="934681"><PRE>
</PRE><A NAME="934731"><PRE> #include “pml.h”
</PRE><A NAME="934682"><PRE>
</PRE><A NAME="934683"><PRE> REAL a1[5][3] = {0, 0, 1, 1, 0, 1, 0, 1, 1,
</PRE><A NAME="934651"><PRE> 0.5, 0.5, 1, 0.7, 0.9, 1};
</PRE><A NAME="934684"><PRE> REAL b1[5][1] = {0.5, 1.4, .93589, 1.16795, 1.52230};
</PRE><A NAME="934688"><PRE>
</PRE><A NAME="934689"><PRE> /* MAIN - a sample program */
</PRE><A NAME="934690"><PRE>
</PRE><A NAME="934691"><PRE> main()
</PRE><A NAME="934692"><PRE> {int i, j, count;
</PRE><A NAME="934693"><PRE> int nrow = 5, ncol = 3;
</PRE><A NAME="934694"><PRE> </a>PM_matrix *m, *t, *b, *a, *c;
</PRE><A NAME="934695"><PRE>
</PRE><A NAME="934696"><PRE> m = </a>PM_create(nrow, ncol);
</PRE><A NAME="934697"><PRE> b = PM_create(nrow, 1);
</PRE><A NAME="934698"><PRE>
</PRE><A NAME="934699"><PRE> count = 0;
</PRE><A NAME="934700"><PRE> for (i = 1; i <= nrow; i++, count++)
</PRE><A NAME="934701"><PRE> {for (j = 1; j <= ncol; j++, count++)
</PRE><A NAME="934702"><PRE> {</a>PM_element(m, i, j) = count;};
</PRE><A NAME="934703"><PRE> PM_element(b, i, 1) = count;};
</PRE><A NAME="934704"><PRE> </a>PM_print(m);
</PRE><A NAME="934705"><PRE>
</PRE><A NAME="934386"><PRE> m->array = (REAL *)a1;
</PRE><A NAME="934706"><PRE> b->array = (REAL *)b1;
</PRE><A NAME="934707"><PRE>
</PRE><A NAME="934708"><PRE> PM_print(m);
</PRE><A NAME="934709"><PRE>
</PRE><A NAME="934710"><PRE> PM_print(b);
</PRE><A NAME="934711"><PRE>
</PRE><A NAME="934712"><PRE> PM_print(t = </a>PM_transpose(m));
</PRE><A NAME="934713"><PRE>
</PRE><A NAME="934714"><PRE> PM_print(a = </a>PM_times(t, m));
</PRE><A NAME="934715"><PRE>
</PRE><A NAME="934716"><PRE> PM_print(c = PM_times(t, b));
</PRE><A NAME="934717"><PRE>
</PRE><A NAME="934718"><PRE> PM_print(</a>PM_solve(a, c));
</PRE><A NAME="934719"><PRE>
</PRE><A NAME="934720"><PRE> PM_print(t = PM_times(PM_transpose(m), m));
</PRE><A NAME="934721"><PRE>
</PRE><A NAME="934722"><PRE> PM_print(c = </a>PM_inverse(t));
</PRE><A NAME="934723"><PRE>
</PRE><A NAME="934724"><PRE> PM_print(PM_times(c, t));
</PRE><A NAME="934685"><PRE>
</PRE><A NAME="934686"><PRE> exit(0);}
</PRE><A NAME="934725"><PRE>
</PRE><a name="934727">
<h2>5.2 </a>Example of </a>ICCG Routine</h2>
</a>
<A NAME="934733"><PRE>
</PRE><A NAME="934734"><PRE> #include “pml.h”
</PRE><A NAME="934735"><PRE>
</PRE><A NAME="934736"><PRE> #define KM 4
</PRE><A NAME="934737"><PRE> #define LM 3
</PRE><A NAME="934738"><PRE> #define KXL 24 /* 2*KM*LM */
</PRE><A NAME="934739"><PRE> #define KL 12 /* KM*LM */
</PRE><A NAME="934740"><PRE> #define EPS 1.0e-6
</PRE><A NAME="934741"><PRE> #define KS 4
</PRE><A NAME="934742"><PRE> #define MAXIT 100
</PRE><A NAME="934743"><PRE>
</PRE><A NAME="934746"><PRE> /* MAIN - solve a Laplace equation */
</PRE><A NAME="934747"><PRE>
</PRE><A NAME="934748"><PRE> main()
</PRE><A NAME="934749"><PRE> {int i, j, k;
</PRE><A NAME="934750"><PRE> REAL ret;
</PRE><A NAME="934751"><PRE> REAL a0[KXL], a1[KXL], b0[KXL], b1[KXL], bm1[KXL],
</PRE><A NAME="934687"><PRE> x[KXL], y[KXL];
</PRE><A NAME="934752"><PRE>
</PRE><A NAME="934753"><PRE> /* solve a Laplacian in the l direction */
</PRE><A NAME="934745"><PRE> for (i = 0; i < KXL; i++)
</PRE><A NAME="934754"><PRE> {a0[i] = 0.0;
</PRE><A NAME="934755"><PRE> a1[i] = 0.0;
</PRE><A NAME="934756"><PRE> b0[i] = 0.0;
</PRE><A NAME="934757"><PRE> b1[i] = 0.0;
</PRE><A NAME="934758"><PRE> bm1[i] = 0.0;
</PRE><A NAME="934759"><PRE> x[i] = 0.0;
</PRE><A NAME="934760"><PRE> y[i] = 0.0;};
</PRE><A NAME="934761"><PRE>
</PRE><A NAME="934762"><PRE> for (i = KM; i < 2*KM; i++)
</PRE><A NAME="934763"><PRE> {a0[i] = 2.0;
</PRE><A NAME="934764"><PRE> a1[i] = -1.0;
</PRE><A NAME="934765"><PRE> b0[i] = 0.0;
</PRE><A NAME="934766"><PRE> b1[i] = 0.0;
</PRE><A NAME="934767"><PRE> bm1[i] = 0.0;};
</PRE><A NAME="934768"><PRE>
</PRE><A NAME="934769"><PRE> a1[2*KM-1] = 0.0;
</PRE><A NAME="934770"><PRE> for (i = 0; i < KM; i++)
</PRE><A NAME="934771"><PRE> {k = KM + i;
</PRE><A NAME="934772"><PRE> y[k] = (REAL) (i+1);};
</PRE><A NAME="934773"><PRE>
</PRE><A NAME="934774"><PRE> printf(“\nProblem y :\n”);
</PRE><A NAME="934775"><PRE> for (j = 0; j < LM; j++)
</PRE><A NAME="934776"><PRE> {printf(“\nRow #%d: \n”, j+1);
</PRE><A NAME="934777"><PRE> for (i = 0; i < KM; i++)
</PRE><A NAME="934778"><PRE> {k = j*KM + i;
</PRE><A NAME="934779"><PRE> printf(“ y(%2d, %2d) = %g “, j+1, i+1, y[k]);};};
</PRE><A NAME="934780"><PRE> printf(“\n”);
</PRE><A NAME="934781"><PRE>
</PRE><A NAME="934782"><PRE> ret = </a>PM_iccg(KM, LM, EPS, KS, MAXIT,
</PRE><A NAME="934726"><PRE> a0, a1, b0, b1, bm1, x, y);
</PRE><A NAME="934783"><PRE>
</PRE><A NAME="934784"><PRE> printf(“\nSolution x (should be (4, 7, 8, 6)) : %g\n”, ret);
</PRE><A NAME="934785"><PRE> for (j = 0; j < LM; j++)
</PRE><A NAME="934786"><PRE> {printf(“\nRow #%d: \n”, j+1);
</PRE><A NAME="934787"><PRE> for (i = 0; i < KM; i++)
</PRE><A NAME="934788"><PRE> {k = j*KM + i;
</PRE><A NAME="934728"><PRE> printf(“ x(%2d, %2d) = %g “,j+1, i+1, x[k]);};};
</PRE><A NAME="934789"><PRE> printf(“\n\n”);
</PRE><A NAME="934790"><PRE>
</PRE><A NAME="934791"><PRE> SC_pause();
</PRE><A NAME="934792"><PRE>
</PRE><A NAME="934793"><PRE> /* solve a Laplacian in the k direction */
</PRE><A NAME="934744"><PRE> for (i = 0; i < KXL; i++)
</PRE><A NAME="934794"><PRE> {a0[i] = 0.0;
</PRE><A NAME="934795"><PRE> a1[i] = 0.0;
</PRE><A NAME="934796"><PRE> b0[i] = 0.0;
</PRE><A NAME="934797"><PRE> b1[i] = 0.0;
</PRE><A NAME="934798"><PRE> bm1[i] = 0.0;
</PRE><A NAME="934799"><PRE> x[i] = 0.0;
</PRE><A NAME="934800"><PRE> y[i] = 0.0;};
</PRE><A NAME="934801"><PRE>
</PRE><A NAME="934802"><PRE> for (i = 1; i < KXL; i += LM)
</PRE><A NAME="934803"><PRE> {a0[i] = 2.0;
</PRE><A NAME="934804"><PRE> a1[i] = 0.0;
</PRE><A NAME="934805"><PRE> b0[i] = -1.0;
</PRE><A NAME="934806"><PRE> b1[i] = 0.0;
</PRE><A NAME="934807"><PRE> bm1[i] = 0.0;};
</PRE><A NAME="934808"><PRE>
</PRE><A NAME="934809"><PRE> for (k = 1, i = 1; i < KXL; i += LM, k++)
</PRE><A NAME="934810"><PRE> y[i] = (REAL) k;
</PRE><A NAME="934811"><PRE>
</PRE><A NAME="934812"><PRE> printf(“\nProblem y :\n”);
</PRE><A NAME="934813"><PRE> for (j = 0; j < KM; j++)
</PRE><A NAME="934814"><PRE> {printf(“\nRow #%d: \n”, j+1);
</PRE><A NAME="934815"><PRE> for (i = 0; i < LM; i++)
</PRE><A NAME="934816"><PRE> {k = j*LM + i;
</PRE><A NAME="934817"><PRE> printf(“ y(%2d, %2d) = %g “, j+1, i+1, y[k]);};};
</PRE><A NAME="934818"><PRE> printf(“\n”);
</PRE><A NAME="934819"><PRE>
</PRE><A NAME="934820"><PRE> ret = </a>PM_iccg(LM, KM, EPS, KS, MAXIT,
</PRE><A NAME="934732"><PRE> a0, a1, b0, b1, bm1, x, y);
</PRE><A NAME="934821"><PRE>
</PRE><A NAME="934822"><PRE> printf(“\nSolution x (should be (4, 7, 8, 6)) : %g\n”, ret);
</PRE><A NAME="934823"><PRE> for (j = 0; j < KM; j++)
</PRE><A NAME="934824"><PRE> {printf(“\nRow #%d: \n”, j+1);
</PRE><A NAME="934825"><PRE> for (i = 0; i < LM; i++)
</PRE><A NAME="934826"><PRE> {k = j*LM + i;
</PRE><A NAME="934827"><PRE> printf(“ x(%2d, %2d) = %g “, j+1, i+1, x[k]);};};
</PRE><A NAME="934828"><PRE> printf(“\n\n”);
</PRE><A NAME="934829"><PRE>
</PRE><A NAME="934831"><PRE> exit(0);}
</PRE><a name="934832">
<h2>5.3 </a>Example of </a>Topological </a>Sort Routine</h2>
</a>
<A NAME="934830"><PRE>
</PRE><A NAME="934835"><PRE> #include “pml.h”
</PRE><A NAME="934836"><PRE>
</PRE><A NAME="934837"><PRE> #define NDEP 10
</PRE><A NAME="934838"><PRE> #define NPTS 9
</PRE><A NAME="934839"><PRE>
</PRE><A NAME="934842"><PRE> main()
</PRE><A NAME="934843"><PRE> {int i, rel[20], *ord;
</PRE><A NAME="934844"><PRE>
</PRE><A NAME="934845"><PRE> rel[0] = 9;
</PRE><A NAME="934846"><PRE> rel[1] = 2;
</PRE><A NAME="934847"><PRE>
</PRE><A NAME="934848"><PRE> rel[2] = 3;
</PRE><A NAME="934849"><PRE> rel[3] = 7;
</PRE><A NAME="934850"><PRE>
</PRE><A NAME="934851"><PRE> rel[4] = 7;
</PRE><A NAME="934852"><PRE> rel[5] = 5;
</PRE><A NAME="934853"><PRE>
</PRE><A NAME="934854"><PRE> rel[6] = 5;
</PRE><A NAME="934855"><PRE> rel[7] = 8;
</PRE><A NAME="934856"><PRE>
</PRE><A NAME="934857"><PRE> rel[8] = 8;
</PRE><A NAME="934858"><PRE> rel[9] = 6;
</PRE><A NAME="934859"><PRE>
</PRE><A NAME="934860"><PRE> rel[10] = 4;
</PRE><A NAME="934861"><PRE> rel[11] = 6;
</PRE><A NAME="934862"><PRE>
</PRE><A NAME="934863"><PRE> rel[12] = 1;
</PRE><A NAME="934864"><PRE> rel[13] = 3;
</PRE><A NAME="934865"><PRE>
</PRE><A NAME="934866"><PRE> rel[14] = 7;
</PRE><A NAME="934867"><PRE> rel[15] = 4;
</PRE><A NAME="934868"><PRE>
</PRE><A NAME="934869"><PRE> rel[16] = 9;
</PRE><A NAME="934870"><PRE> rel[17] = 5;
</PRE><A NAME="934871"><PRE>
</PRE><A NAME="934872"><PRE> rel[18] = 2;
</PRE><A NAME="934873"><PRE> rel[19] = 8;
</PRE><A NAME="934874"><PRE>
</PRE><A NAME="934875"><PRE> printf(“\nRelations List:\n”);
</PRE><A NAME="934876"><PRE> for (i = 0; i < 2*NDEP; i += 2)
</PRE><A NAME="934877"><PRE> printf(“%2d < %2d\n”, rel[i], rel[i+1]);
</PRE><A NAME="934878"><PRE> printf(“\n”);
</PRE><A NAME="934879"><PRE>
</PRE><A NAME="934880"><PRE> ord = </a>PM_t_sort(rel, NDEP, NPTS, NULL);
</PRE><A NAME="934881"><PRE>
</PRE><A NAME="934882"><PRE> printf(“Partially ordered List:\n”);
</PRE><A NAME="934883"><PRE> for (i = 0; i < NPTS; i++)
</PRE><A NAME="934884"><PRE> printf(“A(%2d) = %2d\n”, i, ord[i]);
</PRE><A NAME="934885"><PRE>
</PRE><A NAME="934886"><PRE> printf(“\nThe correct order is:\n 1 9 3 2 7 4 5 8 6\n\n”);
</PRE><A NAME="934887"><PRE>
</PRE><A NAME="934231"><PRE> exit(0);}
</PRE><a name="934311">
<h1>6.0 </a>Related </a>Documentation</h1>
</a>
<a name="934275">
</a>PML is a part of the </a>PACT </a>portable </a>code development and </a>visualization tool set. Interested readers may with to consult other PACT documents.<p>
</a>
<a name="934312">
<p>
</a>
<a name="934338">
The list of PACT Documents is:<p>
</a>
<A NAME="934348"><PRE> PACT User’s Guide, UCRL-MA-112087
</PRE><A NAME="934351"><PRE> SCORE User’s Manual, UCRL-MA-108976 Rev.1
</PRE><A NAME="934360"><PRE> PPC User’s Manual UCRL-MA-108964 Rev.1
</PRE><A NAME="934365"><PRE> PML User’s Manual, UCRL-MA-108965 Rev.1 (this document)
</PRE><A NAME="934369"><PRE> PDBLib User’s Manual, M-270 Rev.2
</PRE><A NAME="934372"><PRE> PGS User’s Manual, UCRL-MA-108966 Rev.1
</PRE><A NAME="934373"><PRE> PANACEA User’s Manual, M-276 Rev.2
</PRE><A NAME="934376"><PRE> ULTRA II User’s Manual, UCRL-MA-108967 Rev.1
</PRE><A NAME="934378"><PRE> PDBDiff User’s Manual, UCRL-MA-108975 Rev.1
</PRE><A NAME="934381"><PRE> PDBView User’s Manual, UCRL-MA-108968 Rev.1
</PRE><A NAME="934340"><PRE> SX User’s Manual, UCRL-MA-112315
</PRE><a name="934228">
<p>
</a>
<a name="934363">
Additional reading:<p>
</a>
<A NAME="934364"><PRE> Abramowitz, Stegun, Handbook of Mathematical Functions, Dover.
</PRE><A NAME="934387"><PRE> Knuth, D.E. The Art of Computer Programming, Vol I - III, Addison-Wesley.
</PRE><A NAME="934383"><PRE> Press, Flannery, Teukolsky, and Vetterling, Numerical Recipes in C, Cambridge Press.
</PRE>
<p><hr>
</body></html>
|
