Actual source code: mlApply.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: mlApply.c,v 1.5 2000/01/10 03:20:33 knepley Exp $";
  3: #endif
  4: /*
  5:    Defines the application of the multilevel preconditioner
  6: */
 7:  #include src/sles/pc/pcimpl.h
 8:  #include ml.h

 10: #undef  __FUNCT__
 12: int DQMV(char *trans, int numRows, PetscReal *Q, int ldQ, PetscReal *Tau, int numCols, PetscReal *x, PetscReal *y)
 13: {
 14:   /* Matvec with orthogonal matrix Q (elementary reflections)
 15:      Assumes LAPACK representation : Q is lower triangular with diagonal entries = 1 (which are not stored)
 16:   */
 17:   PetscTruth isnormal, istrans;
 18:   PetscReal  dot;
 19:   int        start, end, step;
 20:   int        i, j;
 21:   int        ierr;

 24:   if (numCols > numRows) SETERRQ(PETSC_ERR_ARG_WRONG, "Number of reflectors cannot exceed the size of Q");

 26:   PetscMemcpy(y, x, numRows * sizeof(double));
 27:   PetscStrcasecmp(trans, "N", &isnormal);
 28:   PetscStrcasecmp(trans, "T", &istrans);
 29:   if (isnormal == PETSC_TRUE) {
 30:     start = numCols-1;
 31:     end   = -1;
 32:     step  = -1;
 33:   } else if (istrans == PETSC_TRUE) {
 34:     start = 0;
 35:     end   = numCols;
 36:     step  = 1;
 37:   } else {
 38:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Application type must be 'N' or 'T', not %s", trans);
 39:   }

 41:   for(j = start; j != end; j += step)
 42:   {
 43:     if (Tau[j] != 0.0)
 44:     {
 45:       /* dot = v^T_j y */
 46:       dot = y[j];
 47:       for(i = j+1; i < numRows; i++)
 48:         dot  += Q[j*ldQ+i]*y[i];
 49:       /* dot = tau v^T_j y */
 50:       dot  *= Tau[j];
 51:       /* y = (I - tau v_j v^T_j) y */
 52:       y[j] -= dot;
 53:       for(i = j+1; i < numRows; i++)
 54:         y[i] -= Q[j*ldQ+i]*dot;
 55:     }
 56:   }
 57:   PetscLogFlops(numCols*(2 + numRows*2));
 58:   return(0);
 59: }

 61: #undef  __FUNCT__
 63: /*@C PCMultiLevelApplyGradient
 64:         This function applies the gradient to a vector.

 66:   Input Parameters:
 67: + pc - The preconditioner context
 68: - x  - The input vector

 70:   Output Parameter:
 71: . y  - The output vector

 73:   Level: intermediate

 75: .keywords multilevel
 76: .seealso PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans, PCMultiLevelApplyV,
 77:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
 78: @*/
 79: int PCMultiLevelApplyGradient(PC pc, GVec x, GVec y)
 80: {
 81:   PC_Multilevel *ml;
 82:   int            size, rows, cols;
 83:   int            ierr;

 89:   if (pc->setupcalled < 2) {
 90:     PCSetUp(pc);
 91:   }

 93:   ml = (PC_Multilevel *) pc->data;
 94: #ifdef PETSC_USE_BOPT_g
 95:   if (ml->useMath == PETSC_FALSE) {
 96:     PCValidQ_Multilevel(pc);
 97:   }
 98: #endif
 99:   VecGetSize(y, &size);
100:   if (ml->B != PETSC_NULL) {
101:     MatGetSize(ml->B, &rows, &cols);
102:   }
103:   if ((ml->B != PETSC_NULL) && (size == rows)) {
104:     MatMult(ml->B, x, y);
105:   } else {
106:     GVecEvaluateOperatorGalerkinRectangular(y, x, ml->sOrder, ml->sLocOrder, ml->tOrder, ml->tLocOrder,
107:                                                    ml->gradOp, ml->gradAlpha, PETSC_NULL);
108: 
109:   }
110:   return(0);
111: }

113: #undef  __FUNCT__
115: /*@C PCMultiLevelApplyGradientTrans
116:         This function applies the transpose of the
117:   gradient to a vector.

119:   Input Parameters:
120: + pc - The preconditioner context
121: - x  - The input vector

123:   Output Parameter:
124: . y  - The output vector

126:   Level: intermediate

128: .keywords multilevel
129: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyP, PCMultiLevelApplyPTrans, PCMultiLevelApplyV,
130:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
131: @*/
132: int PCMultiLevelApplyGradientTrans(PC pc, GVec x, GVec y)
133: {
134:   PC_Multilevel *ml;
135:   int            size, rows, cols;
136:   int            ierr;

142:   if (pc->setupcalled < 2) {
143:     PCSetUp(pc);
144:   }

146:   ml = (PC_Multilevel *) pc->data;
147: #ifdef PETSC_USE_BOPT_g
148:   if (ml->useMath == PETSC_FALSE) {
149:     PCValidQ_Multilevel(pc);
150:   }
151: #endif
152:   VecGetSize(x, &size);
153:   if (ml->B != PETSC_NULL) {
154:     MatGetSize(ml->B, &rows, &cols);
155:   }
156:   if ((ml->B != PETSC_NULL) && (size == rows)) {
157:     MatMultTranspose(ml->B, x, y);
158:   } else {
159:     GVecEvaluateOperatorGalerkinRectangular(y, x, ml->tOrder, ml->tLocOrder, ml->sOrder, ml->sLocOrder,
160:                                                    ml->divOp, ml->gradAlpha, PETSC_NULL);
161: 
162:   }
163:   return(0);
164: }

166: /* This is for just applying the interior portions of D^{-1} */
167: int PCMultiLevelApplyDInvLoc_Private(PC pc, GVec x, GVec y)
168: {
169: #ifdef HAVE_MATHEMATICA
170:   MLINK          link;
171: #endif
172:   PC_Multilevel *ml = (PC_Multilevel *) pc->data;
173:   double         zeroTol = ml->zeroTol;
174:   PetscScalar   *yArray;
175:   PetscReal     *invSingVals;
176:   int           *colIndices;
177:   int            level, part, col;
178:   int            ierr;

181:   /* Apply D^{-1} which is block diagonal, so we just split the vector and apply each local matrix */
182:   if (ml->useMath == PETSC_FALSE)
183:   {
184: #ifdef PETSC_USE_BOPT_g
185:     PCValidQ_Multilevel(pc);
186: #endif
187:     if (x != y)
188:       {VecCopy(x, y);                                                                              }
189:     /* Apply D^{-1} for each level */
190:     for(level = 0; level < ml->numLevels; level++)
191:     {
192:       /* Apply D^{-1} for each partition */
193:       for(part = 0; part < ml->numPartitions[level]; part++)
194:       {
195:         colIndices  = ml->colPartition[level][part];
196:         invSingVals = ml->factors[level][part][FACT_DINV];
197:         /* Here, null singular values are replaced by 1 instead of zero since these columns
198:            are carried to the next level */
199:         VecGetArray(y, &yArray);
200:         for(col = 0; col < ml->numPartitionCols[level][part]; col++)
201:           if (invSingVals[col] > zeroTol)
202:             yArray[colIndices[col]] *= invSingVals[col];
203:         VecRestoreArray(y, &yArray);
204:         PetscLogFlops(ml->numPartitionCols[level][part]);
205:       }
206:     }
207:   }
208:   else
209:   {
210: #ifdef HAVE_MATHEMATICA
211:     /* The link to Mathematica */
212:     PetscViewerMathematicaGetLink(ml->mathViewer, &link);

214:     /* vec1 = input vector */
215:     PetscViewerMathematicaSetName(ml->mathViewer, "vec1");
216:     VecView(x, ml->mathViewer);

218:     /* vec2 = DInverseApply[mattML,vec] */
219:     MLPutFunction(link, "EvaluatePacket", 1);
220:       MLPutFunction(link, "Set", 2);
221:         MLPutSymbol(link, "vec2");
222:         MLPutFunction(link, "DInverseApply", 2);
223:           MLPutSymbol(link, "mattML");
224:           MLPutSymbol(link, "vec1");
225:     MLEndPacket(link);
226:     /* Skip packets until ReturnPacket */
227:     PetscViewerMathematicaSkipPackets(ml->mathViewer, RETURNPKT);
228:     /* Skip ReturnPacket */
229:     MLNewPacket(link);

231:     /* y = vec2 */
232:     PetscViewerMathematicaSetName(ml->mathViewer, "vec2");
233:     PetscViewerMathematicaGetVector(ml->mathViewer, y);
234:     PetscViewerMathematicaClearName(ml->mathViewer);
235: #endif
236:   }
237:   return(0);
238: }

240: #undef  __FUNCT__
242: /*@C PCMultiLevelApplyDInv
243:         This function applies the inverse of D to a vector.

245:   Input Parameters:
246: + pc - The preconditioner context
247: - x  - The input vector

249:   Output Parameter:
250: . y  - The output vector

252:   Level: intermediate

254: .keywords multilevel
255: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
256:          PCMultiLevelApplyV, PCMultiLevelApplyVTrans, PCMultiLevelApplyDInvTrans
257: @*/
258: int PCMultiLevelApplyDInv(PC pc, GVec x, GVec y)
259: {
260:   PC_Multilevel *ml;
261:   PetscScalar   *rhsArray;
262:   int            numProcs, rank;
263: #ifdef PETSC_HAVE_PLAPACK
264:   PetscScalar    zero = 0.0;
265:   double         one  = 1.0;
266:   PLA_Obj        globalR = PETSC_NULL;
267: #else
268:   int            numRhs;
269: #endif
270:   int            ierr;

273:   /* Setup the PC context */
277:   if (pc->setupcalled < 2)
278:   {
279:     PCSetUp(pc);
280:   }

282:   /* Scatter in interface vector now since D^{-1} zeros out null space rows */
283:   ml = (PC_Multilevel *) pc->data;
284:   MPI_Comm_size(pc->comm, &numProcs);
285:   MPI_Comm_rank(pc->comm, &rank);
286:   if (numProcs > 1)
287:   {
288:     VecScatterBegin(x, ml->interfaceColRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceColScatter);
289:     VecScatterEnd(x, ml->interfaceColRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceColScatter);
290:   }

292:   /* Apply interior portion of D^{-1} */
293:   PCMultiLevelApplyDInvLoc_Private(pc, x, y);

295:   /* Apply R^{-1} from the QR of the interface matrix */
296: #ifdef PETSC_HAVE_PLAPACK
297:   if (numProcs > 1)
298:   {
299:     /* Put result in y */
300:     VecGetArray(ml->interfaceColRhs, &rhsArray);
301:     PLA_Obj_set_to_zero(ml->PLArhsD);
302:     PLA_API_begin();
303:     PLA_Obj_API_open(ml->PLArhsD);
304:     PLA_API_axpy_vector_to_global(ml->numLocNullCols, &one, rhsArray, 1, ml->PLArhsD, ml->firstNullCol[rank]);
305: 
306:     PLA_Obj_API_close(ml->PLArhsD);
307:     PLA_API_end();
308:     /* Solve y <-- R^{-1} y */
309:     PLA_Obj_horz_split_2(ml->interfaceQR, ml->numNullCols, &globalR, PLA_DUMMY);
310:     PLA_Trsv(PLA_UPPER_TRIANGULAR, PLA_NO_TRANSPOSE, PLA_NONUNIT_DIAG, globalR, ml->PLArhsD);
311:     PLA_Obj_free(&globalR);
312:     /* Get result from y */
313:     VecSet(&zero, ml->interfaceColRhs);
314:     PLA_API_begin();
315:     PLA_Obj_API_open(ml->PLArhsD);
316:     PLA_API_axpy_global_to_vector(ml->numLocNullCols, &one, ml->PLArhsD, ml->firstNullCol[rank], rhsArray, 1);
317: 
318:     PLA_Obj_API_close(ml->PLArhsD);
319:     PLA_API_end();
320:     VecRestoreArray(ml->interfaceColRhs, &rhsArray);
321:   }
322: #else
323:   if ((numProcs > 1) && (rank == 0))
324:   {
325:     numRhs = 1;
326:     VecGetArray(ml->interfaceColRhs, &rhsArray);
327:     LAtrtrs_("U", "N", "N", &ml->numNullCols, &numRhs, ml->interfaceQR, &ml->numInterfaceRows, rhsArray, &ml->numNullCols, &ierr);
328: 
329:     VecRestoreArray(ml->interfaceColRhs, &rhsArray);
330:     PetscLogFlops((ml->numNullCols*(ml->numNullCols-1))/2 + 2*ml->numNullCols);
331:   }
332: #endif

334:   /* Scatter out interface vector */
335:   if (numProcs > 1)
336:   {
337:     VecScatterBegin(ml->interfaceColRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceColScatter);
338:     VecScatterEnd(ml->interfaceColRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceColScatter);
339:   }
340:   return(0);
341: }

343: /* This is for just applying the interior portions of D^{-T} */
344: int PCMultiLevelApplyDInvTransLoc_Private(PC pc, GVec x, GVec y)
345: {
346: #ifdef HAVE_MATHEMATICA
347:   MLINK          link;
348: #endif
349:   PC_Multilevel *ml = (PC_Multilevel *) pc->data;
350:   double         zeroTol = ml->zeroTol;
351:   PetscScalar   *yArray;
352:   PetscReal     *invSingVals;
353:   int           *colIndices;
354:   int            level, part, col;
355:   int            ierr;

358:   /* Apply D^{-T} which is block diagonal, so we just split the vector and apply each local matrix */
359:   if (ml->useMath == PETSC_FALSE)
360:   {
361: #ifdef PETSC_USE_BOPT_g
362:     PCValidQ_Multilevel(pc);
363: #endif
364:     if (x != y)
365:       {VecCopy(x, y);                                                                              }
366:     /* Apply D^{-T} for each level */
367:     for(level = ml->numLevels-1; level >= 0; level--)
368:     {
369:       /* Apply D^{-T} for each partition */
370:       for(part = 0; part < ml->numPartitions[level]; part++)
371:       {
372:         colIndices  = ml->colPartition[level][part];
373:         invSingVals = ml->factors[level][part][FACT_DINV];
374:         /* Here, null singular values are replaced by 1 instead of zero since these columns
375:            are carried to the next level */
376:         VecGetArray(y, &yArray);
377:         for(col = 0; col < ml->numPartitionCols[level][part]; col++)
378:           if (invSingVals[col] > zeroTol)
379:             yArray[colIndices[col]] *= invSingVals[col];
380:         VecRestoreArray(y, &yArray);
381:         PetscLogFlops(ml->numPartitionCols[level][part]);
382:       }
383:     }
384:   }
385:   else
386:   {
387: #ifdef HAVE_MATHEMATICA
388:     /* The link to Mathematica */
389:     PetscViewerMathematicaGetLink(ml->mathViewer, &link);

391:     /* vec1 = input vector */
392:     PetscViewerMathematicaSetName(ml->mathViewer, "vec1");
393:     VecView(x, ml->mathViewer);

395:     /* vec2 = DInverseApply[mattML,vec] */
396:     MLPutFunction(link, "EvaluatePacket", 1);
397:       MLPutFunction(link, "Set", 2);
398:         MLPutSymbol(link, "vec2");
399:         MLPutFunction(link, "DInverseTransposeApply", 2);
400:           MLPutSymbol(link, "mattML");
401:           MLPutSymbol(link, "vec1");
402:     MLEndPacket(link);
403:     /* Skip packets until ReturnPacket */
404:     PetscViewerMathematicaSkipPackets(ml->mathViewer, RETURNPKT);
405:     /* Skip ReturnPacket */
406:     MLNewPacket(link);

408:     /* y = vec2 */
409:     PetscViewerMathematicaSetName(ml->mathViewer, "vec2");
410:     PetscViewerMathematicaGetVector(ml->mathViewer, y);
411:     PetscViewerMathematicaClearName(ml->mathViewer);
412: #endif
413:   }
414:   return(0);
415: }

417: #undef  __FUNCT__
419: /*@C PCMultiLevelApplyDInvTrans
420:         This function applies the inverse transpose of D to a vector.

422:   Input Parameters:
423: + pc - The preconditioner context
424: - x  - The input vector

426:   Output Parameter:
427: . y  - The output vector

429:   Level: intermediate

431: .keywords multilevel
432: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
433:          PCMultiLevelApplyV, PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv
434: @*/
435: int PCMultiLevelApplyDInvTrans(PC pc, GVec x, GVec y)
436: {
437:   PC_Multilevel *ml;
438:   PetscScalar   *rhsArray;
439:   int            numProcs, rank;
440: #ifdef PETSC_HAVE_PLAPACK
441:   PetscScalar    zero = 0.0;
442:   double         one  = 1.0;
443:   PLA_Obj        globalR = PETSC_NULL;
444: #else
445:   int            numRhs;
446: #endif
447:   int            ierr;

450:   /* Setup the PC context */
454:   if (pc->setupcalled < 2)
455:   {
456:     PCSetUp(pc);
457:   }

459:   /* Scatter in interface vector now since D^{-T} zeros out null space rows */
460:   ml = (PC_Multilevel *) pc->data;
461:   MPI_Comm_size(pc->comm, &numProcs);
462:   MPI_Comm_rank(pc->comm, &rank);
463:   if (numProcs > 1)
464:   {
465:     VecScatterBegin(x, ml->interfaceColRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceColScatter);
466:     VecScatterEnd(x, ml->interfaceColRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceColScatter);
467:   }

469:   /* Apply interior portion of D^{-T} */
470:   PCMultiLevelApplyDInvTransLoc_Private(pc, x, y);

472:   /* Apply R^{-T} from the QR of the interface matrix */
473: #ifdef PETSC_HAVE_PLAPACK
474:   if (numProcs > 1)
475:   {
476:     /* Put result in y */
477:     VecGetArray(ml->interfaceColRhs, &rhsArray);
478:     PLA_Obj_set_to_zero(ml->PLArhsD);
479:     PLA_API_begin();
480:     PLA_Obj_API_open(ml->PLArhsD);
481:     PLA_API_axpy_vector_to_global(ml->numLocNullCols, &one, rhsArray, 1, ml->PLArhsD, ml->firstNullCol[rank]);
482: 
483:     PLA_Obj_API_close(ml->PLArhsD);
484:     PLA_API_end();
485:     /* Solve y <-- R^{-T} y */
486:     PLA_Obj_horz_split_2(ml->interfaceQR, ml->numNullCols, &globalR, PLA_DUMMY);
487:     PLA_Trsv(PLA_UPPER_TRIANGULAR, PLA_TRANSPOSE, PLA_NONUNIT_DIAG, globalR, ml->PLArhsD);
488:     PLA_Obj_free(&globalR);
489:     /* Get result from y */
490:     VecSet(&zero, ml->interfaceColRhs);
491:     PLA_API_begin();
492:     PLA_Obj_API_open(ml->PLArhsD);
493:     PLA_API_axpy_global_to_vector(ml->numLocNullCols, &one, ml->PLArhsD, ml->firstNullCol[rank], rhsArray, 1);
494: 
495:     PLA_Obj_API_close(ml->PLArhsD);
496:     PLA_API_end();
497:     VecRestoreArray(ml->interfaceColRhs, &rhsArray);
498:   }
499: #else
500:   if ((numProcs > 1) && (rank == 0))
501:   {
502:     numRhs = 1;
503:     VecGetArray(ml->interfaceColRhs, &rhsArray);
504:     LAtrtrs_("U", "T", "N", &ml->numNullCols, &numRhs, ml->interfaceQR, &ml->numInterfaceRows, rhsArray, &ml->numNullCols, &ierr);
505: 
506:     VecRestoreArray(ml->interfaceColRhs, &rhsArray);
507:     PetscLogFlops((ml->numNullCols*(ml->numNullCols-1))/2 + 2*ml->numNullCols);
508:   }
509: #endif

511:   /* Scatter out interface vector */
512:   if (numProcs > 1)
513:   {
514:     VecScatterBegin(ml->interfaceColRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceColScatter);
515:     VecScatterEnd(ml->interfaceColRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceColScatter);
516:   }

518:   return(0);
519: }

521: #undef  __FUNCT__
523: /*@C PCMultiLevelApplyV
524:         This function applies V to a vector.

526:   Input Parameters:
527: + pc - The preconditioner context
528: - x  - The input vector

530:   Output Parameter:
531: . y  - The output vector

533:   Level: intermediate

535: .keywords multilevel
536: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
537:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
538: @*/
539: int PCMultiLevelApplyV(PC pc, GVec x, GVec y)
540: {
541: #ifdef HAVE_MATHEMATICA
542:   MLINK          link;
543: #endif
544:   PC_Multilevel *ml;
545:   PetscScalar   *yArray;
546:   PetscScalar   *colArray;
547:   PetscScalar   *colArray2;
548:   PetscReal     *VArray;
549:   int           *colIndices;
550:   PetscScalar    zero = 0.0;
551:   PetscScalar    one  = 1.0;
552:   int            dummy;
553:   int            level, part, dim, col;
554:   int            ierr;

557:   /* Setup the PC context */
561:   if (pc->setupcalled < 2)
562:   {
563:     PCSetUp(pc);
564:   }

566:   /* Apply V which is block diagonal, so we just split the vector and apply each local matrix */
567:   ml = (PC_Multilevel *) pc->data;
568:   if (ml->useMath == PETSC_FALSE)
569:   {
570: #ifdef PETSC_USE_BOPT_g
571:     PCValidQ_Multilevel(pc);
572: #endif
573:     if (x != y)
574:       {VecCopy(x, y);                                                                              }
575:     /* Apply V for each level */
576:     VecGetArray(y, &yArray);
577:     for(level = ml->numLevels-1; level >= 0; level--)
578:     {
579:       VecGetArray(ml->colReduceVecs[level],  &colArray);
580:       VecGetArray(ml->colReduceVecs2[level], &colArray2);
581:       /* Apply V for each partition */
582:       for(part = 0; part < ml->numPartitions[level]; part++)
583:       {
584:         colIndices = ml->colPartition[level][part];
585:         VArray     = ml->factors[level][part][FACT_V];
586:         dim        = ml->numPartitionCols[level][part];
587:         /* Scatter into work vector */
588:         for(col = 0; col < dim; col++)
589:           colArray[col] = yArray[colIndices[col]];
590:         dummy = 1;
591:         LAgemv_("T", &dim, &dim, &one, VArray, &dim, colArray, &dummy, &zero, colArray2, &dummy);
592:         /* Scatter from work vector */
593:         for(col = 0; col < dim; col++)
594:           yArray[colIndices[col]] = colArray2[col];
595:         PetscLogFlops(2*dim*dim - dim);
596:       }
597:       VecRestoreArray(ml->colReduceVecs[level],  &colArray);
598:       VecRestoreArray(ml->colReduceVecs2[level], &colArray2);
599:     }
600:     VecRestoreArray(y, &yArray);
601:   }
602:   else
603:   {
604: #ifdef HAVE_MATHEMATICA
605:     /* The link to Mathematica */
606:     PetscViewerMathematicaGetLink(ml->mathViewer, &link);

608:     /* vec1 = input vector */
609:     PetscViewerMathematicaSetName(ml->mathViewer, "vec1");
610:     VecView(x, ml->mathViewer);

612:     /* vec2 = VApply[mattML,vec] */
613:     MLPutFunction(link, "EvaluatePacket", 1);
614:       MLPutFunction(link, "Set", 2);
615:         MLPutSymbol(link, "vec2");
616:         MLPutFunction(link, "VTransposeApply", 2);
617:           MLPutSymbol(link, "mattML");
618:           MLPutSymbol(link, "vec1");
619:     MLEndPacket(link);
620:     /* Skip packets until ReturnPacket */
621:     PetscViewerMathematicaSkipPackets(ml->mathViewer, RETURNPKT);
622:     /* Skip ReturnPacket */
623:     MLNewPacket(link);

625:     /* y = vec2 */
626:     PetscViewerMathematicaSetName(ml->mathViewer, "vec2");
627:     PetscViewerMathematicaGetVector(ml->mathViewer, y);
628:     PetscViewerMathematicaClearName(ml->mathViewer);
629: #endif
630:   }
631:   return(0);
632: }

634: #undef  __FUNCT__
636: /*@C PCMultiLevelApplyVTrans
637:         This function applies the transpose of V to a vector.

639:   Input Parameters:
640: + pc - The preconditioner context
641: - x  - The input vector

643:   Output Parameter:
644: . y  - The output vector

646:   Level: intermediate

648: .keywords multilevel
649: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
650:          PCMultiLevelApplyV, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
651: @*/
652: int PCMultiLevelApplyVTrans(PC pc, GVec x, GVec y)
653: {
654: #ifdef HAVE_MATHEMATICA
655:   MLINK          link;
656: #endif
657:   PC_Multilevel *ml;
658:   PetscScalar   *yArray;
659:   PetscScalar   *colArray;
660:   PetscScalar   *colArray2;
661:   PetscReal     *VArray;
662:   int           *colIndices;
663:   PetscScalar    zero = 0.0;
664:   PetscScalar    one  = 1.0;
665:   int            dummy;
666:   int            level, part, dim, col;
667:   int            ierr;

670:   /* Setup the PC context */
674:   if (pc->setupcalled < 2)
675:   {
676:     PCSetUp(pc);
677:   }

679:   /* Apply V^T which is block diagonal, so we just split the vector and apply each local matrix */
680:   ml = (PC_Multilevel *) pc->data;
681:   if (ml->useMath == PETSC_FALSE)
682:   {
683: #ifdef PETSC_USE_BOPT_g
684:     PCValidQ_Multilevel(pc);
685: #endif
686:     if (x != y)
687:       {VecCopy(x, y);                                                                              }
688:     /* Apply V^T for each level */
689:     VecGetArray(y, &yArray);
690:     for(level = 0; level < ml->numLevels; level++)
691:     {
692:       VecGetArray(ml->colReduceVecs[level],  &colArray);
693:       VecGetArray(ml->colReduceVecs2[level], &colArray2);
694:       /* Apply V^T for each partition */
695:       for(part = 0; part < ml->numPartitions[level]; part++)
696:       {
697:         colIndices = ml->colPartition[level][part];
698:         VArray     = ml->factors[level][part][FACT_V];
699:         dim        = ml->numPartitionCols[level][part];
700:         /* Scatter into work vector */
701:         for(col = 0; col < dim; col++)
702:           colArray[col] = yArray[colIndices[col]];
703:         dummy = 1;
704:         LAgemv_("N", &dim, &dim, &one, VArray, &dim, colArray, &dummy, &zero, colArray2, &dummy);
705:         /* Scatter from work vector */
706:         for(col = 0; col < dim; col++)
707:           yArray[colIndices[col]] = colArray2[col];
708:         PetscLogFlops(2*dim*dim - dim);
709:       }
710:       VecRestoreArray(ml->colReduceVecs[level],  &colArray);
711:       VecRestoreArray(ml->colReduceVecs2[level], &colArray2);
712:     }
713:     VecRestoreArray(y, &yArray);
714:   }
715:   else
716:   {
717: #ifdef HAVE_MATHEMATICA
718:     /* The link to Mathematica */
719:     PetscViewerMathematicaGetLink(ml->mathViewer, &link);

721:     /* vec1 = input vector */
722:     PetscViewerMathematicaSetName(ml->mathViewer, "vec1");
723:     VecView(x, ml->mathViewer);

725:     /* vec2 = VApply[mattML,vec] */
726:     MLPutFunction(link, "EvaluatePacket", 1);
727:       MLPutFunction(link, "Set", 2);
728:         MLPutSymbol(link, "vec2");
729:         MLPutFunction(link, "VApply", 2);
730:           MLPutSymbol(link, "mattML");
731:           MLPutSymbol(link, "vec1");
732:     MLEndPacket(link);
733:     /* Skip packets until ReturnPacket */
734:     PetscViewerMathematicaSkipPackets(ml->mathViewer, RETURNPKT);
735:     /* Skip ReturnPacket */
736:     MLNewPacket(link);

738:     /* y = vec2 */
739:     PetscViewerMathematicaSetName(ml->mathViewer, "vec2");
740:     PetscViewerMathematicaGetVector(ml->mathViewer, y);
741:     PetscViewerMathematicaClearName(ml->mathViewer);
742: #endif
743:   }
744:   return(0);
745: }

747: #undef  __FUNCT__
749: /*@C PCMultiLevelApplyP
750:         This function applies P to a vector.

752:   Input Parameters:
753: + pc - The preconditioner context
754: - x  - The input vector

756:   Output Parameter:
757: . y  - The output vector

759:   Level: intermediate

761: .keywords multilevel
762: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyPTrans, PCMultiLevelApplyP1,
763:          PCMultiLevelApplyP1Trans, PCMultiLevelApplyP2, PCMultiLevelApplyP2Trans, PCMultiLevelApplyV
764:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
765: @*/
766: int PCMultiLevelApplyP(PC pc, GVec x, GVec y)
767: {
768: #ifdef HAVE_MATHEMATICA
769:   MLINK          link;
770: #endif
771:   PC_Multilevel *ml;
772:   PetscScalar   *yArray;
773:   PetscScalar   *rhsArray;
774:   PetscScalar   *localWorkArray;
775:   PetscScalar   *interiorArray;
776:   PetscScalar   *interiorArray2;
777:   PetscScalar   *bdArray;
778:   PetscScalar   *colArray;
779:   PetscScalar   *colArray2;
780:   PetscReal     *UArray;
781:   PetscReal     *QRArray;
782:   PetscReal     *TAUArray;
783:   PetscReal     *invSingVals;
784:   PetscReal     *VArray;
785:   int           *rowIndices;
786:   int            numParts, numBdRows, numResRows;
787:   int            partOffset, locColVars;
788:   int            numProcs, rank;
789:   int            nullCol, rangeCol;
790:   PetscScalar    zero =  0.0;
791:   PetscScalar    one  =  1.0;
792:   int            dummy;
793: #ifndef PETSC_HAVE_PLAPACK
794: #ifdef PETSC_USE_DEBUG
795:   int            numInterfaceRows;
796: #endif
797: #endif
798:   int            level, part, dim, col, row;
799:   int            ierr;

802:   /* Setup the PC context */
806:   if (pc->setupcalled < 2)
807:   {
808:     PCSetUp(pc);
809:   }

811:   /* Initialization */
812:   ml   = (PC_Multilevel *) pc->data;
813:   MPI_Comm_size(pc->comm, &numProcs);
814:   MPI_Comm_rank(pc->comm, &rank);

816:   /* Copy x into y if necessary */
817:   if (x != y)
818:     {VecCopy(x, y);                                                                                }

820:   /* Calculate interface values */
821:   if (numProcs > 1)
822:   {
823:     PC_MLLogEventBegin(PC_ML_ApplySymmetricRightParallel, pc, x, y, 0);
824: #ifdef PETSC_HAVE_PLAPACK
825:     /* Get the interface vector and reduce interface columns */
826:     VecScatterBegin(y, ml->interfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceScatter);
827:     VecScatterEnd(y, ml->interfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceScatter);

829:     PC_MLLogEventBegin(PC_ML_ApplyQR, pc, x, y, 0);
830:     VecGetArray(ml->interfaceRhs, &rhsArray);
831:     /* Put result in x */
832:     PLA_Obj_set_to_zero(ml->PLArhsP);
833:     PLA_API_begin();
834:     PLA_Obj_API_open(ml->PLArhsP);
835:     PLA_API_axpy_vector_to_global(ml->numLocInterfaceRows, &one, rhsArray, 1, ml->PLArhsP, ml->firstInterfaceRow[rank]);
836: 
837:     PLA_Obj_API_close(ml->PLArhsP);
838:     PLA_API_end();
839:     /* Apply x <-- Q x */
840:     PLA_Q_solve(PLA_SIDE_LEFT, PLA_TRANS, ml->interfaceQR, ml->interfaceTAU, ml->PLArhsP);
841:     /* Get result from x */
842:     VecSet(&zero, ml->interfaceRhs);
843:     PLA_API_begin();
844:     PLA_Obj_API_open(ml->PLArhsP);
845:     PLA_API_axpy_global_to_vector(ml->numLocInterfaceRows, &one, ml->PLArhsP, ml->firstInterfaceRow[rank], rhsArray, 1);
846: 
847:     PLA_Obj_API_close(ml->PLArhsP);
848:     PLA_API_end();
849:     VecRestoreArray(ml->interfaceRhs, &rhsArray);
850:     PC_MLLogEventEnd(PC_ML_ApplyQR, pc, x, y, 0);

852:     /* Set the interface values */
853:     VecScatterBegin(ml->interfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceScatter);
854:     VecScatterEnd(ml->interfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceScatter);
855: #else
856:     /* Get the interface vector and reduce interface columns */
857:     VecScatterBegin(y, ml->locInterfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->locInterfaceScatter);
858:     VecScatterEnd(y, ml->locInterfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->locInterfaceScatter);
859:     VecGetArray(ml->locInterfaceRhs, &rhsArray);

861:     PC_MLLogEventBegin(PC_ML_ApplyQR, pc, x, y, 0);
862:     if (rank == 0)
863:     {
864:       /* Apply Q from the QR of the interface matrix */
865:       dummy = 1;
866: #ifdef PETSC_USE_DEBUG
867:       VecGetSize(ml->locInterfaceRhs, &numInterfaceRows);
868:       if (numInterfaceRows != ml->numInterfaceRows) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid interface vector");
869: #endif
870:       LAormqr_("L", "N", &ml->numInterfaceRows, &dummy, &ml->numNullCols, ml->interfaceQR, &ml->numInterfaceRows,
871:              ml->interfaceTAU, rhsArray, &ml->numInterfaceRows, ml->work, &ml->workLen, &ierr);
872: 
873:       PetscLogFlops(ml->numNullCols*(2 + ml->numInterfaceRows*2));
874:     }
875:     PC_MLLogEventEnd(PC_ML_ApplyQR, pc, x, y, 0);

877:     /* Set the interface values */
878:     VecRestoreArray(ml->locInterfaceRhs, &rhsArray);
879:     VecScatterBegin(ml->locInterfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->locInterfaceScatter);
880:     VecScatterEnd(ml->locInterfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->locInterfaceScatter);

882:     /* Retrieve the local interface columns */
883:     VecScatterBegin(y, ml->interfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceScatter);
884:     VecScatterEnd(y, ml->interfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interfaceScatter);
885: #endif

887:     /* Multiply x by B^T_I */
888:     MatMultTranspose(ml->interfaceB, ml->interfaceRhs, ml->colWorkVec);

890:     /* Multiply B_p^(IT) x by V^T_p */
891:     PCMultiLevelApplyVTrans(pc, ml->colWorkVec, ml->colWorkVec);

893:     /* Multiply V^T_p B_p^(IT) x by -D_p^{-T} */
894:     PCMultiLevelApplyDInvTransLoc_Private(pc, ml->colWorkVec, ml->colWorkVec);

896:     /* Reduce rows */
897:     VecGetArray(ml->colWorkVec,         &localWorkArray);
898:     VecGetArray(y,                      &yArray);
899:     VarOrderingGetLocalSize(ml->sOrder, &locColVars);
900:     for(col = 0, nullCol = 0, rangeCol = 0; col < locColVars; col++)
901:       if (ml->nullCols[nullCol] == col)
902:         /* yArray[ml->nullCols[nullCol++]] -= localWorkArray[col]; */
903:         nullCol++;
904:       else
905:         yArray[ml->range[rangeCol++]]   -= localWorkArray[col];
906:     PetscLogFlops(rangeCol);
907:     if ((rangeCol != ml->localRank) || (nullCol != ml->numLocNullCols)) {
908:       SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid range space");
909:     }
910:     VecRestoreArray(ml->colWorkVec, &localWorkArray);
911:     VecRestoreArray(y,              &yArray);
912:     /* Here the new interface values are in y, and the old interior values */

914:     /* Scatter in interior values */
915:     VecScatterBegin(y, ml->interiorRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interiorScatter);
916:     VecScatterEnd(y, ml->interiorRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interiorScatter);
917:     PC_MLLogEventEnd(PC_ML_ApplySymmetricRightParallel, pc, x, y, 0);
918:   } else {
919:     ml->interiorRhs = y;
920:   }

922:   /* Apply P */
923:   if (ml->useMath == PETSC_FALSE)
924:   {
925: #ifdef PETSC_USE_BOPT_g
926:     PCValidQ_Multilevel(pc);
927: #endif
928:     /* Apply P for each level */
929:     VecGetArray(ml->interiorRhs, &rhsArray);
930:     interiorArray  = ml->interiorWork;
931:     interiorArray2 = ml->interiorWork2;
932:     for(level = ml->numLevels-1; level >= 0; level--)
933:     {
934:       numParts = ml->numPartitions[level];
935:       /* Scatter in boundary rows */
936:       VecGetArray(ml->bdReduceVecs[level], &bdArray);
937:       rowIndices = ml->rowPartition[level][PART_ROW_BD][0];
938:       numBdRows  = ml->numPartitionRows[level][numParts];
939:       for(row = 0; row < numBdRows; row++)
940:         bdArray[row] = rhsArray[rowIndices[row]];
941:       /* Scatter in residual rows */
942:       rowIndices = ml->rowPartition[level][PART_ROW_RES][0];
943:       numResRows = ml->numPartitionRows[level][numParts+1];
944:       for(row = 0; row < numResRows; row++)
945:         bdArray[row+numBdRows] = rhsArray[rowIndices[row]];
946:       VecRestoreArray(ml->bdReduceVecs[level], &bdArray);
947:       /* Create B^T_Gamma x */
948:       MatMultTranspose(ml->grads[level], ml->bdReduceVecs[level], ml->colReduceVecs[level]);
949:       if (numBdRows+numResRows == 0)
950:       {
951:         /* If ml->grads[level] has no rows, the default behavior is to leave ml->colReduceVecs[level] untouched */
952:         VecSet(&zero, ml->colReduceVecs[level]);
953:       }
954:       /* Reduce interior columns using / I   -D^{-T} V^T B^T_Gamma    and apply U
955:                                         0              I           /                  */
956:       VecGetArray(ml->colReduceVecs[level],  &colArray);
957:       VecGetArray(ml->colReduceVecs2[level], &colArray2);
958:       for(part = 0, partOffset = 0; part < numParts; part++)
959:       {
960:         /* Apply V^T */
961:         VArray = ml->factors[level][part][FACT_V];
962:         dim    = ml->numPartitionCols[level][part];
963:         dummy  = 1;
964:         LAgemv_("N", &dim, &dim, &one, VArray, &dim, &colArray[partOffset], &dummy, &zero, colArray2, &dummy);
965:         PetscLogFlops(2*dim*dim - dim);
966:         partOffset += dim;
967:         /* Apply D^{-T} and reduce, since we take D as rectangular we must watch out for the dimension */
968:         rowIndices  = ml->rowPartition[level][PART_ROW_INT][part];
969:         invSingVals = ml->factors[level][part][FACT_DINV];
970:         dim         = PetscMin(ml->numPartitionCols[level][part], ml->numPartitionRows[level][part]);
971:         for(col = 0; col < dim; col++)
972:           rhsArray[rowIndices[col]] -= colArray2[col] * invSingVals[col];
973:         PetscLogFlops(2*dim);
974:         /* Apply U */
975:         UArray = ml->factors[level][part][FACT_U];
976:         dim    = ml->numPartitionRows[level][part];
977:         if (dim > 0) {
978:           col = ml->numPartitionCols[level][part];
979:           if (PCMultiLevelDoQR_Private(pc, dim, col) == PETSC_TRUE) {
980:             QRArray  = ml->factors[level][part][FACT_QR];
981:             TAUArray = ml->factors[level][part][FACT_TAU];
982:             /* Scatter into work vector */
983:             for(row = 0; row < dim; row++) interiorArray[row] = rhsArray[rowIndices[row]];
984:             dummy = 1;
985:             LAgemv_("N", &col, &col, &one, UArray, &col, interiorArray, &dummy, &zero, interiorArray2, &dummy);
986:             PetscLogFlops(2*col*col - col);
987:             PetscMemcpy(interiorArray2+col, interiorArray+col, (dim - col) * sizeof(double));
988:             DQMV("N", dim, QRArray, dim, TAUArray, col, interiorArray2, interiorArray);
989:             /* Scatter from work vector */
990:             for(row = 0; row < dim; row++)
991:               rhsArray[rowIndices[row]] = interiorArray[row];
992:           }
993:           else
994:           {
995:             /* Scatter into work vector */
996:             for(row = 0; row < dim; row++)
997:               interiorArray[row] = rhsArray[rowIndices[row]];
998:             dummy = 1;
999:             LAgemv_("N", &dim, &dim, &one, UArray, &dim, interiorArray, &dummy, &zero, interiorArray2, &dummy);
1000:             PetscLogFlops(2*dim*dim - dim);
1001:             /* Scatter from work vector */
1002:             for(row = 0; row < dim; row++)
1003:               rhsArray[rowIndices[row]] = interiorArray2[row];
1004:           }
1005:         }
1006:       }
1007:       VecRestoreArray(ml->colReduceVecs[level],  &colArray);
1008:       VecRestoreArray(ml->colReduceVecs2[level], &colArray2);
1009:     }
1010:     VecRestoreArray(ml->interiorRhs, &rhsArray);
1011:   }
1012:   else
1013:   {
1014: #ifdef HAVE_MATHEMATICA
1015:     /* The link to Mathematica */
1016:     PetscViewerMathematicaGetLink(ml->mathViewer, &link);

1018:     /* vec1 = input vector */
1019:     PetscViewerMathematicaSetName(ml->mathViewer, "vec1");
1020:     VecView(ml->interiorRhs, ml->mathViewer);

1022:     /* vec2 = PApply[mattML,vec] */
1023:     MLPutFunction(link, "EvaluatePacket", 1);
1024:       MLPutFunction(link, "Set", 2);
1025:         MLPutSymbol(link, "vec2");
1026:         MLPutFunction(link, "PApply", 2);
1027:           MLPutSymbol(link, "mattML");
1028:           MLPutSymbol(link, "vec1");
1029:     MLEndPacket(link);
1030:     /* Skip packets until ReturnPacket */
1031:     PetscViewerMathematicaSkipPackets(ml->mathViewer, RETURNPKT);
1032:     /* Skip ReturnPacket */
1033:     MLNewPacket(link);

1035:     /* y = vec2 */
1036:     PetscViewerMathematicaSetName(ml->mathViewer, "vec2");
1037:     PetscViewerMathematicaGetVector(ml->mathViewer,  ml->interiorRhs);
1038:     PetscViewerMathematicaClearName(ml->mathViewer);
1039: #endif
1040:   }

1042:   /* Scatter back interior values */
1043:   if (numProcs > 1)
1044:   {
1045:     PC_MLLogEventBegin(PC_ML_ApplySymmetricRightParallel, pc, x, y, 0);
1046:     VecScatterBegin(ml->interiorRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interiorScatter);
1047:     VecScatterEnd(ml->interiorRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interiorScatter);
1048:     PC_MLLogEventEnd(PC_ML_ApplySymmetricRightParallel, pc, x, y, 0);
1049:   }

1051:   /* Scale by the diagonal of A */
1052:   if (ml->diag != PETSC_NULL) {
1053:     VecPointwiseMult(y, ml->diag, y);
1054:   }
1055:   return(0);
1056: }

1058: #undef  __FUNCT__
1060: /*@C PCMultiLevelApplyPTrans
1061:         This function applies the transpose of P a vector.

1063:   Input Parameters:
1064: + pc - The preconditioner context
1065: - x  - The input vector

1067:   Output Parameter:
1068: . y  - The output vector

1070:   Level: intermediate

1072: .keywords multilevel
1073: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyP1,
1074:          PCMultiLevelApplyP1Trans, PCMultiLevelApplyP2, PCMultiLevelApplyP2Trans, PCMultiLevelApplyV
1075:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
1076: @*/
1077: int PCMultiLevelApplyPTrans(PC pc, GVec x, GVec y)
1078: {
1079: #ifdef HAVE_MATHEMATICA
1080:   MLINK          link;
1081: #endif
1082:   PC_Multilevel *ml;
1083:   PetscScalar   *yArray;
1084:   PetscScalar   *rhsArray;
1085:   PetscScalar   *localWorkArray;
1086:   PetscScalar   *interiorArray;
1087:   PetscScalar   *interiorArray2;
1088:   PetscScalar   *bdArray;
1089:   PetscScalar   *colArray;
1090:   PetscScalar   *colArray2;
1091:   PetscReal     *UArray;
1092:   PetscReal     *QRArray;
1093:   PetscReal     *TAUArray;
1094:   PetscReal     *invSingVals;
1095:   PetscReal     *VArray;
1096:   int           *rowIndices;
1097:   int            numParts, numBdRows, numResRows;
1098:   int            partOffset, locColVars;
1099:   int            numProcs, rank;
1100:   int            nullCol, rangeCol;
1101:   PetscScalar    minusOne = -1.0;
1102:   PetscScalar    zero     =  0.0;
1103:   PetscScalar    one      =  1.0;
1104:   int            dummy;
1105: #ifndef PETSC_HAVE_PLAPACK
1106: #ifdef PETSC_USE_DEBUG
1107:   int            numInterfaceRows;
1108: #endif
1109: #endif
1110:   int            level, part, dim, row, col;
1111:   int            ierr;

1114:   /* Setup the PC context */
1118:   if (pc->setupcalled < 2)
1119:   {
1120:     PCSetUp(pc);
1121:   }

1123:   /* Initialization */
1124:   ml   = (PC_Multilevel *) pc->data;
1125:   MPI_Comm_size(pc->comm, &numProcs);
1126:   MPI_Comm_rank(pc->comm, &rank);

1128:   if (ml->diag != PETSC_NULL) {
1129:     /* Scale by the diagonal of A */
1130:     VecPointwiseMult(x, ml->diag, y);
1131:   } else if (x != y) {
1132:     /* Copy x into y if necessary */
1133:     VecCopy(x, y);
1134:   }

1136:   /* Scatter in interior values */
1137:   if (numProcs > 1)
1138:   {
1139:     PC_MLLogEventBegin(PC_ML_ApplySymmetricLeftParallel, pc, x, y, 0);
1140:     VecScatterBegin(y, ml->interiorRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interiorScatter);
1141:     VecScatterEnd(y, ml->interiorRhs, INSERT_VALUES, SCATTER_FORWARD, ml->interiorScatter);
1142:     PC_MLLogEventEnd(PC_ML_ApplySymmetricLeftParallel, pc, x, y, 0);
1143:   }
1144:   else
1145:   {
1146:     ml->interiorRhs = y;
1147:   }

1149:   /* Apply P^T */
1150:   if (ml->useMath == PETSC_FALSE)
1151:   {
1152: #ifdef PETSC_USE_BOPT_g
1153:     PCValidQ_Multilevel(pc);
1154: #endif
1155:     /* Apply P^T for each level */
1156:     VecGetArray(ml->interiorRhs, &rhsArray);
1157:     interiorArray  = ml->interiorWork;
1158:     interiorArray2 = ml->interiorWork2;
1159:     for(level = 0; level < ml->numLevels; level++)
1160:     {
1161:       numParts = ml->numPartitions[level];
1162:       VecGetArray(ml->colReduceVecs[level],  &colArray);
1163:       VecGetArray(ml->colReduceVecs2[level], &colArray2);
1164:       /* Apply U^T for each partition and form V D^{-1} */
1165:       for(part = 0, partOffset = 0; part < numParts; part++)
1166:       {
1167:         /* Apply U^T */
1168:         rowIndices = ml->rowPartition[level][PART_ROW_INT][part];
1169:         UArray     = ml->factors[level][part][FACT_U];
1170:         dim        = ml->numPartitionRows[level][part];
1171:         if (dim  > 0)
1172:         {
1173:           col = ml->numPartitionCols[level][part];
1174:           if (PCMultiLevelDoQR_Private(pc, dim, col) == PETSC_TRUE)
1175:           {
1176:             QRArray  = ml->factors[level][part][FACT_QR];
1177:             TAUArray = ml->factors[level][part][FACT_TAU];
1178:             /* Scatter into work vector */
1179:             for(row = 0; row < dim; row++)
1180:               interiorArray[row] = rhsArray[rowIndices[row]];
1181:             DQMV("T", dim, QRArray, dim, TAUArray, col, interiorArray, interiorArray2);
1182:             dummy = 1;
1183:             LAgemv_("T", &col, &col, &one, UArray, &col, interiorArray2, &dummy, &zero, interiorArray, &dummy);
1184:             PetscLogFlops(2*col*col - col);
1185:             PetscMemcpy(interiorArray+col, interiorArray2+col, (dim - col) * sizeof(double));
1186:             /* Scatter from work vector */
1187:             for(row = 0; row < dim; row++)
1188:               rhsArray[rowIndices[row]] = interiorArray[row];
1189:           }
1190:           else
1191:           {
1192:             /* Scatter into work vector */
1193:             for(row = 0; row < dim; row++)
1194:               interiorArray[row] = rhsArray[rowIndices[row]];
1195:             dummy = 1;
1196:             LAgemv_("T", &dim, &dim, &one, UArray, &dim, interiorArray, &dummy, &zero, interiorArray2, &dummy);
1197:             PetscLogFlops(2*dim*dim - dim);
1198:             /* Scatter from work vector */
1199:             for(row = 0; row < dim; row++)
1200:               rhsArray[rowIndices[row]] = interiorArray2[row];
1201:           }
1202:         }
1203:         /* Apply D^{-1}, since we take D as rectangular we must watch out for the dimension */
1204:         invSingVals = ml->factors[level][part][FACT_DINV];
1205:         dim         = PetscMin(ml->numPartitionCols[level][part], ml->numPartitionRows[level][part]);
1206:         PetscMemzero(colArray2, ml->numPartitionCols[level][part] * sizeof(double));
1207:         for(col = 0; col < dim; col++)
1208:           colArray2[col] = rhsArray[rowIndices[col]] * invSingVals[col];
1209:         PetscLogFlops(dim);
1210:         /* Apply V */
1211:         VArray = ml->factors[level][part][FACT_V];
1212:         dim    = ml->numPartitionCols[level][part];
1213:         dummy = 1;
1214:         LAgemv_("T", &dim, &dim, &one, VArray, &dim, colArray2, &dummy, &zero, &colArray[partOffset], &dummy);
1215:         PetscLogFlops(2*dim*dim - dim);
1216:         partOffset += dim;
1217:       }
1218:       VecRestoreArray(ml->colReduceVecs[level],  &colArray);
1219:       VecRestoreArray(ml->colReduceVecs2[level], &colArray2);

1221:       /* Reduce boundary columns using /          I           0 
1222:                                         -B_Gamma V D^{-1}   I / */
1223:       MatMult(ml->grads[level], ml->colReduceVecs[level], ml->bdReduceVecs[level]);
1224:       VecGetArray(ml->bdReduceVecs[level], &bdArray);
1225:       /* Update boundary rows */
1226:       rowIndices = ml->rowPartition[level][PART_ROW_BD][0];
1227:       numBdRows  = ml->numPartitionRows[level][numParts];
1228:       for(row = 0; row < numBdRows; row++)
1229:         rhsArray[rowIndices[row]] -= bdArray[row];
1230:       /* Update residual rows */
1231:       rowIndices = ml->rowPartition[level][PART_ROW_RES][0];
1232:       numResRows = ml->numPartitionRows[level][numParts+1];
1233:       for(row = 0; row < numResRows; row++)
1234:         rhsArray[rowIndices[row]] -= bdArray[row+numBdRows];
1235:       PetscLogFlops(numBdRows+numResRows);
1236:       VecRestoreArray(ml->bdReduceVecs[level], &bdArray);
1237:     }
1238:     VecRestoreArray(ml->interiorRhs, &rhsArray);
1239:   }
1240:   else
1241:   {
1242: #ifdef HAVE_MATHEMATICA
1243:     /* The link to Mathematica */
1244:     PetscViewerMathematicaGetLink(ml->mathViewer, &link);

1246:     /* vec1 = input vector */
1247:     PetscViewerMathematicaSetName(ml->mathViewer, "vec1");
1248:     VecView(ml->interiorRhs, ml->mathViewer);

1250:     /* vec2 = PApply[mattML,vec] */
1251:     MLPutFunction(link, "EvaluatePacket", 1);
1252:       MLPutFunction(link, "Set", 2);
1253:         MLPutSymbol(link, "vec2");
1254:         MLPutFunction(link, "PTransposeApply", 2);
1255:           MLPutSymbol(link, "mattML");
1256:           MLPutSymbol(link, "vec1");
1257:     MLEndPacket(link);
1258:     /* Skip packets until ReturnPacket */
1259:     PetscViewerMathematicaSkipPackets(ml->mathViewer, RETURNPKT);
1260:     /* Skip ReturnPacket */
1261:     MLNewPacket(link);

1263:     /* y = vec2 */
1264:     PetscViewerMathematicaSetName(ml->mathViewer, "vec2");
1265:     PetscViewerMathematicaGetVector(ml->mathViewer,  ml->interiorRhs);
1266:     PetscViewerMathematicaClearName(ml->mathViewer);
1267: #endif
1268:   }

1270:   if (numProcs > 1)
1271:   {
1272:     PC_MLLogEventBegin(PC_ML_ApplySymmetricLeftParallel, pc, x, y, 0);
1273:     /* Scatter back interior values */
1274:     VecScatterBegin(ml->interiorRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interiorScatter);
1275:     VecScatterEnd(ml->interiorRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interiorScatter);

1277:     /* Calculate interface values */
1278:     /* Apply (I - D^{-1} D): Scatter some of the interior of y into a column work vector */
1279:     VecGetArray(y,                      &yArray);
1280:     VecGetArray(ml->colWorkVec,         &localWorkArray);
1281:     VarOrderingGetLocalSize(ml->sOrder, &locColVars);
1282:     PetscMemzero(localWorkArray, locColVars * sizeof(PetscScalar));
1283:     for(col = 0, nullCol = 0, rangeCol = 0; col < locColVars; col++)
1284:       if (ml->nullCols[nullCol] == col)
1285:         /* localWorkArray[col] = yArray[ml->nullCols[nullCol++]]; */
1286:         nullCol++;
1287:       else
1288:         localWorkArray[col] = yArray[ml->range[rangeCol++]];
1289:     if ((rangeCol != ml->localRank) || (nullCol != ml->numLocNullCols)) {
1290:       SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid range space");
1291:     }
1292:     VecRestoreArray(y,                    &yArray);
1293:     VecRestoreArray(ml->colWorkVec,       &localWorkArray);

1295:     /* Multiply y by D_p^{-1} */
1296:     PCMultiLevelApplyDInvLoc_Private(pc, ml->colWorkVec, ml->colWorkVec);

1298:     /* Multiply (I - D_p^{-1} D_p) y by V_p */
1299:     PCMultiLevelApplyV(pc, ml->colWorkVec, ml->colWorkVec);

1301:     /* Multiply V (I - D^{-1} D) y by -B_I */
1302:     MatMult(ml->interfaceB, ml->colWorkVec, ml->interfaceRhs);
1303:     VecScale(&minusOne, ml->interfaceRhs);

1305: #ifdef PETSC_HAVE_PLAPACK
1306:     /* Reduce the local interface columns */
1307:     VecScatterBegin(y, ml->interfaceRhs, ADD_VALUES, SCATTER_FORWARD, ml->interfaceScatter);
1308:     VecScatterEnd(y, ml->interfaceRhs, ADD_VALUES, SCATTER_FORWARD, ml->interfaceScatter);

1310:     PC_MLLogEventBegin(PC_ML_ApplyQR, pc, x, y, 0);
1311:     VecGetArray(ml->interfaceRhs, &rhsArray);
1312:     /* Put result in x */
1313:     PLA_Obj_set_to_zero(ml->PLArhsP);
1314:     PLA_API_begin();
1315:     PLA_Obj_API_open(ml->PLArhsP);
1316:     PLA_API_axpy_vector_to_global(ml->numLocInterfaceRows, &one, rhsArray, 1, ml->PLArhsP, ml->firstInterfaceRow[rank]);
1317: 
1318:     PLA_Obj_API_close(ml->PLArhsP);
1319:     PLA_API_end();
1320:     /* Apply x <-- Q^T x */
1321:     PLA_Q_solve(PLA_SIDE_LEFT, PLA_NO_TRANS, ml->interfaceQR, ml->interfaceTAU, ml->PLArhsP);
1322:     /* Get result from x */
1323:     VecSet(&zero, ml->interfaceRhs);
1324:     PLA_API_begin();
1325:     PLA_Obj_API_open(ml->PLArhsP);
1326:     PLA_API_axpy_global_to_vector(ml->numLocInterfaceRows, &one, ml->PLArhsP, ml->firstInterfaceRow[rank], rhsArray, 1);
1327: 
1328:     PLA_Obj_API_close(ml->PLArhsP);
1329:     PLA_API_end();
1330:     VecRestoreArray(ml->interfaceRhs, &rhsArray);
1331:     PC_MLLogEventEnd(PC_ML_ApplyQR, pc, x, y, 0);

1333:     /* Set the interface values */
1334:     VecScatterBegin(ml->interfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceScatter);
1335:     VecScatterEnd(ml->interfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->interfaceScatter);
1336: #else
1337:     /* Reduce the local interface columns */
1338:     VecScatterBegin(ml->interfaceRhs, y, ADD_VALUES, SCATTER_REVERSE, ml->interfaceScatter);
1339:     VecScatterEnd(ml->interfaceRhs, y, ADD_VALUES, SCATTER_REVERSE, ml->interfaceScatter);

1341:     /* Get the interface vector */
1342:     VecScatterBegin(y, ml->locInterfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->locInterfaceScatter);
1343:     VecScatterEnd(y, ml->locInterfaceRhs, INSERT_VALUES, SCATTER_FORWARD, ml->locInterfaceScatter);
1344:     VecGetArray(ml->locInterfaceRhs, &rhsArray);

1346:     PC_MLLogEventBegin(PC_ML_ApplyQR, pc, x, y, 0);
1347:     if (rank == 0) {
1348:       /* Apply Q^T from the QR of the interface matrix */
1349:       dummy = 1;
1350: #ifdef PETSC_USE_DEBUG
1351:       VecGetSize(ml->locInterfaceRhs, &numInterfaceRows);
1352:       if (numInterfaceRows != ml->numInterfaceRows) SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid interface vector");
1353: #endif
1354:       LAormqr_("L", "T", &ml->numInterfaceRows, &dummy, &ml->numNullCols, ml->interfaceQR, &ml->numInterfaceRows,
1355:              ml->interfaceTAU, rhsArray, &ml->numInterfaceRows, ml->work, &ml->workLen, &ierr);
1356:       PetscLogFlops(ml->numNullCols*(2 + ml->numInterfaceRows*2));
1357: 
1358:     }
1359:     PC_MLLogEventEnd(PC_ML_ApplyQR, pc, x, y, 0);

1361:     /* Set the interface values */
1362:     VecRestoreArray(ml->locInterfaceRhs, &rhsArray);
1363:     VecScatterBegin(ml->locInterfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->locInterfaceScatter);
1364:     VecScatterEnd(ml->locInterfaceRhs, y, INSERT_VALUES, SCATTER_REVERSE, ml->locInterfaceScatter);
1365: #endif
1366:     PC_MLLogEventEnd(PC_ML_ApplySymmetricLeftParallel, pc, x, y, 0);
1367:   }
1368:   return(0);
1369: }

1371: #undef  __FUNCT__
1373: /*@C PCMultiLevelApplyP1
1374:         This function applies P_1, the projector on the range of B, to a vector.

1376:   Input Parameters:
1377: + pc - The preconditioner context
1378: - x  - The input vector

1380:   Output Parameter:
1381: . y  - The output vector

1383:   Level: intermediate

1385: .keywords multilevel
1386: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
1387:          PCMultiLevelApplyP1Trans, PCMultiLevelApplyP2, PCMultiLevelApplyP2Trans, PCMultiLevelApplyV
1388:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
1389: @*/
1390: int PCMultiLevelApplyP1(PC pc, GVec x, GVec y)
1391: {
1392:   PC_Multilevel *ml;
1393:   PetscScalar    zero = 0.0;
1394:   int            ierr;

1397:   /* Setup the PC context */
1401:   if (pc->setupcalled < 2) {
1402:     PCSetUp(pc);
1403:   }

1405:   /* Scatter the column vector x to the solution vector y */
1406:   ml   = (PC_Multilevel *) pc->data;
1407:   VecSet(&zero, y);
1408:   VecScatterBegin(x, y, INSERT_VALUES, SCATTER_REVERSE, ml->rangeScatter);
1409:   VecScatterEnd(x, y, INSERT_VALUES, SCATTER_REVERSE, ml->rangeScatter);

1411:   /* Apply P */
1412:   PCMultiLevelApplyP(pc, y, y);
1413:   return(0);
1414: }

1416: #undef  __FUNCT__
1418: /*@C PCMultiLevelApplyP1Trans
1419:         This function applies the transpose of P_1, the projector on the range of B, to a vector.

1421:   Input Parameters:
1422: + pc - The preconditioner context
1423: - x  - The input vector

1425:   Output Parameter:
1426: . y  - The output vector

1428:   Level: intermediate

1430: .keywords multilevel
1431: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
1432:          PCMultiLevelApplyP1, PCMultiLevelApplyP2, PCMultiLevelApplyP2Trans, PCMultiLevelApplyV
1433:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
1434: @*/
1435: int PCMultiLevelApplyP1Trans(PC pc, GVec x, GVec y)
1436: {
1437:   PC_Multilevel *ml;
1438:   GVec           z;
1439:   PetscScalar    zero = 0.0;
1440:   int            ierr;

1443:   /* Setup the PC context */
1447:   if (pc->setupcalled < 2)
1448:   {
1449:     PCSetUp(pc);
1450:   }

1452:   /* Apply P^T */
1453:   VecDuplicate(x, &z);
1454:   PCMultiLevelApplyPTrans(pc, x, z);

1456:   /* Scatter the solution vector z to the column vector y */
1457:   ml = (PC_Multilevel *) pc->data;
1458:   VecSet(&zero, y);
1459:   VecScatterBegin(z, y, INSERT_VALUES, SCATTER_FORWARD, ml->rangeScatter);
1460:   VecScatterEnd(z, y, INSERT_VALUES, SCATTER_FORWARD, ml->rangeScatter);
1461:   VecDestroy(z);
1462:   return(0);
1463: }

1465: #undef  __FUNCT__
1467: /*@C PCMultiLevelApplyP2
1468:         This function applies P_2, the projector on the nullspace of B, to a vector.
1469:   Note that this function is setup to take two vectors of the same size. The
1470:   component of x and y in the range of B is zero on output.

1472:   Input Parameters:
1473: + pc - The preconditioner context
1474: - x  - The input vector

1476:   Output Parameter:
1477: . y  - The output vector

1479:   Level: intermediate

1481: .keywords multilevel
1482: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
1483:          PCMultiLevelApplyP1, PCMultiLevelApplyP1Trans, PCMultiLevelApplyP2Trans, PCMultiLevelApplyV
1484:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
1485: @*/
1486: int PCMultiLevelApplyP2(PC pc, GVec x, GVec y)
1487: {
1488:   PC_Multilevel *ml;
1489:   PetscScalar   *xArray;
1490:   int            row;
1491:   int            ierr;

1494:   /* Setup the PC context */
1498:   if (pc->setupcalled < 2)
1499:   {
1500:     PCSetUp(pc);
1501:   }

1503:   /* Zero out range space of B */
1504:   ml = (PC_Multilevel *) pc->data;
1505:   VecGetArray(x, &xArray);
1506:   for(row = 0; row < ml->rank; row++)
1507:   {
1508:     xArray[ml->range[row]] = 0.0;
1509:   }
1510:   VecRestoreArray(x, &xArray);

1512:   /* Apply P */
1513:   PCMultiLevelApplyP(pc, x, y);
1514:   return(0);
1515: }

1517: #undef  __FUNCT__
1519: /*@C PCMultiLevelApplyP2Trans
1520:         This function applies the transpose of P_2, the projector on the nullspace of B, to a vector.
1521:   Note that this function is setup to take two vectors of the same size. The component of y in
1522:   the range of B is zero on output.

1524:   Input Parameters:
1525: + pc - The preconditioner context
1526: - x  - The input vector

1528:   Output Parameter:
1529: . y  - The output vector

1531:   Level: intermediate

1533: .keywords multilevel
1534: .seealso PCMultiLevelApplyGradient, PCMultiLevelApplyGradientTrans, PCMultiLevelApplyP, PCMultiLevelApplyPTrans,
1535:          PCMultiLevelApplyP1, PCMultiLevelApplyP1Trans, PCMultiLevelApplyP2, PCMultiLevelApplyV
1536:          PCMultiLevelApplyVTrans, PCMultiLevelApplyDInv, PCMultiLevelApplyDInvTrans
1537: @*/
1538: int PCMultiLevelApplyP2Trans(PC pc, GVec x, GVec y)
1539: {
1540:   PC_Multilevel *ml;
1541:   PetscScalar   *yArray;
1542:   int            row;
1543:   int            ierr;

1549:   /* Apply P */
1550:   PCMultiLevelApplyPTrans(pc, x, y);

1552:   /* Zero out range space of B */
1553:   ml = (PC_Multilevel *) pc->data;
1554:   VecGetArray(y, &yArray);
1555:   for(row = 0; row < ml->rank; row++)
1556:   {
1557:     yArray[ml->range[row]] = 0.0;
1558:   }
1559:   VecRestoreArray(y, &yArray);
1560:   return(0);
1561: }