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