1 #include <../src/mat/impls/nest/matnestimpl.h> /*I "petscmat.h" I*/
2 #include <../src/mat/impls/aij/seq/aij.h>
3 #include <petscsf.h>
4
5 static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
6 static PetscErrorCode MatCreateVecs_Nest(Mat,Vec*,Vec*);
7 static PetscErrorCode MatReset_Nest(Mat);
8
9 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat,MatType,MatReuse,Mat*);
10
11 /* private functions */
MatNestGetSizes_Private(Mat A,PetscInt * m,PetscInt * n,PetscInt * M,PetscInt * N)12 static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N)
13 {
14 Mat_Nest *bA = (Mat_Nest*)A->data;
15 PetscInt i,j;
16 PetscErrorCode ierr;
17
18 PetscFunctionBegin;
19 *m = *n = *M = *N = 0;
20 for (i=0; i<bA->nr; i++) { /* rows */
21 PetscInt sm,sM;
22 ierr = ISGetLocalSize(bA->isglobal.row[i],&sm);CHKERRQ(ierr);
23 ierr = ISGetSize(bA->isglobal.row[i],&sM);CHKERRQ(ierr);
24 *m += sm;
25 *M += sM;
26 }
27 for (j=0; j<bA->nc; j++) { /* cols */
28 PetscInt sn,sN;
29 ierr = ISGetLocalSize(bA->isglobal.col[j],&sn);CHKERRQ(ierr);
30 ierr = ISGetSize(bA->isglobal.col[j],&sN);CHKERRQ(ierr);
31 *n += sn;
32 *N += sN;
33 }
34 PetscFunctionReturn(0);
35 }
36
37 /* operations */
MatMult_Nest(Mat A,Vec x,Vec y)38 static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
39 {
40 Mat_Nest *bA = (Mat_Nest*)A->data;
41 Vec *bx = bA->right,*by = bA->left;
42 PetscInt i,j,nr = bA->nr,nc = bA->nc;
43 PetscErrorCode ierr;
44
45 PetscFunctionBegin;
46 for (i=0; i<nr; i++) {ierr = VecGetSubVector(y,bA->isglobal.row[i],&by[i]);CHKERRQ(ierr);}
47 for (i=0; i<nc; i++) {ierr = VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);}
48 for (i=0; i<nr; i++) {
49 ierr = VecZeroEntries(by[i]);CHKERRQ(ierr);
50 for (j=0; j<nc; j++) {
51 if (!bA->m[i][j]) continue;
52 /* y[i] <- y[i] + A[i][j] * x[j] */
53 ierr = MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);CHKERRQ(ierr);
54 }
55 }
56 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);CHKERRQ(ierr);}
57 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);}
58 PetscFunctionReturn(0);
59 }
60
MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)61 static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
62 {
63 Mat_Nest *bA = (Mat_Nest*)A->data;
64 Vec *bx = bA->right,*bz = bA->left;
65 PetscInt i,j,nr = bA->nr,nc = bA->nc;
66 PetscErrorCode ierr;
67
68 PetscFunctionBegin;
69 for (i=0; i<nr; i++) {ierr = VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);CHKERRQ(ierr);}
70 for (i=0; i<nc; i++) {ierr = VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);}
71 for (i=0; i<nr; i++) {
72 if (y != z) {
73 Vec by;
74 ierr = VecGetSubVector(y,bA->isglobal.row[i],&by);CHKERRQ(ierr);
75 ierr = VecCopy(by,bz[i]);CHKERRQ(ierr);
76 ierr = VecRestoreSubVector(y,bA->isglobal.row[i],&by);CHKERRQ(ierr);
77 }
78 for (j=0; j<nc; j++) {
79 if (!bA->m[i][j]) continue;
80 /* y[i] <- y[i] + A[i][j] * x[j] */
81 ierr = MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);CHKERRQ(ierr);
82 }
83 }
84 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);CHKERRQ(ierr);}
85 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);}
86 PetscFunctionReturn(0);
87 }
88
89 typedef struct {
90 Mat *workC; /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
91 PetscScalar *tarray; /* buffer for storing all temporary products A[i][j] B[j] */
92 PetscInt *dm,*dn,k; /* displacements and number of submatrices */
93 } Nest_Dense;
94
MatProductNumeric_Nest_Dense(Mat C)95 PETSC_INTERN PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
96 {
97 Mat_Nest *bA;
98 Nest_Dense *contents;
99 Mat viewB,viewC,productB,workC;
100 const PetscScalar *barray;
101 PetscScalar *carray;
102 PetscInt i,j,M,N,nr,nc,ldb,ldc;
103 PetscErrorCode ierr;
104 Mat A,B;
105
106 PetscFunctionBegin;
107 MatCheckProduct(C,3);
108 A = C->product->A;
109 B = C->product->B;
110 ierr = MatGetSize(B,NULL,&N);CHKERRQ(ierr);
111 if (!N) {
112 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
113 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
114 PetscFunctionReturn(0);
115 }
116 contents = (Nest_Dense*)C->product->data;
117 if (!contents) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
118 bA = (Mat_Nest*)A->data;
119 nr = bA->nr;
120 nc = bA->nc;
121 ierr = MatDenseGetLDA(B,&ldb);CHKERRQ(ierr);
122 ierr = MatDenseGetLDA(C,&ldc);CHKERRQ(ierr);
123 ierr = MatZeroEntries(C);CHKERRQ(ierr);
124 ierr = MatDenseGetArrayRead(B,&barray);CHKERRQ(ierr);
125 ierr = MatDenseGetArrayWrite(C,&carray);CHKERRQ(ierr);
126 for (i=0; i<nr; i++) {
127 ierr = ISGetSize(bA->isglobal.row[i],&M);CHKERRQ(ierr);
128 ierr = MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);CHKERRQ(ierr);
129 ierr = MatDenseSetLDA(viewC,ldc);CHKERRQ(ierr);
130 for (j=0; j<nc; j++) {
131 if (!bA->m[i][j]) continue;
132 ierr = ISGetSize(bA->isglobal.col[j],&M);CHKERRQ(ierr);
133 ierr = MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);CHKERRQ(ierr);
134 ierr = MatDenseSetLDA(viewB,ldb);CHKERRQ(ierr);
135
136 /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
137 workC = contents->workC[i*nc + j];
138 productB = workC->product->B;
139 workC->product->B = viewB; /* use newly created dense matrix viewB */
140 ierr = MatProductNumeric(workC);CHKERRQ(ierr);
141 ierr = MatDestroy(&viewB);CHKERRQ(ierr);
142 workC->product->B = productB; /* resume original B */
143
144 /* C[i] <- workC + C[i] */
145 ierr = MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);CHKERRQ(ierr);
146 }
147 ierr = MatDestroy(&viewC);CHKERRQ(ierr);
148 }
149 ierr = MatDenseRestoreArrayWrite(C,&carray);CHKERRQ(ierr);
150 ierr = MatDenseRestoreArrayRead(B,&barray);CHKERRQ(ierr);
151
152 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
153 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
154 PetscFunctionReturn(0);
155 }
156
MatNest_DenseDestroy(void * ctx)157 PetscErrorCode MatNest_DenseDestroy(void *ctx)
158 {
159 Nest_Dense *contents = (Nest_Dense*)ctx;
160 PetscInt i;
161 PetscErrorCode ierr;
162
163 PetscFunctionBegin;
164 ierr = PetscFree(contents->tarray);CHKERRQ(ierr);
165 for (i=0; i<contents->k; i++) {
166 ierr = MatDestroy(contents->workC + i);CHKERRQ(ierr);
167 }
168 ierr = PetscFree3(contents->dm,contents->dn,contents->workC);CHKERRQ(ierr);
169 ierr = PetscFree(contents);CHKERRQ(ierr);
170 PetscFunctionReturn(0);
171 }
172
MatProductSymbolic_Nest_Dense(Mat C)173 PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
174 {
175 Mat_Nest *bA;
176 Mat viewB,workC;
177 const PetscScalar *barray;
178 PetscInt i,j,M,N,m,n,nr,nc,maxm = 0,ldb;
179 Nest_Dense *contents=NULL;
180 PetscBool cisdense;
181 PetscErrorCode ierr;
182 Mat A,B;
183 PetscReal fill;
184
185 PetscFunctionBegin;
186 MatCheckProduct(C,4);
187 if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
188 A = C->product->A;
189 B = C->product->B;
190 fill = C->product->fill;
191 bA = (Mat_Nest*)A->data;
192 nr = bA->nr;
193 nc = bA->nc;
194 ierr = MatGetLocalSize(B,NULL,&n);CHKERRQ(ierr);
195 ierr = MatGetSize(B,NULL,&N);CHKERRQ(ierr);
196 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
197 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
198 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr);
199 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATMPIDENSE,MATSEQDENSECUDA,MATMPIDENSECUDA,"");CHKERRQ(ierr);
200 if (!cisdense) {
201 ierr = MatSetType(C,((PetscObject)B)->type_name);CHKERRQ(ierr);
202 }
203 ierr = MatSetUp(C);CHKERRQ(ierr);
204 if (!N) {
205 C->ops->productnumeric = MatProductNumeric_Nest_Dense;
206 PetscFunctionReturn(0);
207 }
208
209 ierr = PetscNew(&contents);CHKERRQ(ierr);
210 C->product->data = contents;
211 C->product->destroy = MatNest_DenseDestroy;
212 ierr = PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);CHKERRQ(ierr);
213 contents->k = nr*nc;
214 for (i=0; i<nr; i++) {
215 ierr = ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);CHKERRQ(ierr);
216 maxm = PetscMax(maxm,contents->dm[i+1]);
217 contents->dm[i+1] += contents->dm[i];
218 }
219 for (i=0; i<nc; i++) {
220 ierr = ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);CHKERRQ(ierr);
221 contents->dn[i+1] += contents->dn[i];
222 }
223 ierr = PetscMalloc1(maxm*N,&contents->tarray);CHKERRQ(ierr);
224 ierr = MatDenseGetLDA(B,&ldb);CHKERRQ(ierr);
225 ierr = MatGetSize(B,NULL,&N);CHKERRQ(ierr);
226 ierr = MatDenseGetArrayRead(B,&barray);CHKERRQ(ierr);
227 /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
228 for (j=0; j<nc; j++) {
229 ierr = ISGetSize(bA->isglobal.col[j],&M);CHKERRQ(ierr);
230 ierr = MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);CHKERRQ(ierr);
231 ierr = MatDenseSetLDA(viewB,ldb);CHKERRQ(ierr);
232 for (i=0; i<nr; i++) {
233 if (!bA->m[i][j]) continue;
234 /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */
235
236 ierr = MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);CHKERRQ(ierr);
237 workC = contents->workC[i*nc + j];
238 ierr = MatProductSetType(workC,MATPRODUCT_AB);CHKERRQ(ierr);
239 ierr = MatProductSetAlgorithm(workC,"default");CHKERRQ(ierr);
240 ierr = MatProductSetFill(workC,fill);CHKERRQ(ierr);
241 ierr = MatProductSetFromOptions(workC);CHKERRQ(ierr);
242 ierr = MatProductSymbolic(workC);CHKERRQ(ierr);
243
244 /* since tarray will be shared by all Mat */
245 ierr = MatSeqDenseSetPreallocation(workC,contents->tarray);CHKERRQ(ierr);
246 ierr = MatMPIDenseSetPreallocation(workC,contents->tarray);CHKERRQ(ierr);
247 }
248 ierr = MatDestroy(&viewB);CHKERRQ(ierr);
249 }
250 ierr = MatDenseRestoreArrayRead(B,&barray);CHKERRQ(ierr);
251
252 C->ops->productnumeric = MatProductNumeric_Nest_Dense;
253 PetscFunctionReturn(0);
254 }
255
256 /* --------------------------------------------------------- */
MatProductSetFromOptions_Nest_Dense_AB(Mat C)257 static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
258 {
259 PetscFunctionBegin;
260 C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
261 PetscFunctionReturn(0);
262 }
263
MatProductSetFromOptions_Nest_Dense(Mat C)264 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
265 {
266 PetscErrorCode ierr;
267 Mat_Product *product = C->product;
268
269 PetscFunctionBegin;
270 if (product->type == MATPRODUCT_AB) {
271 ierr = MatProductSetFromOptions_Nest_Dense_AB(C);CHKERRQ(ierr);
272 }
273 PetscFunctionReturn(0);
274 }
275 /* --------------------------------------------------------- */
276
MatMultTranspose_Nest(Mat A,Vec x,Vec y)277 static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
278 {
279 Mat_Nest *bA = (Mat_Nest*)A->data;
280 Vec *bx = bA->left,*by = bA->right;
281 PetscInt i,j,nr = bA->nr,nc = bA->nc;
282 PetscErrorCode ierr;
283
284 PetscFunctionBegin;
285 for (i=0; i<nr; i++) {ierr = VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);}
286 for (i=0; i<nc; i++) {ierr = VecGetSubVector(y,bA->isglobal.col[i],&by[i]);CHKERRQ(ierr);}
287 for (j=0; j<nc; j++) {
288 ierr = VecZeroEntries(by[j]);CHKERRQ(ierr);
289 for (i=0; i<nr; i++) {
290 if (!bA->m[i][j]) continue;
291 /* y[j] <- y[j] + (A[i][j])^T * x[i] */
292 ierr = MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);CHKERRQ(ierr);
293 }
294 }
295 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);}
296 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);CHKERRQ(ierr);}
297 PetscFunctionReturn(0);
298 }
299
MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)300 static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
301 {
302 Mat_Nest *bA = (Mat_Nest*)A->data;
303 Vec *bx = bA->left,*bz = bA->right;
304 PetscInt i,j,nr = bA->nr,nc = bA->nc;
305 PetscErrorCode ierr;
306
307 PetscFunctionBegin;
308 for (i=0; i<nr; i++) {ierr = VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);}
309 for (i=0; i<nc; i++) {ierr = VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);CHKERRQ(ierr);}
310 for (j=0; j<nc; j++) {
311 if (y != z) {
312 Vec by;
313 ierr = VecGetSubVector(y,bA->isglobal.col[j],&by);CHKERRQ(ierr);
314 ierr = VecCopy(by,bz[j]);CHKERRQ(ierr);
315 ierr = VecRestoreSubVector(y,bA->isglobal.col[j],&by);CHKERRQ(ierr);
316 }
317 for (i=0; i<nr; i++) {
318 if (!bA->m[i][j]) continue;
319 /* z[j] <- y[j] + (A[i][j])^T * x[i] */
320 ierr = MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);CHKERRQ(ierr);
321 }
322 }
323 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);}
324 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);CHKERRQ(ierr);}
325 PetscFunctionReturn(0);
326 }
327
MatTranspose_Nest(Mat A,MatReuse reuse,Mat * B)328 static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
329 {
330 Mat_Nest *bA = (Mat_Nest*)A->data, *bC;
331 Mat C;
332 PetscInt i,j,nr = bA->nr,nc = bA->nc;
333 PetscErrorCode ierr;
334
335 PetscFunctionBegin;
336 if (reuse == MAT_INPLACE_MATRIX && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place");
337
338 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
339 Mat *subs;
340 IS *is_row,*is_col;
341
342 ierr = PetscCalloc1(nr * nc,&subs);CHKERRQ(ierr);
343 ierr = PetscMalloc2(nr,&is_row,nc,&is_col);CHKERRQ(ierr);
344 ierr = MatNestGetISs(A,is_row,is_col);CHKERRQ(ierr);
345 if (reuse == MAT_INPLACE_MATRIX) {
346 for (i=0; i<nr; i++) {
347 for (j=0; j<nc; j++) {
348 subs[i + nr * j] = bA->m[i][j];
349 }
350 }
351 }
352
353 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);CHKERRQ(ierr);
354 ierr = PetscFree(subs);CHKERRQ(ierr);
355 ierr = PetscFree2(is_row,is_col);CHKERRQ(ierr);
356 } else {
357 C = *B;
358 }
359
360 bC = (Mat_Nest*)C->data;
361 for (i=0; i<nr; i++) {
362 for (j=0; j<nc; j++) {
363 if (bA->m[i][j]) {
364 ierr = MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));CHKERRQ(ierr);
365 } else {
366 bC->m[j][i] = NULL;
367 }
368 }
369 }
370
371 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
372 *B = C;
373 } else {
374 ierr = MatHeaderMerge(A, &C);CHKERRQ(ierr);
375 }
376 PetscFunctionReturn(0);
377 }
378
MatNestDestroyISList(PetscInt n,IS ** list)379 static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
380 {
381 PetscErrorCode ierr;
382 IS *lst = *list;
383 PetscInt i;
384
385 PetscFunctionBegin;
386 if (!lst) PetscFunctionReturn(0);
387 for (i=0; i<n; i++) if (lst[i]) {ierr = ISDestroy(&lst[i]);CHKERRQ(ierr);}
388 ierr = PetscFree(lst);CHKERRQ(ierr);
389 *list = NULL;
390 PetscFunctionReturn(0);
391 }
392
MatReset_Nest(Mat A)393 static PetscErrorCode MatReset_Nest(Mat A)
394 {
395 Mat_Nest *vs = (Mat_Nest*)A->data;
396 PetscInt i,j;
397 PetscErrorCode ierr;
398
399 PetscFunctionBegin;
400 /* release the matrices and the place holders */
401 ierr = MatNestDestroyISList(vs->nr,&vs->isglobal.row);CHKERRQ(ierr);
402 ierr = MatNestDestroyISList(vs->nc,&vs->isglobal.col);CHKERRQ(ierr);
403 ierr = MatNestDestroyISList(vs->nr,&vs->islocal.row);CHKERRQ(ierr);
404 ierr = MatNestDestroyISList(vs->nc,&vs->islocal.col);CHKERRQ(ierr);
405
406 ierr = PetscFree(vs->row_len);CHKERRQ(ierr);
407 ierr = PetscFree(vs->col_len);CHKERRQ(ierr);
408 ierr = PetscFree(vs->nnzstate);CHKERRQ(ierr);
409
410 ierr = PetscFree2(vs->left,vs->right);CHKERRQ(ierr);
411
412 /* release the matrices and the place holders */
413 if (vs->m) {
414 for (i=0; i<vs->nr; i++) {
415 for (j=0; j<vs->nc; j++) {
416 ierr = MatDestroy(&vs->m[i][j]);CHKERRQ(ierr);
417 }
418 ierr = PetscFree(vs->m[i]);CHKERRQ(ierr);
419 }
420 ierr = PetscFree(vs->m);CHKERRQ(ierr);
421 }
422
423 /* restore defaults */
424 vs->nr = 0;
425 vs->nc = 0;
426 vs->splitassembly = PETSC_FALSE;
427 PetscFunctionReturn(0);
428 }
429
MatDestroy_Nest(Mat A)430 static PetscErrorCode MatDestroy_Nest(Mat A)
431 {
432 PetscErrorCode ierr;
433
434 ierr = MatReset_Nest(A);CHKERRQ(ierr);
435 ierr = PetscFree(A->data);CHKERRQ(ierr);
436
437 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);CHKERRQ(ierr);
438 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);CHKERRQ(ierr);
439 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);CHKERRQ(ierr);
440 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);CHKERRQ(ierr);
441 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);CHKERRQ(ierr);
442 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);CHKERRQ(ierr);
443 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);CHKERRQ(ierr);
444 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);CHKERRQ(ierr);
445 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);CHKERRQ(ierr);
446 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);CHKERRQ(ierr);
447 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);CHKERRQ(ierr);
448 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);CHKERRQ(ierr);
449 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);CHKERRQ(ierr);
450 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);CHKERRQ(ierr);
451 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);CHKERRQ(ierr);
452 PetscFunctionReturn(0);
453 }
454
MatMissingDiagonal_Nest(Mat mat,PetscBool * missing,PetscInt * dd)455 static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
456 {
457 Mat_Nest *vs = (Mat_Nest*)mat->data;
458 PetscInt i;
459 PetscErrorCode ierr;
460
461 PetscFunctionBegin;
462 if (dd) *dd = 0;
463 if (!vs->nr) {
464 *missing = PETSC_TRUE;
465 PetscFunctionReturn(0);
466 }
467 *missing = PETSC_FALSE;
468 for (i = 0; i < vs->nr && !(*missing); i++) {
469 *missing = PETSC_TRUE;
470 if (vs->m[i][i]) {
471 ierr = MatMissingDiagonal(vs->m[i][i],missing,NULL);CHKERRQ(ierr);
472 if (*missing && dd) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"First missing entry not yet implemented");
473 }
474 }
475 PetscFunctionReturn(0);
476 }
477
MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)478 static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
479 {
480 Mat_Nest *vs = (Mat_Nest*)A->data;
481 PetscInt i,j;
482 PetscErrorCode ierr;
483 PetscBool nnzstate = PETSC_FALSE;
484
485 PetscFunctionBegin;
486 for (i=0; i<vs->nr; i++) {
487 for (j=0; j<vs->nc; j++) {
488 PetscObjectState subnnzstate = 0;
489 if (vs->m[i][j]) {
490 ierr = MatAssemblyBegin(vs->m[i][j],type);CHKERRQ(ierr);
491 if (!vs->splitassembly) {
492 /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
493 * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
494 * already performing an assembly, but the result would by more complicated and appears to offer less
495 * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
496 * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
497 */
498 ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr);
499 ierr = MatGetNonzeroState(vs->m[i][j],&subnnzstate);CHKERRQ(ierr);
500 }
501 }
502 nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
503 vs->nnzstate[i*vs->nc+j] = subnnzstate;
504 }
505 }
506 if (nnzstate) A->nonzerostate++;
507 PetscFunctionReturn(0);
508 }
509
MatAssemblyEnd_Nest(Mat A,MatAssemblyType type)510 static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
511 {
512 Mat_Nest *vs = (Mat_Nest*)A->data;
513 PetscInt i,j;
514 PetscErrorCode ierr;
515
516 PetscFunctionBegin;
517 for (i=0; i<vs->nr; i++) {
518 for (j=0; j<vs->nc; j++) {
519 if (vs->m[i][j]) {
520 if (vs->splitassembly) {
521 ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr);
522 }
523 }
524 }
525 }
526 PetscFunctionReturn(0);
527 }
528
MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat * B)529 static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
530 {
531 PetscErrorCode ierr;
532 Mat_Nest *vs = (Mat_Nest*)A->data;
533 PetscInt j;
534 Mat sub;
535
536 PetscFunctionBegin;
537 sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
538 for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
539 if (sub) {ierr = MatSetUp(sub);CHKERRQ(ierr);} /* Ensure that the sizes are available */
540 *B = sub;
541 PetscFunctionReturn(0);
542 }
543
MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat * B)544 static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
545 {
546 PetscErrorCode ierr;
547 Mat_Nest *vs = (Mat_Nest*)A->data;
548 PetscInt i;
549 Mat sub;
550
551 PetscFunctionBegin;
552 sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
553 for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
554 if (sub) {ierr = MatSetUp(sub);CHKERRQ(ierr);} /* Ensure that the sizes are available */
555 *B = sub;
556 PetscFunctionReturn(0);
557 }
558
MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt * found)559 static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
560 {
561 PetscErrorCode ierr;
562 PetscInt i;
563 PetscBool flg;
564
565 PetscFunctionBegin;
566 PetscValidPointer(list,3);
567 PetscValidHeaderSpecific(is,IS_CLASSID,4);
568 PetscValidIntPointer(found,5);
569 *found = -1;
570 for (i=0; i<n; i++) {
571 if (!list[i]) continue;
572 ierr = ISEqualUnsorted(list[i],is,&flg);CHKERRQ(ierr);
573 if (flg) {
574 *found = i;
575 PetscFunctionReturn(0);
576 }
577 }
578 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
579 PetscFunctionReturn(0);
580 }
581
582 /* Get a block row as a new MatNest */
MatNestGetRow(Mat A,PetscInt row,Mat * B)583 static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
584 {
585 Mat_Nest *vs = (Mat_Nest*)A->data;
586 char keyname[256];
587 PetscErrorCode ierr;
588
589 PetscFunctionBegin;
590 *B = NULL;
591 ierr = PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);CHKERRQ(ierr);
592 ierr = PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);CHKERRQ(ierr);
593 if (*B) PetscFunctionReturn(0);
594
595 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);CHKERRQ(ierr);
596
597 (*B)->assembled = A->assembled;
598
599 ierr = PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);CHKERRQ(ierr);
600 ierr = PetscObjectDereference((PetscObject)*B);CHKERRQ(ierr); /* Leave the only remaining reference in the composition */
601 PetscFunctionReturn(0);
602 }
603
MatNestFindSubMat(Mat A,struct MatNestISPair * is,IS isrow,IS iscol,Mat * B)604 static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
605 {
606 Mat_Nest *vs = (Mat_Nest*)A->data;
607 PetscErrorCode ierr;
608 PetscInt row,col;
609 PetscBool same,isFullCol,isFullColGlobal;
610
611 PetscFunctionBegin;
612 /* Check if full column space. This is a hack */
613 isFullCol = PETSC_FALSE;
614 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);CHKERRQ(ierr);
615 if (same) {
616 PetscInt n,first,step,i,an,am,afirst,astep;
617 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr);
618 ierr = ISGetLocalSize(iscol,&n);CHKERRQ(ierr);
619 isFullCol = PETSC_TRUE;
620 for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
621 ierr = PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);CHKERRQ(ierr);
622 ierr = ISGetLocalSize(is->col[i],&am);CHKERRQ(ierr);
623 if (same) {
624 ierr = ISStrideGetInfo(is->col[i],&afirst,&astep);CHKERRQ(ierr);
625 if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
626 } else isFullCol = PETSC_FALSE;
627 an += am;
628 }
629 if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
630 }
631 ierr = MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));CHKERRQ(ierr);
632
633 if (isFullColGlobal && vs->nc > 1) {
634 PetscInt row;
635 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr);
636 ierr = MatNestGetRow(A,row,B);CHKERRQ(ierr);
637 } else {
638 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr);
639 ierr = MatNestFindIS(A,vs->nc,is->col,iscol,&col);CHKERRQ(ierr);
640 if (!vs->m[row][col]) {
641 PetscInt lr,lc;
642
643 ierr = MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);CHKERRQ(ierr);
644 ierr = ISGetLocalSize(vs->isglobal.row[row],&lr);CHKERRQ(ierr);
645 ierr = ISGetLocalSize(vs->isglobal.col[col],&lc);CHKERRQ(ierr);
646 ierr = MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
647 ierr = MatSetType(vs->m[row][col],MATAIJ);CHKERRQ(ierr);
648 ierr = MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);CHKERRQ(ierr);
649 ierr = MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);CHKERRQ(ierr);
650 ierr = MatSetUp(vs->m[row][col]);CHKERRQ(ierr);
651 ierr = MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
652 ierr = MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
653 }
654 *B = vs->m[row][col];
655 }
656 PetscFunctionReturn(0);
657 }
658
659 /*
660 TODO: This does not actually returns a submatrix we can modify
661 */
MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat * B)662 static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
663 {
664 PetscErrorCode ierr;
665 Mat_Nest *vs = (Mat_Nest*)A->data;
666 Mat sub;
667
668 PetscFunctionBegin;
669 ierr = MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);CHKERRQ(ierr);
670 switch (reuse) {
671 case MAT_INITIAL_MATRIX:
672 if (sub) { ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr); }
673 *B = sub;
674 break;
675 case MAT_REUSE_MATRIX:
676 if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
677 break;
678 case MAT_IGNORE_MATRIX: /* Nothing to do */
679 break;
680 case MAT_INPLACE_MATRIX: /* Nothing to do */
681 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
682 break;
683 }
684 PetscFunctionReturn(0);
685 }
686
MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat * B)687 PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
688 {
689 PetscErrorCode ierr;
690 Mat_Nest *vs = (Mat_Nest*)A->data;
691 Mat sub;
692
693 PetscFunctionBegin;
694 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr);
695 /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
696 if (sub) {ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr);}
697 *B = sub;
698 PetscFunctionReturn(0);
699 }
700
MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat * B)701 static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
702 {
703 PetscErrorCode ierr;
704 Mat_Nest *vs = (Mat_Nest*)A->data;
705 Mat sub;
706
707 PetscFunctionBegin;
708 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr);
709 if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
710 if (sub) {
711 if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
712 ierr = MatDestroy(B);CHKERRQ(ierr);
713 }
714 PetscFunctionReturn(0);
715 }
716
MatGetDiagonal_Nest(Mat A,Vec v)717 static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
718 {
719 Mat_Nest *bA = (Mat_Nest*)A->data;
720 PetscInt i;
721 PetscErrorCode ierr;
722
723 PetscFunctionBegin;
724 for (i=0; i<bA->nr; i++) {
725 Vec bv;
726 ierr = VecGetSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr);
727 if (bA->m[i][i]) {
728 ierr = MatGetDiagonal(bA->m[i][i],bv);CHKERRQ(ierr);
729 } else {
730 ierr = VecSet(bv,0.0);CHKERRQ(ierr);
731 }
732 ierr = VecRestoreSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr);
733 }
734 PetscFunctionReturn(0);
735 }
736
MatDiagonalScale_Nest(Mat A,Vec l,Vec r)737 static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
738 {
739 Mat_Nest *bA = (Mat_Nest*)A->data;
740 Vec bl,*br;
741 PetscInt i,j;
742 PetscErrorCode ierr;
743
744 PetscFunctionBegin;
745 ierr = PetscCalloc1(bA->nc,&br);CHKERRQ(ierr);
746 if (r) {
747 for (j=0; j<bA->nc; j++) {ierr = VecGetSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);}
748 }
749 bl = NULL;
750 for (i=0; i<bA->nr; i++) {
751 if (l) {
752 ierr = VecGetSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr);
753 }
754 for (j=0; j<bA->nc; j++) {
755 if (bA->m[i][j]) {
756 ierr = MatDiagonalScale(bA->m[i][j],bl,br[j]);CHKERRQ(ierr);
757 }
758 }
759 if (l) {
760 ierr = VecRestoreSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr);
761 }
762 }
763 if (r) {
764 for (j=0; j<bA->nc; j++) {ierr = VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);}
765 }
766 ierr = PetscFree(br);CHKERRQ(ierr);
767 PetscFunctionReturn(0);
768 }
769
MatScale_Nest(Mat A,PetscScalar a)770 static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
771 {
772 Mat_Nest *bA = (Mat_Nest*)A->data;
773 PetscInt i,j;
774 PetscErrorCode ierr;
775
776 PetscFunctionBegin;
777 for (i=0; i<bA->nr; i++) {
778 for (j=0; j<bA->nc; j++) {
779 if (bA->m[i][j]) {
780 ierr = MatScale(bA->m[i][j],a);CHKERRQ(ierr);
781 }
782 }
783 }
784 PetscFunctionReturn(0);
785 }
786
MatShift_Nest(Mat A,PetscScalar a)787 static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
788 {
789 Mat_Nest *bA = (Mat_Nest*)A->data;
790 PetscInt i;
791 PetscErrorCode ierr;
792 PetscBool nnzstate = PETSC_FALSE;
793
794 PetscFunctionBegin;
795 for (i=0; i<bA->nr; i++) {
796 PetscObjectState subnnzstate = 0;
797 if (!bA->m[i][i]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for shifting an empty diagonal block, insert a matrix in block (%D,%D)",i,i);
798 ierr = MatShift(bA->m[i][i],a);CHKERRQ(ierr);
799 ierr = MatGetNonzeroState(bA->m[i][i],&subnnzstate);CHKERRQ(ierr);
800 nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
801 bA->nnzstate[i*bA->nc+i] = subnnzstate;
802 }
803 if (nnzstate) A->nonzerostate++;
804 PetscFunctionReturn(0);
805 }
806
MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)807 static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
808 {
809 Mat_Nest *bA = (Mat_Nest*)A->data;
810 PetscInt i;
811 PetscErrorCode ierr;
812 PetscBool nnzstate = PETSC_FALSE;
813
814 PetscFunctionBegin;
815 for (i=0; i<bA->nr; i++) {
816 PetscObjectState subnnzstate = 0;
817 Vec bv;
818 ierr = VecGetSubVector(D,bA->isglobal.row[i],&bv);CHKERRQ(ierr);
819 if (bA->m[i][i]) {
820 ierr = MatDiagonalSet(bA->m[i][i],bv,is);CHKERRQ(ierr);
821 ierr = MatGetNonzeroState(bA->m[i][i],&subnnzstate);CHKERRQ(ierr);
822 }
823 ierr = VecRestoreSubVector(D,bA->isglobal.row[i],&bv);CHKERRQ(ierr);
824 nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
825 bA->nnzstate[i*bA->nc+i] = subnnzstate;
826 }
827 if (nnzstate) A->nonzerostate++;
828 PetscFunctionReturn(0);
829 }
830
MatSetRandom_Nest(Mat A,PetscRandom rctx)831 static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
832 {
833 Mat_Nest *bA = (Mat_Nest*)A->data;
834 PetscInt i,j;
835 PetscErrorCode ierr;
836
837 PetscFunctionBegin;
838 for (i=0; i<bA->nr; i++) {
839 for (j=0; j<bA->nc; j++) {
840 if (bA->m[i][j]) {
841 ierr = MatSetRandom(bA->m[i][j],rctx);CHKERRQ(ierr);
842 }
843 }
844 }
845 PetscFunctionReturn(0);
846 }
847
MatCreateVecs_Nest(Mat A,Vec * right,Vec * left)848 static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
849 {
850 Mat_Nest *bA = (Mat_Nest*)A->data;
851 Vec *L,*R;
852 MPI_Comm comm;
853 PetscInt i,j;
854 PetscErrorCode ierr;
855
856 PetscFunctionBegin;
857 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
858 if (right) {
859 /* allocate R */
860 ierr = PetscMalloc1(bA->nc, &R);CHKERRQ(ierr);
861 /* Create the right vectors */
862 for (j=0; j<bA->nc; j++) {
863 for (i=0; i<bA->nr; i++) {
864 if (bA->m[i][j]) {
865 ierr = MatCreateVecs(bA->m[i][j],&R[j],NULL);CHKERRQ(ierr);
866 break;
867 }
868 }
869 if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
870 }
871 ierr = VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);CHKERRQ(ierr);
872 /* hand back control to the nest vector */
873 for (j=0; j<bA->nc; j++) {
874 ierr = VecDestroy(&R[j]);CHKERRQ(ierr);
875 }
876 ierr = PetscFree(R);CHKERRQ(ierr);
877 }
878
879 if (left) {
880 /* allocate L */
881 ierr = PetscMalloc1(bA->nr, &L);CHKERRQ(ierr);
882 /* Create the left vectors */
883 for (i=0; i<bA->nr; i++) {
884 for (j=0; j<bA->nc; j++) {
885 if (bA->m[i][j]) {
886 ierr = MatCreateVecs(bA->m[i][j],NULL,&L[i]);CHKERRQ(ierr);
887 break;
888 }
889 }
890 if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
891 }
892
893 ierr = VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);CHKERRQ(ierr);
894 for (i=0; i<bA->nr; i++) {
895 ierr = VecDestroy(&L[i]);CHKERRQ(ierr);
896 }
897
898 ierr = PetscFree(L);CHKERRQ(ierr);
899 }
900 PetscFunctionReturn(0);
901 }
902
MatView_Nest(Mat A,PetscViewer viewer)903 static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
904 {
905 Mat_Nest *bA = (Mat_Nest*)A->data;
906 PetscBool isascii,viewSub = PETSC_FALSE;
907 PetscInt i,j;
908 PetscErrorCode ierr;
909
910 PetscFunctionBegin;
911 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
912 if (isascii) {
913
914 ierr = PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);CHKERRQ(ierr);
915 ierr = PetscViewerASCIIPrintf(viewer,"Matrix object: \n");CHKERRQ(ierr);
916 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
917 ierr = PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);CHKERRQ(ierr);
918
919 ierr = PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");CHKERRQ(ierr);
920 for (i=0; i<bA->nr; i++) {
921 for (j=0; j<bA->nc; j++) {
922 MatType type;
923 char name[256] = "",prefix[256] = "";
924 PetscInt NR,NC;
925 PetscBool isNest = PETSC_FALSE;
926
927 if (!bA->m[i][j]) {
928 CHKERRQ(ierr);PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);CHKERRQ(ierr);
929 continue;
930 }
931 ierr = MatGetSize(bA->m[i][j],&NR,&NC);CHKERRQ(ierr);
932 ierr = MatGetType(bA->m[i][j], &type);CHKERRQ(ierr);
933 if (((PetscObject)bA->m[i][j])->name) {ierr = PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);CHKERRQ(ierr);}
934 if (((PetscObject)bA->m[i][j])->prefix) {ierr = PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);CHKERRQ(ierr);}
935 ierr = PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);CHKERRQ(ierr);
936
937 ierr = PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);CHKERRQ(ierr);
938
939 if (isNest || viewSub) {
940 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); /* push1 */
941 ierr = MatView(bA->m[i][j],viewer);CHKERRQ(ierr);
942 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop1 */
943 }
944 }
945 }
946 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop0 */
947 }
948 PetscFunctionReturn(0);
949 }
950
MatZeroEntries_Nest(Mat A)951 static PetscErrorCode MatZeroEntries_Nest(Mat A)
952 {
953 Mat_Nest *bA = (Mat_Nest*)A->data;
954 PetscInt i,j;
955 PetscErrorCode ierr;
956
957 PetscFunctionBegin;
958 for (i=0; i<bA->nr; i++) {
959 for (j=0; j<bA->nc; j++) {
960 if (!bA->m[i][j]) continue;
961 ierr = MatZeroEntries(bA->m[i][j]);CHKERRQ(ierr);
962 }
963 }
964 PetscFunctionReturn(0);
965 }
966
MatCopy_Nest(Mat A,Mat B,MatStructure str)967 static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
968 {
969 Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
970 PetscInt i,j,nr = bA->nr,nc = bA->nc;
971 PetscErrorCode ierr;
972 PetscBool nnzstate = PETSC_FALSE;
973
974 PetscFunctionBegin;
975 if (nr != bB->nr || nc != bB->nc) SETERRQ4(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot copy a Mat_Nest of block size (%D,%D) to a Mat_Nest of block size (%D,%D)",bB->nr,bB->nc,nr,nc);
976 for (i=0; i<nr; i++) {
977 for (j=0; j<nc; j++) {
978 PetscObjectState subnnzstate = 0;
979 if (bA->m[i][j] && bB->m[i][j]) {
980 ierr = MatCopy(bA->m[i][j],bB->m[i][j],str);CHKERRQ(ierr);
981 } else if (bA->m[i][j] || bB->m[i][j]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D",i,j);
982 ierr = MatGetNonzeroState(bB->m[i][j],&subnnzstate);CHKERRQ(ierr);
983 nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
984 bB->nnzstate[i*nc+j] = subnnzstate;
985 }
986 }
987 if (nnzstate) B->nonzerostate++;
988 PetscFunctionReturn(0);
989 }
990
MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)991 static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
992 {
993 Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
994 PetscInt i,j,nr = bY->nr,nc = bY->nc;
995 PetscErrorCode ierr;
996 PetscBool nnzstate = PETSC_FALSE;
997
998 PetscFunctionBegin;
999 if (nr != bX->nr || nc != bX->nc) SETERRQ4(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Cannot AXPY a MatNest of block size (%D,%D) with a MatNest of block size (%D,%D)",bX->nr,bX->nc,nr,nc);
1000 for (i=0; i<nr; i++) {
1001 for (j=0; j<nc; j++) {
1002 PetscObjectState subnnzstate = 0;
1003 if (bY->m[i][j] && bX->m[i][j]) {
1004 ierr = MatAXPY(bY->m[i][j],a,bX->m[i][j],str);CHKERRQ(ierr);
1005 } else if (bX->m[i][j]) {
1006 Mat M;
1007
1008 if (str != DIFFERENT_NONZERO_PATTERN) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D. Use DIFFERENT_NONZERO_PATTERN",i,j);
1009 ierr = MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);CHKERRQ(ierr);
1010 ierr = MatNestSetSubMat(Y,i,j,M);CHKERRQ(ierr);
1011 ierr = MatDestroy(&M);CHKERRQ(ierr);
1012 }
1013 if (bY->m[i][j]) { ierr = MatGetNonzeroState(bY->m[i][j],&subnnzstate);CHKERRQ(ierr); }
1014 nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
1015 bY->nnzstate[i*nc+j] = subnnzstate;
1016 }
1017 }
1018 if (nnzstate) Y->nonzerostate++;
1019 PetscFunctionReturn(0);
1020 }
1021
MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat * B)1022 static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1023 {
1024 Mat_Nest *bA = (Mat_Nest*)A->data;
1025 Mat *b;
1026 PetscInt i,j,nr = bA->nr,nc = bA->nc;
1027 PetscErrorCode ierr;
1028
1029 PetscFunctionBegin;
1030 ierr = PetscMalloc1(nr*nc,&b);CHKERRQ(ierr);
1031 for (i=0; i<nr; i++) {
1032 for (j=0; j<nc; j++) {
1033 if (bA->m[i][j]) {
1034 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr);
1035 } else {
1036 b[i*nc+j] = NULL;
1037 }
1038 }
1039 }
1040 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr);
1041 /* Give the new MatNest exclusive ownership */
1042 for (i=0; i<nr*nc; i++) {
1043 ierr = MatDestroy(&b[i]);CHKERRQ(ierr);
1044 }
1045 ierr = PetscFree(b);CHKERRQ(ierr);
1046
1047 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1048 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1049 PetscFunctionReturn(0);
1050 }
1051
1052 /* nest api */
MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat * mat)1053 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1054 {
1055 Mat_Nest *bA = (Mat_Nest*)A->data;
1056
1057 PetscFunctionBegin;
1058 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1059 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1060 *mat = bA->m[idxm][jdxm];
1061 PetscFunctionReturn(0);
1062 }
1063
1064 /*@
1065 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.
1066
1067 Not collective
1068
1069 Input Parameters:
1070 + A - nest matrix
1071 . idxm - index of the matrix within the nest matrix
1072 - jdxm - index of the matrix within the nest matrix
1073
1074 Output Parameter:
1075 . sub - matrix at index idxm,jdxm within the nest matrix
1076
1077 Level: developer
1078
1079 .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1080 MatNestGetLocalISs(), MatNestGetISs()
1081 @*/
MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat * sub)1082 PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1083 {
1084 PetscErrorCode ierr;
1085
1086 PetscFunctionBegin;
1087 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr);
1088 PetscFunctionReturn(0);
1089 }
1090
MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)1091 PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1092 {
1093 Mat_Nest *bA = (Mat_Nest*)A->data;
1094 PetscInt m,n,M,N,mi,ni,Mi,Ni;
1095 PetscErrorCode ierr;
1096
1097 PetscFunctionBegin;
1098 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1099 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1100 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr);
1101 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
1102 ierr = ISGetLocalSize(bA->isglobal.row[idxm],&mi);CHKERRQ(ierr);
1103 ierr = ISGetSize(bA->isglobal.row[idxm],&Mi);CHKERRQ(ierr);
1104 ierr = ISGetLocalSize(bA->isglobal.col[jdxm],&ni);CHKERRQ(ierr);
1105 ierr = ISGetSize(bA->isglobal.col[jdxm],&Ni);CHKERRQ(ierr);
1106 if (M != Mi || N != Ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix dimension (%D,%D) incompatible with nest block (%D,%D)",M,N,Mi,Ni);
1107 if (m != mi || n != ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix local dimension (%D,%D) incompatible with nest block (%D,%D)",m,n,mi,ni);
1108
1109 /* do not increase object state */
1110 if (mat == bA->m[idxm][jdxm]) PetscFunctionReturn(0);
1111
1112 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr);
1113 ierr = MatDestroy(&bA->m[idxm][jdxm]);CHKERRQ(ierr);
1114 bA->m[idxm][jdxm] = mat;
1115 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr);
1116 ierr = MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);CHKERRQ(ierr);
1117 A->nonzerostate++;
1118 PetscFunctionReturn(0);
1119 }
1120
1121 /*@
1122 MatNestSetSubMat - Set a single submatrix in the nest matrix.
1123
1124 Logically collective on the submatrix communicator
1125
1126 Input Parameters:
1127 + A - nest matrix
1128 . idxm - index of the matrix within the nest matrix
1129 . jdxm - index of the matrix within the nest matrix
1130 - sub - matrix at index idxm,jdxm within the nest matrix
1131
1132 Notes:
1133 The new submatrix must have the same size and communicator as that block of the nest.
1134
1135 This increments the reference count of the submatrix.
1136
1137 Level: developer
1138
1139 .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1140 MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1141 @*/
MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)1142 PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1143 {
1144 PetscErrorCode ierr;
1145
1146 PetscFunctionBegin;
1147 ierr = PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));CHKERRQ(ierr);
1148 PetscFunctionReturn(0);
1149 }
1150
MatNestGetSubMats_Nest(Mat A,PetscInt * M,PetscInt * N,Mat *** mat)1151 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1152 {
1153 Mat_Nest *bA = (Mat_Nest*)A->data;
1154
1155 PetscFunctionBegin;
1156 if (M) *M = bA->nr;
1157 if (N) *N = bA->nc;
1158 if (mat) *mat = bA->m;
1159 PetscFunctionReturn(0);
1160 }
1161
1162 /*@C
1163 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.
1164
1165 Not collective
1166
1167 Input Parameters:
1168 . A - nest matrix
1169
1170 Output Parameter:
1171 + M - number of rows in the nest matrix
1172 . N - number of cols in the nest matrix
1173 - mat - 2d array of matrices
1174
1175 Notes:
1176
1177 The user should not free the array mat.
1178
1179 In Fortran, this routine has a calling sequence
1180 $ call MatNestGetSubMats(A, M, N, mat, ierr)
1181 where the space allocated for the optional argument mat is assumed large enough (if provided).
1182
1183 Level: developer
1184
1185 .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1186 MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1187 @*/
MatNestGetSubMats(Mat A,PetscInt * M,PetscInt * N,Mat *** mat)1188 PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1189 {
1190 PetscErrorCode ierr;
1191
1192 PetscFunctionBegin;
1193 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr);
1194 PetscFunctionReturn(0);
1195 }
1196
MatNestGetSize_Nest(Mat A,PetscInt * M,PetscInt * N)1197 PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1198 {
1199 Mat_Nest *bA = (Mat_Nest*)A->data;
1200
1201 PetscFunctionBegin;
1202 if (M) *M = bA->nr;
1203 if (N) *N = bA->nc;
1204 PetscFunctionReturn(0);
1205 }
1206
1207 /*@
1208 MatNestGetSize - Returns the size of the nest matrix.
1209
1210 Not collective
1211
1212 Input Parameters:
1213 . A - nest matrix
1214
1215 Output Parameter:
1216 + M - number of rows in the nested mat
1217 - N - number of cols in the nested mat
1218
1219 Notes:
1220
1221 Level: developer
1222
1223 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1224 MatNestGetISs()
1225 @*/
MatNestGetSize(Mat A,PetscInt * M,PetscInt * N)1226 PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1227 {
1228 PetscErrorCode ierr;
1229
1230 PetscFunctionBegin;
1231 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr);
1232 PetscFunctionReturn(0);
1233 }
1234
MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])1235 static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1236 {
1237 Mat_Nest *vs = (Mat_Nest*)A->data;
1238 PetscInt i;
1239
1240 PetscFunctionBegin;
1241 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1242 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1243 PetscFunctionReturn(0);
1244 }
1245
1246 /*@C
1247 MatNestGetISs - Returns the index sets partitioning the row and column spaces
1248
1249 Not collective
1250
1251 Input Parameters:
1252 . A - nest matrix
1253
1254 Output Parameter:
1255 + rows - array of row index sets
1256 - cols - array of column index sets
1257
1258 Level: advanced
1259
1260 Notes:
1261 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1262
1263 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1264 MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1265 @*/
MatNestGetISs(Mat A,IS rows[],IS cols[])1266 PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[])
1267 {
1268 PetscErrorCode ierr;
1269
1270 PetscFunctionBegin;
1271 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1272 ierr = PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr);
1273 PetscFunctionReturn(0);
1274 }
1275
MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])1276 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1277 {
1278 Mat_Nest *vs = (Mat_Nest*)A->data;
1279 PetscInt i;
1280
1281 PetscFunctionBegin;
1282 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1283 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1284 PetscFunctionReturn(0);
1285 }
1286
1287 /*@C
1288 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces
1289
1290 Not collective
1291
1292 Input Parameters:
1293 . A - nest matrix
1294
1295 Output Parameter:
1296 + rows - array of row index sets (or NULL to ignore)
1297 - cols - array of column index sets (or NULL to ignore)
1298
1299 Level: advanced
1300
1301 Notes:
1302 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.
1303
1304 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1305 MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1306 @*/
MatNestGetLocalISs(Mat A,IS rows[],IS cols[])1307 PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1308 {
1309 PetscErrorCode ierr;
1310
1311 PetscFunctionBegin;
1312 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1313 ierr = PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr);
1314 PetscFunctionReturn(0);
1315 }
1316
MatNestSetVecType_Nest(Mat A,VecType vtype)1317 PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype)
1318 {
1319 PetscErrorCode ierr;
1320 PetscBool flg;
1321
1322 PetscFunctionBegin;
1323 ierr = PetscStrcmp(vtype,VECNEST,&flg);CHKERRQ(ierr);
1324 /* In reality, this only distinguishes VECNEST and "other" */
1325 if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1326 else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1327 PetscFunctionReturn(0);
1328 }
1329
1330 /*@C
1331 MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()
1332
1333 Not collective
1334
1335 Input Parameters:
1336 + A - nest matrix
1337 - vtype - type to use for creating vectors
1338
1339 Notes:
1340
1341 Level: developer
1342
1343 .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1344 @*/
MatNestSetVecType(Mat A,VecType vtype)1345 PetscErrorCode MatNestSetVecType(Mat A,VecType vtype)
1346 {
1347 PetscErrorCode ierr;
1348
1349 PetscFunctionBegin;
1350 ierr = PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));CHKERRQ(ierr);
1351 PetscFunctionReturn(0);
1352 }
1353
MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])1354 PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1355 {
1356 Mat_Nest *s = (Mat_Nest*)A->data;
1357 PetscInt i,j,m,n,M,N;
1358 PetscErrorCode ierr;
1359 PetscBool cong;
1360
1361 PetscFunctionBegin;
1362 ierr = MatReset_Nest(A);CHKERRQ(ierr);
1363
1364 s->nr = nr;
1365 s->nc = nc;
1366
1367 /* Create space for submatrices */
1368 ierr = PetscMalloc1(nr,&s->m);CHKERRQ(ierr);
1369 for (i=0; i<nr; i++) {
1370 ierr = PetscMalloc1(nc,&s->m[i]);CHKERRQ(ierr);
1371 }
1372 for (i=0; i<nr; i++) {
1373 for (j=0; j<nc; j++) {
1374 s->m[i][j] = a[i*nc+j];
1375 if (a[i*nc+j]) {
1376 ierr = PetscObjectReference((PetscObject)a[i*nc+j]);CHKERRQ(ierr);
1377 }
1378 }
1379 }
1380
1381 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr);
1382
1383 ierr = PetscMalloc1(nr,&s->row_len);CHKERRQ(ierr);
1384 ierr = PetscMalloc1(nc,&s->col_len);CHKERRQ(ierr);
1385 for (i=0; i<nr; i++) s->row_len[i]=-1;
1386 for (j=0; j<nc; j++) s->col_len[j]=-1;
1387
1388 ierr = PetscCalloc1(nr*nc,&s->nnzstate);CHKERRQ(ierr);
1389 for (i=0; i<nr; i++) {
1390 for (j=0; j<nc; j++) {
1391 if (s->m[i][j]) {
1392 ierr = MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);CHKERRQ(ierr);
1393 }
1394 }
1395 }
1396
1397 ierr = MatNestGetSizes_Private(A,&m,&n,&M,&N);CHKERRQ(ierr);
1398
1399 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr);
1400 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr);
1401 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr);
1402 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr);
1403
1404 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
1405 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
1406
1407 /* disable operations that are not supported for non-square matrices,
1408 or matrices for which is_row != is_col */
1409 ierr = MatHasCongruentLayouts(A,&cong);CHKERRQ(ierr);
1410 if (cong && nr != nc) cong = PETSC_FALSE;
1411 if (cong) {
1412 for (i = 0; cong && i < nr; i++) {
1413 ierr = ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);CHKERRQ(ierr);
1414 }
1415 }
1416 if (!cong) {
1417 A->ops->missingdiagonal = NULL;
1418 A->ops->getdiagonal = NULL;
1419 A->ops->shift = NULL;
1420 A->ops->diagonalset = NULL;
1421 }
1422
1423 ierr = PetscCalloc2(nr,&s->left,nc,&s->right);CHKERRQ(ierr);
1424 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr);
1425 A->nonzerostate++;
1426 PetscFunctionReturn(0);
1427 }
1428
1429 /*@
1430 MatNestSetSubMats - Sets the nested submatrices
1431
1432 Collective on Mat
1433
1434 Input Parameter:
1435 + A - nested matrix
1436 . nr - number of nested row blocks
1437 . is_row - index sets for each nested row block, or NULL to make contiguous
1438 . nc - number of nested column blocks
1439 . is_col - index sets for each nested column block, or NULL to make contiguous
1440 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1441
1442 Notes: this always resets any submatrix information previously set
1443
1444 Level: advanced
1445
1446 .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1447 @*/
MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])1448 PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1449 {
1450 PetscErrorCode ierr;
1451 PetscInt i;
1452
1453 PetscFunctionBegin;
1454 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1455 if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1456 if (nr && is_row) {
1457 PetscValidPointer(is_row,3);
1458 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3);
1459 }
1460 if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1461 if (nc && is_col) {
1462 PetscValidPointer(is_col,5);
1463 for (i=0; i<nc; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5);
1464 }
1465 if (nr*nc > 0) PetscValidPointer(a,6);
1466 ierr = PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));CHKERRQ(ierr);
1467 PetscFunctionReturn(0);
1468 }
1469
MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping * ltog)1470 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1471 {
1472 PetscErrorCode ierr;
1473 PetscBool flg;
1474 PetscInt i,j,m,mi,*ix;
1475
1476 PetscFunctionBegin;
1477 *ltog = NULL;
1478 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1479 if (islocal[i]) {
1480 ierr = ISGetLocalSize(islocal[i],&mi);CHKERRQ(ierr);
1481 flg = PETSC_TRUE; /* We found a non-trivial entry */
1482 } else {
1483 ierr = ISGetLocalSize(isglobal[i],&mi);CHKERRQ(ierr);
1484 }
1485 m += mi;
1486 }
1487 if (!flg) PetscFunctionReturn(0);
1488
1489 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr);
1490 for (i=0,m=0; i<n; i++) {
1491 ISLocalToGlobalMapping smap = NULL;
1492 Mat sub = NULL;
1493 PetscSF sf;
1494 PetscLayout map;
1495 const PetscInt *ix2;
1496
1497 if (!colflg) {
1498 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr);
1499 } else {
1500 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr);
1501 }
1502 if (sub) {
1503 if (!colflg) {
1504 ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr);
1505 } else {
1506 ierr = MatGetLocalToGlobalMapping(sub,NULL,&smap);CHKERRQ(ierr);
1507 }
1508 }
1509 /*
1510 Now we need to extract the monolithic global indices that correspond to the given split global indices.
1511 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1512 */
1513 ierr = ISGetIndices(isglobal[i],&ix2);CHKERRQ(ierr);
1514 if (islocal[i]) {
1515 PetscInt *ilocal,*iremote;
1516 PetscInt mil,nleaves;
1517
1518 ierr = ISGetLocalSize(islocal[i],&mi);CHKERRQ(ierr);
1519 if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map");
1520 for (j=0; j<mi; j++) ix[m+j] = j;
1521 ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr);
1522
1523 /* PetscSFSetGraphLayout does not like negative indices */
1524 ierr = PetscMalloc2(mi,&ilocal,mi,&iremote);CHKERRQ(ierr);
1525 for (j=0, nleaves = 0; j<mi; j++) {
1526 if (ix[m+j] < 0) continue;
1527 ilocal[nleaves] = j;
1528 iremote[nleaves] = ix[m+j];
1529 nleaves++;
1530 }
1531 ierr = ISGetLocalSize(isglobal[i],&mil);CHKERRQ(ierr);
1532 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);CHKERRQ(ierr);
1533 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);CHKERRQ(ierr);
1534 ierr = PetscLayoutSetLocalSize(map,mil);CHKERRQ(ierr);
1535 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr);
1536 ierr = PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);CHKERRQ(ierr);
1537 ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr);
1538 ierr = PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr);
1539 ierr = PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr);
1540 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);
1541 ierr = PetscFree2(ilocal,iremote);CHKERRQ(ierr);
1542 } else {
1543 ierr = ISGetLocalSize(isglobal[i],&mi);CHKERRQ(ierr);
1544 for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1545 }
1546 ierr = ISRestoreIndices(isglobal[i],&ix2);CHKERRQ(ierr);
1547 m += mi;
1548 }
1549 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr);
1550 PetscFunctionReturn(0);
1551 }
1552
1553
1554 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1555 /*
1556 nprocessors = NP
1557 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1558 proc 0: => (g_0,h_0,)
1559 proc 1: => (g_1,h_1,)
1560 ...
1561 proc nprocs-1: => (g_NP-1,h_NP-1,)
1562
1563 proc 0: proc 1: proc nprocs-1:
1564 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1))
1565
1566 proc 0:
1567 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1568 proc 1:
1569 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)
1570
1571 proc NP-1:
1572 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1573 */
MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])1574 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1575 {
1576 Mat_Nest *vs = (Mat_Nest*)A->data;
1577 PetscInt i,j,offset,n,nsum,bs;
1578 PetscErrorCode ierr;
1579 Mat sub = NULL;
1580
1581 PetscFunctionBegin;
1582 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr);
1583 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr);
1584 if (is_row) { /* valid IS is passed in */
1585 /* refs on is[] are incremeneted */
1586 for (i=0; i<vs->nr; i++) {
1587 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr);
1588
1589 vs->isglobal.row[i] = is_row[i];
1590 }
1591 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1592 nsum = 0;
1593 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1594 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr);
1595 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1596 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr);
1597 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1598 nsum += n;
1599 }
1600 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr);
1601 offset -= nsum;
1602 for (i=0; i<vs->nr; i++) {
1603 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr);
1604 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr);
1605 ierr = MatGetBlockSizes(sub,&bs,NULL);CHKERRQ(ierr);
1606 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr);
1607 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr);
1608 offset += n;
1609 }
1610 }
1611
1612 if (is_col) { /* valid IS is passed in */
1613 /* refs on is[] are incremeneted */
1614 for (j=0; j<vs->nc; j++) {
1615 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr);
1616
1617 vs->isglobal.col[j] = is_col[j];
1618 }
1619 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1620 offset = A->cmap->rstart;
1621 nsum = 0;
1622 for (j=0; j<vs->nc; j++) {
1623 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr);
1624 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1625 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr);
1626 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1627 nsum += n;
1628 }
1629 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr);
1630 offset -= nsum;
1631 for (j=0; j<vs->nc; j++) {
1632 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr);
1633 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr);
1634 ierr = MatGetBlockSizes(sub,NULL,&bs);CHKERRQ(ierr);
1635 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr);
1636 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr);
1637 offset += n;
1638 }
1639 }
1640
1641 /* Set up the local ISs */
1642 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr);
1643 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr);
1644 for (i=0,offset=0; i<vs->nr; i++) {
1645 IS isloc;
1646 ISLocalToGlobalMapping rmap = NULL;
1647 PetscInt nlocal,bs;
1648 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr);
1649 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);}
1650 if (rmap) {
1651 ierr = MatGetBlockSizes(sub,&bs,NULL);CHKERRQ(ierr);
1652 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr);
1653 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr);
1654 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr);
1655 } else {
1656 nlocal = 0;
1657 isloc = NULL;
1658 }
1659 vs->islocal.row[i] = isloc;
1660 offset += nlocal;
1661 }
1662 for (i=0,offset=0; i<vs->nc; i++) {
1663 IS isloc;
1664 ISLocalToGlobalMapping cmap = NULL;
1665 PetscInt nlocal,bs;
1666 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr);
1667 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);}
1668 if (cmap) {
1669 ierr = MatGetBlockSizes(sub,NULL,&bs);CHKERRQ(ierr);
1670 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr);
1671 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr);
1672 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr);
1673 } else {
1674 nlocal = 0;
1675 isloc = NULL;
1676 }
1677 vs->islocal.col[i] = isloc;
1678 offset += nlocal;
1679 }
1680
1681 /* Set up the aggregate ISLocalToGlobalMapping */
1682 {
1683 ISLocalToGlobalMapping rmap,cmap;
1684 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr);
1685 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr);
1686 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);}
1687 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr);
1688 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr);
1689 }
1690
1691 if (PetscDefined(USE_DEBUG)) {
1692 for (i=0; i<vs->nr; i++) {
1693 for (j=0; j<vs->nc; j++) {
1694 PetscInt m,n,M,N,mi,ni,Mi,Ni;
1695 Mat B = vs->m[i][j];
1696 if (!B) continue;
1697 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr);
1698 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr);
1699 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr);
1700 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr);
1701 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr);
1702 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr);
1703 if (M != Mi || N != Ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Global sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",M,N,i,j,Mi,Ni);
1704 if (m != mi || n != ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Local sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",m,n,i,j,mi,ni);
1705 }
1706 }
1707 }
1708
1709 /* Set A->assembled if all non-null blocks are currently assembled */
1710 for (i=0; i<vs->nr; i++) {
1711 for (j=0; j<vs->nc; j++) {
1712 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0);
1713 }
1714 }
1715 A->assembled = PETSC_TRUE;
1716 PetscFunctionReturn(0);
1717 }
1718
1719 /*@C
1720 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately
1721
1722 Collective on Mat
1723
1724 Input Parameter:
1725 + comm - Communicator for the new Mat
1726 . nr - number of nested row blocks
1727 . is_row - index sets for each nested row block, or NULL to make contiguous
1728 . nc - number of nested column blocks
1729 . is_col - index sets for each nested column block, or NULL to make contiguous
1730 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL
1731
1732 Output Parameter:
1733 . B - new matrix
1734
1735 Level: advanced
1736
1737 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1738 MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1739 MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1740 @*/
MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat * B)1741 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1742 {
1743 Mat A;
1744 PetscErrorCode ierr;
1745
1746 PetscFunctionBegin;
1747 *B = NULL;
1748 ierr = MatCreate(comm,&A);CHKERRQ(ierr);
1749 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr);
1750 A->preallocated = PETSC_TRUE;
1751 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr);
1752 *B = A;
1753 PetscFunctionReturn(0);
1754 }
1755
MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)1756 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1757 {
1758 Mat_Nest *nest = (Mat_Nest*)A->data;
1759 Mat *trans;
1760 PetscScalar **avv;
1761 PetscScalar *vv;
1762 PetscInt **aii,**ajj;
1763 PetscInt *ii,*jj,*ci;
1764 PetscInt nr,nc,nnz,i,j;
1765 PetscBool done;
1766 PetscErrorCode ierr;
1767
1768 PetscFunctionBegin;
1769 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr);
1770 if (reuse == MAT_REUSE_MATRIX) {
1771 PetscInt rnr;
1772
1773 ierr = MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);CHKERRQ(ierr);
1774 if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1775 if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1776 ierr = MatSeqAIJGetArray(*newmat,&vv);CHKERRQ(ierr);
1777 }
1778 /* extract CSR for nested SeqAIJ matrices */
1779 nnz = 0;
1780 ierr = PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);CHKERRQ(ierr);
1781 for (i=0; i<nest->nr; ++i) {
1782 for (j=0; j<nest->nc; ++j) {
1783 Mat B = nest->m[i][j];
1784 if (B) {
1785 PetscScalar *naa;
1786 PetscInt *nii,*njj,nnr;
1787 PetscBool istrans;
1788
1789 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr);
1790 if (istrans) {
1791 Mat Bt;
1792
1793 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr);
1794 ierr = MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);CHKERRQ(ierr);
1795 B = trans[i*nest->nc+j];
1796 }
1797 ierr = MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);CHKERRQ(ierr);
1798 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1799 ierr = MatSeqAIJGetArray(B,&naa);CHKERRQ(ierr);
1800 nnz += nii[nnr];
1801
1802 aii[i*nest->nc+j] = nii;
1803 ajj[i*nest->nc+j] = njj;
1804 avv[i*nest->nc+j] = naa;
1805 }
1806 }
1807 }
1808 if (reuse != MAT_REUSE_MATRIX) {
1809 ierr = PetscMalloc1(nr+1,&ii);CHKERRQ(ierr);
1810 ierr = PetscMalloc1(nnz,&jj);CHKERRQ(ierr);
1811 ierr = PetscMalloc1(nnz,&vv);CHKERRQ(ierr);
1812 } else {
1813 if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1814 }
1815
1816 /* new row pointer */
1817 ierr = PetscArrayzero(ii,nr+1);CHKERRQ(ierr);
1818 for (i=0; i<nest->nr; ++i) {
1819 PetscInt ncr,rst;
1820
1821 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr);
1822 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr);
1823 for (j=0; j<nest->nc; ++j) {
1824 if (aii[i*nest->nc+j]) {
1825 PetscInt *nii = aii[i*nest->nc+j];
1826 PetscInt ir;
1827
1828 for (ir=rst; ir<ncr+rst; ++ir) {
1829 ii[ir+1] += nii[1]-nii[0];
1830 nii++;
1831 }
1832 }
1833 }
1834 }
1835 for (i=0; i<nr; i++) ii[i+1] += ii[i];
1836
1837 /* construct CSR for the new matrix */
1838 ierr = PetscCalloc1(nr,&ci);CHKERRQ(ierr);
1839 for (i=0; i<nest->nr; ++i) {
1840 PetscInt ncr,rst;
1841
1842 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr);
1843 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr);
1844 for (j=0; j<nest->nc; ++j) {
1845 if (aii[i*nest->nc+j]) {
1846 PetscScalar *nvv = avv[i*nest->nc+j];
1847 PetscInt *nii = aii[i*nest->nc+j];
1848 PetscInt *njj = ajj[i*nest->nc+j];
1849 PetscInt ir,cst;
1850
1851 ierr = ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);CHKERRQ(ierr);
1852 for (ir=rst; ir<ncr+rst; ++ir) {
1853 PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];
1854
1855 for (ij=0;ij<rsize;ij++) {
1856 jj[ist+ij] = *njj+cst;
1857 vv[ist+ij] = *nvv;
1858 njj++;
1859 nvv++;
1860 }
1861 ci[ir] += rsize;
1862 nii++;
1863 }
1864 }
1865 }
1866 }
1867 ierr = PetscFree(ci);CHKERRQ(ierr);
1868
1869 /* restore info */
1870 for (i=0; i<nest->nr; ++i) {
1871 for (j=0; j<nest->nc; ++j) {
1872 Mat B = nest->m[i][j];
1873 if (B) {
1874 PetscInt nnr = 0, k = i*nest->nc+j;
1875
1876 B = (trans[k] ? trans[k] : B);
1877 ierr = MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);CHKERRQ(ierr);
1878 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1879 ierr = MatSeqAIJRestoreArray(B,&avv[k]);CHKERRQ(ierr);
1880 ierr = MatDestroy(&trans[k]);CHKERRQ(ierr);
1881 }
1882 }
1883 }
1884 ierr = PetscFree4(aii,ajj,avv,trans);CHKERRQ(ierr);
1885
1886 /* finalize newmat */
1887 if (reuse == MAT_INITIAL_MATRIX) {
1888 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);CHKERRQ(ierr);
1889 } else if (reuse == MAT_INPLACE_MATRIX) {
1890 Mat B;
1891
1892 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);CHKERRQ(ierr);
1893 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr);
1894 }
1895 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1896 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1897 {
1898 Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1899 a->free_a = PETSC_TRUE;
1900 a->free_ij = PETSC_TRUE;
1901 }
1902 PetscFunctionReturn(0);
1903 }
1904
MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)1905 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1906 {
1907 PetscErrorCode ierr;
1908 Mat_Nest *nest = (Mat_Nest*)A->data;
1909 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1910 PetscInt cstart,cend;
1911 PetscMPIInt size;
1912 Mat C;
1913
1914 PetscFunctionBegin;
1915 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
1916 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1917 PetscInt nf;
1918 PetscBool fast;
1919
1920 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr);
1921 if (!fast) {
1922 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr);
1923 }
1924 for (i=0; i<nest->nr && fast; ++i) {
1925 for (j=0; j<nest->nc && fast; ++j) {
1926 Mat B = nest->m[i][j];
1927 if (B) {
1928 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr);
1929 if (!fast) {
1930 PetscBool istrans;
1931
1932 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr);
1933 if (istrans) {
1934 Mat Bt;
1935
1936 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr);
1937 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);CHKERRQ(ierr);
1938 }
1939 }
1940 }
1941 }
1942 }
1943 for (i=0, nf=0; i<nest->nr && fast; ++i) {
1944 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr);
1945 if (fast) {
1946 PetscInt f,s;
1947
1948 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr);
1949 if (f != nf || s != 1) { fast = PETSC_FALSE; }
1950 else {
1951 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr);
1952 nf += f;
1953 }
1954 }
1955 }
1956 for (i=0, nf=0; i<nest->nc && fast; ++i) {
1957 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr);
1958 if (fast) {
1959 PetscInt f,s;
1960
1961 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr);
1962 if (f != nf || s != 1) { fast = PETSC_FALSE; }
1963 else {
1964 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr);
1965 nf += f;
1966 }
1967 }
1968 }
1969 if (fast) {
1970 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr);
1971 PetscFunctionReturn(0);
1972 }
1973 }
1974 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
1975 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
1976 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr);
1977 switch (reuse) {
1978 case MAT_INITIAL_MATRIX:
1979 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr);
1980 ierr = MatSetType(C,newtype);CHKERRQ(ierr);
1981 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr);
1982 *newmat = C;
1983 break;
1984 case MAT_REUSE_MATRIX:
1985 C = *newmat;
1986 break;
1987 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1988 }
1989 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr);
1990 onnz = dnnz + m;
1991 for (k=0; k<m; k++) {
1992 dnnz[k] = 0;
1993 onnz[k] = 0;
1994 }
1995 for (j=0; j<nest->nc; ++j) {
1996 IS bNis;
1997 PetscInt bN;
1998 const PetscInt *bNindices;
1999 /* Using global column indices and ISAllGather() is not scalable. */
2000 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr);
2001 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr);
2002 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr);
2003 for (i=0; i<nest->nr; ++i) {
2004 PetscSF bmsf;
2005 PetscSFNode *iremote;
2006 Mat B;
2007 PetscInt bm, *sub_dnnz,*sub_onnz, br;
2008 const PetscInt *bmindices;
2009 B = nest->m[i][j];
2010 if (!B) continue;
2011 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr);
2012 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr);
2013 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr);
2014 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr);
2015 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr);
2016 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr);
2017 for (k = 0; k < bm; ++k){
2018 sub_dnnz[k] = 0;
2019 sub_onnz[k] = 0;
2020 }
2021 /*
2022 Locate the owners for all of the locally-owned global row indices for this row block.
2023 These determine the roots of PetscSF used to communicate preallocation data to row owners.
2024 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2025 */
2026 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr);
2027 for (br = 0; br < bm; ++br) {
2028 PetscInt row = bmindices[br], brncols, col;
2029 const PetscInt *brcols;
2030 PetscInt rowrel = 0; /* row's relative index on its owner rank */
2031 PetscMPIInt rowowner = 0;
2032 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr);
2033 /* how many roots */
2034 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2035 /* get nonzero pattern */
2036 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr);
2037 for (k=0; k<brncols; k++) {
2038 col = bNindices[brcols[k]];
2039 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2040 sub_dnnz[br]++;
2041 } else {
2042 sub_onnz[br]++;
2043 }
2044 }
2045 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr);
2046 }
2047 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr);
2048 /* bsf will have to take care of disposing of bedges. */
2049 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr);
2050 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr);
2051 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr);
2052 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr);
2053 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr);
2054 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr);
2055 ierr = PetscFree(sub_onnz);CHKERRQ(ierr);
2056 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr);
2057 }
2058 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr);
2059 ierr = ISDestroy(&bNis);CHKERRQ(ierr);
2060 }
2061 /* Resize preallocation if overestimated */
2062 for (i=0;i<m;i++) {
2063 dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2064 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2065 }
2066 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr);
2067 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr);
2068 ierr = PetscFree(dnnz);CHKERRQ(ierr);
2069
2070 /* Fill by row */
2071 for (j=0; j<nest->nc; ++j) {
2072 /* Using global column indices and ISAllGather() is not scalable. */
2073 IS bNis;
2074 PetscInt bN;
2075 const PetscInt *bNindices;
2076 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr);
2077 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr);
2078 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr);
2079 for (i=0; i<nest->nr; ++i) {
2080 Mat B;
2081 PetscInt bm, br;
2082 const PetscInt *bmindices;
2083 B = nest->m[i][j];
2084 if (!B) continue;
2085 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr);
2086 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr);
2087 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr);
2088 for (br = 0; br < bm; ++br) {
2089 PetscInt row = bmindices[br], brncols, *cols;
2090 const PetscInt *brcols;
2091 const PetscScalar *brcoldata;
2092 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr);
2093 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr);
2094 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2095 /*
2096 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2097 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2098 */
2099 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr);
2100 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr);
2101 ierr = PetscFree(cols);CHKERRQ(ierr);
2102 }
2103 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr);
2104 }
2105 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr);
2106 ierr = ISDestroy(&bNis);CHKERRQ(ierr);
2107 }
2108 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2109 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2110 PetscFunctionReturn(0);
2111 }
2112
MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool * has)2113 PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2114 {
2115 Mat_Nest *bA = (Mat_Nest*)mat->data;
2116 MatOperation opAdd;
2117 PetscInt i,j,nr = bA->nr,nc = bA->nc;
2118 PetscBool flg;
2119 PetscErrorCode ierr;
2120 PetscFunctionBegin;
2121
2122 *has = PETSC_FALSE;
2123 if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2124 opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2125 for (j=0; j<nc; j++) {
2126 for (i=0; i<nr; i++) {
2127 if (!bA->m[i][j]) continue;
2128 ierr = MatHasOperation(bA->m[i][j],opAdd,&flg);CHKERRQ(ierr);
2129 if (!flg) PetscFunctionReturn(0);
2130 }
2131 }
2132 }
2133 if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2134 PetscFunctionReturn(0);
2135 }
2136
2137 /*MC
2138 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.
2139
2140 Level: intermediate
2141
2142 Notes:
2143 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
2144 It allows the use of symmetric and block formats for parts of multi-physics simulations.
2145 It is usually used with DMComposite and DMCreateMatrix()
2146
2147 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
2148 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
2149 than the nest matrix.
2150
2151 .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2152 VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2153 MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2154 M*/
MatCreate_Nest(Mat A)2155 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2156 {
2157 Mat_Nest *s;
2158 PetscErrorCode ierr;
2159
2160 PetscFunctionBegin;
2161 ierr = PetscNewLog(A,&s);CHKERRQ(ierr);
2162 A->data = (void*)s;
2163
2164 s->nr = -1;
2165 s->nc = -1;
2166 s->m = NULL;
2167 s->splitassembly = PETSC_FALSE;
2168
2169 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr);
2170
2171 A->ops->mult = MatMult_Nest;
2172 A->ops->multadd = MatMultAdd_Nest;
2173 A->ops->multtranspose = MatMultTranspose_Nest;
2174 A->ops->multtransposeadd = MatMultTransposeAdd_Nest;
2175 A->ops->transpose = MatTranspose_Nest;
2176 A->ops->assemblybegin = MatAssemblyBegin_Nest;
2177 A->ops->assemblyend = MatAssemblyEnd_Nest;
2178 A->ops->zeroentries = MatZeroEntries_Nest;
2179 A->ops->copy = MatCopy_Nest;
2180 A->ops->axpy = MatAXPY_Nest;
2181 A->ops->duplicate = MatDuplicate_Nest;
2182 A->ops->createsubmatrix = MatCreateSubMatrix_Nest;
2183 A->ops->destroy = MatDestroy_Nest;
2184 A->ops->view = MatView_Nest;
2185 A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2186 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest;
2187 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2188 A->ops->getdiagonal = MatGetDiagonal_Nest;
2189 A->ops->diagonalscale = MatDiagonalScale_Nest;
2190 A->ops->scale = MatScale_Nest;
2191 A->ops->shift = MatShift_Nest;
2192 A->ops->diagonalset = MatDiagonalSet_Nest;
2193 A->ops->setrandom = MatSetRandom_Nest;
2194 A->ops->hasoperation = MatHasOperation_Nest;
2195 A->ops->missingdiagonal = MatMissingDiagonal_Nest;
2196
2197 A->spptr = NULL;
2198 A->assembled = PETSC_FALSE;
2199
2200 /* expose Nest api's */
2201 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr);
2202 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr);
2203 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr);
2204 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr);
2205 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr);
2206 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr);
2207 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr);
2208 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr);
2209 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr);
2210 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr);
2211 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr);
2212 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr);
2213 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);CHKERRQ(ierr);
2214 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);CHKERRQ(ierr);
2215 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);CHKERRQ(ierr);
2216
2217 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr);
2218 PetscFunctionReturn(0);
2219 }
2220