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