Actual source code: ml.h
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: ml.h,v 1.6 2001/01/27 21:32:36 knepley Exp $";
3: #endif
4: /*
5: Private data structure for Multilevel preconditioner.
6: */
10: #include gvec.h
11: #ifdef PETSC_HAVE_PLAPACK
12: #include "PLA.h"
13: #endif
14: #include petscblaslapack.h
16: enum {PC_ML_SetUpInit, PC_ML_SetUpConstrained, PC_ML_SetUpConstrainedBd, PC_ML_SetUpParallel,
17: PC_ML_ReducePartitionMesh, PC_ML_ReducePartitionRowCol, PC_ML_ReduceFactor, PC_ML_ReduceBdGrad, PC_ML_ReduceBdGradExtract,
18: PC_ML_ReduceBdGradRowPartLocalToGlobal, PC_ML_ReduceShrinkMesh, PC_ML_ApplySymmetricLeftParallel,
19: PC_ML_ApplySymmetricRightParallel, PC_ML_QRFactorization, PC_ML_ApplyQR, PC_ML_CreateBdGrad, PC_ML_MAX_EVENTS};
20: extern int PC_MLEvents[PC_ML_MAX_EVENTS];
21: #define PC_MLLogEventBegin(e,o1,o2,o3,o4) PetscLogEventBegin(PC_MLEvents[e],o1,o2,o3,o4)
22: #define PC_MLLogEventEnd(e,o1,o2,o3,o4) PetscLogEventEnd(PC_MLEvents[e],o1,o2,o3,o4)
24: enum {MESH_OFFSETS, MESH_ADJ, MESH_ELEM, NUM_MESH_DIV};
25: enum {PART_ROW_INT, PART_ROW_BD, PART_ROW_RES, NUM_PART_ROW_DIV};
26: enum {FACT_U, FACT_DINV, FACT_V, FACT_QR, FACT_TAU, NUM_FACT_DIV};
28: typedef struct {
29: Grid grid; /* The problem grid */
30: PetscTruth useMath; /* Use Mathematica functions instead */
31: PetscViewer mathViewer; /* The link to Mathematica */
32: PetscViewer factorizationViewer; /* The visualization of the factorization process */
34: /* Gradient variabels */
35: int gradOp; /* The gradient operator */
36: int divOp; /* The divergence operator */
37: int sField; /* The shape function field */
38: int tField; /* The test function field */
39: VarOrdering sOrder; /* The global variable ordering for the shape function in B */
40: LocalVarOrdering sLocOrder; /* The local variable ordering for the shape function in B */
41: VarOrdering tOrder; /* The global variable ordering for the test function in B */
42: LocalVarOrdering tLocOrder; /* The local variable ordering for the test function in B */
43: PetscScalar gradAlpha; /* The scalar multiplier for the gradient operator */
44: GMat B; /* The gradient matrix */
45: Mat locB; /* The local gradient matrix */
46: GVec reduceSol; /* The solution in the range of B */
47: GVec nonlinearIterate; /* The current iterate in a Newton solve */
48: GVec diag; /* The diagonal of the jacobian (approximately) */
50: /* Mesh variables */
51: int numMeshes; /* The number of meshes corresponding to the factorization hierarchy */
52: int partSize; /* The target size of column partitions */
53: int *numElements; /* The number of elements in each mesh */
54: int *numVertices; /* The number of vertices in each mesh */
55: double *vertices; /* The coordinates for each vertex */
56: int numEdges; /* The number of edges in the original mesh */
57: int ***meshes; /* The mesh hierarchy in CSR format */
58:
59: /* Numbering variables */
60: int *numPartitions; /* The number of partitions for each level */
61: int **numPartitionCols; /* The number of columns in each partition and the final QR */
62: int ***colPartition; /* The column or node partition */
63: int **numPartitionRows; /* The number of rows in each partition followed by the boundary and residual */
64: int ****rowPartition; /* The row or edge partition */
65: int rank; /* The rank of B^p which is the length of range */
66: int *range; /* The rows of P in the range of B ordered so that D is diagonal */
67: VecScatter rangeScatter; /* The map from a solution vector to a column vector */
68: int numLocNullCols; /* The number of columns corresponding to zero singular values */
69: int *nullCols; /* The columns corresponding to zero singular values */
71: /* Factorization variables */
72: int numRows; /* The number of rows in B */
73: int numCols; /* The number of columns in B */
74: int numLevels; /* The number of levels in the factorization */
75: int maxLevels; /* The threshold for new allocation */
76: int QRthresh; /* The node threshold for performing a final QR in this domain */
77: PetscReal DoQRFactor; /* The multiplier m for row > cols*m which determines when to QR U */
78: double zeroTol; /* The numeric tolerance used for determining zero singular values */
79: PetscReal ****factors; /* The factor list */
80: Mat *grads; /* The boundary gradients for each level */
81: Vec *bdReduceVecs; /* The work vectors for the boundary+residual rows at each level */
82: Vec *colReduceVecs; /* The work vectors for the remaining columns at each level */
83: Vec *colReduceVecs2; /* A duplicate of colReduceVecs[] */
84: PetscReal *interiorWork; /* A work vector the size of the maximum rows in any partition */
85: PetscReal *interiorWork2; /* A duplicate of interiorWork[] */
86: int interiorWorkLen; /* The size of interiorWork[] */
88: /* Parallel variables */
89: int globalRank; /* The rank of B, or total number of range space vectors */
90: int localRank; /* The rank of B_p */
91: int numNullCols; /* The total number of columns corresponding to zero singular values */
92: int *firstNullCol; /* The first column with a zero singular value in each domain */
93: int numLocInterfaceRows; /* The number of rows with nonzeros in another domain */
94: int numInterfaceRows; /* The global number of interface rows */
95: int *interfaceRows; /* The row numbers of interface rows in the original matrix */
96: int *firstInterfaceRow; /* The first interface row in each domain */
97: int numInteriorRows; /* The number of rows which are not interface rows in this domain */
98: int *interiorRows; /* The row numbers of interior rows in the original matrix */
99: #ifdef PETSC_HAVE_PLAPACK
100: PLA_Template PLAlayout; /* Layout of PLA structures across processors */
101: PLA_Obj interfaceQR; /* The QR factorization of the interface matrix B_I */
102: PLA_Obj interfaceTAU; /* The elementary reflection multipliers for interfaceQR */
103: PLA_Obj PLArhsD; /* The work vector for applying D size matrices */
104: PLA_Obj PLArhsP; /* The work vector for applying P size matrices */
105: #else
106: PetscReal *interfaceQR; /* The QR factorization of the interface matrix B_I */
107: PetscReal *interfaceTAU; /* The elementary reflection multipliers for interfaceQR */
108: #ifdef CHOLESKY_SCHEME
109: PetscReal *interfaceGrad; /* The interface matrix B_I */
110: #endif
111: #endif
112: GVec *interfaceV; /* The column of V^i corresponding to the zero singular value in that domain */
113: Mat interfaceB; /* The interface rows of B for each domain */
114: Vec interiorRhs; /* The rhs for solves in the interior */
115: VecScatter interiorScatter; /* The map from a global solution vector to an interior vector */
116: #ifndef PETSC_HAVE_PLAPACK
117: Vec locInterfaceRhs; /* The rhs for reduction of the local interface */
118: VecScatter locInterfaceScatter; /* The map from a global solution vector to a local interface vector */
119: #endif
120: Vec interfaceRhs; /* The rhs for solves on the interface: only on proc 1 */
121: VecScatter interfaceScatter; /* The map from a global solution vector to an interface vector: only on proc 1 */
122: Vec interfaceColRhs; /* The rhs for solves on the column interface: only on proc 1 */
123: VecScatter interfaceColScatter; /* The map from a global column vector to a column interface vector: only on proc 1 */
124: VecScatter reduceScatter; /* The map from a parallel column vector to a serial column vector on each proc */
125: GVec colWorkVec; /* A Work vector in the column space of B */
126: Vec globalColWorkVec; /* A work vector the size of colWorkVec, but duplicated on each processor */
128: /* Work space for LAPACK routines */
129: int workLen;
130: PetscScalar *work;
131: int svdWorkLen;
132: PetscScalar *svdWork;
133: } PC_Multilevel;
135: /* mlSerial prototypes */
136: extern PetscTruth PCMultiLevelDoQR_Private(PC, int, int);
137: extern int PCMultiLevelCalcVertexCoords_Private(PC, int, int, int *);
138: extern int PCMultiLevelFilterVertexCoords_Tutte(PC, int, double *, int **);
139: extern int PCMultiLevelReduceView(PC, int, int, int *, int, int *);
140: extern int PCMultiLevelReduce(PC);
141: extern int PCMultiLevelView_Private(PC, PetscViewer, int, int, int, int **);
143: /* mlApply prototypes */
144: extern int PCMultiLevelApplyDInv(PC, GVec, GVec);
145: extern int PCMultiLevelApplyDInvTrans(PC, GVec, GVec);
146: extern int PCMultiLevelApplyV(PC, GVec, GVec);
147: extern int PCMultiLevelApplyVTrans(PC, GVec, GVec);
148: extern int PCMultiLevelApplyP(PC, GVec, GVec);
149: extern int PCMultiLevelApplyPTrans(PC, GVec, GVec);
150: extern int PCMultiLevelApplyP1(PC, GVec, GVec);
151: extern int PCMultiLevelApplyP1Trans(PC, GVec, GVec);
152: extern int PCMultiLevelApplyP2(PC, GVec, GVec);
153: extern int PCMultiLevelApplyP2Trans(PC, GVec, GVec);
155: /* mlCheck prototypes */
156: extern int PCValidQ_Multilevel(PC);
157: extern int PCDebug_Multilevel(PC);
158: extern int PCDebugPrintMat_Multilevel(PC, const char []);
160: #endif