1 /* mkNDlinsysQR.c */
2
3 #include "../../FrontMtx.h"
4 #include "../../EGraph.h"
5 #include "../../SymbFac.h"
6 #include "../../Drand.h"
7 #include "../../misc.h"
8 #include "../../timings.h"
9
10 /*--------------------------------------------------------------------*/
11 /*
12 ---------------------------------------------------------------------
13 purpose -- to create an overdetermined linear system A*X = B
14 for a nested dissection ordering of a 2-d or 3-d regular grid
15
16 input --
17
18 n1 -- # of nodes in first direction
19 n2 -- # of nodes in second direction
20 n3 -- # of nodes in third direction
21 type -- type of entries
22 SPOOLES_REAL or SPOOLES_COMPLEX
23 nrhs -- number of right hand sides
24 seed -- seed for random number generator
25 msglvl -- message level
26 msgFile -- message file
27
28 output --
29
30 pfrontETree -- to be filled with address of front tree
31 psymbfacIVL -- to be filled with address of symbolic factorization
32 pmtxA -- to be filled with address of matrix object A
33 pmtxX -- to be filled with address of matrix object X
34 pmtxB -- to be filled with address of matrix object B
35
36 created -- 98may29, cca
37 ---------------------------------------------------------------------
38 */
39 void
mkNDlinsysQR(int n1,int n2,int n3,int type,int nrhs,int seed,int msglvl,FILE * msgFile,ETree ** pfrontETree,IVL ** psymbfacIVL,InpMtx ** pmtxA,DenseMtx ** pmtxX,DenseMtx ** pmtxB)40 mkNDlinsysQR (
41 int n1,
42 int n2,
43 int n3,
44 int type,
45 int nrhs,
46 int seed,
47 int msglvl,
48 FILE *msgFile,
49 ETree **pfrontETree,
50 IVL **psymbfacIVL,
51 InpMtx **pmtxA,
52 DenseMtx **pmtxX,
53 DenseMtx **pmtxB
54 ) {
55 DenseMtx *mtxB, *mtxX ;
56 double one[2] = {1.0, 0.0} ;
57 Drand drand ;
58 double t1, t2 ;
59 double *entries ;
60 EGraph *egraph ;
61 ETree *etree, *frontETree ;
62 Graph *graph ;
63 InpMtx *mtxA ;
64 int ielem, ii, irow, jrow, mrow, ndim, nelem, nentA,
65 neqns, nfront, nrowA, size ;
66 int *indices, *newToOld, *oldToNew ;
67 IV *nzerosIV, *oldToNewIV ;
68 IVL *adjIVL, *symbfacIVL ;
69 /*
70 --------------------------------------
71 initialize the random number generator
72 --------------------------------------
73 */
74 Drand_setDefaultFields(&drand) ;
75 Drand_init(&drand) ;
76 Drand_setSeed(&drand, seed) ;
77 Drand_setNormal(&drand, 0.0, 1.0) ;
78 /*
79 -----------------------------------------------------
80 get the grid adjacency structure and set up the graph
81 -----------------------------------------------------
82 */
83 neqns = n1*n2*n3 ;
84 MARKTIME(t1) ;
85 if ( (n1 == 1 && n2 == 1)
86 || (n1 == 1 && n3 == 1)
87 || (n1 == 2 && n3 == 1) ) {
88 fprintf(stderr, "\n fatal error in mkNDlinsysQR"
89 "\n n1 = %d, n2 = %d, n3 = %d\n", n1, n2, n3) ;
90 exit(-1) ;
91 }
92 if ( n1 == 1 ) {
93 adjIVL = IVL_make9P(n2, n3, 1) ;
94 ndim = 2 ;
95 } else if ( n2 == 1 ) {
96 adjIVL = IVL_make9P(n1, n3, 1) ;
97 ndim = 2 ;
98 } else if ( n3 == 1 ) {
99 adjIVL = IVL_make9P(n1, n2, 1) ;
100 ndim = 2 ;
101 } else {
102 adjIVL = IVL_make27P(n1, n2, n3, 1) ;
103 ndim = 3 ;
104 }
105 graph = Graph_new() ;
106 Graph_init2(graph, 0, neqns, 0, adjIVL->tsize, neqns,
107 adjIVL->tsize, adjIVL, NULL, NULL) ;
108 MARKTIME(t2) ;
109 fprintf(msgFile, "\n CPU %8.3f : create the grid graph",
110 t2 - t1) ;
111 if ( msglvl > 3 ) {
112 fprintf(msgFile, "\n\n grid graph") ;
113 Graph_writeForHumanEye(graph, msgFile) ;
114 }
115 /*
116 ---------------------------------------------
117 make the nested dissection permutation vector
118 ---------------------------------------------
119 */
120 MARKTIME(t1) ;
121 newToOld = IVinit(neqns, -1) ;
122 oldToNew = IVinit(neqns, -1) ;
123 mkNDperm(n1, n2, n3, newToOld, 0, n1-1, 0, n2-1, 0, n3-1) ;
124 for ( ii = 0 ; ii < neqns ; ii++ ) {
125 oldToNew[newToOld[ii]] = ii ;
126 }
127 MARKTIME(t2) ;
128 fprintf(msgFile, "\n CPU %8.3f : make the nested dissection ordering",
129 t2 - t1) ;
130 if ( msglvl > 3 ) {
131 fprintf(msgFile, "\n\n oldToNew") ;
132 IVfprintf(msgFile, neqns, oldToNew) ;
133 }
134 /*
135 ----------------------------------
136 create the elimination tree object
137 ----------------------------------
138 */
139 MARKTIME(t1) ;
140 etree = ETree_new() ;
141 ETree_initFromGraphWithPerms(etree, graph, newToOld, oldToNew) ;
142 if ( msglvl > 3 ) {
143 fprintf(msgFile, "\n\n elimination tree") ;
144 ETree_writeForHumanEye(etree, msgFile) ;
145 }
146 nzerosIV = IV_new() ;
147 IV_init(nzerosIV, neqns, NULL) ;
148 IV_fill(nzerosIV, 0) ;
149 frontETree = ETree_mergeFrontsOne(etree, 0, nzerosIV) ;
150 IV_free(nzerosIV) ;
151 IVfree(newToOld) ;
152 IVfree(oldToNew) ;
153 ETree_free(etree) ;
154 nfront = frontETree->nfront ;
155 if ( msglvl > 3 ) {
156 fprintf(msgFile, "\n\n front tree") ;
157 ETree_writeForHumanEye(frontETree, msgFile) ;
158 }
159 MARKTIME(t2) ;
160 fprintf(msgFile, "\n CPU %8.3f : create the front tree",
161 t2 - t1) ;
162 /*
163 --------------------------------------------
164 create the symbolic factorization IVL object
165 --------------------------------------------
166 */
167 MARKTIME(t1) ;
168 symbfacIVL = SymbFac_initFromGraph(frontETree, graph) ;
169 MARKTIME(t2) ;
170 fprintf(msgFile, "\n CPU %8.3f : compute the symbolic factorization",
171 t2 - t1) ;
172 if ( msglvl > 1 ) {
173 fprintf(msgFile,
174 "\n\n symbolic factorization IVL object in original ordering") ;
175 if ( msglvl == 2 ) {
176 IVL_writeStats(symbfacIVL, msgFile) ;
177 } else {
178 IVL_writeForHumanEye(symbfacIVL, msgFile) ;
179 }
180 fflush(msgFile) ;
181 }
182 /*
183 --------------------------------------
184 permute the vertices in the front tree
185 --------------------------------------
186 */
187 MARKTIME(t1) ;
188 oldToNewIV = ETree_oldToNewVtxPerm(frontETree) ;
189 oldToNew = IV_entries(oldToNewIV) ;
190 ETree_permuteVertices(frontETree, oldToNewIV) ;
191 MARKTIME(t2) ;
192 fprintf(msgFile, "\n CPU %8.3f : permute the front tree",
193 t2 - t1) ;
194 /*
195 ------------------------------------------------------
196 convert the symbolic factorization to the new ordering
197 ------------------------------------------------------
198 */
199 IVL_overwrite(symbfacIVL, oldToNewIV) ;
200 if ( msglvl > 1 ) {
201 fprintf(msgFile,
202 "\n\n overwritten symbolic factorization IVL object ") ;
203 if ( msglvl == 2 ) {
204 IVL_writeStats(symbfacIVL, msgFile) ;
205 } else {
206 IVL_writeForHumanEye(symbfacIVL, msgFile) ;
207 }
208 fflush(msgFile) ;
209 }
210 IVL_sortUp(symbfacIVL) ;
211 if ( msglvl > 1 ) {
212 fprintf(msgFile,
213 "\n\n sorted symbolic factorization IVL object ") ;
214 if ( msglvl == 2 ) {
215 IVL_writeStats(symbfacIVL, msgFile) ;
216 } else {
217 IVL_writeForHumanEye(symbfacIVL, msgFile) ;
218 }
219 fflush(msgFile) ;
220 }
221 /*
222 --------------------------------
223 create the natural factor matrix
224 --------------------------------
225 */
226 MARKTIME(t1) ;
227 if ( ndim == 2 ) {
228 if ( n1 == 1 ) {
229 egraph = EGraph_make9P(n2, n3, 1) ;
230 nelem = (n2-1)*(n3-1) ;
231 } else if ( n2 == 1 ) {
232 egraph = EGraph_make9P(n1, n3, 1) ;
233 nelem = (n1-1)*(n3-1) ;
234 } else if ( n3 == 1 ) {
235 egraph = EGraph_make9P(n1, n2, 1) ;
236 nelem = (n1-1)*(n2-1) ;
237 }
238 mrow = 4 ;
239 } else {
240 egraph = EGraph_make27P(n1, n2, n3, 1) ;
241 mrow = 8 ;
242 nelem = (n1-1)*(n2-1)*(n3-1) ;
243 }
244 nrowA = mrow*nelem ;
245 if ( type == SPOOLES_REAL ) {
246 entries = DVinit(mrow, 0.0) ;
247 } else if ( type == SPOOLES_COMPLEX ) {
248 entries = DVinit(2*mrow, 0.0) ;
249 }
250 nentA = mrow*nrowA ;
251 MARKTIME(t2) ;
252 fprintf(msgFile, "\n CPU %8.3f : create egraph ", t2 - t1) ;
253 if ( msglvl > 2 ) {
254 EGraph_writeForHumanEye(egraph, msgFile) ;
255 fflush(msgFile) ;
256 }
257 MARKTIME(t1) ;
258 mtxA = InpMtx_new() ;
259 InpMtx_init(mtxA, INPMTX_BY_ROWS, type, nentA, nrowA) ;
260 MARKTIME(t2) ;
261 fprintf(msgFile, "\n CPU %8.3f : initialize the InpMtx object",
262 t2 - t1) ;
263 if ( msglvl > 1 ) {
264 fprintf(msgFile, "\n InpMtx after initialization") ;
265 InpMtx_writeForHumanEye(mtxA, msgFile) ;
266 }
267 for ( ielem = 0, jrow = 0 ; ielem < nelem ; ielem++ ) {
268 IVL_listAndSize(egraph->adjIVL, ielem, &size, &indices) ;
269 for ( irow = 0 ; irow < mrow ; irow++, jrow++ ) {
270 if ( type == SPOOLES_REAL ) {
271 Drand_fillDvector(&drand, size, entries) ;
272 InpMtx_inputRealRow(mtxA, jrow, size, indices, entries) ;
273 } else if ( type == SPOOLES_COMPLEX ) {
274 Drand_fillDvector(&drand, 2*size, entries) ;
275 InpMtx_inputComplexRow(mtxA, jrow, size, indices, entries) ;
276 }
277 }
278 }
279 DVfree(entries) ;
280 EGraph_free(egraph) ;
281 if ( msglvl > 1 ) {
282 fprintf(msgFile, "\n\n natural factor matrix") ;
283 if ( msglvl == 2 ) {
284 InpMtx_writeStats(mtxA, msgFile) ;
285 } else {
286 InpMtx_writeForHumanEye(mtxA, msgFile) ;
287 }
288 fflush(msgFile) ;
289 }
290 /*
291 -------------------------------------------------------------
292 permute the InpMtx object into the nested dissection ordering
293 -------------------------------------------------------------
294 */
295 MARKTIME(t1) ;
296 InpMtx_permute(mtxA, NULL, IV_entries(oldToNewIV)) ;
297 InpMtx_changeCoordType(mtxA, INPMTX_BY_ROWS) ;
298 InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
299 MARKTIME(t2) ;
300 fprintf(msgFile, "\n CPU %8.3f : permute inpmtx ", t2 - t1) ;
301 if ( msglvl > 1 ) {
302 fprintf(msgFile, "\n\n after permute InpMtx object ") ;
303 if ( msglvl == 2 ) {
304 InpMtx_writeStats(mtxA, msgFile) ;
305 } else {
306 InpMtx_writeForHumanEye(mtxA, msgFile) ;
307 }
308 }
309 if ( msglvl > 5 ) {
310 fprintf(msgFile, "\n A = zeros(%d,%d) ;", nrowA, neqns) ;
311 InpMtx_writeForMatlab(mtxA, "A", msgFile) ;
312 }
313 /*
314 --------------------------------------------------------
315 generate the linear system
316 1. generate solution matrix and fill with random numbers
317 2. generate rhs matrix and fill with zeros
318 3. compute matrix-matrix multiply
319 --------------------------------------------------------
320 */
321 mtxX = DenseMtx_new() ;
322 DenseMtx_init(mtxX, type, 0, -1, neqns, nrhs, 1, neqns) ;
323 DenseMtx_fillRandomEntries(mtxX, &drand) ;
324 mtxB = DenseMtx_new() ;
325 DenseMtx_init(mtxB, type, 1, -1, nrowA, nrhs, 1, nrowA) ;
326 DenseMtx_zero(mtxB) ;
327 InpMtx_nonsym_mmm(mtxA, mtxB, one, mtxX) ;
328 if ( msglvl > 3 ) {
329 fprintf(msgFile, "\n\n mtxX") ;
330 DenseMtx_writeForHumanEye(mtxX, msgFile) ;
331 fprintf(msgFile, "\n\n mtxB") ;
332 DenseMtx_writeForHumanEye(mtxB, msgFile) ;
333 fflush(msgFile) ;
334 }
335 if ( msglvl > 5 ) {
336 fprintf(msgFile, "\n X = zeros(%d,%d) ;", mtxX->nrow, mtxX->ncol) ;
337 DenseMtx_writeForMatlab(mtxX, "X", msgFile) ;
338 fprintf(msgFile, "\n B = zeros(%d,%d) ;", mtxB->nrow, mtxB->ncol) ;
339 DenseMtx_writeForMatlab(mtxB, "B", msgFile) ;
340 }
341 /*
342 -----------------------
343 set the output pointers
344 -----------------------
345 */
346 *pfrontETree = frontETree ;
347 *psymbfacIVL = symbfacIVL ;
348 *pmtxA = mtxA ;
349 *pmtxX = mtxX ;
350 *pmtxB = mtxB ;
351 /*
352 ------------------------
353 free the working storage
354 ------------------------
355 */
356 Graph_free(graph) ;
357 IV_free(oldToNewIV) ;
358
359 return ; }
360
361 /*--------------------------------------------------------------------*/
362