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