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: }