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