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