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 */ 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 */ 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 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 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 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 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 /* --------------------------------------------------------- */ 257 static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C) 258 { 259 PetscFunctionBegin; 260 C->ops->productsymbolic = MatProductSymbolic_Nest_Dense; 261 PetscFunctionReturn(0); 262 } 263 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 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 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 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 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 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 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 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 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 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 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 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 559 static PetscErrorCode MatNestFindISRange(Mat A,PetscInt n,const IS list[],IS is,PetscInt *begin,PetscInt *end) 560 { 561 PetscInt i,j,size,m; 562 PetscBool flg; 563 IS out,concatenate[2]; 564 PetscErrorCode ierr; 565 566 PetscFunctionBegin; 567 PetscValidPointer(list,3); 568 PetscValidHeaderSpecific(is,IS_CLASSID,4); 569 if (begin) { 570 PetscValidIntPointer(begin,5); 571 *begin = -1; 572 } 573 if (end) { 574 PetscValidIntPointer(end,6); 575 *end = -1; 576 } 577 for (i=0; i<n; i++) { 578 if (!list[i]) continue; 579 ierr = ISEqualUnsorted(list[i],is,&flg);CHKERRQ(ierr); 580 if (flg) { 581 if (begin) *begin = i; 582 if (end) *end = i+1; 583 PetscFunctionReturn(0); 584 } 585 } 586 ierr = ISGetSize(is,&size);CHKERRQ(ierr); 587 for (i=0; i<n-1; i++) { 588 if (!list[i]) continue; 589 m = 0; 590 ierr = ISConcatenate(PetscObjectComm((PetscObject)A),2,list+i,&out);CHKERRQ(ierr); 591 ierr = ISGetSize(out,&m);CHKERRQ(ierr); 592 for (j=i+2; j<n && m<size; j++) { 593 if (list[j]) { 594 concatenate[0] = out; 595 concatenate[1] = list[j]; 596 ierr = ISConcatenate(PetscObjectComm((PetscObject)A),2,concatenate,&out);CHKERRQ(ierr); 597 ierr = ISDestroy(concatenate);CHKERRQ(ierr); 598 ierr = ISGetSize(out,&m);CHKERRQ(ierr); 599 } 600 } 601 if (m == size) { 602 ierr = ISEqualUnsorted(out,is,&flg);CHKERRQ(ierr); 603 if (flg) { 604 if (begin) *begin = i; 605 if (end) *end = j; 606 ierr = ISDestroy(&out);CHKERRQ(ierr); 607 PetscFunctionReturn(0); 608 } 609 } 610 ierr = ISDestroy(&out);CHKERRQ(ierr); 611 } 612 PetscFunctionReturn(0); 613 } 614 615 616 static PetscErrorCode MatNestFillEmptyMat_Private(Mat A,PetscInt i,PetscInt j,Mat *B) 617 { 618 Mat_Nest *vs = (Mat_Nest*)A->data; 619 PetscInt lr,lc; 620 PetscErrorCode ierr; 621 622 PetscFunctionBegin; 623 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 624 ierr = ISGetLocalSize(vs->isglobal.row[i],&lr);CHKERRQ(ierr); 625 ierr = ISGetLocalSize(vs->isglobal.col[j],&lc);CHKERRQ(ierr); 626 ierr = MatSetSizes(*B,lr,lc,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 627 ierr = MatSetType(*B,MATAIJ);CHKERRQ(ierr); 628 ierr = MatSeqAIJSetPreallocation(*B,0,NULL);CHKERRQ(ierr); 629 ierr = MatMPIAIJSetPreallocation(*B,0,NULL,0,NULL);CHKERRQ(ierr); 630 ierr = MatSetUp(*B);CHKERRQ(ierr); 631 ierr = MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 632 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 633 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 634 PetscFunctionReturn(0); 635 } 636 637 static PetscErrorCode MatNestGetBlock_Private(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *B) 638 { 639 Mat_Nest *vs = (Mat_Nest*)A->data; 640 Mat *a; 641 PetscInt i,j,k,l,nr=rend-rbegin,nc=cend-cbegin; 642 char keyname[256]; 643 PetscBool *b; 644 PetscBool flg; 645 PetscErrorCode ierr; 646 647 PetscFunctionBegin; 648 *B = NULL; 649 ierr = PetscSNPrintf(keyname,sizeof(keyname),"NestBlock_%D-%Dx%D-%D",rbegin,rend,cbegin,cend);CHKERRQ(ierr); 650 ierr = PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);CHKERRQ(ierr); 651 if (*B) PetscFunctionReturn(0); 652 653 ierr = PetscMalloc2(nr*nc,&a,nr*nc,&b);CHKERRQ(ierr); 654 for (i=0; i<nr; i++) { 655 for (j=0; j<nc; j++) { 656 a[i*nc + j] = vs->m[rbegin+i][cbegin+j]; 657 b[i*nc + j] = PETSC_FALSE; 658 } 659 } 660 if (nc!=vs->nc&&nr!=vs->nr) { 661 for (i=0; i<nr; i++) { 662 for (j=0; j<nc; j++) { 663 flg = PETSC_FALSE; 664 for (k=0; (k<nr&&!flg); k++) { 665 if (a[j + k*nc]) flg = PETSC_TRUE; 666 } 667 if (flg) { 668 flg = PETSC_FALSE; 669 for (l=0; (l<nc&&!flg); l++) { 670 if (a[i*nc + l]) flg = PETSC_TRUE; 671 } 672 } 673 if (!flg) { 674 b[i*nc + j] = PETSC_TRUE; 675 ierr = MatNestFillEmptyMat_Private(A,rbegin+i,cbegin+j,a + i*nc + j);CHKERRQ(ierr); 676 } 677 } 678 } 679 } 680 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nr,nr!=vs->nr?NULL:vs->isglobal.row,nc,nc!=vs->nc?NULL:vs->isglobal.col,a,B);CHKERRQ(ierr); 681 for (i=0; i<nr; i++) { 682 for (j=0; j<nc; j++) { 683 if (b[i*nc + j]) { 684 ierr = MatDestroy(a + i*nc + j);CHKERRQ(ierr); 685 } 686 } 687 } 688 ierr = PetscFree2(a,b);CHKERRQ(ierr); 689 (*B)->assembled = A->assembled; 690 ierr = PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);CHKERRQ(ierr); 691 ierr = PetscObjectDereference((PetscObject)*B);CHKERRQ(ierr); /* Leave the only remaining reference in the composition */ 692 PetscFunctionReturn(0); 693 } 694 695 static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B) 696 { 697 Mat_Nest *vs = (Mat_Nest*)A->data; 698 PetscInt rbegin,rend,cbegin,cend; 699 PetscErrorCode ierr; 700 701 PetscFunctionBegin; 702 ierr = MatNestFindISRange(A,vs->nr,is->row,isrow,&rbegin,&rend);CHKERRQ(ierr); 703 ierr = MatNestFindISRange(A,vs->nc,is->col,iscol,&cbegin,&cend);CHKERRQ(ierr); 704 if (rend == rbegin + 1 && cend == cbegin + 1) { 705 if (!vs->m[rbegin][cbegin]) { 706 ierr = MatNestFillEmptyMat_Private(A,rbegin,cbegin,vs->m[rbegin] + cbegin);CHKERRQ(ierr); 707 } 708 *B = vs->m[rbegin][cbegin]; 709 } else if (rbegin != -1 && cbegin != -1) { 710 ierr = MatNestGetBlock_Private(A,rbegin,rend,cbegin,cend,B);CHKERRQ(ierr); 711 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set"); 712 PetscFunctionReturn(0); 713 } 714 715 /* 716 TODO: This does not actually returns a submatrix we can modify 717 */ 718 static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B) 719 { 720 PetscErrorCode ierr; 721 Mat_Nest *vs = (Mat_Nest*)A->data; 722 Mat sub; 723 724 PetscFunctionBegin; 725 ierr = MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);CHKERRQ(ierr); 726 switch (reuse) { 727 case MAT_INITIAL_MATRIX: 728 if (sub) { ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr); } 729 *B = sub; 730 break; 731 case MAT_REUSE_MATRIX: 732 if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call"); 733 break; 734 case MAT_IGNORE_MATRIX: /* Nothing to do */ 735 break; 736 case MAT_INPLACE_MATRIX: /* Nothing to do */ 737 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet"); 738 } 739 PetscFunctionReturn(0); 740 } 741 742 PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 743 { 744 PetscErrorCode ierr; 745 Mat_Nest *vs = (Mat_Nest*)A->data; 746 Mat sub; 747 748 PetscFunctionBegin; 749 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 750 /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */ 751 if (sub) {ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr);} 752 *B = sub; 753 PetscFunctionReturn(0); 754 } 755 756 static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 757 { 758 PetscErrorCode ierr; 759 Mat_Nest *vs = (Mat_Nest*)A->data; 760 Mat sub; 761 762 PetscFunctionBegin; 763 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 764 if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten"); 765 if (sub) { 766 if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times"); 767 ierr = MatDestroy(B);CHKERRQ(ierr); 768 } 769 PetscFunctionReturn(0); 770 } 771 772 static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v) 773 { 774 Mat_Nest *bA = (Mat_Nest*)A->data; 775 PetscInt i; 776 PetscErrorCode ierr; 777 778 PetscFunctionBegin; 779 for (i=0; i<bA->nr; i++) { 780 Vec bv; 781 ierr = VecGetSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 782 if (bA->m[i][i]) { 783 ierr = MatGetDiagonal(bA->m[i][i],bv);CHKERRQ(ierr); 784 } else { 785 ierr = VecSet(bv,0.0);CHKERRQ(ierr); 786 } 787 ierr = VecRestoreSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 788 } 789 PetscFunctionReturn(0); 790 } 791 792 static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r) 793 { 794 Mat_Nest *bA = (Mat_Nest*)A->data; 795 Vec bl,*br; 796 PetscInt i,j; 797 PetscErrorCode ierr; 798 799 PetscFunctionBegin; 800 ierr = PetscCalloc1(bA->nc,&br);CHKERRQ(ierr); 801 if (r) { 802 for (j=0; j<bA->nc; j++) {ierr = VecGetSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);} 803 } 804 bl = NULL; 805 for (i=0; i<bA->nr; i++) { 806 if (l) { 807 ierr = VecGetSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr); 808 } 809 for (j=0; j<bA->nc; j++) { 810 if (bA->m[i][j]) { 811 ierr = MatDiagonalScale(bA->m[i][j],bl,br[j]);CHKERRQ(ierr); 812 } 813 } 814 if (l) { 815 ierr = VecRestoreSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr); 816 } 817 } 818 if (r) { 819 for (j=0; j<bA->nc; j++) {ierr = VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);} 820 } 821 ierr = PetscFree(br);CHKERRQ(ierr); 822 PetscFunctionReturn(0); 823 } 824 825 static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a) 826 { 827 Mat_Nest *bA = (Mat_Nest*)A->data; 828 PetscInt i,j; 829 PetscErrorCode ierr; 830 831 PetscFunctionBegin; 832 for (i=0; i<bA->nr; i++) { 833 for (j=0; j<bA->nc; j++) { 834 if (bA->m[i][j]) { 835 ierr = MatScale(bA->m[i][j],a);CHKERRQ(ierr); 836 } 837 } 838 } 839 PetscFunctionReturn(0); 840 } 841 842 static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a) 843 { 844 Mat_Nest *bA = (Mat_Nest*)A->data; 845 PetscInt i; 846 PetscErrorCode ierr; 847 PetscBool nnzstate = PETSC_FALSE; 848 849 PetscFunctionBegin; 850 for (i=0; i<bA->nr; i++) { 851 PetscObjectState subnnzstate = 0; 852 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); 853 ierr = MatShift(bA->m[i][i],a);CHKERRQ(ierr); 854 ierr = MatGetNonzeroState(bA->m[i][i],&subnnzstate);CHKERRQ(ierr); 855 nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate); 856 bA->nnzstate[i*bA->nc+i] = subnnzstate; 857 } 858 if (nnzstate) A->nonzerostate++; 859 PetscFunctionReturn(0); 860 } 861 862 static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is) 863 { 864 Mat_Nest *bA = (Mat_Nest*)A->data; 865 PetscInt i; 866 PetscErrorCode ierr; 867 PetscBool nnzstate = PETSC_FALSE; 868 869 PetscFunctionBegin; 870 for (i=0; i<bA->nr; i++) { 871 PetscObjectState subnnzstate = 0; 872 Vec bv; 873 ierr = VecGetSubVector(D,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 874 if (bA->m[i][i]) { 875 ierr = MatDiagonalSet(bA->m[i][i],bv,is);CHKERRQ(ierr); 876 ierr = MatGetNonzeroState(bA->m[i][i],&subnnzstate);CHKERRQ(ierr); 877 } 878 ierr = VecRestoreSubVector(D,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 879 nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate); 880 bA->nnzstate[i*bA->nc+i] = subnnzstate; 881 } 882 if (nnzstate) A->nonzerostate++; 883 PetscFunctionReturn(0); 884 } 885 886 static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx) 887 { 888 Mat_Nest *bA = (Mat_Nest*)A->data; 889 PetscInt i,j; 890 PetscErrorCode ierr; 891 892 PetscFunctionBegin; 893 for (i=0; i<bA->nr; i++) { 894 for (j=0; j<bA->nc; j++) { 895 if (bA->m[i][j]) { 896 ierr = MatSetRandom(bA->m[i][j],rctx);CHKERRQ(ierr); 897 } 898 } 899 } 900 PetscFunctionReturn(0); 901 } 902 903 static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left) 904 { 905 Mat_Nest *bA = (Mat_Nest*)A->data; 906 Vec *L,*R; 907 MPI_Comm comm; 908 PetscInt i,j; 909 PetscErrorCode ierr; 910 911 PetscFunctionBegin; 912 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 913 if (right) { 914 /* allocate R */ 915 ierr = PetscMalloc1(bA->nc, &R);CHKERRQ(ierr); 916 /* Create the right vectors */ 917 for (j=0; j<bA->nc; j++) { 918 for (i=0; i<bA->nr; i++) { 919 if (bA->m[i][j]) { 920 ierr = MatCreateVecs(bA->m[i][j],&R[j],NULL);CHKERRQ(ierr); 921 break; 922 } 923 } 924 if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 925 } 926 ierr = VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);CHKERRQ(ierr); 927 /* hand back control to the nest vector */ 928 for (j=0; j<bA->nc; j++) { 929 ierr = VecDestroy(&R[j]);CHKERRQ(ierr); 930 } 931 ierr = PetscFree(R);CHKERRQ(ierr); 932 } 933 934 if (left) { 935 /* allocate L */ 936 ierr = PetscMalloc1(bA->nr, &L);CHKERRQ(ierr); 937 /* Create the left vectors */ 938 for (i=0; i<bA->nr; i++) { 939 for (j=0; j<bA->nc; j++) { 940 if (bA->m[i][j]) { 941 ierr = MatCreateVecs(bA->m[i][j],NULL,&L[i]);CHKERRQ(ierr); 942 break; 943 } 944 } 945 if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 946 } 947 948 ierr = VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);CHKERRQ(ierr); 949 for (i=0; i<bA->nr; i++) { 950 ierr = VecDestroy(&L[i]);CHKERRQ(ierr); 951 } 952 953 ierr = PetscFree(L);CHKERRQ(ierr); 954 } 955 PetscFunctionReturn(0); 956 } 957 958 static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer) 959 { 960 Mat_Nest *bA = (Mat_Nest*)A->data; 961 PetscBool isascii,viewSub = PETSC_FALSE; 962 PetscInt i,j; 963 PetscErrorCode ierr; 964 965 PetscFunctionBegin; 966 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 967 if (isascii) { 968 969 ierr = PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);CHKERRQ(ierr); 970 ierr = PetscViewerASCIIPrintf(viewer,"Matrix object: \n");CHKERRQ(ierr); 971 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 972 ierr = PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);CHKERRQ(ierr); 973 974 ierr = PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");CHKERRQ(ierr); 975 for (i=0; i<bA->nr; i++) { 976 for (j=0; j<bA->nc; j++) { 977 MatType type; 978 char name[256] = "",prefix[256] = ""; 979 PetscInt NR,NC; 980 PetscBool isNest = PETSC_FALSE; 981 982 if (!bA->m[i][j]) { 983 ierr = PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);CHKERRQ(ierr); 984 continue; 985 } 986 ierr = MatGetSize(bA->m[i][j],&NR,&NC);CHKERRQ(ierr); 987 ierr = MatGetType(bA->m[i][j], &type);CHKERRQ(ierr); 988 if (((PetscObject)bA->m[i][j])->name) {ierr = PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);CHKERRQ(ierr);} 989 if (((PetscObject)bA->m[i][j])->prefix) {ierr = PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);CHKERRQ(ierr);} 990 ierr = PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);CHKERRQ(ierr); 991 992 ierr = PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);CHKERRQ(ierr); 993 994 if (isNest || viewSub) { 995 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); /* push1 */ 996 ierr = MatView(bA->m[i][j],viewer);CHKERRQ(ierr); 997 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop1 */ 998 } 999 } 1000 } 1001 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop0 */ 1002 } 1003 PetscFunctionReturn(0); 1004 } 1005 1006 static PetscErrorCode MatZeroEntries_Nest(Mat A) 1007 { 1008 Mat_Nest *bA = (Mat_Nest*)A->data; 1009 PetscInt i,j; 1010 PetscErrorCode ierr; 1011 1012 PetscFunctionBegin; 1013 for (i=0; i<bA->nr; i++) { 1014 for (j=0; j<bA->nc; j++) { 1015 if (!bA->m[i][j]) continue; 1016 ierr = MatZeroEntries(bA->m[i][j]);CHKERRQ(ierr); 1017 } 1018 } 1019 PetscFunctionReturn(0); 1020 } 1021 1022 static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str) 1023 { 1024 Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data; 1025 PetscInt i,j,nr = bA->nr,nc = bA->nc; 1026 PetscErrorCode ierr; 1027 PetscBool nnzstate = PETSC_FALSE; 1028 1029 PetscFunctionBegin; 1030 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); 1031 for (i=0; i<nr; i++) { 1032 for (j=0; j<nc; j++) { 1033 PetscObjectState subnnzstate = 0; 1034 if (bA->m[i][j] && bB->m[i][j]) { 1035 ierr = MatCopy(bA->m[i][j],bB->m[i][j],str);CHKERRQ(ierr); 1036 } 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); 1037 ierr = MatGetNonzeroState(bB->m[i][j],&subnnzstate);CHKERRQ(ierr); 1038 nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate); 1039 bB->nnzstate[i*nc+j] = subnnzstate; 1040 } 1041 } 1042 if (nnzstate) B->nonzerostate++; 1043 PetscFunctionReturn(0); 1044 } 1045 1046 static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str) 1047 { 1048 Mat_Nest *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data; 1049 PetscInt i,j,nr = bY->nr,nc = bY->nc; 1050 PetscErrorCode ierr; 1051 PetscBool nnzstate = PETSC_FALSE; 1052 1053 PetscFunctionBegin; 1054 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); 1055 for (i=0; i<nr; i++) { 1056 for (j=0; j<nc; j++) { 1057 PetscObjectState subnnzstate = 0; 1058 if (bY->m[i][j] && bX->m[i][j]) { 1059 ierr = MatAXPY(bY->m[i][j],a,bX->m[i][j],str);CHKERRQ(ierr); 1060 } else if (bX->m[i][j]) { 1061 Mat M; 1062 1063 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); 1064 ierr = MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);CHKERRQ(ierr); 1065 ierr = MatNestSetSubMat(Y,i,j,M);CHKERRQ(ierr); 1066 ierr = MatDestroy(&M);CHKERRQ(ierr); 1067 } 1068 if (bY->m[i][j]) { ierr = MatGetNonzeroState(bY->m[i][j],&subnnzstate);CHKERRQ(ierr); } 1069 nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate); 1070 bY->nnzstate[i*nc+j] = subnnzstate; 1071 } 1072 } 1073 if (nnzstate) Y->nonzerostate++; 1074 PetscFunctionReturn(0); 1075 } 1076 1077 static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B) 1078 { 1079 Mat_Nest *bA = (Mat_Nest*)A->data; 1080 Mat *b; 1081 PetscInt i,j,nr = bA->nr,nc = bA->nc; 1082 PetscErrorCode ierr; 1083 1084 PetscFunctionBegin; 1085 ierr = PetscMalloc1(nr*nc,&b);CHKERRQ(ierr); 1086 for (i=0; i<nr; i++) { 1087 for (j=0; j<nc; j++) { 1088 if (bA->m[i][j]) { 1089 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr); 1090 } else { 1091 b[i*nc+j] = NULL; 1092 } 1093 } 1094 } 1095 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr); 1096 /* Give the new MatNest exclusive ownership */ 1097 for (i=0; i<nr*nc; i++) { 1098 ierr = MatDestroy(&b[i]);CHKERRQ(ierr); 1099 } 1100 ierr = PetscFree(b);CHKERRQ(ierr); 1101 1102 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1103 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 /* nest api */ 1108 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 1109 { 1110 Mat_Nest *bA = (Mat_Nest*)A->data; 1111 1112 PetscFunctionBegin; 1113 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 1114 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 1115 *mat = bA->m[idxm][jdxm]; 1116 PetscFunctionReturn(0); 1117 } 1118 1119 /*@ 1120 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 1121 1122 Not collective 1123 1124 Input Parameters: 1125 + A - nest matrix 1126 . idxm - index of the matrix within the nest matrix 1127 - jdxm - index of the matrix within the nest matrix 1128 1129 Output Parameter: 1130 . sub - matrix at index idxm,jdxm within the nest matrix 1131 1132 Level: developer 1133 1134 .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(), 1135 MatNestGetLocalISs(), MatNestGetISs() 1136 @*/ 1137 PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 1138 { 1139 PetscErrorCode ierr; 1140 1141 PetscFunctionBegin; 1142 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr); 1143 PetscFunctionReturn(0); 1144 } 1145 1146 PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat) 1147 { 1148 Mat_Nest *bA = (Mat_Nest*)A->data; 1149 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1150 PetscErrorCode ierr; 1151 1152 PetscFunctionBegin; 1153 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 1154 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 1155 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1156 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1157 ierr = ISGetLocalSize(bA->isglobal.row[idxm],&mi);CHKERRQ(ierr); 1158 ierr = ISGetSize(bA->isglobal.row[idxm],&Mi);CHKERRQ(ierr); 1159 ierr = ISGetLocalSize(bA->isglobal.col[jdxm],&ni);CHKERRQ(ierr); 1160 ierr = ISGetSize(bA->isglobal.col[jdxm],&Ni);CHKERRQ(ierr); 1161 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); 1162 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); 1163 1164 /* do not increase object state */ 1165 if (mat == bA->m[idxm][jdxm]) PetscFunctionReturn(0); 1166 1167 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 1168 ierr = MatDestroy(&bA->m[idxm][jdxm]);CHKERRQ(ierr); 1169 bA->m[idxm][jdxm] = mat; 1170 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1171 ierr = MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);CHKERRQ(ierr); 1172 A->nonzerostate++; 1173 PetscFunctionReturn(0); 1174 } 1175 1176 /*@ 1177 MatNestSetSubMat - Set a single submatrix in the nest matrix. 1178 1179 Logically collective on the submatrix communicator 1180 1181 Input Parameters: 1182 + A - nest matrix 1183 . idxm - index of the matrix within the nest matrix 1184 . jdxm - index of the matrix within the nest matrix 1185 - sub - matrix at index idxm,jdxm within the nest matrix 1186 1187 Notes: 1188 The new submatrix must have the same size and communicator as that block of the nest. 1189 1190 This increments the reference count of the submatrix. 1191 1192 Level: developer 1193 1194 .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(), 1195 MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize() 1196 @*/ 1197 PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub) 1198 { 1199 PetscErrorCode ierr; 1200 1201 PetscFunctionBegin; 1202 ierr = PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));CHKERRQ(ierr); 1203 PetscFunctionReturn(0); 1204 } 1205 1206 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 1207 { 1208 Mat_Nest *bA = (Mat_Nest*)A->data; 1209 1210 PetscFunctionBegin; 1211 if (M) *M = bA->nr; 1212 if (N) *N = bA->nc; 1213 if (mat) *mat = bA->m; 1214 PetscFunctionReturn(0); 1215 } 1216 1217 /*@C 1218 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 1219 1220 Not collective 1221 1222 Input Parameters: 1223 . A - nest matrix 1224 1225 Output Parameter: 1226 + M - number of rows in the nest matrix 1227 . N - number of cols in the nest matrix 1228 - mat - 2d array of matrices 1229 1230 Notes: 1231 1232 The user should not free the array mat. 1233 1234 In Fortran, this routine has a calling sequence 1235 $ call MatNestGetSubMats(A, M, N, mat, ierr) 1236 where the space allocated for the optional argument mat is assumed large enough (if provided). 1237 1238 Level: developer 1239 1240 .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(), 1241 MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat() 1242 @*/ 1243 PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 1244 { 1245 PetscErrorCode ierr; 1246 1247 PetscFunctionBegin; 1248 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr); 1249 PetscFunctionReturn(0); 1250 } 1251 1252 PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 1253 { 1254 Mat_Nest *bA = (Mat_Nest*)A->data; 1255 1256 PetscFunctionBegin; 1257 if (M) *M = bA->nr; 1258 if (N) *N = bA->nc; 1259 PetscFunctionReturn(0); 1260 } 1261 1262 /*@ 1263 MatNestGetSize - Returns the size of the nest matrix. 1264 1265 Not collective 1266 1267 Input Parameters: 1268 . A - nest matrix 1269 1270 Output Parameter: 1271 + M - number of rows in the nested mat 1272 - N - number of cols in the nested mat 1273 1274 Notes: 1275 1276 Level: developer 1277 1278 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(), 1279 MatNestGetISs() 1280 @*/ 1281 PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 1282 { 1283 PetscErrorCode ierr; 1284 1285 PetscFunctionBegin; 1286 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr); 1287 PetscFunctionReturn(0); 1288 } 1289 1290 static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[]) 1291 { 1292 Mat_Nest *vs = (Mat_Nest*)A->data; 1293 PetscInt i; 1294 1295 PetscFunctionBegin; 1296 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i]; 1297 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i]; 1298 PetscFunctionReturn(0); 1299 } 1300 1301 /*@C 1302 MatNestGetISs - Returns the index sets partitioning the row and column spaces 1303 1304 Not collective 1305 1306 Input Parameters: 1307 . A - nest matrix 1308 1309 Output Parameter: 1310 + rows - array of row index sets 1311 - cols - array of column index sets 1312 1313 Level: advanced 1314 1315 Notes: 1316 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 1317 1318 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST, 1319 MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats() 1320 @*/ 1321 PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[]) 1322 { 1323 PetscErrorCode ierr; 1324 1325 PetscFunctionBegin; 1326 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1327 ierr = PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 1328 PetscFunctionReturn(0); 1329 } 1330 1331 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[]) 1332 { 1333 Mat_Nest *vs = (Mat_Nest*)A->data; 1334 PetscInt i; 1335 1336 PetscFunctionBegin; 1337 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i]; 1338 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i]; 1339 PetscFunctionReturn(0); 1340 } 1341 1342 /*@C 1343 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces 1344 1345 Not collective 1346 1347 Input Parameters: 1348 . A - nest matrix 1349 1350 Output Parameter: 1351 + rows - array of row index sets (or NULL to ignore) 1352 - cols - array of column index sets (or NULL to ignore) 1353 1354 Level: advanced 1355 1356 Notes: 1357 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 1358 1359 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(), 1360 MATNEST, MatNestSetSubMats(), MatNestSetSubMat() 1361 @*/ 1362 PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[]) 1363 { 1364 PetscErrorCode ierr; 1365 1366 PetscFunctionBegin; 1367 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1368 ierr = PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 1369 PetscFunctionReturn(0); 1370 } 1371 1372 PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype) 1373 { 1374 PetscErrorCode ierr; 1375 PetscBool flg; 1376 1377 PetscFunctionBegin; 1378 ierr = PetscStrcmp(vtype,VECNEST,&flg);CHKERRQ(ierr); 1379 /* In reality, this only distinguishes VECNEST and "other" */ 1380 if (flg) A->ops->getvecs = MatCreateVecs_Nest; 1381 else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0; 1382 PetscFunctionReturn(0); 1383 } 1384 1385 /*@C 1386 MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs() 1387 1388 Not collective 1389 1390 Input Parameters: 1391 + A - nest matrix 1392 - vtype - type to use for creating vectors 1393 1394 Notes: 1395 1396 Level: developer 1397 1398 .seealso: MatCreateVecs(), MATNEST, MatCreateNest() 1399 @*/ 1400 PetscErrorCode MatNestSetVecType(Mat A,VecType vtype) 1401 { 1402 PetscErrorCode ierr; 1403 1404 PetscFunctionBegin; 1405 ierr = PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));CHKERRQ(ierr); 1406 PetscFunctionReturn(0); 1407 } 1408 1409 PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1410 { 1411 Mat_Nest *s = (Mat_Nest*)A->data; 1412 PetscInt i,j,m,n,M,N; 1413 PetscErrorCode ierr; 1414 PetscBool cong; 1415 1416 PetscFunctionBegin; 1417 ierr = MatReset_Nest(A);CHKERRQ(ierr); 1418 1419 s->nr = nr; 1420 s->nc = nc; 1421 1422 /* Create space for submatrices */ 1423 ierr = PetscMalloc1(nr,&s->m);CHKERRQ(ierr); 1424 for (i=0; i<nr; i++) { 1425 ierr = PetscMalloc1(nc,&s->m[i]);CHKERRQ(ierr); 1426 } 1427 for (i=0; i<nr; i++) { 1428 for (j=0; j<nc; j++) { 1429 s->m[i][j] = a[i*nc+j]; 1430 if (a[i*nc+j]) { 1431 ierr = PetscObjectReference((PetscObject)a[i*nc+j]);CHKERRQ(ierr); 1432 } 1433 } 1434 } 1435 1436 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr); 1437 1438 ierr = PetscMalloc1(nr,&s->row_len);CHKERRQ(ierr); 1439 ierr = PetscMalloc1(nc,&s->col_len);CHKERRQ(ierr); 1440 for (i=0; i<nr; i++) s->row_len[i]=-1; 1441 for (j=0; j<nc; j++) s->col_len[j]=-1; 1442 1443 ierr = PetscCalloc1(nr*nc,&s->nnzstate);CHKERRQ(ierr); 1444 for (i=0; i<nr; i++) { 1445 for (j=0; j<nc; j++) { 1446 if (s->m[i][j]) { 1447 ierr = MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);CHKERRQ(ierr); 1448 } 1449 } 1450 } 1451 1452 ierr = MatNestGetSizes_Private(A,&m,&n,&M,&N);CHKERRQ(ierr); 1453 1454 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1455 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1456 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1457 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1458 1459 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1460 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1461 1462 /* disable operations that are not supported for non-square matrices, 1463 or matrices for which is_row != is_col */ 1464 ierr = MatHasCongruentLayouts(A,&cong);CHKERRQ(ierr); 1465 if (cong && nr != nc) cong = PETSC_FALSE; 1466 if (cong) { 1467 for (i = 0; cong && i < nr; i++) { 1468 ierr = ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);CHKERRQ(ierr); 1469 } 1470 } 1471 if (!cong) { 1472 A->ops->missingdiagonal = NULL; 1473 A->ops->getdiagonal = NULL; 1474 A->ops->shift = NULL; 1475 A->ops->diagonalset = NULL; 1476 } 1477 1478 ierr = PetscCalloc2(nr,&s->left,nc,&s->right);CHKERRQ(ierr); 1479 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1480 A->nonzerostate++; 1481 PetscFunctionReturn(0); 1482 } 1483 1484 /*@ 1485 MatNestSetSubMats - Sets the nested submatrices 1486 1487 Collective on Mat 1488 1489 Input Parameter: 1490 + A - nested matrix 1491 . nr - number of nested row blocks 1492 . is_row - index sets for each nested row block, or NULL to make contiguous 1493 . nc - number of nested column blocks 1494 . is_col - index sets for each nested column block, or NULL to make contiguous 1495 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1496 1497 Notes: this always resets any submatrix information previously set 1498 1499 Level: advanced 1500 1501 .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats() 1502 @*/ 1503 PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1504 { 1505 PetscErrorCode ierr; 1506 PetscInt i; 1507 1508 PetscFunctionBegin; 1509 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1510 if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative"); 1511 if (nr && is_row) { 1512 PetscValidPointer(is_row,3); 1513 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3); 1514 } 1515 if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative"); 1516 if (nc && is_col) { 1517 PetscValidPointer(is_col,5); 1518 for (i=0; i<nc; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5); 1519 } 1520 if (nr*nc > 0) PetscValidPointer(a,6); 1521 ierr = PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));CHKERRQ(ierr); 1522 PetscFunctionReturn(0); 1523 } 1524 1525 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog) 1526 { 1527 PetscErrorCode ierr; 1528 PetscBool flg; 1529 PetscInt i,j,m,mi,*ix; 1530 1531 PetscFunctionBegin; 1532 *ltog = NULL; 1533 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1534 if (islocal[i]) { 1535 ierr = ISGetLocalSize(islocal[i],&mi);CHKERRQ(ierr); 1536 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1537 } else { 1538 ierr = ISGetLocalSize(isglobal[i],&mi);CHKERRQ(ierr); 1539 } 1540 m += mi; 1541 } 1542 if (!flg) PetscFunctionReturn(0); 1543 1544 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1545 for (i=0,m=0; i<n; i++) { 1546 ISLocalToGlobalMapping smap = NULL; 1547 Mat sub = NULL; 1548 PetscSF sf; 1549 PetscLayout map; 1550 const PetscInt *ix2; 1551 1552 if (!colflg) { 1553 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1554 } else { 1555 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1556 } 1557 if (sub) { 1558 if (!colflg) { 1559 ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr); 1560 } else { 1561 ierr = MatGetLocalToGlobalMapping(sub,NULL,&smap);CHKERRQ(ierr); 1562 } 1563 } 1564 /* 1565 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1566 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1567 */ 1568 ierr = ISGetIndices(isglobal[i],&ix2);CHKERRQ(ierr); 1569 if (islocal[i]) { 1570 PetscInt *ilocal,*iremote; 1571 PetscInt mil,nleaves; 1572 1573 ierr = ISGetLocalSize(islocal[i],&mi);CHKERRQ(ierr); 1574 if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map"); 1575 for (j=0; j<mi; j++) ix[m+j] = j; 1576 ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr); 1577 1578 /* PetscSFSetGraphLayout does not like negative indices */ 1579 ierr = PetscMalloc2(mi,&ilocal,mi,&iremote);CHKERRQ(ierr); 1580 for (j=0, nleaves = 0; j<mi; j++) { 1581 if (ix[m+j] < 0) continue; 1582 ilocal[nleaves] = j; 1583 iremote[nleaves] = ix[m+j]; 1584 nleaves++; 1585 } 1586 ierr = ISGetLocalSize(isglobal[i],&mil);CHKERRQ(ierr); 1587 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);CHKERRQ(ierr); 1588 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);CHKERRQ(ierr); 1589 ierr = PetscLayoutSetLocalSize(map,mil);CHKERRQ(ierr); 1590 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 1591 ierr = PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);CHKERRQ(ierr); 1592 ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr); 1593 ierr = PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr); 1594 ierr = PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr); 1595 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1596 ierr = PetscFree2(ilocal,iremote);CHKERRQ(ierr); 1597 } else { 1598 ierr = ISGetLocalSize(isglobal[i],&mi);CHKERRQ(ierr); 1599 for (j=0; j<mi; j++) ix[m+j] = ix2[i]; 1600 } 1601 ierr = ISRestoreIndices(isglobal[i],&ix2);CHKERRQ(ierr); 1602 m += mi; 1603 } 1604 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1605 PetscFunctionReturn(0); 1606 } 1607 1608 1609 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1610 /* 1611 nprocessors = NP 1612 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1613 proc 0: => (g_0,h_0,) 1614 proc 1: => (g_1,h_1,) 1615 ... 1616 proc nprocs-1: => (g_NP-1,h_NP-1,) 1617 1618 proc 0: proc 1: proc nprocs-1: 1619 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1620 1621 proc 0: 1622 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1623 proc 1: 1624 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1625 1626 proc NP-1: 1627 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1628 */ 1629 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1630 { 1631 Mat_Nest *vs = (Mat_Nest*)A->data; 1632 PetscInt i,j,offset,n,nsum,bs; 1633 PetscErrorCode ierr; 1634 Mat sub = NULL; 1635 1636 PetscFunctionBegin; 1637 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr); 1638 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr); 1639 if (is_row) { /* valid IS is passed in */ 1640 /* refs on is[] are incremeneted */ 1641 for (i=0; i<vs->nr; i++) { 1642 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1643 1644 vs->isglobal.row[i] = is_row[i]; 1645 } 1646 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1647 nsum = 0; 1648 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1649 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1650 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1651 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1652 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1653 nsum += n; 1654 } 1655 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRMPI(ierr); 1656 offset -= nsum; 1657 for (i=0; i<vs->nr; i++) { 1658 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1659 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1660 ierr = MatGetBlockSizes(sub,&bs,NULL);CHKERRQ(ierr); 1661 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1662 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1663 offset += n; 1664 } 1665 } 1666 1667 if (is_col) { /* valid IS is passed in */ 1668 /* refs on is[] are incremeneted */ 1669 for (j=0; j<vs->nc; j++) { 1670 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1671 1672 vs->isglobal.col[j] = is_col[j]; 1673 } 1674 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1675 offset = A->cmap->rstart; 1676 nsum = 0; 1677 for (j=0; j<vs->nc; j++) { 1678 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1679 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1680 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1681 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1682 nsum += n; 1683 } 1684 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRMPI(ierr); 1685 offset -= nsum; 1686 for (j=0; j<vs->nc; j++) { 1687 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1688 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1689 ierr = MatGetBlockSizes(sub,NULL,&bs);CHKERRQ(ierr); 1690 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1691 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1692 offset += n; 1693 } 1694 } 1695 1696 /* Set up the local ISs */ 1697 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1698 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1699 for (i=0,offset=0; i<vs->nr; i++) { 1700 IS isloc; 1701 ISLocalToGlobalMapping rmap = NULL; 1702 PetscInt nlocal,bs; 1703 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1704 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1705 if (rmap) { 1706 ierr = MatGetBlockSizes(sub,&bs,NULL);CHKERRQ(ierr); 1707 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1708 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1709 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1710 } else { 1711 nlocal = 0; 1712 isloc = NULL; 1713 } 1714 vs->islocal.row[i] = isloc; 1715 offset += nlocal; 1716 } 1717 for (i=0,offset=0; i<vs->nc; i++) { 1718 IS isloc; 1719 ISLocalToGlobalMapping cmap = NULL; 1720 PetscInt nlocal,bs; 1721 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1722 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1723 if (cmap) { 1724 ierr = MatGetBlockSizes(sub,NULL,&bs);CHKERRQ(ierr); 1725 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1726 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1727 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1728 } else { 1729 nlocal = 0; 1730 isloc = NULL; 1731 } 1732 vs->islocal.col[i] = isloc; 1733 offset += nlocal; 1734 } 1735 1736 /* Set up the aggregate ISLocalToGlobalMapping */ 1737 { 1738 ISLocalToGlobalMapping rmap,cmap; 1739 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr); 1740 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr); 1741 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1742 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1743 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1744 } 1745 1746 if (PetscDefined(USE_DEBUG)) { 1747 for (i=0; i<vs->nr; i++) { 1748 for (j=0; j<vs->nc; j++) { 1749 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1750 Mat B = vs->m[i][j]; 1751 if (!B) continue; 1752 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1753 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1754 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1755 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1756 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1757 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1758 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); 1759 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); 1760 } 1761 } 1762 } 1763 1764 /* Set A->assembled if all non-null blocks are currently assembled */ 1765 for (i=0; i<vs->nr; i++) { 1766 for (j=0; j<vs->nc; j++) { 1767 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1768 } 1769 } 1770 A->assembled = PETSC_TRUE; 1771 PetscFunctionReturn(0); 1772 } 1773 1774 /*@C 1775 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1776 1777 Collective on Mat 1778 1779 Input Parameter: 1780 + comm - Communicator for the new Mat 1781 . nr - number of nested row blocks 1782 . is_row - index sets for each nested row block, or NULL to make contiguous 1783 . nc - number of nested column blocks 1784 . is_col - index sets for each nested column block, or NULL to make contiguous 1785 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1786 1787 Output Parameter: 1788 . B - new matrix 1789 1790 Level: advanced 1791 1792 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(), 1793 MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(), 1794 MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats() 1795 @*/ 1796 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1797 { 1798 Mat A; 1799 PetscErrorCode ierr; 1800 1801 PetscFunctionBegin; 1802 *B = NULL; 1803 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1804 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1805 A->preallocated = PETSC_TRUE; 1806 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1807 *B = A; 1808 PetscFunctionReturn(0); 1809 } 1810 1811 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1812 { 1813 Mat_Nest *nest = (Mat_Nest*)A->data; 1814 Mat *trans; 1815 PetscScalar **avv; 1816 PetscScalar *vv; 1817 PetscInt **aii,**ajj; 1818 PetscInt *ii,*jj,*ci; 1819 PetscInt nr,nc,nnz,i,j; 1820 PetscBool done; 1821 PetscErrorCode ierr; 1822 1823 PetscFunctionBegin; 1824 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 1825 if (reuse == MAT_REUSE_MATRIX) { 1826 PetscInt rnr; 1827 1828 ierr = MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);CHKERRQ(ierr); 1829 if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ"); 1830 if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows"); 1831 ierr = MatSeqAIJGetArray(*newmat,&vv);CHKERRQ(ierr); 1832 } 1833 /* extract CSR for nested SeqAIJ matrices */ 1834 nnz = 0; 1835 ierr = PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);CHKERRQ(ierr); 1836 for (i=0; i<nest->nr; ++i) { 1837 for (j=0; j<nest->nc; ++j) { 1838 Mat B = nest->m[i][j]; 1839 if (B) { 1840 PetscScalar *naa; 1841 PetscInt *nii,*njj,nnr; 1842 PetscBool istrans; 1843 1844 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1845 if (istrans) { 1846 Mat Bt; 1847 1848 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1849 ierr = MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);CHKERRQ(ierr); 1850 B = trans[i*nest->nc+j]; 1851 } 1852 ierr = MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);CHKERRQ(ierr); 1853 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ"); 1854 ierr = MatSeqAIJGetArray(B,&naa);CHKERRQ(ierr); 1855 nnz += nii[nnr]; 1856 1857 aii[i*nest->nc+j] = nii; 1858 ajj[i*nest->nc+j] = njj; 1859 avv[i*nest->nc+j] = naa; 1860 } 1861 } 1862 } 1863 if (reuse != MAT_REUSE_MATRIX) { 1864 ierr = PetscMalloc1(nr+1,&ii);CHKERRQ(ierr); 1865 ierr = PetscMalloc1(nnz,&jj);CHKERRQ(ierr); 1866 ierr = PetscMalloc1(nnz,&vv);CHKERRQ(ierr); 1867 } else { 1868 if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros"); 1869 } 1870 1871 /* new row pointer */ 1872 ierr = PetscArrayzero(ii,nr+1);CHKERRQ(ierr); 1873 for (i=0; i<nest->nr; ++i) { 1874 PetscInt ncr,rst; 1875 1876 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1877 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1878 for (j=0; j<nest->nc; ++j) { 1879 if (aii[i*nest->nc+j]) { 1880 PetscInt *nii = aii[i*nest->nc+j]; 1881 PetscInt ir; 1882 1883 for (ir=rst; ir<ncr+rst; ++ir) { 1884 ii[ir+1] += nii[1]-nii[0]; 1885 nii++; 1886 } 1887 } 1888 } 1889 } 1890 for (i=0; i<nr; i++) ii[i+1] += ii[i]; 1891 1892 /* construct CSR for the new matrix */ 1893 ierr = PetscCalloc1(nr,&ci);CHKERRQ(ierr); 1894 for (i=0; i<nest->nr; ++i) { 1895 PetscInt ncr,rst; 1896 1897 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1898 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1899 for (j=0; j<nest->nc; ++j) { 1900 if (aii[i*nest->nc+j]) { 1901 PetscScalar *nvv = avv[i*nest->nc+j]; 1902 PetscInt *nii = aii[i*nest->nc+j]; 1903 PetscInt *njj = ajj[i*nest->nc+j]; 1904 PetscInt ir,cst; 1905 1906 ierr = ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);CHKERRQ(ierr); 1907 for (ir=rst; ir<ncr+rst; ++ir) { 1908 PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir]; 1909 1910 for (ij=0;ij<rsize;ij++) { 1911 jj[ist+ij] = *njj+cst; 1912 vv[ist+ij] = *nvv; 1913 njj++; 1914 nvv++; 1915 } 1916 ci[ir] += rsize; 1917 nii++; 1918 } 1919 } 1920 } 1921 } 1922 ierr = PetscFree(ci);CHKERRQ(ierr); 1923 1924 /* restore info */ 1925 for (i=0; i<nest->nr; ++i) { 1926 for (j=0; j<nest->nc; ++j) { 1927 Mat B = nest->m[i][j]; 1928 if (B) { 1929 PetscInt nnr = 0, k = i*nest->nc+j; 1930 1931 B = (trans[k] ? trans[k] : B); 1932 ierr = MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);CHKERRQ(ierr); 1933 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ"); 1934 ierr = MatSeqAIJRestoreArray(B,&avv[k]);CHKERRQ(ierr); 1935 ierr = MatDestroy(&trans[k]);CHKERRQ(ierr); 1936 } 1937 } 1938 } 1939 ierr = PetscFree4(aii,ajj,avv,trans);CHKERRQ(ierr); 1940 1941 /* finalize newmat */ 1942 if (reuse == MAT_INITIAL_MATRIX) { 1943 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);CHKERRQ(ierr); 1944 } else if (reuse == MAT_INPLACE_MATRIX) { 1945 Mat B; 1946 1947 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);CHKERRQ(ierr); 1948 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1949 } 1950 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1951 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1952 { 1953 Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data); 1954 a->free_a = PETSC_TRUE; 1955 a->free_ij = PETSC_TRUE; 1956 } 1957 PetscFunctionReturn(0); 1958 } 1959 1960 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1961 { 1962 PetscErrorCode ierr; 1963 Mat_Nest *nest = (Mat_Nest*)A->data; 1964 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1965 PetscInt cstart,cend; 1966 PetscMPIInt size; 1967 Mat C; 1968 1969 PetscFunctionBegin; 1970 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRMPI(ierr); 1971 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1972 PetscInt nf; 1973 PetscBool fast; 1974 1975 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr); 1976 if (!fast) { 1977 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr); 1978 } 1979 for (i=0; i<nest->nr && fast; ++i) { 1980 for (j=0; j<nest->nc && fast; ++j) { 1981 Mat B = nest->m[i][j]; 1982 if (B) { 1983 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr); 1984 if (!fast) { 1985 PetscBool istrans; 1986 1987 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1988 if (istrans) { 1989 Mat Bt; 1990 1991 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1992 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);CHKERRQ(ierr); 1993 } 1994 } 1995 } 1996 } 1997 } 1998 for (i=0, nf=0; i<nest->nr && fast; ++i) { 1999 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr); 2000 if (fast) { 2001 PetscInt f,s; 2002 2003 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr); 2004 if (f != nf || s != 1) { fast = PETSC_FALSE; } 2005 else { 2006 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr); 2007 nf += f; 2008 } 2009 } 2010 } 2011 for (i=0, nf=0; i<nest->nc && fast; ++i) { 2012 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr); 2013 if (fast) { 2014 PetscInt f,s; 2015 2016 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr); 2017 if (f != nf || s != 1) { fast = PETSC_FALSE; } 2018 else { 2019 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr); 2020 nf += f; 2021 } 2022 } 2023 } 2024 if (fast) { 2025 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr); 2026 PetscFunctionReturn(0); 2027 } 2028 } 2029 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 2030 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 2031 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 2032 switch (reuse) { 2033 case MAT_INITIAL_MATRIX: 2034 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2035 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 2036 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 2037 *newmat = C; 2038 break; 2039 case MAT_REUSE_MATRIX: 2040 C = *newmat; 2041 break; 2042 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 2043 } 2044 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 2045 onnz = dnnz + m; 2046 for (k=0; k<m; k++) { 2047 dnnz[k] = 0; 2048 onnz[k] = 0; 2049 } 2050 for (j=0; j<nest->nc; ++j) { 2051 IS bNis; 2052 PetscInt bN; 2053 const PetscInt *bNindices; 2054 /* Using global column indices and ISAllGather() is not scalable. */ 2055 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 2056 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 2057 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 2058 for (i=0; i<nest->nr; ++i) { 2059 PetscSF bmsf; 2060 PetscSFNode *iremote; 2061 Mat B; 2062 PetscInt bm, *sub_dnnz,*sub_onnz, br; 2063 const PetscInt *bmindices; 2064 B = nest->m[i][j]; 2065 if (!B) continue; 2066 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 2067 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 2068 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 2069 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 2070 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 2071 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 2072 for (k = 0; k < bm; ++k){ 2073 sub_dnnz[k] = 0; 2074 sub_onnz[k] = 0; 2075 } 2076 /* 2077 Locate the owners for all of the locally-owned global row indices for this row block. 2078 These determine the roots of PetscSF used to communicate preallocation data to row owners. 2079 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 2080 */ 2081 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 2082 for (br = 0; br < bm; ++br) { 2083 PetscInt row = bmindices[br], brncols, col; 2084 const PetscInt *brcols; 2085 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 2086 PetscMPIInt rowowner = 0; 2087 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 2088 /* how many roots */ 2089 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 2090 /* get nonzero pattern */ 2091 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 2092 for (k=0; k<brncols; k++) { 2093 col = bNindices[brcols[k]]; 2094 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 2095 sub_dnnz[br]++; 2096 } else { 2097 sub_onnz[br]++; 2098 } 2099 } 2100 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 2101 } 2102 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 2103 /* bsf will have to take care of disposing of bedges. */ 2104 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 2105 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 2106 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 2107 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 2108 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 2109 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 2110 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 2111 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 2112 } 2113 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 2114 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 2115 } 2116 /* Resize preallocation if overestimated */ 2117 for (i=0;i<m;i++) { 2118 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 2119 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 2120 } 2121 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 2122 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 2123 ierr = PetscFree(dnnz);CHKERRQ(ierr); 2124 2125 /* Fill by row */ 2126 for (j=0; j<nest->nc; ++j) { 2127 /* Using global column indices and ISAllGather() is not scalable. */ 2128 IS bNis; 2129 PetscInt bN; 2130 const PetscInt *bNindices; 2131 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 2132 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 2133 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 2134 for (i=0; i<nest->nr; ++i) { 2135 Mat B; 2136 PetscInt bm, br; 2137 const PetscInt *bmindices; 2138 B = nest->m[i][j]; 2139 if (!B) continue; 2140 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 2141 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 2142 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 2143 for (br = 0; br < bm; ++br) { 2144 PetscInt row = bmindices[br], brncols, *cols; 2145 const PetscInt *brcols; 2146 const PetscScalar *brcoldata; 2147 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 2148 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 2149 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 2150 /* 2151 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 2152 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 2153 */ 2154 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 2155 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 2156 ierr = PetscFree(cols);CHKERRQ(ierr); 2157 } 2158 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 2159 } 2160 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 2161 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 2162 } 2163 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2164 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2165 PetscFunctionReturn(0); 2166 } 2167 2168 PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has) 2169 { 2170 Mat_Nest *bA = (Mat_Nest*)mat->data; 2171 MatOperation opAdd; 2172 PetscInt i,j,nr = bA->nr,nc = bA->nc; 2173 PetscBool flg; 2174 PetscErrorCode ierr; 2175 PetscFunctionBegin; 2176 2177 *has = PETSC_FALSE; 2178 if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 2179 opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD); 2180 for (j=0; j<nc; j++) { 2181 for (i=0; i<nr; i++) { 2182 if (!bA->m[i][j]) continue; 2183 ierr = MatHasOperation(bA->m[i][j],opAdd,&flg);CHKERRQ(ierr); 2184 if (!flg) PetscFunctionReturn(0); 2185 } 2186 } 2187 } 2188 if (((void**)mat->ops)[op]) *has = PETSC_TRUE; 2189 PetscFunctionReturn(0); 2190 } 2191 2192 /*MC 2193 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 2194 2195 Level: intermediate 2196 2197 Notes: 2198 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 2199 It allows the use of symmetric and block formats for parts of multi-physics simulations. 2200 It is usually used with DMComposite and DMCreateMatrix() 2201 2202 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 2203 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 2204 than the nest matrix. 2205 2206 .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(), 2207 VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(), 2208 MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats() 2209 M*/ 2210 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 2211 { 2212 Mat_Nest *s; 2213 PetscErrorCode ierr; 2214 2215 PetscFunctionBegin; 2216 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 2217 A->data = (void*)s; 2218 2219 s->nr = -1; 2220 s->nc = -1; 2221 s->m = NULL; 2222 s->splitassembly = PETSC_FALSE; 2223 2224 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 2225 2226 A->ops->mult = MatMult_Nest; 2227 A->ops->multadd = MatMultAdd_Nest; 2228 A->ops->multtranspose = MatMultTranspose_Nest; 2229 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 2230 A->ops->transpose = MatTranspose_Nest; 2231 A->ops->assemblybegin = MatAssemblyBegin_Nest; 2232 A->ops->assemblyend = MatAssemblyEnd_Nest; 2233 A->ops->zeroentries = MatZeroEntries_Nest; 2234 A->ops->copy = MatCopy_Nest; 2235 A->ops->axpy = MatAXPY_Nest; 2236 A->ops->duplicate = MatDuplicate_Nest; 2237 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 2238 A->ops->destroy = MatDestroy_Nest; 2239 A->ops->view = MatView_Nest; 2240 A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 2241 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 2242 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 2243 A->ops->getdiagonal = MatGetDiagonal_Nest; 2244 A->ops->diagonalscale = MatDiagonalScale_Nest; 2245 A->ops->scale = MatScale_Nest; 2246 A->ops->shift = MatShift_Nest; 2247 A->ops->diagonalset = MatDiagonalSet_Nest; 2248 A->ops->setrandom = MatSetRandom_Nest; 2249 A->ops->hasoperation = MatHasOperation_Nest; 2250 A->ops->missingdiagonal = MatMissingDiagonal_Nest; 2251 2252 A->spptr = NULL; 2253 A->assembled = PETSC_FALSE; 2254 2255 /* expose Nest api's */ 2256 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 2257 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 2258 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 2259 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 2260 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 2261 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 2262 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 2263 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 2264 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr); 2265 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr); 2266 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr); 2267 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr); 2268 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);CHKERRQ(ierr); 2269 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);CHKERRQ(ierr); 2270 ierr = PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);CHKERRQ(ierr); 2271 2272 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 2273 PetscFunctionReturn(0); 2274 } 2275