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