1 /* factor.c */
2
3 #include "../FrontMtx.h"
4 #include "../../timings.h"
5
6 /*--------------------------------------------------------------------*/
7 /*
8 -------------------------------------------------------------------
9 compute an (U^T + I)D(I + U) or (L + I)D(I + L) factorization of A.
10 this is a wrapper method around FrontMtx_factorPencil().
11
12 input --
13
14 frontmtx -- pointer to the FrontMtx object that will hold
15 the factorization
16 pencil -- pointer to the Pencil object that holds A + sigma*B
17 tau -- upper bound on entries in L and U,
18 used only when pivoting enabled
19 droptol -- lower bound on entries in L and U,
20 used only when sparsity enabled
21 perror -- error flag, on return
22 *perror >= 0 --> factorization failed at front *perror
23 *perror < 0 --> factorization succeeded
24 cpus[] -- timing array
25 cpus[0] -- initialize fronts
26 cpus[1] -- load original entries
27 cpus[2] -- get updates from descendents
28 cpus[3] -- assembled postponed data
29 cpus[4] -- factor the front
30 cpus[5] -- extract postponed data
31 cpus[6] -- store factor entries
32 cpus[7] -- miscellaneous time
33 cpus[8] -- total time
34 stats[] -- statistics array
35 stats[0] -- # of pivots
36 stats[1] -- # of pivot tests
37 stats[2] -- # of delayed rows and columns
38 stats[3] -- # of entries in D
39 stats[4] -- # of entries in L
40 stats[5] -- # of entries in U
41 msglvl -- message level
42 msgFile -- message file
43
44 created -- 98mar27, cca
45 modified -- 98mar27, cca
46 perror added to argument list
47 -------------------------------------------------------------------
48 */
49 Chv *
FrontMtx_factorInpMtx(FrontMtx * frontmtx,InpMtx * inpmtx,double tau,double droptol,ChvManager * chvmanager,int * perror,double cpus[],int stats[],int msglvl,FILE * msgFile)50 FrontMtx_factorInpMtx (
51 FrontMtx *frontmtx,
52 InpMtx *inpmtx,
53 double tau,
54 double droptol,
55 ChvManager *chvmanager,
56 int *perror,
57 double cpus[],
58 int stats[],
59 int msglvl,
60 FILE *msgFile
61 ) {
62 double zero[2] = {0.0, 0.0} ;
63 Chv *rootchv ;
64 Pencil pencil ;
65
66 Pencil_setDefaultFields(&pencil) ;
67 Pencil_init(&pencil, frontmtx->type, frontmtx->symmetryflag,
68 inpmtx, zero, NULL) ;
69 rootchv = FrontMtx_factorPencil(frontmtx, &pencil, tau, droptol,
70 chvmanager, perror, cpus, stats, msglvl, msgFile) ;
71
72 return(rootchv) ; }
73
74 /*--------------------------------------------------------------------*/
75 /*
76 -------------------------------------------------------------------
77 compute an (U^T + I)D(I + U) or (L + I)D(I + L)
78 factorization of A + sigma*B.
79
80 input --
81
82 frontmtx -- pointer to the FrontMtx object that will hold
83 the factorization
84 pencil -- pointer to the Pencil object that holds A + sigma*B
85 tau -- upper bound on entries in L and U,
86 used only when pivoting enabled
87 droptol -- lower bound on entries in L and U,
88 used only when sparsity enabled
89 perror -- error flag, on return
90 *perror >= 0 --> factorization failed at front *perror
91 *perror < 0 --> factorization succeeded
92 cpus[] -- timing array
93 cpus[0] -- initialize fronts
94 cpus[1] -- load original entries
95 cpus[2] -- get updates from descendents
96 cpus[3] -- assembled postponed data
97 cpus[4] -- factor the front
98 cpus[5] -- extract postponed data
99 cpus[6] -- store factor entries
100 cpus[7] -- miscellaneous time
101 cpus[8] -- total time
102 stats[] -- statistics array
103 stats[0] -- # of pivots
104 stats[1] -- # of pivot tests
105 stats[2] -- # of delayed rows and columns
106 stats[3] -- # of entries in D
107 stats[4] -- # of entries in L
108 stats[5] -- # of entries in U
109 msglvl -- message level
110 msgFile -- message file
111
112 created -- 98mar27, cca
113 modified -- 98mar27, cca
114 perror added to argument list
115 -------------------------------------------------------------------
116 */
117 Chv *
FrontMtx_factorPencil(FrontMtx * frontmtx,Pencil * pencil,double tau,double droptol,ChvManager * chvmanager,int * perror,double cpus[],int stats[],int msglvl,FILE * msgFile)118 FrontMtx_factorPencil (
119 FrontMtx *frontmtx,
120 Pencil *pencil,
121 double tau,
122 double droptol,
123 ChvManager *chvmanager,
124 int *perror,
125 double cpus[],
126 int stats[],
127 int msglvl,
128 FILE *msgFile
129 ) {
130 char *status ;
131 ChvList *postList ;
132 Chv *rootchv ;
133 Chv **fronts ;
134 double t0, t3 ;
135 DV tmpDV ;
136 ETree *frontETree ;
137 int J, K, ndelayed, nfront, npivots, ntests ;
138 int *par ;
139 IP **heads ;
140 IV pivotsizesIV ;
141 Tree *tree ;
142 /*
143 ---------------
144 check the input
145 ---------------
146 */
147 MARKTIME(t0) ;
148 if ( frontmtx == NULL || pencil == NULL || cpus == NULL || stats == NULL
149 || (msglvl > 0 && msgFile == NULL) ) {
150 fprintf(stderr, "\n fatal error in FrontMtx_factorPencil()"
151 "\n frontmtx = %p, pencil = %p"
152 "\n tau = %e, droptol = %e, cpus = %p"
153 "\n msglvl = %d, msgFile = %p"
154 "\n bad input\n",
155 frontmtx, pencil, tau, droptol, cpus, msglvl, msgFile) ;
156 exit(-1) ;
157 }
158 if ( msglvl > 2 ) {
159 fprintf(msgFile, "\n\n INSIDE FrontMtx_factorPencil()") ;
160 fflush(msgFile) ;
161 }
162 /*
163 -------------------------------
164 extract pointers and dimensions
165 -------------------------------
166 */
167 frontETree = frontmtx->frontETree ;
168 nfront = ETree_nfront(frontETree) ;
169 tree = ETree_tree(frontETree) ;
170 par = ETree_par(frontETree) ;
171 if ( msglvl > 2 ) {
172 fprintf(msgFile, "\n got pointers and dimensions") ;
173 fflush(msgFile) ;
174 }
175 /*
176 ---------------------------------------
177 allocate and initialize working storage
178 ---------------------------------------
179 */
180 heads = FrontMtx_factorSetup(frontmtx, NULL, 0, msglvl, msgFile) ;
181 status = CVinit(nfront, 'W') ;
182 ALLOCATE(fronts, Chv *, nfront) ;
183 for ( J = 0 ; J < nfront ; J++ ) {
184 fronts[J] = NULL ;
185 }
186 DV_setDefaultFields(&tmpDV) ;
187 IV_setDefaultFields(&pivotsizesIV) ;
188 if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
189 postList = ChvList_new() ;
190 ChvList_init(postList, nfront+1, NULL, 0, NULL) ;
191 } else {
192 postList = NULL ;
193 }
194 if ( msglvl > 1 ) {
195 fprintf(msgFile, "\n allocated working storage") ;
196 fflush(msgFile) ;
197 }
198 /*
199 --------------------------------------------
200 loop over the tree in a post-order traversal
201 --------------------------------------------
202 */
203 *perror = -1 ;
204 npivots = ndelayed = ntests = 0 ;
205 rootchv = NULL ;
206 for ( J = Tree_postOTfirst(tree) ;
207 J != -1 ;
208 J = Tree_postOTnext(tree, J) ) {
209 K = par[J] ;
210 if ( msglvl > 1 ) {
211 fprintf(msgFile, "\n\n ##### working on front %d, parent %d",
212 J, K) ;
213 fflush(msgFile) ;
214 }
215 FrontMtx_factorVisit(frontmtx, pencil, J, 0, NULL, fronts, 0,
216 tau, droptol, status, heads, &pivotsizesIV,
217 &tmpDV, par, NULL, postList,
218 chvmanager, stats, cpus, msglvl, msgFile) ;
219 if ( status[J] == 'E' ) {
220 /*
221 -------------------------------
222 error found, unable to continue
223 -------------------------------
224 */
225 *perror = J ;
226 break ;
227 } else if ( status[J] != 'F' ) {
228 fprintf(stderr, "\n fatal error, return %c from front %d",
229 status[J], J) ;
230 exit(-1) ;
231 }
232 }
233 /*
234 ---------------------------------
235 get a pointer to the root chevron
236 ---------------------------------
237 */
238 if ( postList == NULL ) {
239 rootchv = NULL ;
240 } else {
241 rootchv = ChvList_getList(postList, nfront) ;
242 }
243 /*
244 ------------------
245 set the statistics
246 ------------------
247 */
248 stats[3] = frontmtx->nentD ;
249 stats[4] = frontmtx->nentL ;
250 stats[5] = frontmtx->nentU ;
251 /*
252 ------------------------
253 free the working storage
254 ------------------------
255 */
256 IP_free(heads[nfront+1]) ;
257 FREE(heads) ;
258 DV_clearData(&tmpDV) ;
259 IV_clearData(&pivotsizesIV) ;
260 CVfree(status) ;
261 FREE(fronts) ;
262 if ( postList != NULL ) {
263 ChvList_free(postList) ;
264 }
265 /*
266 --------------------------------
267 set final and miscellaneous time
268 --------------------------------
269 */
270 MARKTIME(t3) ;
271 cpus[8] = t3 - t0 ;
272 cpus[7] = cpus[8] - cpus[0] - cpus[1] - cpus[2]
273 - cpus[3] - cpus[4] - cpus[5] - cpus[6] ;
274
275 return(rootchv) ; }
276
277 /*--------------------------------------------------------------------*/
278