1 /*  testOptPart.c  */
2 
3 #include "../../SymbFac.h"
4 #include "../../timings.h"
5 
6 /*--------------------------------------------------------------------*/
7 int
main(int argc,char * argv[])8 main ( int argc, char *argv[] )
9 /*
10    ------------------------------------------------------
11    (1) read in an ETree object.
12    (2) read in an Graph object.
13    (3) find the optimal domain/schur complement partition
14        for a semi-implicit factorization
15 
16    created -- 96oct03, cca
17    ------------------------------------------------------
18 */
19 {
20 char     *inETreeFileName, *inGraphFileName, *outIVfileName ;
21 double   alpha, nA21, nfent1, nfops1, nL11, nL22, nPhi, nV, t1, t2 ;
22 Graph    *graph ;
23 int      ii, inside, J, K, msglvl, nfind1, nfront, nJ, nleaves1,
24          nnode1, nvtx, rc, sizeJ, totalgain, vsize, v, w ;
25 int      *adjJ, *compids, *nodwghts, *vadj, *vtxToFront, *vwghts ;
26 IV       *compidsIV ;
27 IVL      *symbfacIVL ;
28 ETree    *etree ;
29 FILE     *msgFile ;
30 Tree     *tree ;
31 
32 if ( argc != 7 ) {
33    fprintf(stdout,
34 "\n\n usage : %s msglvl msgFile inETreeFile inGraphFile alpha"
35 "\n         outIVfile "
36 "\n    msglvl       -- message level"
37 "\n    msgFile      -- message file"
38 "\n    inETreeFile  -- input file, must be *.etreef or *.etreeb"
39 "\n    inGraphFile  -- input file, must be *.graphf or *.graphb"
40 "\n    alpha        -- weight parameter"
41 "\n       alpha = 0 --> minimize storage"
42 "\n       alpha = 1 --> minimize solve ops"
43 "\n    outIVfile    -- output file for oldToNew vector,"
44 "\n                    must be *.ivf or *.ivb"
45 "\n", argv[0]) ;
46    return(0) ;
47 }
48 msglvl = atoi(argv[1]) ;
49 if ( strcmp(argv[2], "stdout") == 0 ) {
50    msgFile = stdout ;
51 } else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
52    fprintf(stderr, "\n fatal error in %s"
53            "\n unable to open file %s\n",
54            argv[0], argv[2]) ;
55    return(-1) ;
56 }
57 inETreeFileName  = argv[3] ;
58 inGraphFileName  = argv[4] ;
59 alpha            = atof(argv[5]) ;
60 outIVfileName    = argv[6] ;
61 fprintf(msgFile,
62         "\n %s "
63         "\n msglvl        -- %d"
64         "\n msgFile       -- %s"
65         "\n inETreeFile   -- %s"
66         "\n inGraphFile   -- %s"
67         "\n alpha         -- %f"
68         "\n outIVfile     -- %s"
69         "\n",
70         argv[0], msglvl, argv[2],
71         inETreeFileName, inGraphFileName, alpha, outIVfileName) ;
72 fflush(msgFile) ;
73 /*
74    ------------------------
75    read in the ETree object
76    ------------------------
77 */
78 if ( strcmp(inETreeFileName, "none") == 0 ) {
79    fprintf(msgFile, "\n no file to read from") ;
80    exit(0) ;
81 }
82 etree = ETree_new() ;
83 MARKTIME(t1) ;
84 rc = ETree_readFromFile(etree, inETreeFileName) ;
85 MARKTIME(t2) ;
86 fprintf(msgFile, "\n CPU %9.5f : read in etree from file %s",
87         t2 - t1, inETreeFileName) ;
88 if ( rc != 1 ) {
89    fprintf(msgFile, "\n return value %d from ETree_readFromFile(%p,%s)",
90            rc, etree, inETreeFileName) ;
91    exit(-1) ;
92 }
93 ETree_leftJustify(etree) ;
94 fprintf(msgFile, "\n\n after reading ETree object from file %s",
95         inETreeFileName) ;
96 if ( msglvl > 2 ) {
97    ETree_writeForHumanEye(etree, msgFile) ;
98 } else {
99    ETree_writeStats(etree, msgFile) ;
100 }
101 fflush(msgFile) ;
102 nfront     = ETree_nfront(etree) ;
103 tree       = ETree_tree(etree) ;
104 nodwghts   = ETree_nodwghts(etree) ;
105 vtxToFront = ETree_vtxToFront(etree) ;
106 /*
107    ------------------------
108    read in the Graph object
109    ------------------------
110 */
111 if ( strcmp(inGraphFileName, "none") == 0 ) {
112    fprintf(msgFile, "\n no file to read from") ;
113    exit(0) ;
114 }
115 graph = Graph_new() ;
116 MARKTIME(t1) ;
117 rc = Graph_readFromFile(graph, inGraphFileName) ;
118 nvtx = graph->nvtx ;
119 vwghts = graph->vwghts ;
120 MARKTIME(t2) ;
121 fprintf(msgFile, "\n CPU %9.5f : read in graph from file %s",
122         t2 - t1, inGraphFileName) ;
123 if ( rc != 1 ) {
124    fprintf(msgFile, "\n return value %d from Graph_readFromFile(%p,%s)",
125            rc, graph, inGraphFileName) ;
126    exit(-1) ;
127 }
128 fprintf(msgFile, "\n\n after reading Graph object from file %s",
129         inGraphFileName) ;
130 if ( msglvl > 2 ) {
131    Graph_writeForHumanEye(graph, msgFile) ;
132 } else {
133    Graph_writeStats(graph, msgFile) ;
134 }
135 fflush(msgFile) ;
136 /*
137    ----------------------
138    compute the statistics
139    ----------------------
140 */
141 nnode1 = etree->tree->n ;
142 nfind1 = ETree_nFactorIndices(etree) ;
143 nfent1 = ETree_nFactorEntries(etree, SPOOLES_SYMMETRIC) ;
144 nfops1 = ETree_nFactorOps(etree, SPOOLES_REAL, SPOOLES_SYMMETRIC) ;
145 nleaves1 = Tree_nleaves(etree->tree) ;
146 fprintf(stdout, "\n root front %d has %d vertices",
147         etree->tree->root,
148         etree->nodwghtsIV->vec[etree->tree->root]) ;
149 /*
150    ---------------------------------
151    create the symbolic factorization
152    ---------------------------------
153 */
154 symbfacIVL = SymbFac_initFromGraph(etree, graph) ;
155 if ( msglvl > 2 ) {
156    IVL_writeForHumanEye(symbfacIVL, msgFile) ;
157 } else {
158    IVL_writeStats(symbfacIVL, msgFile) ;
159 }
160 fflush(msgFile) ;
161 /*
162    --------------------------
163    find the optimal partition
164    --------------------------
165 */
166 compidsIV = ETree_optPart(etree, graph, symbfacIVL, alpha,
167                           &totalgain, msglvl, msgFile) ;
168 if ( msglvl > 2 ) {
169    IV_writeForHumanEye(compidsIV, msgFile) ;
170 } else {
171    IV_writeStats(compidsIV, msgFile) ;
172 }
173 fflush(msgFile) ;
174 compids = IV_entries(compidsIV) ;
175 /*
176    ------------------------------------------------------
177    compute the number of vertices in the schur complement
178    ------------------------------------------------------
179 */
180 for ( J = 0, nPhi = nV = 0. ; J < nfront ; J++ ) {
181    if ( compids[J] == 0 ) {
182       nPhi += nodwghts[J] ;
183    }
184    nV += nodwghts[J] ;
185 }
186 /*
187    --------------------------------------------
188    compute the number of entries in L11 and L22
189    --------------------------------------------
190 */
191 nL11 = nL22 = 0 ;
192 for ( J = Tree_postOTfirst(tree) ;
193       J != -1 ;
194       J = Tree_postOTnext(tree, J) ) {
195    nJ = nodwghts[J] ;
196    if ( msglvl > 3 ) {
197       fprintf(msgFile, "\n\n front %d, nJ = %d", J, nJ) ;
198    }
199    IVL_listAndSize(symbfacIVL, J, &sizeJ, &adjJ) ;
200    for ( ii = 0, inside = 0 ; ii < sizeJ ; ii++ ) {
201       w = adjJ[ii] ;
202       K = vtxToFront[w] ;
203       if ( msglvl > 3 ) {
204          fprintf(msgFile, "\n    w = %d, K = %d", w, K) ;
205       }
206       if ( K > J && compids[K] == compids[J] ) {
207          inside += (vwghts == NULL) ? 1 : vwghts[w] ;
208          if ( msglvl > 3 ) {
209             fprintf(msgFile, ", inside") ;
210          }
211       }
212    }
213    if ( compids[J] != 0 ) {
214       if ( msglvl > 3 ) {
215          fprintf(msgFile, "\n    inside = %d, adding %d to L11",
216                  inside, nJ*nJ + 2*nJ*inside) ;
217       }
218       nL11 += (nJ*(nJ+1))/2 + nJ*inside ;
219    } else {
220       if ( msglvl > 3 ) {
221          fprintf(msgFile, "\n    inside = %d, adding %d to L22",
222                  inside, (nJ*(nJ+1))/2 + nJ*inside) ;
223       }
224       nL22 += (nJ*(nJ+1))/2 + nJ*inside ;
225    }
226 }
227 if ( msglvl > 0 ) {
228    fprintf(msgFile, "\n |L| = %.0f, |L11| = %.0f, |L22| = %.0f",
229            nfent1, nL11, nL22) ;
230 }
231 /*
232    ------------------------------------
233    compute the number of entries in A21
234    ------------------------------------
235 */
236 nA21 = 0 ;
237 if ( vwghts != NULL ) {
238    for ( v = 0 ; v < nvtx ; v++ ) {
239       J = vtxToFront[v] ;
240       if ( compids[J] != 0 ) {
241          Graph_adjAndSize(graph, v, &vsize, &vadj) ;
242          for ( ii = 0 ; ii < vsize ; ii++ ) {
243             w = vadj[ii] ;
244             K = vtxToFront[w] ;
245             if ( compids[K] == 0 ) {
246                if ( msglvl > 3 ) {
247                   fprintf(msgFile, "\n A21 : v = %d, w = %d", v, w) ;
248                }
249                nA21 += vwghts[v] * vwghts[w] ;
250             }
251          }
252       }
253    }
254 } else {
255    for ( v = 0 ; v < nvtx ; v++ ) {
256       J = vtxToFront[v] ;
257       if ( compids[J] != 0 ) {
258          Graph_adjAndSize(graph, v, &vsize, &vadj) ;
259          for ( ii = 0 ; ii < vsize ; ii++ ) {
260             w = vadj[ii] ;
261             K = vtxToFront[w] ;
262             if ( compids[K] == 0 ) {
263                if ( msglvl > 3 ) {
264                   fprintf(msgFile, "\n A21 : v = %d, w = %d", v, w) ;
265                }
266                nA21++ ;
267             }
268          }
269       }
270    }
271 }
272 if ( msglvl > 0 ) {
273    fprintf(msgFile,
274            "\n |L| = %.0f, |L11| = %.0f, |L22| = %.0f, |A21| = %.0f",
275            nfent1, nL11, nL22, nA21) ;
276    fprintf(msgFile,
277       "\n storage: explicit = %.0f, semi-implicit = %.0f, ratio = %.3f"
278       "\n opcount: explicit = %.0f, semi-implicit = %.0f, ratio = %.3f",
279       nfent1, nL11 + nA21 + nL22,
280       nfent1/(nL11 + nA21 + nL22),
281       2*nfent1, 4*nL11 + 2*nA21 + 2*nL22,
282       2*nfent1/(4*nL11 + 2*nA21 + 2*nL22)) ;
283    fprintf(msgFile, "\n ratios %8.3f %8.3f %8.3f",
284            nPhi/nV,
285            nfent1/(nL11 + nA21 + nL22),
286            2*nfent1/(4*nL11 + 2*nA21 + 2*nL22)) ;
287 }
288 /*
289    ----------------
290    free the objects
291    ----------------
292 */
293 ETree_free(etree) ;
294 Graph_free(graph) ;
295 IVL_free(symbfacIVL) ;
296 
297 fprintf(msgFile, "\n") ;
298 fclose(msgFile) ;
299 
300 return(1) ; }
301 
302 /*--------------------------------------------------------------------*/
303