Actual source code: mlSerial.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: mlSerial.c,v 1.12 2000/10/22 20:27:52 knepley Exp $";
3: #endif
4: /*
5: Defines the serial multilevel preconditioner routines
6: */
7: #include src/sles/pc/pcimpl.h
8: #include ml.h
10: #undef __FUNCT__
12: PetscTruth PCMultiLevelDoQR_Private(PC pc, int numRows, int numCols)
13: {
14: /* This functions decides whether to QR factorize the matrix U of an SVD.
15: This is a memory saving device when B = UDV is long and skinny.
16: */
17: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
20: if (numRows > (int) numCols*ml->DoQRFactor) {
21: PetscFunctionReturn(PETSC_TRUE);
22: } else {
23: PetscFunctionReturn(PETSC_FALSE);
24: }
25: }
27: #undef __FUNCT__
29: int PCMultiLevelPartitionMesh_Private(PC pc, int level, int **mesh, int *numParts, int *colPartition)
30: {
31: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
32: PetscTruth change;
33: int numCols, freeCols;
34: int numAdj, newNumAdj, numFound;
35: int *offsets;
36: int *adj;
37: int *sizes;
38: int part, col, newCol, adjCol, newAdjCol;
39: int ierr;
42: PC_MLLogEventBegin(PC_ML_ReducePartitionMesh, pc, 0, 0, 0);
43: if (ml->useMath == PETSC_FALSE) {
44: offsets = mesh[MESH_OFFSETS];
45: adj = mesh[MESH_ADJ];
46: /* Set initial number of partitions */
47: numCols = ml->numVertices[level];
48: *numParts = numCols/ml->partSize;
49: /* If the graph is small enough, just return the whole thing */
50: if (*numParts == 0) {
51: *numParts = 1;
52: PetscMemzero(colPartition, numCols * sizeof(int));
53: PC_MLLogEventEnd(PC_ML_ReducePartitionMesh, pc, 0, 0, 0);
54: return(0);
55: }
56: /* Get connectivity information */
57: /* Initialize partition and partition sizes */
58: PetscMalloc((*numParts) * sizeof(int), &sizes);
59: for(col = 0; col < numCols; col++) colPartition[col] = -1;
60: /* Pick numParts starting points -- Conforms to the Mathematica */
61: for(part = 1; part <= *numParts; part++) {
62: colPartition[ml->partSize*part-1] = part-1;
63: sizes[part-1] = 1;
64: }
65: /* Start enlarging partitions */
66: freeCols = numCols - (*numParts);
67: change = PETSC_TRUE;
68: while(freeCols > 0)
69: {
70: /* Quit if nothing more can be done */
71: if (change == PETSC_FALSE)
72: break;
73: change = PETSC_FALSE;
74: /* Add a node to each partition if possible */
75: for(part = 0; part < (*numParts); part++)
76: {
77: /* Look for the column adjacent to the most columns in this partition */
78: newCol = -1;
79: numAdj = 0;
80: for(col = 0, numFound = 0; (col < numCols) && (numFound < sizes[part]); col++)
81: {
82: if (colPartition[col] == part)
83: {
84: numFound++;
85: /* Look for adjacent columns ignoring self edge */
86: for(adjCol = offsets[col]+1; adjCol < offsets[col+1]; adjCol++)
87: {
88: /* Ignore columns already in a partition */
89: if (colPartition[adj[adjCol]] >= 0)
90: continue;
91: newNumAdj = 0;
92: /* Count the number of columns adjacent to this one that are already in the partition */
93: for(newAdjCol = offsets[adj[adjCol]]+1; newAdjCol < offsets[adj[adjCol]+1]; newAdjCol++)
94: if (colPartition[adj[newAdjCol]] == part)
95: newNumAdj++;
96: if ((newNumAdj > numAdj) || ((newNumAdj == numAdj) && (adj[adjCol] < newCol)))
97: {
98: newCol = adj[adjCol];
99: numAdj = newNumAdj;
100: }
101: }
102: }
103: }
104: if (newCol >= 0)
105: {
106: colPartition[newCol] = part;
107: sizes[part]++;
108: freeCols--;
109: change = PETSC_TRUE;
110: }
111: }
112: }
113: /* Take care of extra columns */
114: if (change == PETSC_FALSE) {
115: for(col = 0; col < numCols; col++)
116: if (colPartition[col] < 0)
117: colPartition[col] = (*numParts)++;
118: }
119: PetscFree(sizes);
120: } else {
121: #ifdef PETSC_HAVE_MATHEMATICA
122: ViewerMathematicaPartitionMesh(ml->mathViewer, ml->numElements[level], ml->numVertices[level],
123: ml->vertices, mesh, numParts, colPartition);
124:
125: #else
126: SETERRQ(PETSC_ERR_SUP, " ");
127: #endif
128: }
129: PC_MLLogEventEnd(PC_ML_ReducePartitionMesh, pc, 0, 0, 0);
131: return(0);
132: }
134: #undef __FUNCT__
136: int PCMultiLevelPartitionCols_Private(PC pc, int level, int numParts, int numCols, int *colPartition, int *globalCols)
137: {
138: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
139: int numLocCols;
140: int part, col;
141: int ierr;
144: PetscMalloc(numParts * sizeof(int), &ml->numPartitionCols[level]);
145: PetscLogObjectMemory(pc, numParts * sizeof(int));
146: PetscMemzero(ml->numPartitionCols[level], numParts * sizeof(int));
147: for(col = 0; col < numCols; col++) {
148: ml->numPartitionCols[level][colPartition[col]]++;
149: }
150: PetscMalloc(numParts * sizeof(int *), &ml->colPartition[level]);
151: PetscLogObjectMemory(pc, numParts * sizeof(int *));
152: for(part = 0; part < numParts; part++) {
153: numLocCols = ml->numPartitionCols[level][part];
154: PetscMalloc(numLocCols * sizeof(int), &ml->colPartition[level][part]);
155: PetscLogObjectMemory(pc, numLocCols * sizeof(int));
156: ml->numPartitionCols[level][part] = 0;
157: }
158: for(col = 0; col < numCols; col++) {
159: ml->colPartition[level][colPartition[col]][ml->numPartitionCols[level][colPartition[col]]++] = globalCols[col];
160: }
161: for(part = 0, col = 0; part < numParts; part++) {
162: col += ml->numPartitionCols[level][part];
163: }
164: if (col != numCols) {
165: int rank;
167: MPI_Comm_rank(pc->comm, &rank);
168: PetscPrintf(PETSC_COMM_SELF, "[%d] %d does not equal the number of columns %dn", rank, col, numCols);
169: for(col = 0; col < numCols; col++) {
170: PetscPrintf(PETSC_COMM_SELF, "[%d] part: %d col: %dn", rank, part, colPartition[col]);
171: }
172: SETERRQ(PETSC_ERR_PLIB, "Invalid column partition");
173: }
175: return(0);
176: }
178: #undef __FUNCT__
180: int PCMultiLevelPartitionRows_Private(PC pc, int level, int numParts, int numCols, int numRows, int *colPartition,
181: int *rowPartition, int *globalCols, Mat grad)
182: {
183: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
184: int numLocRows, numLocCols, ncols, numSlots;
185: int *cols;
186: PetscScalar *vals;
187: int part, col, locCol, row;
188: int ierr;
191: /* Allocation */
192: numSlots = numParts + NUM_PART_ROW_DIV - 1;
193: PetscMalloc(numSlots * sizeof(int), &ml->numPartitionRows[level]);
194: PetscMalloc(NUM_PART_ROW_DIV * sizeof(int **), &ml->rowPartition[level]);
195: PetscMalloc(numParts * sizeof(int *), &ml->rowPartition[level][PART_ROW_INT]);
196: PetscMalloc(1 * sizeof(int *), &ml->rowPartition[level][PART_ROW_BD]);
197: PetscMalloc(1 * sizeof(int *), &ml->rowPartition[level][PART_ROW_RES]);
198: PetscLogObjectMemory(pc, numSlots*sizeof(int) + NUM_PART_ROW_DIV*sizeof(int **) + numSlots*sizeof(int *));
200: /* Create row partition */
201: for(row = 0; row < numRows; row++) {
202: /* Find partitions adjacent to this row */
203: MatGetRow(grad, row, &ncols, &cols, &vals);
204: PetscMemzero(ml->numPartitionRows[level], numParts * sizeof(int));
205: for(col = 0; col < ncols; col++) {
206: /* Could invert this mapping explicitly or do bisection */
207: for(locCol = 0; locCol < numCols; locCol++)
208: if (globalCols[locCol] == cols[col])
209: break;
210: if (locCol == numCols) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid boundary gradient");
211: ml->numPartitionRows[level][colPartition[locCol]] = 1;
212: }
213: MatRestoreRow(grad, row, &ncols, &cols, &vals);
215: /* Classify row */
216: if (ncols == 0)
217: {
218: /* This is caused by rows with nonzeros only in full rank partitions, so that
219: B_Gamma (I - D^{-1} D) has a null row. */
220: rowPartition[row] = 0;
221: continue;
222: }
223: for(part = 0, numLocCols = 0; part < numParts; part++) {
224: if (ml->numPartitionRows[level][part] > 0) {
225: numLocCols++;
226: rowPartition[row] = part;
227: }
228: }
229: if (numLocCols == 2) {
230: rowPartition[row] = numParts + PART_ROW_BD - 1;
231: } else if (numLocCols > 2) {
232: rowPartition[row] = numParts + PART_ROW_RES - 1;
233: } else if (numLocCols < 1) {
234: SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid boundary gradient");
235: }
236: }
238: /* Copy rows */
239: PetscMemzero(ml->numPartitionRows[level], (numParts + NUM_PART_ROW_DIV - 1) * sizeof(int));
240: for(row = 0; row < numRows; row++) {
241: ml->numPartitionRows[level][rowPartition[row]]++;
242: }
243: for(part = 0; part < numParts; part++) {
244: numLocRows = ml->numPartitionRows[level][part];
245: if (numLocRows) {
246: PetscMalloc(numLocRows * sizeof(int), &ml->rowPartition[level][PART_ROW_INT][part]);
247: PetscLogObjectMemory(pc, numLocRows * sizeof(int));
248: }
249: ml->numPartitionRows[level][part] = 0;
250: }
251: numLocRows = ml->numPartitionRows[level][numParts];
252: if (numLocRows) {
253: PetscMalloc(numLocRows * sizeof(int), &ml->rowPartition[level][PART_ROW_BD][0]);
254: PetscLogObjectMemory(pc, numLocRows * sizeof(int));
255: }
256: ml->numPartitionRows[level][numParts] = 0;
257: numLocRows = ml->numPartitionRows[level][numParts+1];
258: if (numLocRows) {
259: PetscMalloc(numLocRows * sizeof(int), &ml->rowPartition[level][PART_ROW_RES][0]);
260: PetscLogObjectMemory(pc, numLocRows * sizeof(int));
261: }
262: ml->numPartitionRows[level][numParts+1] = 0;
263: for(row = 0; row < numRows; row++) {
264: if (rowPartition[row] < numParts) {
265: ml->rowPartition[level][PART_ROW_INT][rowPartition[row]][ml->numPartitionRows[level][rowPartition[row]]++] = row;
266: } else if (rowPartition[row] == numParts) {
267: ml->rowPartition[level][PART_ROW_BD][0][ml->numPartitionRows[level][rowPartition[row]]++] = row;
268: } else if (rowPartition[row] > numParts) {
269: ml->rowPartition[level][PART_ROW_RES][0][ml->numPartitionRows[level][rowPartition[row]]++] = row;
270: }
271: }
273: return(0);
274: }
276: #undef __FUNCT__
278: int PCMultiLevelFactorPartitions_Private(PC pc, int level, int numParts, Mat grad)
279: {
280: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
281: int ncols;
282: int *cols;
283: PetscScalar *vals;
284: PetscReal *B;
285: int maxSize;
286: int numRows, numCols;
287: int part, row, col, locCol;
288: int ierr;
291: PC_MLLogEventBegin(PC_ML_ReduceFactor, pc, grad, 0, 0);
292: /* Allocation */
293: for(part = 0, maxSize = 0; part < numParts; part++) {
294: numRows = ml->numPartitionRows[level][part];
295: numCols = ml->numPartitionCols[level][part];
296: maxSize = PetscMax(maxSize, numRows*numCols);
297: }
298: if (maxSize == 0) {
299: PetscFunctionReturn(2);
300: /* SETERRQ(1, "Empty interior (increase partition size)"); */
301: }
302: PetscMalloc(numParts * sizeof(PetscReal **), &ml->factors[level]);
303: PetscLogObjectMemory(pc, numParts * sizeof(PetscReal **));
304: PetscMalloc(maxSize * sizeof(PetscReal), &B);
305:
306: for(part = 0; part < numParts; part++) {
307: /* Allocation */
308: numRows = ml->numPartitionRows[level][part];
309: numCols = ml->numPartitionCols[level][part];
310: PetscMalloc(NUM_FACT_DIV * sizeof(PetscReal *), &ml->factors[level][part]);
311: PetscMalloc(numCols * sizeof(PetscReal), &ml->factors[level][part][FACT_DINV]);
312: PetscMalloc(numCols*numCols * sizeof(PetscReal), &ml->factors[level][part][FACT_V]);
313: PetscLogObjectMemory(pc, NUM_FACT_DIV * sizeof(double *) + (numCols + numCols*numCols)*sizeof(PetscReal));
315: if (numRows > 0) {
316: /* If numCols < numRows, LAPACK will leave a bogus values in the end of D^{-1} */
317: if (numCols > numRows) {
318: PetscMemzero(ml->factors[level][part][FACT_DINV], numCols * sizeof(PetscReal));
319: }
321: if (PCMultiLevelDoQR_Private(pc, numRows, numCols) == PETSC_TRUE) {
322: /* Allocation */
323: PetscMalloc(numRows*numCols * sizeof(PetscReal), &ml->factors[level][part][FACT_QR]);
324: PetscMalloc(numCols * sizeof(PetscReal), &ml->factors[level][part][FACT_TAU]);
325: PetscMalloc(numCols*numCols * sizeof(PetscReal), &ml->factors[level][part][FACT_U]);
326: PetscLogObjectMemory(pc, (numRows*numCols + numCols + numCols*numCols)*sizeof(PetscReal));
328: /* Extract the interior gradient for this partition */
329: PetscMemzero(ml->factors[level][part][FACT_QR], numRows*numCols * sizeof(PetscReal));
330: for(row = 0; row < numRows; row++) {
331: MatGetRow(grad, ml->rowPartition[level][PART_ROW_INT][part][row], &ncols, &cols, &vals);
332: for(col = 0; col < ncols; col++) {
333: for(locCol = 0; locCol < numCols; locCol++) {
334: if (cols[col] == ml->colPartition[level][part][locCol]) break;
335: }
336: if (locCol < numCols) {
337: ml->factors[level][part][FACT_QR][numRows*locCol+row] = vals[col];
338: }
339: }
340: MatRestoreRow(grad, ml->rowPartition[level][PART_ROW_INT][part][row], &ncols, &cols, &vals);
341: }
343: /* Do QR of B */
344: LAgeqrf_(&numRows, &numCols, ml->factors[level][part][FACT_QR],
345: &numRows, ml->factors[level][part][FACT_TAU], ml->svdWork, &ml->svdWorkLen, &ierr);
346:
348: /* Put R in B */
349: PetscMemzero(B, numCols*numCols * sizeof(PetscReal));
350: for(col = 0; col < numCols; col++)
351: for(row = 0; row <= col; row++)
352: B[col*numCols+row] = ml->factors[level][part][FACT_QR][col*numRows+row];
354: /* Do SVD of R */
355: #if defined(PETSC_MISSING_LAPACK_GESVD)
356: SETERRQ(PETSC_ERR_SUP,"GESVD - Lapack routine is unavilable.");
357: #else
358: LAgesvd_("A", "A", &numCols, &numCols, B, &numCols,
359: ml->factors[level][part][FACT_DINV], ml->factors[level][part][FACT_U], &numCols,
360: ml->factors[level][part][FACT_V], &numCols, ml->svdWork, &ml->svdWorkLen, &ierr);
361: #endif
362: } else {
363: /* Allocation */
364: PetscMalloc(numRows*numRows * sizeof(PetscReal), &ml->factors[level][part][FACT_U]);
365: PetscLogObjectMemory(pc, numRows*numRows * sizeof(PetscReal));
367: /* Extract the interior gradient for this partition */
368: PetscMemzero(B, numRows*numCols * sizeof(PetscReal));
369: for(row = 0; row < numRows; row++) {
370: MatGetRow(grad, ml->rowPartition[level][PART_ROW_INT][part][row], &ncols, &cols, &vals);
371: for(col = 0; col < ncols; col++) {
372: for(locCol = 0; locCol < numCols; locCol++) {
373: if (cols[col] == ml->colPartition[level][part][locCol]) break;
374: }
375: if (locCol < numCols) {
376: B[numRows*locCol+row] = vals[col];
377: }
378: }
379: MatRestoreRow(grad, ml->rowPartition[level][PART_ROW_INT][part][row], &ncols, &cols, &vals);
380: }
382: /* Factor the gradient */
383: #if defined(PETSC_MISSING_LAPACK_GESVD)
384: SETERRQ(PETSC_ERR_SUP,"GESVD - Lapack routine is unavilable.");
385: #else
386: LAgesvd_("A", "A", &numRows, &numCols, B, &numRows, ml->factors[level][part][FACT_DINV],
387: ml->factors[level][part][FACT_U], &numRows, ml->factors[level][part][FACT_V], &numCols,
388: ml->svdWork, &ml->svdWorkLen, &ierr);
389: #endif
390:
391: }
393: /* Invert the singular values */
394: for(col = 0; col < numCols; col++) {
395: if (ml->factors[level][part][FACT_DINV][col] > ml->zeroTol)
396: ml->factors[level][part][FACT_DINV][col] = 1.0/ml->factors[level][part][FACT_DINV][col];
397: }
398: } else {
399: /* Create null singular values for D^{-1} and the identity for V */
400: PetscMemzero(ml->factors[level][part][FACT_DINV], numCols * sizeof(PetscReal));
401: PetscMemzero(ml->factors[level][part][FACT_V], numCols*numCols * sizeof(PetscReal));
402: for(col = 0; col < numCols; col++)
403: ml->factors[level][part][FACT_V][col*numCols+col] = 1.0;
404: }
405: }
407: /* Cleanup */
408: PetscFree(B);
409: PC_MLLogEventEnd(PC_ML_ReduceFactor, pc, grad, 0, 0);
410: return(0);
411: }
413: #undef __FUNCT__
415: int PCMultiLevelExtractBoundaryGradient_Private(PC pc, int level, int numParts, int colsRemain, int *newCols,
416: PetscScalar *newVals, Mat grad)
417: {
418: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
419: int ncols;
420: int *cols;
421: PetscScalar *vals;
422: PetscScalar *temp;
423: int *rowLens;
424: int numBdRows, numResRows, rowsRemain, numCols, numVals, prevCols;
425: int part, row, newRow, col, locCol = 0;
426: int ierr;
429: PC_MLLogEventBegin(PC_ML_ReduceBdGradExtract, pc, grad, 0, 0);
430: /* Initialization */
431: numBdRows = ml->numPartitionRows[level][numParts];
432: numResRows = ml->numPartitionRows[level][numParts+1];
433: rowsRemain = numBdRows + numResRows;
435: /* Get the nonzero structure */
436: VecGetArray(ml->bdReduceVecs[level], &temp);
437: rowLens = (int *) temp;
438: for(row = 0; row < numBdRows; row++)
439: {
440: MatGetRow(grad, ml->rowPartition[level][PART_ROW_BD][0][row], &ncols, &cols, PETSC_NULL);
441: for(col = 0, numVals = 0; col < ncols; col++)
442: {
443: for(part = 0, prevCols = 0; part < numParts; part++)
444: {
445: numCols = ml->numPartitionCols[level][part];
446: for(locCol = 0; locCol < numCols; locCol++)
447: if (cols[col] == ml->colPartition[level][part][locCol])
448: break;
449: if (locCol < numCols)
450: break;
451: prevCols += numCols;
452: }
453: if (part < numParts)
454: numVals++;
455: }
456: if (numVals == 0) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Null boundary row in boundary gradient");
457: rowLens[row] = numVals;
458: MatRestoreRow(grad, ml->rowPartition[level][PART_ROW_BD][0][row], &ncols, &cols, PETSC_NULL);
459: }
460: for(row = 0; row < numResRows; row++)
461: {
462: MatGetRow(grad, ml->rowPartition[level][PART_ROW_RES][0][row], &ncols, &cols, PETSC_NULL);
463: for(col = 0, numVals = 0; col < ncols; col++)
464: {
465: for(part = 0, prevCols = 0; part < numParts; part++)
466: {
467: numCols = ml->numPartitionCols[level][part];
468: for(locCol = 0; locCol < numCols; locCol++)
469: if (cols[col] == ml->colPartition[level][part][locCol])
470: break;
471: if (locCol < numCols)
472: break;
473: prevCols += numCols;
474: }
475: if (part < numParts)
476: numVals++;
477: }
478: if (numVals == 0) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Null residual row in boundary gradient");
479: rowLens[row+numBdRows] = numVals;
480: MatRestoreRow(grad, ml->rowPartition[level][PART_ROW_RES][0][row], &ncols, &cols, PETSC_NULL);
481: }
483: /* Create the boundary gradient B_Gamma */
484: MatCreateSeqAIJ(PETSC_COMM_SELF, rowsRemain, colsRemain, 0, rowLens, &ml->grads[level]);
485: MatSetOption(ml->grads[level], MAT_NEW_NONZERO_ALLOCATION_ERR);
486: VecRestoreArray(ml->bdReduceVecs[level], &temp);
488: /* Extract boundary rows */
489: for(row = 0; row < numBdRows; row++) {
490: MatGetRow(grad, ml->rowPartition[level][PART_ROW_BD][0][row], &ncols, &cols, &vals);
491: for(col = 0, numVals = 0; col < ncols; col++)
492: {
493: for(part = 0, prevCols = 0; part < numParts; part++)
494: {
495: numCols = ml->numPartitionCols[level][part];
496: for(locCol = 0; locCol < numCols; locCol++)
497: if (cols[col] == ml->colPartition[level][part][locCol])
498: break;
499: if (locCol < numCols)
500: break;
501: prevCols += numCols;
502: }
503: if (part < numParts)
504: {
505: newCols[numVals] = locCol + prevCols;
506: newVals[numVals] = vals[col];
507: numVals++;
508: }
509: }
510: MatSetValues(ml->grads[level], 1, &row, numVals, newCols, newVals, INSERT_VALUES);
511: MatRestoreRow(grad, ml->rowPartition[level][PART_ROW_BD][0][row], &ncols, &cols, &vals);
512: }
514: /* Extract residual rows */
515: for(row = 0, newRow = numBdRows; row < numResRows; row++, newRow++)
516: {
517: MatGetRow(grad, ml->rowPartition[level][PART_ROW_RES][0][row], &ncols, &cols, &vals);
518: for(col = 0, numVals = 0; col < ncols; col++)
519: {
520: for(part = 0, prevCols = 0; part < numParts; part++)
521: {
522: numCols = ml->numPartitionCols[level][part];
523: for(locCol = 0; locCol < numCols; locCol++)
524: if (cols[col] == ml->colPartition[level][part][locCol])
525: break;
526: if (locCol < numCols)
527: break;
528: prevCols += numCols;
529: }
530: if (part < numParts)
531: {
532: newCols[numVals] = locCol + prevCols;
533: newVals[numVals] = vals[col];
534: numVals++;
535: }
536: }
537: MatSetValues(ml->grads[level], 1, &newRow, numVals, newCols, newVals, INSERT_VALUES);
538: MatRestoreRow(grad, ml->rowPartition[level][PART_ROW_RES][0][row], &ncols, &cols, &vals);
539: }
541: MatAssemblyBegin(ml->grads[level], MAT_FINAL_ASSEMBLY);
542: MatAssemblyEnd(ml->grads[level], MAT_FINAL_ASSEMBLY);
543: PC_MLLogEventEnd(PC_ML_ReduceBdGradExtract, pc, grad, 0, 0);
544: return(0);
545: }
547: #undef __FUNCT__
549: int PCMultiLevelShrinkMesh_Private(PC pc, int level, int numNewCols, int *colPartition, int numElements, int **oldMesh,
550: int *numNewElements, int ***newMesh)
551: {
552: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
553: int numCols = ml->numVertices[level];
554: int *adjCols;
555: int maxSize;
556: int col, oldCol, newCol, adjCol, count;
557: int elem;
558: int ierr;
561: PC_MLLogEventBegin(PC_ML_ReduceShrinkMesh, pc, 0, 0, 0);
562: ml->numMeshes++;
564: /* Allocation */
565: maxSize = PetscMax(PetscMin(oldMesh[MESH_OFFSETS][numCols], numNewCols*numNewCols + 1), 1);
566: PetscMalloc(NUM_MESH_DIV * sizeof(int *), newMesh);
567: PetscMalloc((numNewCols+1) * sizeof(int), &(*newMesh)[MESH_OFFSETS]);
568: PetscMalloc(maxSize * sizeof(int), &(*newMesh)[MESH_ADJ]);
569: PetscMalloc((numNewCols+1) * sizeof(int), &adjCols);
570: PetscLogObjectMemory(pc, NUM_MESH_DIV * sizeof(int *) + (numNewCols+1 + maxSize)*sizeof(int));
571: /* Shrink adjacency lists */
572: (*newMesh)[MESH_OFFSETS][0] = 0;
573: for(col = 0, count = 0; col < numNewCols; col++) {
574: PetscMemzero(adjCols, numNewCols * sizeof(int));
575: /* Put in self edge */
576: (*newMesh)[MESH_ADJ][count++] = col;
577: adjCols[col] = 1;
578: /* Get all columns in this partition */
579: for(oldCol = 0; oldCol < numCols; oldCol++) {
580: if (colPartition[oldCol] == col) {
581: for(adjCol = oldMesh[MESH_OFFSETS][oldCol]; adjCol < oldMesh[MESH_OFFSETS][oldCol+1]; adjCol++) {
582: /* Check for adjacent column */
583: newCol = colPartition[oldMesh[MESH_ADJ][adjCol]];
584: /* Null partitions have number -1 */
585: if ((newCol < 0) || (adjCols[newCol]))
586: continue;
587: (*newMesh)[MESH_ADJ][count++] = newCol;
588: adjCols[newCol] = 1;
589: }
590: }
591: }
592: (*newMesh)[MESH_OFFSETS][col+1] = count;
593: }
594: PetscFree(adjCols);
596: /* Eliminate redundant elements */
597: maxSize = PetscMax(PetscMin(numElements, PetscMax(numCols-3, 0)*2 + 1), 1);
598: PetscMalloc(maxSize*3 * sizeof(int), &(*newMesh)[MESH_ELEM]);
599: PetscLogObjectMemory(pc, maxSize*3 * sizeof(int));
600: for(elem = 0, *numNewElements = 0; elem < numElements; elem++) {
601: if ((colPartition[oldMesh[MESH_ELEM][elem*3]-1] != colPartition[oldMesh[MESH_ELEM][elem*3+1]-1]) &&
602: (colPartition[oldMesh[MESH_ELEM][elem*3+1]-1] != colPartition[oldMesh[MESH_ELEM][elem*3+2]-1]) &&
603: (colPartition[oldMesh[MESH_ELEM][elem*3+2]-1] != colPartition[oldMesh[MESH_ELEM][elem*3]-1]) &&
604: (colPartition[oldMesh[MESH_ELEM][elem*3]-1] >= 0) &&
605: (colPartition[oldMesh[MESH_ELEM][elem*3+1]-1] >= 0) &&
606: (colPartition[oldMesh[MESH_ELEM][elem*3+2]-1] >= 0))
607: {
608: (*newMesh)[MESH_ELEM][(*numNewElements)*3] = colPartition[oldMesh[MESH_ELEM][elem*3]-1] + 1;
609: (*newMesh)[MESH_ELEM][(*numNewElements)*3+1] = colPartition[oldMesh[MESH_ELEM][elem*3+1]-1] + 1;
610: (*newMesh)[MESH_ELEM][(*numNewElements)*3+2] = colPartition[oldMesh[MESH_ELEM][elem*3+2]-1] + 1;
611: (*numNewElements)++;
612: }
613: }
614: PC_MLLogEventEnd(PC_ML_ReduceShrinkMesh, pc, 0, 0, 0);
616: return(0);
617: }
619: #undef __FUNCT__
621: int PCMultiLevelRowPartitionLocalToGlobal_Private(PC pc, int level)
622: {
623: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
624: PetscScalar *bdArray;
625: int *prevIndices;
626: int *rowIndices;
627: int numParts, numIntRows, numBdRows, numResRows;
628: int part, row;
629: int ierr;
632: if (level <= 0)
633: return(0);
634: PC_MLLogEventBegin(PC_ML_ReduceBdGradRowPartLocalToGlobal, pc, 0, 0, 0);
635: VecGetArray(ml->bdReduceVecs[level-1], &bdArray);
637: /* Load unreduced rows from last level into a work vector */
638: prevIndices = (int *) bdArray;
639: numParts = ml->numPartitions[level-1];
640: numBdRows = ml->numPartitionRows[level-1][numParts];
641: PetscMemcpy(prevIndices, ml->rowPartition[level-1][PART_ROW_BD][0], numBdRows * sizeof(int));
642: numResRows = ml->numPartitionRows[level-1][numParts+1];
643: PetscMemcpy(prevIndices+numBdRows, ml->rowPartition[level-1][PART_ROW_RES][0], numResRows * sizeof(int));
644: /* Interior edges */
645: numParts = ml->numPartitions[level];
646: for(part = 0; part < numParts; part++)
647: {
648: rowIndices = ml->rowPartition[level][PART_ROW_INT][part];
649: numIntRows = ml->numPartitionRows[level][part];
650: for(row = 0; row < numIntRows; row++)
651: rowIndices[row] = prevIndices[rowIndices[row]];
652: }
653: /* Boundary edges */
654: rowIndices = ml->rowPartition[level][PART_ROW_BD][0];
655: numBdRows = ml->numPartitionRows[level][numParts];
656: for(row = 0; row < numBdRows; row++)
657: rowIndices[row] = prevIndices[rowIndices[row]];
658: /* Residual edges */
659: rowIndices = ml->rowPartition[level][PART_ROW_RES][0];
660: numResRows = ml->numPartitionRows[level][numParts+1];
661: for(row = 0; row < numResRows; row++)
662: rowIndices[row] = prevIndices[rowIndices[row]];
664: VecRestoreArray(ml->bdReduceVecs[level-1], &bdArray);
665: PC_MLLogEventEnd(PC_ML_ReduceBdGradRowPartLocalToGlobal, pc, 0, 0, 0);
667: return(0);
668: }
670: #undef __FUNCT__
672: /*@C PCMultiLevelReduce
673: This function creates a multilevel factorization of the gradient matrix
674: which may be used for preconditioning a linear solve. This is called by
675: PCSetUp() when the PC_MULTILEVEL type is selected.
677: Input Parameters:
678: . pc - The preconditioner context
680: Level: intermediate
682: .keywords multilevel
683: .seealso
684: @*/
685: int PCMultiLevelReduce(PC pc)
686: {
687: PC_Multilevel *ml;
688: PetscTruth adj;
689: Mat grad;
690: int ncols;
691: int *cols;
692: int *newCols;
693: PetscScalar *vals;
694: PetscScalar *newVals;
695: PetscScalar *colReduceArray;
696: PetscScalar *colReduceArray2;
697: int *colPartition;
698: int *rowPartition;
699: int *globalCols;
700: int *localCols;
701: int *globalPart;
702: int rowsRemain, colsRemain, newColsRemain;
703: int numParts, numNewParts, numCols, numNullCols, numVals;
704: int numBdRows, numResRows;
705: int rankDef;
706: PetscScalar val;
707: int size, sizeMax;
708: int invalidNull;
709: int level, part, newPart, row, col, col2, locCol, nullCol, vCol, loop, anyAgain;
710: #define NEW_BD_GRAD
711: #ifdef NEW_BD_GRAD
712: int *firstCol, *colPart;
713: #ifdef NEW_BD_GRAD_CHECK
714: int oldPart, oldLocCol;
715: #endif
716: #endif
717: int ierr, ierr2;
721: ml = (PC_Multilevel *) pc->data;
723: /* Create P and Z such that
725: P^T B Z = / D
726: 0 /
728: which is almost and SVD
729: */
730: if (ml->useMath == PETSC_FALSE) {
731: /* Initialization */
732: /* MatConvert(ml->locB, MATSAME, &grad); */
733: grad = ml->locB;
734: colsRemain = ml->numCols;
735: rowsRemain = ml->numRows;
736: ml->rank = 0;
737: sizeMax = PetscMax(ml->numRows, ml->numCols);
738: ml->svdWorkLen = PetscMax(3*PetscMin(ml->numRows,ml->numCols) + PetscMax(ml->numRows,ml->numCols),
739: 5*PetscMin(ml->numRows,ml->numCols) - 4);
741: /* Allocation */
742: PetscMalloc(ml->maxLevels * sizeof(int), &ml->numPartitions);
743: PetscMalloc(ml->maxLevels * sizeof(int *), &ml->numPartitionCols);
744: PetscMalloc(ml->maxLevels * sizeof(int **), &ml->colPartition);
745: PetscMalloc(ml->maxLevels * sizeof(int *), &ml->numPartitionRows);
746: PetscMalloc(ml->maxLevels * sizeof(int ***), &ml->rowPartition);
747: PetscMalloc(ml->maxLevels * sizeof(double ***), &ml->factors);
748: PetscMalloc(ml->maxLevels * sizeof(Mat), &ml->grads);
749: PetscMalloc(ml->maxLevels * sizeof(Vec), &ml->bdReduceVecs);
750: PetscMalloc(ml->maxLevels * sizeof(Vec), &ml->colReduceVecs);
751: PetscMalloc(ml->maxLevels * sizeof(Vec), &ml->colReduceVecs2);
752: PetscMalloc(ml->svdWorkLen * sizeof(double), &ml->svdWork);
753: PetscMalloc(ml->numCols * sizeof(int), &ml->range);
754: PetscLogObjectMemory(pc, ml->maxLevels*(sizeof(int) + sizeof(int *) + sizeof(int **) + sizeof(int *) + sizeof(int ***) + sizeof(double ***) +
755: sizeof(Mat) + sizeof(Vec)*3) + ml->svdWorkLen * sizeof(double) + ml->numCols * sizeof(int));
756: PetscMalloc(ml->numCols * sizeof(int), &colPartition);
757: PetscMalloc(ml->numCols * sizeof(int), &globalCols);
758: PetscMalloc(sizeMax * sizeof(int), &rowPartition);
760: for(col = 0; col < ml->numCols; col++) {
761: ml->range[col] = -1;
762: globalCols[col] = col;
763: }
765: /* Loop until the graph has ml->QRthresh nodes or all rows have been reduced */
766: for(level = 0, anyAgain = 0; (colsRemain > ml->QRthresh) && (rowsRemain > 0); level++) {
767: if (level >= ml->maxLevels) SETERRQ(PETSC_ERR_ARG_SIZ, "Exceeded maximum level");
769: loop = 2;
770: while(loop) {
771: /* Partition the columns (nodes): Give a list with the partition number of each column */
772: PCMultiLevelPartitionMesh_Private(pc, level, ml->meshes[level], &numParts, colPartition);
773: ml->numPartitions[level] = numParts;
775: /* Get the global columns in each partition */
776: PC_MLLogEventBegin(PC_ML_ReducePartitionRowCol, pc, 0, 0, 0);
777: PCMultiLevelPartitionCols_Private(pc, level, numParts, colsRemain, colPartition, globalCols);
778:
780: /* Partition rows (edges): Give the local rows in each partition, then boundary and residual rows */
781: PCMultiLevelPartitionRows_Private(pc, level, numParts, colsRemain, rowsRemain, colPartition,
782: rowPartition, globalCols, grad);
783:
784: PC_MLLogEventEnd(PC_ML_ReducePartitionRowCol, pc, 0, 0, 0);
786: /* Factor the interior of each partition */
787: loop = PCMultiLevelFactorPartitions_Private(pc, level, numParts, grad);
789: if (loop == 2) {
790: /* Free memory from bad partition */
791: for(part = 0; part < numParts; part++) {
792: if (ml->numPartitionCols[level][part] > 0)
793: {PetscFree(ml->colPartition[level][part]); }
794: if (ml->numPartitionRows[level][part] > 0)
795: {PetscFree(ml->rowPartition[level][PART_ROW_INT][part]); }
796: }
797: if (ml->numPartitionRows[level][ml->numPartitions[level]] > 0)
798: {PetscFree(ml->rowPartition[level][PART_ROW_BD][0]); }
799: if (ml->numPartitionRows[level][ml->numPartitions[level]+1] > 0)
800: {PetscFree(ml->rowPartition[level][PART_ROW_RES][0]); }
801: PetscFree(ml->rowPartition[level][PART_ROW_INT]);
802: PetscFree(ml->rowPartition[level][PART_ROW_BD]);
803: PetscFree(ml->rowPartition[level][PART_ROW_RES]);
804: PetscFree(ml->colPartition[level]);
805: PetscFree(ml->rowPartition[level]);
806: PetscFree(ml->numPartitionCols[level]);
807: PetscFree(ml->numPartitionRows[level]);
808: /* Increase partition size */
809: PetscLogInfo(pc, "PCMultiLevelReduce: Increasing partition size from %d to %dn", ml->partSize, ml->partSize*2);
810: if (ml->partSize > ml->numCols) SETERRQ(PETSC_ERR_PLIB, "Partition has grown too large.");
811: ml->partSize *= 2;
812: } else if (loop != 0) {
813: CHKERRQ(loop);
814: }
815: }
817: /* QR factor boundary matrices */
819: /* Create the boundary gradient B_Gamma and workspace for applying it */
820: PC_MLLogEventBegin(PC_ML_ReduceBdGrad, pc, 0, 0, 0);
821: numBdRows = ml->numPartitionRows[level][numParts];
822: numResRows = ml->numPartitionRows[level][numParts+1];
823: rowsRemain = numBdRows + numResRows;
824: VecCreateSeq(PETSC_COMM_SELF, rowsRemain, &ml->bdReduceVecs[level]);
825: VecCreateSeq(PETSC_COMM_SELF, colsRemain, &ml->colReduceVecs[level]);
826: VecDuplicate(ml->colReduceVecs[level], &ml->colReduceVecs2[level]);
828: /* Workspace for matrix reduction */
829: VecGetArray(ml->colReduceVecs[0], &colReduceArray);
830: VecGetArray(ml->colReduceVecs2[0], &colReduceArray2);
831: newCols = (int *) colReduceArray;
832: newVals = ml->svdWork;
834: /* Extract the boundary gradient B_Gamma in CSR format */
835: PCMultiLevelExtractBoundaryGradient_Private(pc, level, numParts, colsRemain, newCols, newVals, grad);
836:
838: /* Convert the local (per level) numbering to the global (original) numbering */
839: PCMultiLevelRowPartitionLocalToGlobal_Private(pc, level);
841: /* Determine the columns active at the next level and fix up colPartition[] */
842: globalPart = (int *) colReduceArray2;
843: localCols = rowPartition;
844: for(part = 0; part < numParts; part++) globalCols[part] = -1;
845: for(part = 0, newPart = 0, newColsRemain = 0, numNewParts = numParts; part < numParts; part++, newPart++) {
846: for(col = 0, rankDef = 0; col < ml->numPartitionCols[level][part]; col++) {
847: if (ml->factors[level][part][FACT_DINV][col] < ml->zeroTol) {
848: newColsRemain++;
849: rankDef++;
851: if (rankDef == 1) {
852: globalCols[newPart] = ml->colPartition[level][part][col];
853: globalPart[newPart] = part;
854: localCols[newPart] = col;
855: } else {
856: globalCols[numNewParts] = ml->colPartition[level][part][col];
857: globalPart[numNewParts] = part;
858: localCols[numNewParts] = col;
860: /* Put the column in its own partition */
861: numCols = ml->numPartitionCols[level][part] - (col+1);
862: for(col2 = colsRemain-1; numCols >= 0; col2--) {
863: if (colPartition[col2] == newPart) {
864: numCols--;
865: if (numCols < 0)
866: colPartition[col2] = numNewParts++;
867: }
868: }
869: }
870: } else {
871: ml->rank++;
872: #ifdef PETSC_USE_BOPT_g
873: if (ml->range[ml->colPartition[level][part][col]] != -1) {
874: PetscLogInfo(pc, "Level %d partition %d col %d: column %d has already been reducedn",
875: level, part, col, ml->colPartition[level][part][col]);
876: SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid reduction");
877: }
878: if ((ml->rowPartition[level][PART_ROW_INT][part][col] < 0) ||
879: (ml->rowPartition[level][PART_ROW_INT][part][col] > ml->numRows)) {
880: SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid row index in range");
881: }
882: #endif
883: ml->range[ml->colPartition[level][part][col]] = ml->rowPartition[level][PART_ROW_INT][part][col];
884: }
885: }
886: #ifdef PETSC_USE_BOPT_g
887: PetscLogInfo(pc, "Lvl: %d Partition: %d has rank deficiency %dn", level, part, rankDef);
888: #endif
889: if (rankDef == 0) {
890: /* Eliminate this partition */
891: for(col = 0; col < colsRemain; col++) {
892: if (colPartition[col] == newPart) {
893: colPartition[col] = -1;
894: } else if (colPartition[col] > newPart) {
895: colPartition[col]--;
896: }
897: }
898: /* Shift everything down */
899: size = PetscMax(numNewParts, numParts) - part;
900: PetscMemmove(&globalCols[newPart], &globalCols[newPart+1], size * sizeof(int));
901: PetscMemmove(&globalPart[newPart], &globalPart[newPart+1], size * sizeof(int));
902: PetscMemmove(&localCols[newPart], &localCols[newPart+1], size * sizeof(int));
903: newPart--;
904: numNewParts--;
905: }
906: }
908: #ifdef PETSC_USE_BOPT_g
909: for(col = 0, row = 0; col < ml->numCols; col++) {
910: if (ml->range[col] == -1) {
911: row++;
912: } else if (ml->range[col] < 0) {
913: SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid range space");
914: }
915: }
916: if (row != newColsRemain) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid range space");
918: for(vCol = 0; vCol < newColsRemain; vCol++) {
919: if (globalCols[vCol] < 0) {
920: PetscPrintf(PETSC_COMM_SELF, "Invalid column in reduced boundary gradientn");
921: PetscPrintf(PETSC_COMM_SELF, "Lvl: %d newRows: %d newCols: %dn", level, rowsRemain, newColsRemain);
922: for(col = 0; col < newColsRemain; col++) {
923: PetscPrintf(PETSC_COMM_SELF, " globalCols[%d]: %dn", col, globalCols[col]);
924: }
925: SETERRQ(PETSC_ERR_ARG_CORRUPT, "Negative column");
926: }
927: }
928: #endif
929: PC_MLLogEventEnd(PC_ML_ReduceBdGrad, pc, 0, 0, 0);
931: /* Transform B_Gamma to B_Gamma V (I - D^{-1} D) */
932: PC_MLLogEventBegin(PC_ML_CreateBdGrad, pc, 0, 0, 0);
933: if (level > 0) {
934: MatDestroy(grad);
935: }
936: MatCreateSeqAIJ(PETSC_COMM_SELF, rowsRemain, ml->numCols, newColsRemain, PETSC_NULL, &grad);
937: if (ierr) {
938: PetscPrintf(PETSC_COMM_SELF, "Interface Gradient failure:");
939: PetscPrintf(PETSC_COMM_SELF, " level: %d", level);
940: PetscPrintf(PETSC_COMM_SELF, " rows: %d cols: %dn", rowsRemain, ml->numCols);
941: PetscPrintf(PETSC_COMM_SELF, " nonzeroes per row: %dn", newColsRemain);
942:
943: }
944: #ifdef NEW_BD_GRAD
945: PetscMalloc((numParts+1) * sizeof(int), &firstCol);
946: PetscMalloc(ml->numCols * sizeof(int), &colPart);
947: for(part = 1, firstCol[0] = 0; part <= numParts; part++) {
948: firstCol[part] = firstCol[part-1] + ml->numPartitionCols[level][part-1];
949: for(col = 0; col < ml->numPartitionCols[level][part-1]; col++) {
950: colPart[col+firstCol[part-1]] = part-1;
951: }
952: }
953: #endif
954: for(row = 0; row < rowsRemain; row++) {
955: MatGetRow(ml->grads[level], row, &ncols, &cols, &vals);
956: /* Sparse VecDot with the active columns of V */
957: for(vCol = 0, numVals = 0; vCol < newColsRemain; vCol++) {
958: /* Calculate the dot product */
959: for(col = 0, adj = PETSC_FALSE, val = 0.0; col < ncols; col++) {
960: /* Get the local column number */
961: #ifndef NEW_BD_GRAD
962: locCol = cols[col];
963: for(part = 0; part < PetscMin(numParts, globalPart[vCol]+2); part++) {
964: numCols = ml->numPartitionCols[level][part];
965: if (locCol >= numCols) {
966: locCol -= numCols;
967: continue;
968: }
969: break;
970: }
971: #else
972: part = colPart[cols[col]];
973: locCol = cols[col] - firstCol[part];
974: numCols = ml->numPartitionCols[level][part];
975: #ifdef NEW_BD_GRAD_CHECK
976: /***** CHECK THE RESULT ***/
977: oldLocCol = cols[col];
978: for(oldPart = 0; oldPart < PetscMin(numParts, globalPart[vCol]+2); oldPart++) {
979: numCols = ml->numPartitionCols[level][oldPart];
980: if (oldLocCol >= numCols) {
981: oldLocCol -= numCols;
982: continue;
983: }
984: break;
985: }
986: if ((oldPart == globalPart[vCol]) && (oldPart != part)) SETERRQ2(PETSC_ERR_PLIB, "Invalid partition %d should be %d", part, oldPart);
987: if ((oldPart == globalPart[vCol]) && (oldLocCol != locCol)) SETERRQ2(PETSC_ERR_PLIB, "Invalid local column %d should be %d", locCol, oldLocCol);
988: #endif
989: #endif
990: /* Multiply by the correct element of V (I - D^{-1} D) */
991: if (part == globalPart[vCol]) {
992: val += vals[col]*ml->factors[level][part][FACT_V][locCol*numCols+localCols[vCol]];
993: adj = PETSC_TRUE;
994: }
995: }
996: /* Store the dot product in B_Gamma V (I - D^{-1} D) */
997: if (adj == PETSC_TRUE) {
998: newCols[numVals] = globalCols[vCol];
999: newVals[numVals] = val;
1000: numVals++;
1001: }
1002: }
1003: /* This is an overestimate */
1004: PetscLogFlops(numVals*ncols);
1005: #if 0
1006: #ifdef PETSC_USE_BOPT_g
1007: if (numVals == 0) PetscLogInfo(pc, "Null row %d in reduced boundary gradient matrixn", row);
1008: #endif
1009: #endif
1010: MatSetValues(grad, 1, &row, numVals, newCols, newVals, INSERT_VALUES);
1011: MatRestoreRow(ml->grads[level], row, &ncols, &cols, &vals);
1012: }
1013: MatAssemblyBegin(grad, MAT_FINAL_ASSEMBLY);
1014: MatAssemblyEnd(grad, MAT_FINAL_ASSEMBLY);
1015: #ifdef NEW_BD_GRAD
1016: PetscFree(firstCol);
1017: PetscFree(colPart);
1018: #endif
1019: PC_MLLogEventEnd(PC_ML_CreateBdGrad, pc, 0, 0, 0);
1021: /* Construct the coarse graph */
1022: ml->numVertices[level+1] = newColsRemain;
1023: PCMultiLevelShrinkMesh_Private(pc, level, newColsRemain, colPartition, ml->numElements[level],
1024: ml->meshes[level], &ml->numElements[level+1], &ml->meshes[level+1]);
1025:
1027: /* Update variables */
1028: colsRemain = newColsRemain;
1029: VecRestoreArray(ml->colReduceVecs[0], &colReduceArray);
1030: VecRestoreArray(ml->colReduceVecs2[0], &colReduceArray2);
1032: /* Visualize factorization */
1033: PCMultiLevelReduceView(pc, level, colsRemain, colPartition,
1034: ((colsRemain > ml->QRthresh) && (rowsRemain > 0)), &anyAgain);
1035:
1036: }
1037: while(anyAgain) {
1038: PCMultiLevelReduceView(pc, -level, 0, colPartition, 0, &anyAgain);
1039: }
1040: ml->numLevels = level;
1042: /* Cleanup */
1043: PetscFree(colPartition);
1044: PetscFree(rowPartition);
1045: PetscFree(globalCols);
1046: if (ml->numLevels > 0) {
1047: MatDestroy(grad);
1048: }
1050: /* Workspace allocation */
1051: ml->interiorWorkLen = 1;
1052: for(level = 0; level < ml->numLevels; level++)
1053: for(part = 0; part < ml->numPartitions[level]; part++)
1054: ml->interiorWorkLen = PetscMax(ml->interiorWorkLen, ml->numPartitionRows[level][part]);
1055: PetscMalloc(ml->interiorWorkLen * sizeof(double), &ml->interiorWork);
1056: PetscMalloc(ml->interiorWorkLen * sizeof(double), &ml->interiorWork2);
1057: PetscLogObjectMemory(pc, ml->interiorWorkLen*2 * sizeof(double));
1059: /* Get the null space rows and compress the range */
1060: ml->globalRank = ml->rank;
1061: ml->numLocNullCols = ml->numCols - ml->rank;
1062: PetscMalloc(PetscMax(ml->numLocNullCols, 1) * sizeof(int), &ml->nullCols);
1063: PetscLogObjectMemory(pc, PetscMax(ml->numLocNullCols, 1) * sizeof(int));
1064: ml->nullCols[0] = -1;
1065: for(col = 0, nullCol = 0, numNullCols = 0, ierr2 = 0; nullCol < ml->numCols; col++, nullCol++) {
1066: if (ml->range[col] < 0) {
1067: if (numNullCols == ml->numLocNullCols) {
1068: int rank;
1070: MPI_Comm_rank(pc->comm, &rank);
1071: PetscPrintf(PETSC_COMM_SELF, "[%d]Range %d, Null %d:n", rank, ml->rank, ml->numLocNullCols);
1072: for(col = 0; col < ml->numCols; col++) PetscPrintf(PETSC_COMM_SELF, " %d", ml->range[col]);
1073: PetscPrintf(PETSC_COMM_SELF, "n");
1074: ierr2 = 1;
1075: break;
1076: }
1077: ml->nullCols[numNullCols++] = nullCol;
1078: PetscMemmove(&ml->range[col], &ml->range[col+1], (ml->numCols - (col+1)) * sizeof(int));
1079: col--;
1080: }
1081: }
1082: MPI_Allreduce(&ierr2, &invalidNull, 1, MPI_INT, MPI_LOR, pc->comm);
1083: if (invalidNull) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid null space");
1084: if (numNullCols != ml->numLocNullCols) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid number of null space columns");
1085: if (ml->numLocNullCols + ml->rank != ml->numCols) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid space decomposition");
1087: #ifdef PETSC_USE_BOPT_g
1088: PCValidQ_Multilevel(pc);
1089: #endif
1090: } else {
1091: Mesh mesh;
1093: /* Cleanup */
1094: GridGetMesh(ml->grid, &mesh);
1095: MeshDestroyLocalCSR(mesh, ml->meshes[0][MESH_OFFSETS], ml->meshes[0][MESH_ADJ], ml->vertices);
1096: #ifdef PETSC_HAVE_MATHEMATICA
1097: ViewerMathematicaReduce(ml->mathViewer, pc, ml->QRthresh);
1098: PetscFree(ml->meshes[0][MESH_ELEM]);
1099: PetscFree(ml->meshes[0]);
1100: PetscFree(ml->meshes);
1101: PetscFree(ml->numElements);
1102: PetscFree(ml->numVertices);
1103: /* PCConvert_Multilevel(pc); */
1104: #else
1105: SETERRQ(PETSC_ERR_SUP, " ");
1106: #endif
1107: }
1108: return(0);
1109: }
1111: #undef __FUNCT__
1113: int PCMultiLevelFilterVertexCoords_Tutte(PC pc, int numCols, double *vertices, int **mesh)
1114: {
1115: PetscRandom gen;
1116: double residual = 1.0;
1117: double tol = 1.0e-8;
1118: int iter = 0;
1119: int maxIter = 10000;
1120: int xminNode, xmaxNode, yminNode, ymaxNode;
1121: double xmin, xmax, ymin, ymax;
1122: PetscReal x, y;
1123: int node, neighbor;
1124: int ierr;
1127: if (numCols < 3) return(0);
1128: /* Fix boundary */
1129: PetscRandomCreate(PETSC_COMM_SELF, RANDOM_DEFAULT, &gen);
1130: PetscRandomSetInterval(gen, 0, numCols);
1131: xmin = xmax = vertices[0];
1132: ymin = ymax = vertices[1];
1133: xminNode = xmaxNode = yminNode = ymaxNode = 0;
1134: for(node = 0; node < numCols; node++) {
1135: if (vertices[node*2] < xmin) {
1136: xmin = vertices[node*2];
1137: xminNode = node;
1138: }
1139: if (vertices[node*2] > xmax) {
1140: xmax = vertices[node*2];
1141: xmaxNode = node;
1142: }
1143: if (vertices[node*2+1] < ymin) {
1144: ymin = vertices[node*2+1];
1145: yminNode = node;
1146: }
1147: if (vertices[node*2+1] > ymax) {
1148: ymax = vertices[node*2+1];
1149: ymaxNode = node;
1150: }
1151: }
1152: if ((xmaxNode == yminNode) && (xminNode == ymaxNode)) {
1153: while((xmaxNode == yminNode) || (xmaxNode == xminNode)) {
1154: PetscRandomGetValue(gen, &x);
1155: xmaxNode = (int) floor(x);
1156: }
1157: }
1158: if ((xmaxNode == ymaxNode) && (xminNode == yminNode)) {
1159: while((xmaxNode == ymaxNode) || (xmaxNode == xminNode)) {
1160: PetscRandomGetValue(gen, &x);
1161: xmaxNode = (int) floor(x);
1162: }
1163: }
1164: PetscRandomDestroy(gen);
1165: /* First pass: Jacobi Relaxation */
1166: while((residual > tol) && (iter < maxIter)) {
1167: for(node = 0; node < numCols; node++) {
1168: if ((node == xminNode) || (node == xmaxNode) || (node == yminNode) || (node == ymaxNode)) continue;
1169: vertices[node*2] = 0.0;
1170: vertices[node*2+1] = 0.0;
1171: for(neighbor = mesh[MESH_OFFSETS][node]; neighbor < mesh[MESH_OFFSETS][node+1]; neighbor++) {
1172: if (mesh[MESH_ADJ][neighbor] != node) {
1173: vertices[node*2] += vertices[mesh[MESH_ADJ][neighbor]*2];
1174: vertices[node*2+1] += vertices[mesh[MESH_ADJ][neighbor]*2+1];
1175: }
1176: }
1177: /* Remember to get rid of self edge */
1178: vertices[node*2] /= mesh[MESH_OFFSETS][node+1] - mesh[MESH_OFFSETS][node] - 1;
1179: vertices[node*2+1] /= mesh[MESH_OFFSETS][node+1] - mesh[MESH_OFFSETS][node] - 1;
1180: }
1181: for(node = 0, residual = 0.0; node < numCols; node++) {
1182: if ((node == xminNode) || (node == xmaxNode) || (node == yminNode) || (node == ymaxNode)) continue;
1183: x = 0.0;
1184: y = 0.0;
1185: for(neighbor = mesh[MESH_OFFSETS][node]; neighbor < mesh[MESH_OFFSETS][node+1]; neighbor++) {
1186: if (mesh[MESH_ADJ][neighbor] != node) {
1187: x += vertices[mesh[MESH_ADJ][neighbor]*2];
1188: y += vertices[mesh[MESH_ADJ][neighbor]*2+1];
1189: }
1190: }
1191: x /= mesh[MESH_OFFSETS][node+1] - mesh[MESH_OFFSETS][node] - 1;
1192: y /= mesh[MESH_OFFSETS][node+1] - mesh[MESH_OFFSETS][node] - 1;
1193: residual += PetscAbsReal(x - vertices[node*2]) + PetscAbsReal(y - vertices[node*2+1]);
1194: }
1195: residual /= numCols;
1196: iter++;
1197: }
1198: return(0);
1199: }
1201: #undef __FUNCT__
1203: int PCMultiLevelCalcVertexCoords_Private(PC pc, int numCols, int numNewCols, int *colPartition)
1204: {
1205: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
1206: double *vertices;
1207: int *counts;
1208: int node;
1209: int ierr;
1212: if (numNewCols > 0) {
1213: PetscMalloc(numNewCols*2 * sizeof(double), &vertices);
1214: PetscMalloc(numNewCols * sizeof(int), &counts);
1215: PetscMemzero(vertices, numNewCols*2 * sizeof(double));
1216: PetscMemzero(counts, numNewCols * sizeof(int));
1217: for(node = 0; node < numCols; node++) {
1218: if (colPartition[node] >= 0) {
1219: vertices[colPartition[node]*2] += ml->vertices[node*2];
1220: vertices[colPartition[node]*2+1] += ml->vertices[node*2+1];
1221: counts[colPartition[node]]++;
1222: }
1223: }
1224: for(node = 0; node < numNewCols; node++) {
1225: vertices[node*2] /= counts[node];
1226: vertices[node*2+1] /= counts[node];
1227: }
1228: PetscFree(ml->vertices);
1229: PetscFree(counts);
1230: ml->vertices = vertices;
1231: }
1232: return(0);
1233: }
1235: #undef __FUNCT__
1237: int PCMultiLevelView_Mesh(PC pc, int numCols, int numElements, int **mesh)
1238: {
1239: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
1240: PetscViewer viewer = ml->factorizationViewer;
1241: Mesh smallMesh;
1242: int *faces;
1243: int face;
1244: int ierr;
1247: if (numElements > 0) {
1248: PetscMalloc(numElements*3 * sizeof(int), &faces);
1249: PetscMemcpy(faces, mesh[MESH_ELEM], numElements*3 * sizeof(int));
1250: for(face = 0; face < numElements*3; face++) {
1251: faces[face]--;
1252: }
1253: }
1254: MeshCreateTriangular2DCSR(pc->comm, numCols, numElements, ml->vertices, mesh[MESH_OFFSETS], mesh[MESH_ADJ], faces, &smallMesh);
1255:
1256: if (numElements > 0) {
1257: PetscFree(faces);
1258: }
1259: MeshView(smallMesh, viewer);
1260: MeshDestroy(smallMesh);
1261: return(0);
1262: }
1264: #undef __FUNCT__
1266: int PCMultiLevelView_Private(PC pc, PetscViewer viewer, int level, int numCols, int numElements, int **mesh)
1267: {
1268: PetscDraw draw;
1269: char title[256];
1270: int maxLevel;
1271: int ierr;
1274: PetscViewerDrawClear(viewer);
1275: PetscViewerDrawGetDraw(viewer, 0, &draw);
1276: MPI_Allreduce(&level, &maxLevel, 1, MPI_INT, MPI_MAX, pc->comm);
1277: sprintf(title, "ML Mesh, Level %d", maxLevel);
1278: PetscDrawSetTitle(draw, title);
1279: PCMultiLevelView_Mesh(pc, numCols, numElements, mesh);
1280: return(0);
1281: }
1283: #undef __FUNCT__
1285: int PCMultiLevelReduceView(PC pc, int level, int numNewCols, int *colPartition, int again, int *anyAgain)
1286: {
1287: PC_Multilevel *ml = (PC_Multilevel *) pc->data;
1288: PetscViewer viewer = ml->factorizationViewer;
1289: int **newMesh;
1290: int numCols, numNewElements;
1291: int ierr;
1294: if (viewer != PETSC_NULL) {
1295: if (level < 0) {
1296: PCMultiLevelView_Private(pc, viewer, -1, numNewCols, 0, ml->meshes[-level]);
1297: } else {
1298: numCols = ml->numVertices[level];
1299: numNewElements = ml->numElements[level+1];
1300: newMesh = ml->meshes[level+1];
1301: PCMultiLevelCalcVertexCoords_Private(pc, numCols, numNewCols, colPartition);
1302: PCMultiLevelFilterVertexCoords_Tutte(pc, numNewCols, ml->vertices, newMesh);
1303: PCMultiLevelView_Private(pc, viewer, level+1, numNewCols, numNewElements, newMesh);
1304: }
1305: PetscViewerFlush(viewer);
1306: MPI_Allreduce(&again, anyAgain, 1, MPI_INT, MPI_LOR, pc->comm);
1307: }
1308: return(0);
1309: }