Actual source code: mscatter.c

petsc-3.3-p7 2013-05-11
  2: /*
  3:    This provides a matrix that applies a VecScatter to a vector.
  4: */

  6: #include <petsc-private/matimpl.h>        /*I "petscmat.h" I*/
  7: #include <petsc-private/vecimpl.h>  

  9: typedef struct {
 10:   VecScatter scatter;
 11: } Mat_Scatter;

 15: /*@
 16:     MatScatterGetVecScatter - Returns the user-provided scatter set with MatScatterSetVecScatter()

 18:     Not Collective, but not cannot use scatter if not used collectively on Mat

 20:     Input Parameter:
 21: .   mat - the matrix, should have been created with MatCreateScatter() or have type MATSCATTER

 23:     Output Parameter:
 24: .   scatter - the scatter context

 26:     Level: intermediate

 28: .keywords: matrix, scatter, get

 30: .seealso: MatCreateScatter(), MatScatterSetVecScatter(), MATSCATTER
 31: @*/
 32: PetscErrorCode  MatScatterGetVecScatter(Mat mat,VecScatter *scatter)
 33: {
 34:   Mat_Scatter    *mscatter;

 39:   mscatter = (Mat_Scatter*)mat->data;
 40:   *scatter = mscatter->scatter;
 41:   return(0);
 42: }

 46: PetscErrorCode MatDestroy_Scatter(Mat mat)
 47: {
 49:   Mat_Scatter    *scatter = (Mat_Scatter*)mat->data;

 52:   VecScatterDestroy(&scatter->scatter);
 53:   PetscFree(mat->data);
 54:   return(0);
 55: }

 59: PetscErrorCode MatMult_Scatter(Mat A,Vec x,Vec y)
 60: {
 61:   Mat_Scatter    *scatter = (Mat_Scatter*)A->data;

 65:   if (!scatter->scatter) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to first call MatScatterSetScatter()");
 66:   VecZeroEntries(y);
 67:   VecScatterBegin(scatter->scatter,x,y,ADD_VALUES,SCATTER_FORWARD);
 68:   VecScatterEnd(scatter->scatter,x,y,ADD_VALUES,SCATTER_FORWARD);
 69:   return(0);
 70: }

 74: PetscErrorCode MatMultAdd_Scatter(Mat A,Vec x,Vec y,Vec z)
 75: {
 76:   Mat_Scatter    *scatter = (Mat_Scatter*)A->data;

 80:   if (!scatter->scatter) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to first call MatScatterSetScatter()");
 81:   if (z != y) {VecCopy(y,z);}
 82:   VecScatterBegin(scatter->scatter,x,z,ADD_VALUES,SCATTER_FORWARD);
 83:   VecScatterEnd(scatter->scatter,x,z,ADD_VALUES,SCATTER_FORWARD);
 84:   return(0);
 85: }

 89: PetscErrorCode MatMultTranspose_Scatter(Mat A,Vec x,Vec y)
 90: {
 91:   Mat_Scatter    *scatter = (Mat_Scatter*)A->data;

 95:   if (!scatter->scatter) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to first call MatScatterSetScatter()");
 96:   VecZeroEntries(y);
 97:   VecScatterBegin(scatter->scatter,x,y,ADD_VALUES,SCATTER_REVERSE);
 98:   VecScatterEnd(scatter->scatter,x,y,ADD_VALUES,SCATTER_REVERSE);
 99:   return(0);
100: }

104: PetscErrorCode MatMultTransposeAdd_Scatter(Mat A,Vec x,Vec y,Vec z)
105: {
106:   Mat_Scatter    *scatter = (Mat_Scatter*)A->data;

110:   if (!scatter->scatter) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to first call MatScatterSetScatter()");
111:   if (z != y) {VecCopy(y,z);}
112:   VecScatterBegin(scatter->scatter,x,z,ADD_VALUES,SCATTER_REVERSE);
113:   VecScatterEnd(scatter->scatter,x,z,ADD_VALUES,SCATTER_REVERSE);
114:   return(0);
115: }

117: static struct _MatOps MatOps_Values = {0,
118:        0,
119:        0,
120:        MatMult_Scatter,
121: /* 4*/ MatMultAdd_Scatter,
122:        MatMultTranspose_Scatter,
123:        MatMultTransposeAdd_Scatter,
124:        0,
125:        0,
126:        0,
127: /*10*/ 0,
128:        0,
129:        0,
130:        0,
131:        0,
132: /*15*/ 0,
133:        0,
134:        0,
135:        0,
136:        0,
137: /*20*/ 0,
138:        0,
139:        0,
140:        0,
141: /*24*/ 0,
142:        0,
143:        0,
144:        0,
145:        0,
146: /*29*/ 0,
147:        0,
148:        0,
149:        0,
150:        0,
151: /*34*/ 0,
152:        0,
153:        0,
154:        0,
155:        0,
156: /*39*/ 0,
157:        0,
158:        0,
159:        0,
160:        0,
161: /*44*/ 0,
162:        0,
163:        0,
164:        0,
165:        0,
166: /*49*/ 0,
167:        0,
168:        0,
169:        0,
170:        0,
171: /*54*/ 0,
172:        0,
173:        0,
174:        0,
175:        0,
176: /*59*/ 0,
177:        MatDestroy_Scatter,
178:        0,
179:        0,
180:        0,
181: /*64*/ 0,
182:        0,
183:        0,
184:        0,
185:        0,
186: /*69*/ 0,
187:        0,
188:        0,
189:        0,
190:        0,
191: /*74*/ 0,
192:        0,
193:        0,
194:        0,
195:        0,
196: /*79*/ 0,
197:        0,
198:        0,
199:        0,
200:        0,
201: /*84*/ 0,
202:        0,
203:        0,
204:        0,
205:        0,
206: /*89*/ 0,
207:        0,
208:        0,
209:        0,
210:        0,
211: /*94*/ 0,
212:        0,
213:        0,
214:        0};

216: /*MC
217:    MATSCATTER - MATSCATTER = "scatter" - A matrix type that simply applies a VecScatterBegin/End()

219:   Level: advanced

221: .seealso: MatCreateScatter(), MatScatterSetVecScatter(), MatScatterGetVecScatter()

223: M*/

225: EXTERN_C_BEGIN
228: PetscErrorCode  MatCreate_Scatter(Mat A)
229: {
230:   Mat_Scatter    *b;

234:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
235:   PetscNewLog(A,Mat_Scatter,&b);

237:   A->data = (void*)b;

239:   PetscLayoutSetUp(A->rmap);
240:   PetscLayoutSetUp(A->cmap);

242:   A->assembled     = PETSC_TRUE;
243:   A->preallocated  = PETSC_FALSE;

245:   PetscObjectChangeTypeName((PetscObject)A,MATSCATTER);
246:   return(0);
247: }
248: EXTERN_C_END

252: /*@C
253:    MatCreateScatter - Creates a new matrix based on a VecScatter

255:   Collective on MPI_Comm

257:    Input Parameters:
258: +  comm - MPI communicator
259: -  scatter - a VecScatterContext

261:    Output Parameter:
262: .  A - the matrix

264:    Level: intermediate

266:    PETSc requires that matrices and vectors being used for certain
267:    operations are partitioned accordingly.  For example, when
268:    creating a scatter matrix, A, that supports parallel matrix-vector
269:    products using MatMult(A,x,y) the user should set the number
270:    of local matrix rows to be the number of local elements of the
271:    corresponding result vector, y. Note that this is information is
272:    required for use of the matrix interface routines, even though
273:    the scatter matrix may not actually be physically partitioned.

275: .keywords: matrix, scatter, create

277: .seealso: MatScatterSetVecScatter(), MatScatterGetVecScatter(), MATSCATTER
278: @*/
279: PetscErrorCode  MatCreateScatter(MPI_Comm comm,VecScatter scatter,Mat *A)
280: {

284:   MatCreate(comm,A);
285:   MatSetSizes(*A,scatter->to_n,scatter->from_n,PETSC_DETERMINE,PETSC_DETERMINE);
286:   MatSetType(*A,MATSCATTER);
287:   MatScatterSetVecScatter(*A,scatter);
288:   MatSetUp(*A);
289:   return(0);
290: }

294: /*@
295:     MatScatterSetVecScatter - sets that scatter that the matrix is to apply as its linear operator

297:    Collective on Mat

299:     Input Parameters:
300: +   mat - the scatter matrix
301: -   scatter - the scatter context create with VecScatterCreate()

303:    Level: advanced


306: .seealso: MatCreateScatter(), MATSCATTER
307: @*/
308: PetscErrorCode  MatScatterSetVecScatter(Mat mat,VecScatter scatter)
309: {
310:   Mat_Scatter    *mscatter = (Mat_Scatter*)mat->data;

317:   if (mat->rmap->n != scatter->to_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number of local rows in matrix %D not equal local scatter size %D",mat->rmap->n,scatter->to_n);
318:   if (mat->cmap->n != scatter->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number of local columns in matrix %D not equal local scatter size %D",mat->cmap->n,scatter->from_n);

320:   PetscObjectReference((PetscObject)scatter);
321:   VecScatterDestroy(&mscatter->scatter);
322:   mscatter->scatter = scatter;
323:   return(0);
324: }