1 /* ms.c */
2
3 #include "../ETree.h"
4
5 #define MYDEBUG 0
6
7 /*--------------------------------------------------------------------*/
8 /*
9 ------------------------------------------------
10 returns a compidsIV IV object that maps the
11 vertices to a domain (compids[v] > 1)
12 or to the multisector (compids[v] = 0).
13 the vertices in the multisector is specified
14 by their depth of their front in the tree.
15
16 created -- 96jan04, cca
17 ------------------------------------------------
18 */
19 IV *
ETree_msByDepth(ETree * etree,int depth)20 ETree_msByDepth (
21 ETree *etree,
22 int depth
23 ) {
24 int front, nfront, nvtx, v ;
25 int *compids, *dmetric, *vtxToFront ;
26 IV *compidsIV, *vmetricIV, *dmetricIV ;
27 /*
28 ---------------
29 check the input
30 ---------------
31 */
32 if ( etree == NULL
33 || (nfront = etree->nfront) <= 0
34 || (nvtx = etree->nvtx) <= 0
35 || depth <= 0 ) {
36 fprintf(stderr, "\n fatal error in ETree_msByDepth(%p,%d)"
37 "\n bad input\n", etree, depth) ;
38 exit(-1) ;
39 }
40 vtxToFront = IV_entries(etree->vtxToFrontIV) ;
41 /*
42 --------------------
43 get the depth metric
44 --------------------
45 */
46 vmetricIV = IV_new() ;
47 IV_init(vmetricIV, nfront, NULL) ;
48 IV_fill(vmetricIV, 1) ;
49 dmetricIV = Tree_setDepthImetric(etree->tree, vmetricIV) ;
50 dmetric = IV_entries(dmetricIV) ;
51 #if MYDEBUG > 0
52 { int ierr ;
53 fprintf(stdout, "\n ETree_msByDepth") ;
54 fprintf(stdout, "\n vmetric") ;
55 IV_writeForHumanEye(vmetricIV, stdout) ;
56 fprintf(stdout, "\n dmetric") ;
57 IV_writeForHumanEye(dmetricIV, stdout) ;
58 }
59 #endif
60 IV_free(vmetricIV) ;
61 /*
62 --------------------------------------------------
63 fill the IV object with the ids in the multisector
64 --------------------------------------------------
65 */
66 compidsIV = IV_new() ;
67 IV_init(compidsIV, nvtx, NULL) ;
68 compids = IV_entries(compidsIV) ;
69 for ( v = 0 ; v < nvtx ; v++ ) {
70 front = vtxToFront[v] ;
71 if ( dmetric[front] <= depth ) {
72 compids[v] = 0 ;
73 } else {
74 compids[v] = 1 ;
75 }
76 }
77 IV_free(dmetricIV) ;
78
79 return(compidsIV) ; }
80
81 /*--------------------------------------------------------------------*/
82 /*
83 ----------------------------------------------------------------
84 construct a multisector based on vertices found in a subtree.
85
86 created -- 96jan04, cca
87 ----------------------------------------------------------------
88 */
89 IV *
ETree_msByNvtxCutoff(ETree * etree,double cutoff)90 ETree_msByNvtxCutoff (
91 ETree *etree,
92 double cutoff
93 ) {
94 int front, nfront, nvtx, v ;
95 int *compids, *tmetric, *vtxToFront ;
96 IV *compidsIV, *tmetricIV, *vmetricIV ;
97 /*
98 ---------------
99 check the input
100 ---------------
101 */
102 if ( etree == NULL
103 || (nfront = etree->nfront) <= 0
104 || (nvtx = etree->nvtx) <= 0 ) {
105 fprintf(stderr, "\n fatal error in ETree_msByCutoff(%p,%f)"
106 "\n bad input\n", etree, cutoff) ;
107 exit(-1) ;
108 }
109 vtxToFront = IV_entries(etree->vtxToFrontIV) ;
110 /*
111 ----------------------
112 get the subtree metric
113 ----------------------
114 */
115 vmetricIV = ETree_nvtxMetric(etree) ;
116 tmetricIV = Tree_setSubtreeImetric(etree->tree, vmetricIV) ;
117 #if MYDEBUG > 0
118 fprintf(stdout, "\n ETree_msByNvtxCutoff") ;
119 fprintf(stdout, "\n vmetric") ;
120 IV_writeForHumanEye(vmetricIV, stdout) ;
121 fprintf(stdout, "\n tmetric") ;
122 IV_writeForHumanEye(tmetricIV, stdout) ;
123 fflush(stdout) ;
124 #endif
125 IV_free(vmetricIV) ;
126 cutoff = cutoff * IV_max(tmetricIV) ;
127 /*
128 --------------------------------------------------
129 fill the IV object with the ids in the multisector
130 --------------------------------------------------
131 */
132 compidsIV = IV_new() ;
133 IV_init(compidsIV, nvtx, NULL) ;
134 compids = IV_entries(compidsIV) ;
135 tmetric = IV_entries(tmetricIV) ;
136 for ( v = 0 ; v < nvtx ; v++ ) {
137 front = vtxToFront[v] ;
138 if ( tmetric[front] >= cutoff ) {
139 compids[v] = 0 ;
140 } else {
141 compids[v] = 1 ;
142 }
143 }
144 IV_free(tmetricIV) ;
145
146 return(compidsIV) ; }
147
148 /*--------------------------------------------------------------------*/
149 /*
150 --------------------------------------------------
151 construct a multisector based on the number
152 of factor entries found in a subtree.
153
154 symflag -- symmetry flag, one of SPOOLES_SYMMETRIC
155 SPOOLES_HERMITIAN or SPOOLES_NONSYMMETRIC
156
157 created -- 96jan04, cca
158 --------------------------------------------------
159 */
160 IV *
ETree_msByNentCutoff(ETree * etree,double cutoff,int symflag)161 ETree_msByNentCutoff (
162 ETree *etree,
163 double cutoff,
164 int symflag
165 ) {
166 int front, nfront, nvtx, v ;
167 int *compids, *tmetric, *vtxToFront ;
168 IV *compidsIV, *tmetricIV, *vmetricIV ;
169 /*
170 ---------------
171 check the input
172 ---------------
173 */
174 if ( etree == NULL
175 || (nfront = etree->nfront) <= 0
176 || (nvtx = etree->nvtx) <= 0 ) {
177 fprintf(stderr, "\n fatal error in ETree_msByCutoff(%p,%f,%d)"
178 "\n bad input\n", etree, cutoff, symflag) ;
179 exit(-1) ;
180 }
181 vtxToFront = IV_entries(etree->vtxToFrontIV) ;
182 /*
183 ----------------------
184 get the subtree metric
185 ----------------------
186 */
187 vmetricIV = ETree_nentMetric(etree, symflag) ;
188 tmetricIV = Tree_setSubtreeImetric(etree->tree, vmetricIV) ;
189 #if MYDEBUG > 0
190 fprintf(stdout, "\n ETree_msByNentCutoff") ;
191 fprintf(stdout, "\n vmetric") ;
192 IV_writeForHumanEye(vmetricIV, stdout) ;
193 fprintf(stdout, "\n tmetric") ;
194 IV_writeForHumanEye(tmetricIV, stdout) ;
195 fflush(stdout) ;
196 #endif
197 IV_free(vmetricIV) ;
198 cutoff = cutoff * IV_max(tmetricIV) ;
199 /*
200 --------------------------------------------------
201 fill the IV object with the ids in the multisector
202 --------------------------------------------------
203 */
204 compidsIV = IV_new() ;
205 IV_init(compidsIV, nvtx, NULL) ;
206 compids = IV_entries(compidsIV) ;
207 tmetric = IV_entries(tmetricIV) ;
208 for ( v = 0 ; v < nvtx ; v++ ) {
209 front = vtxToFront[v] ;
210 if ( tmetric[front] >= cutoff ) {
211 compids[v] = 0 ;
212 } else {
213 compids[v] = 1 ;
214 }
215 }
216 IV_free(tmetricIV) ;
217
218 return(compidsIV) ; }
219
220 /*--------------------------------------------------------------------*/
221 /*
222 --------------------------------------------------
223 construct a multisector based on the number
224 of factor operations found in a subtree.
225
226 type -- type of entries,
227 SPOOLES_REAL or SPOOLES_COMPLEX
228
229 symflag -- symmetry flag, one of SPOOLES_SYMMETRIC
230 SPOOLES_HERMITIAN or SPOOLES_NONSYMMETRIC
231
232 created -- 96jan04, cca
233 --------------------------------------------------
234 */
235 IV *
ETree_msByNopsCutoff(ETree * etree,double cutoff,int type,int symflag)236 ETree_msByNopsCutoff (
237 ETree *etree,
238 double cutoff,
239 int type,
240 int symflag
241 ) {
242 double *tmetric ;
243 DV *tmetricDV, *vmetricDV ;
244 int front, nfront, nvtx, v ;
245 int *compids, *vtxToFront ;
246 IV *compidsIV ;
247 /*
248 ---------------
249 check the input
250 ---------------
251 */
252 if ( etree == NULL
253 || (nfront = etree->nfront) <= 0
254 || (nvtx = etree->nvtx) <= 0 ) {
255 fprintf(stderr, "\n fatal error in ETree_msByCutoff(%p,%f,%d)"
256 "\n bad input\n", etree, cutoff, symflag) ;
257 exit(-1) ;
258 }
259 vtxToFront = IV_entries(etree->vtxToFrontIV) ;
260 /*
261 ----------------------
262 get the subtree metric
263 ----------------------
264 */
265 vmetricDV = ETree_nopsMetric(etree, type, symflag) ;
266 tmetricDV = Tree_setSubtreeDmetric(etree->tree, vmetricDV) ;
267 #if MYDEBUG >= 0
268 fprintf(stdout, "\n ETree_msByNopsCutoff") ;
269 fprintf(stdout, "\n vmetric") ;
270 DV_writeForHumanEye(vmetricDV, stdout) ;
271 fprintf(stdout, "\n tmetric") ;
272 DV_writeForHumanEye(tmetricDV, stdout) ;
273 fflush(stdout) ;
274 fprintf(stdout, "\n max(tmetricDV) = %.0f, sum(vmetricDV) = %.0f",
275 DV_max(tmetricDV), DV_sum(vmetricDV)) ;
276 fprintf(stdout, "\n cutoff = %.0f", cutoff * DV_max(tmetricDV)) ;
277 #endif
278 cutoff = cutoff * DV_max(tmetricDV) ;
279 /*
280 --------------------------------------------------
281 fill the IV object with the ids in the multisector
282 --------------------------------------------------
283 */
284 compidsIV = IV_new() ;
285 IV_init(compidsIV, nvtx, NULL) ;
286 compids = IV_entries(compidsIV) ;
287 tmetric = DV_entries(tmetricDV) ;
288 for ( v = 0 ; v < nvtx ; v++ ) {
289 front = vtxToFront[v] ;
290 if ( tmetric[front] >= cutoff ) {
291 compids[v] = 0 ;
292 } else {
293 compids[v] = 1 ;
294 }
295 }
296 {
297 double domops, schurops ;
298 double *vmetric ;
299 int J ;
300
301 vmetric = DV_entries(vmetricDV) ;
302 domops = schurops = 0.0 ;
303 for ( J = 0 ; J < nfront ; J++ ) {
304 if ( tmetric[J] >= cutoff ) {
305 schurops += vmetric[J] ;
306 } else {
307 domops += vmetric[J] ;
308 }
309 }
310 fprintf(stdout, "\n domops = %.0f, schurops = %.0f, total = %.0f",
311 domops, schurops, domops + schurops) ;
312 }
313 /*
314 ------------------------
315 free the working storage
316 ------------------------
317 */
318 DV_free(vmetricDV) ;
319 DV_free(tmetricDV) ;
320
321 return(compidsIV) ; }
322
323 /*--------------------------------------------------------------------*/
324 /*
325 --------------------------------------------------------------
326 purpose -- given a front tree and a multisector map vector,
327 fill the map vector with domain ids and the three statistics
328 arrays with domain and schur complement statistics.
329
330 frontETree -- front tree object, unchanged on output
331 msIV -- map from fronts to domains or schur complement
332 on input, ms[v] = 0 --> v is in the schur complement
333 ms[v] = 1 --> v is not in the schur complement
334 on output, ms[v] = 0 --> v is in the schur complement
335 ms[v] != 0 --> v is in domain ms[v]
336 on output
337 nvtxIV -- nvtx[ireg] = # of dof in region ireg
338 nzfIV -- nzf[ireg] = # of factor entries in region ireg
339 opsIV -- ops[ireg] = # of factor ops in region ireg
340
341 type -- type of entries, SPOOLES_REAL or SPOOLES_COMPLEX
342
343 symflag -- symmetry flag, one of SPOOLES_SYMMETRIC
344 SPOOLES_HERMITIAN or SPOOLES_NONSYMMETRIC
345
346 created -- 98jan30, cca
347 --------------------------------------------------------------
348 */
349 void
ETree_msStats(ETree * frontETree,IV * msIV,IV * nvtxIV,IV * nzfIV,DV * opsDV,int type,int symflag)350 ETree_msStats (
351 ETree *frontETree,
352 IV *msIV,
353 IV *nvtxIV,
354 IV *nzfIV,
355 DV *opsDV,
356 int type,
357 int symflag
358 ) {
359 double *opsreg, *opsvec ;
360 DV *opsvecDV ;
361 int J, K, ndom, nfront, nvtx, reg, v ;
362 int *map, *ms, *nodwghts, *nvtxreg,
363 *nzfreg, *nzfvec, *par, *vtxToFront ;
364 IV *nzfvecIV ;
365 Tree *tree ;
366 /*
367 ---------------
368 check the input
369 ---------------
370 */
371 if ( frontETree == NULL || msIV == NULL || nvtxIV == NULL
372 || nzfIV == NULL || opsDV == NULL ) {
373 fprintf(stderr, "\n fatal error in ETree_msStats()"
374 "\n frontETree = %p, msIV = %p, nvtxIV = %p"
375 "\n nzfIV = %p, opsDV = %p, symflag = %d\n",
376 frontETree, msIV, nvtxIV, nzfIV, opsDV, symflag) ;
377 exit(-1) ;
378 }
379 nfront = ETree_nfront(frontETree) ;
380 nvtx = ETree_nvtx(frontETree) ;
381 tree = ETree_tree(frontETree) ;
382 par = ETree_par(frontETree) ;
383 vtxToFront = ETree_vtxToFront(frontETree) ;
384 ms = IV_entries(msIV) ;
385 /*
386 ----------------------------------------
387 if J is not in the schur complement then
388 fill ms[J] with its domain id
389 endif
390 ----------------------------------------
391 */
392 map = IVinit(nfront, -1) ;
393 for ( v = 0 ; v < nvtx ; v++ ) {
394 J = vtxToFront[v] ;
395 map[J] = ms[v] ;
396 /*
397 fprintf(stdout, "\n vertex %d, front %d, map %d", v, J, map[J]) ;
398 */
399 }
400 ndom = 0 ;
401 for ( J = Tree_preOTfirst(tree) ;
402 J != -1 ;
403 J = Tree_preOTnext(tree, J) ) {
404 /*
405 fprintf(stdout, "\n J = %d", J) ;
406 */
407 if ( map[J] != 0 ) {
408 if ( (K = par[J]) != -1 ) {
409 if ( map[K] == 0 ) {
410 map[J] = ++ndom ;
411 } else {
412 map[J] = map[K] ;
413 }
414 } else {
415 map[J] = ++ndom ;
416 }
417 /*
418 fprintf(stdout, ", in domain %d", map[J]) ;
419 } else {
420 fprintf(stdout, ", schur complement front") ;
421 */
422 }
423 }
424 for ( v = 0 ; v < nvtx ; v++ ) {
425 J = vtxToFront[v] ;
426 ms[v] = map[J] ;
427 /*
428 fprintf(stdout, "\n vertex %d, front %d, region %d", v, J, map[J]) ;
429 */
430 }
431 /*
432 --------------------------------------------------
433 set sizes of the nvtxIV, nzfIV and opsV vectors
434 to hold the domain and schur complement statistics
435 --------------------------------------------------
436 */
437 IV_setSize(nvtxIV, ndom+1) ;
438 IV_setSize(nzfIV, ndom+1) ;
439 DV_setSize(opsDV, ndom+1) ;
440 nvtxreg = IV_entries(nvtxIV) ;
441 nzfreg = IV_entries(nzfIV) ;
442 opsreg = DV_entries(opsDV) ;
443 IVzero(ndom+1, nvtxreg) ;
444 IVzero(ndom+1, nzfreg) ;
445 DVzero(ndom+1, opsreg) ;
446 /*
447 ---------------------------------
448 get the statistics for the fronts
449 ---------------------------------
450 */
451 nodwghts = ETree_nodwghts(frontETree) ;
452 nzfvecIV = ETree_factorEntriesIV(frontETree, symflag) ;
453 nzfvec = IV_entries(nzfvecIV) ;
454 opsvecDV = ETree_forwardOps(frontETree, type, symflag) ;
455 opsvec = DV_entries(opsvecDV) ;
456 /*
457 ----------------------------------
458 fill the region statistics vectors
459 ----------------------------------
460 */
461 for ( J = 0 ; J < nfront ; J++ ) {
462 reg = map[J] ;
463 nvtxreg[reg] += nodwghts[J] ;
464 nzfreg[reg] += nzfvec[J] ;
465 opsreg[reg] += opsvec[J] ;
466 }
467 /*
468 ------------------------
469 free the working storage
470 ------------------------
471 */
472 IV_free(nzfvecIV) ;
473 DV_free(opsvecDV) ;
474 IVfree(map) ;
475
476 return ; }
477
478 /*--------------------------------------------------------------------*/
479