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 CHKERRQ(ISGetLocalSize(bA->isglobal.row[i],&sm)); 22 CHKERRQ(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 CHKERRQ(ISGetLocalSize(bA->isglobal.col[j],&sn)); 29 CHKERRQ(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++) CHKERRQ(VecGetSubVector(y,bA->isglobal.row[i],&by[i])); 45 for (i=0; i<nc; i++) CHKERRQ(VecGetSubVector(x,bA->isglobal.col[i],&bx[i])); 46 for (i=0; i<nr; i++) { 47 CHKERRQ(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 CHKERRQ(MatMultAdd(bA->m[i][j],bx[j],by[i],by[i])); 52 } 53 } 54 for (i=0; i<nr; i++) CHKERRQ(VecRestoreSubVector(y,bA->isglobal.row[i],&by[i])); 55 for (i=0; i<nc; i++) CHKERRQ(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++) CHKERRQ(VecGetSubVector(z,bA->isglobal.row[i],&bz[i])); 67 for (i=0; i<nc; i++) CHKERRQ(VecGetSubVector(x,bA->isglobal.col[i],&bx[i])); 68 for (i=0; i<nr; i++) { 69 if (y != z) { 70 Vec by; 71 CHKERRQ(VecGetSubVector(y,bA->isglobal.row[i],&by)); 72 CHKERRQ(VecCopy(by,bz[i])); 73 CHKERRQ(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 CHKERRQ(MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i])); 79 } 80 } 81 for (i=0; i<nr; i++) CHKERRQ(VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i])); 82 for (i=0; i<nc; i++) CHKERRQ(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 CHKERRQ(MatGetSize(B,NULL,&N)); 107 if (!N) { 108 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 109 CHKERRQ(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 CHKERRQ(MatDenseGetLDA(B,&ldb)); 118 CHKERRQ(MatDenseGetLDA(C,&ldc)); 119 CHKERRQ(MatZeroEntries(C)); 120 CHKERRQ(MatDenseGetArrayRead(B,&barray)); 121 CHKERRQ(MatDenseGetArray(C,&carray)); 122 for (i=0; i<nr; i++) { 123 CHKERRQ(ISGetSize(bA->isglobal.row[i],&M)); 124 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC)); 125 CHKERRQ(MatDenseSetLDA(viewC,ldc)); 126 for (j=0; j<nc; j++) { 127 if (!bA->m[i][j]) continue; 128 CHKERRQ(ISGetSize(bA->isglobal.col[j],&M)); 129 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB)); 130 CHKERRQ(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 CHKERRQ(MatProductNumeric(workC)); 137 CHKERRQ(MatDestroy(&viewB)); 138 workC->product->B = productB; /* resume original B */ 139 140 /* C[i] <- workC + C[i] */ 141 CHKERRQ(MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN)); 142 } 143 CHKERRQ(MatDestroy(&viewC)); 144 } 145 CHKERRQ(MatDenseRestoreArray(C,&carray)); 146 CHKERRQ(MatDenseRestoreArrayRead(B,&barray)); 147 148 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 149 CHKERRQ(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 CHKERRQ(PetscFree(contents->tarray)); 160 for (i=0; i<contents->k; i++) { 161 CHKERRQ(MatDestroy(contents->workC + i)); 162 } 163 CHKERRQ(PetscFree3(contents->dm,contents->dn,contents->workC)); 164 CHKERRQ(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 CHKERRQ(MatGetLocalSize(C,&m,&n)); 189 CHKERRQ(MatGetSize(C,&M,&N)); 190 if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) { 191 CHKERRQ(MatGetLocalSize(B,NULL,&n)); 192 CHKERRQ(MatGetSize(B,NULL,&N)); 193 CHKERRQ(MatGetLocalSize(A,&m,NULL)); 194 CHKERRQ(MatGetSize(A,&M,NULL)); 195 CHKERRQ(MatSetSizes(C,m,n,M,N)); 196 } 197 CHKERRQ(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATMPIDENSE,MATSEQDENSECUDA,MATMPIDENSECUDA,"")); 198 if (!cisdense) { 199 CHKERRQ(MatSetType(C,((PetscObject)B)->type_name)); 200 } 201 CHKERRQ(MatSetUp(C)); 202 if (!N) { 203 C->ops->productnumeric = MatProductNumeric_Nest_Dense; 204 PetscFunctionReturn(0); 205 } 206 207 CHKERRQ(PetscNew(&contents)); 208 C->product->data = contents; 209 C->product->destroy = MatNest_DenseDestroy; 210 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1)); 219 contents->dn[i+1] += contents->dn[i]; 220 } 221 CHKERRQ(PetscMalloc1(maxm*N,&contents->tarray)); 222 CHKERRQ(MatDenseGetLDA(B,&ldb)); 223 CHKERRQ(MatGetSize(B,NULL,&N)); 224 CHKERRQ(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 CHKERRQ(ISGetSize(bA->isglobal.col[j],&M)); 228 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB)); 229 CHKERRQ(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 CHKERRQ(MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j])); 235 workC = contents->workC[i*nc + j]; 236 CHKERRQ(MatProductSetType(workC,MATPRODUCT_AB)); 237 CHKERRQ(MatProductSetAlgorithm(workC,"default")); 238 CHKERRQ(MatProductSetFill(workC,fill)); 239 CHKERRQ(MatProductSetFromOptions(workC)); 240 CHKERRQ(MatProductSymbolic(workC)); 241 242 /* since tarray will be shared by all Mat */ 243 CHKERRQ(MatSeqDenseSetPreallocation(workC,contents->tarray)); 244 CHKERRQ(MatMPIDenseSetPreallocation(workC,contents->tarray)); 245 } 246 CHKERRQ(MatDestroy(&viewB)); 247 } 248 CHKERRQ(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 CHKERRQ(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++) CHKERRQ(VecGetSubVector(x,bA->isglobal.row[i],&bx[i])); 282 for (i=0; i<nc; i++) CHKERRQ(VecGetSubVector(y,bA->isglobal.col[i],&by[i])); 283 for (j=0; j<nc; j++) { 284 CHKERRQ(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 CHKERRQ(MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j])); 289 } 290 } 291 for (i=0; i<nr; i++) CHKERRQ(VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i])); 292 for (i=0; i<nc; i++) CHKERRQ(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++) CHKERRQ(VecGetSubVector(x,bA->isglobal.row[i],&bx[i])); 304 for (i=0; i<nc; i++) CHKERRQ(VecGetSubVector(z,bA->isglobal.col[i],&bz[i])); 305 for (j=0; j<nc; j++) { 306 if (y != z) { 307 Vec by; 308 CHKERRQ(VecGetSubVector(y,bA->isglobal.col[j],&by)); 309 CHKERRQ(VecCopy(by,bz[j])); 310 CHKERRQ(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 CHKERRQ(MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j])); 316 } 317 } 318 for (i=0; i<nr; i++) CHKERRQ(VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i])); 319 for (i=0; i<nc; i++) CHKERRQ(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 PetscCheckFalse(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 CHKERRQ(PetscCalloc1(nr * nc,&subs)); 337 CHKERRQ(PetscMalloc2(nr,&is_row,nc,&is_col)); 338 CHKERRQ(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 CHKERRQ(MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C)); 348 CHKERRQ(PetscFree(subs)); 349 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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]) CHKERRQ(ISDestroy(&lst[i])); 381 CHKERRQ(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 CHKERRQ(MatNestDestroyISList(vs->nr,&vs->isglobal.row)); 394 CHKERRQ(MatNestDestroyISList(vs->nc,&vs->isglobal.col)); 395 CHKERRQ(MatNestDestroyISList(vs->nr,&vs->islocal.row)); 396 CHKERRQ(MatNestDestroyISList(vs->nc,&vs->islocal.col)); 397 398 CHKERRQ(PetscFree(vs->row_len)); 399 CHKERRQ(PetscFree(vs->col_len)); 400 CHKERRQ(PetscFree(vs->nnzstate)); 401 402 CHKERRQ(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 CHKERRQ(MatDestroy(&vs->m[i][j])); 409 } 410 CHKERRQ(PetscFree(vs->m[i])); 411 } 412 CHKERRQ(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 CHKERRQ(MatReset_Nest(A)); 426 CHKERRQ(PetscFree(A->data)); 427 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",NULL)); 428 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",NULL)); 429 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",NULL)); 430 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",NULL)); 431 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",NULL)); 432 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",NULL)); 433 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",NULL)); 434 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",NULL)); 435 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",NULL)); 436 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",NULL)); 437 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",NULL)); 438 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",NULL)); 439 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",NULL)); 440 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",NULL)); 441 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL)); 442 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL)); 443 CHKERRQ(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 CHKERRQ(MatMissingDiagonal(vs->m[i][i],missing,NULL)); 463 PetscCheckFalse(*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 CHKERRQ(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 CHKERRQ(MatAssemblyEnd(vs->m[i][j],type)); 489 CHKERRQ(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 CHKERRQ(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) CHKERRQ(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) CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISGetSize(is,&size)); 573 for (i=0; i<n-1; i++) { 574 if (!list[i]) continue; 575 m = 0; 576 CHKERRQ(ISConcatenate(PetscObjectComm((PetscObject)A),2,list+i,&out)); 577 CHKERRQ(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 CHKERRQ(ISConcatenate(PetscObjectComm((PetscObject)A),2,concatenate,&out)); 583 CHKERRQ(ISDestroy(concatenate)); 584 CHKERRQ(ISGetSize(out,&m)); 585 } 586 } 587 if (m == size) { 588 CHKERRQ(ISEqualUnsorted(out,is,&flg)); 589 if (flg) { 590 if (begin) *begin = i; 591 if (end) *end = j; 592 CHKERRQ(ISDestroy(&out)); 593 PetscFunctionReturn(0); 594 } 595 } 596 CHKERRQ(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 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),B)); 608 CHKERRQ(ISGetLocalSize(vs->isglobal.row[i],&lr)); 609 CHKERRQ(ISGetLocalSize(vs->isglobal.col[j],&lc)); 610 CHKERRQ(MatSetSizes(*B,lr,lc,PETSC_DECIDE,PETSC_DECIDE)); 611 CHKERRQ(MatSetType(*B,MATAIJ)); 612 CHKERRQ(MatSeqAIJSetPreallocation(*B,0,NULL)); 613 CHKERRQ(MatMPIAIJSetPreallocation(*B,0,NULL,0,NULL)); 614 CHKERRQ(MatSetUp(*B)); 615 CHKERRQ(MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE)); 616 CHKERRQ(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 617 CHKERRQ(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 CHKERRQ(PetscSNPrintf(keyname,sizeof(keyname),"NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT,rbegin,rend,cbegin,cend)); 633 CHKERRQ(PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B)); 634 if (*B) PetscFunctionReturn(0); 635 636 CHKERRQ(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 CHKERRQ(MatNestFillEmptyMat_Private(A,rbegin+i,cbegin+j,a + i*nc + j)); 659 } 660 } 661 } 662 } 663 CHKERRQ(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 CHKERRQ(MatDestroy(a + i*nc + j)); 668 } 669 } 670 } 671 CHKERRQ(PetscFree2(a,b)); 672 (*B)->assembled = A->assembled; 673 CHKERRQ(PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B)); 674 CHKERRQ(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 CHKERRQ(MatNestFindISRange(A,vs->nr,is->row,isrow,&rbegin,&rend)); 685 CHKERRQ(MatNestFindISRange(A,vs->nc,is->col,iscol,&cbegin,&cend)); 686 if (rend == rbegin + 1 && cend == cbegin + 1) { 687 if (!vs->m[rbegin][cbegin]) { 688 CHKERRQ(MatNestFillEmptyMat_Private(A,rbegin,cbegin,vs->m[rbegin] + cbegin)); 689 } 690 *B = vs->m[rbegin][cbegin]; 691 } else if (rbegin != -1 && cbegin != -1) { 692 CHKERRQ(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 CHKERRQ(MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub)); 707 switch (reuse) { 708 case MAT_INITIAL_MATRIX: 709 if (sub) CHKERRQ(PetscObjectReference((PetscObject)sub)); 710 *B = sub; 711 break; 712 case MAT_REUSE_MATRIX: 713 PetscCheckFalse(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 CHKERRQ(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) CHKERRQ(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 CHKERRQ(MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub)); 743 PetscCheckFalse(*B != sub,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten"); 744 if (sub) { 745 PetscCheckFalse(((PetscObject)sub)->refct <= 1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times"); 746 CHKERRQ(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 CHKERRQ(VecGetSubVector(v,bA->isglobal.row[i],&bv)); 760 if (bA->m[i][i]) { 761 CHKERRQ(MatGetDiagonal(bA->m[i][i],bv)); 762 } else { 763 CHKERRQ(VecSet(bv,0.0)); 764 } 765 CHKERRQ(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 CHKERRQ(PetscCalloc1(bA->nc,&br)); 778 if (r) { 779 for (j=0; j<bA->nc; j++) CHKERRQ(VecGetSubVector(r,bA->isglobal.col[j],&br[j])); 780 } 781 bl = NULL; 782 for (i=0; i<bA->nr; i++) { 783 if (l) { 784 CHKERRQ(VecGetSubVector(l,bA->isglobal.row[i],&bl)); 785 } 786 for (j=0; j<bA->nc; j++) { 787 if (bA->m[i][j]) { 788 CHKERRQ(MatDiagonalScale(bA->m[i][j],bl,br[j])); 789 } 790 } 791 if (l) { 792 CHKERRQ(VecRestoreSubVector(l,bA->isglobal.row[i],&bl)); 793 } 794 } 795 if (r) { 796 for (j=0; j<bA->nc; j++) CHKERRQ(VecRestoreSubVector(r,bA->isglobal.col[j],&br[j])); 797 } 798 CHKERRQ(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 CHKERRQ(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 PetscCheckFalse(!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 CHKERRQ(MatShift(bA->m[i][i],a)); 829 CHKERRQ(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 CHKERRQ(VecGetSubVector(D,bA->isglobal.row[i],&bv)); 848 if (bA->m[i][i]) { 849 CHKERRQ(MatDiagonalSet(bA->m[i][i],bv,is)); 850 CHKERRQ(MatGetNonzeroState(bA->m[i][i],&subnnzstate)); 851 } 852 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)A,&comm)); 885 if (right) { 886 /* allocate R */ 887 CHKERRQ(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 CHKERRQ(MatCreateVecs(bA->m[i][j],&R[j],NULL)); 893 break; 894 } 895 } 896 PetscCheckFalse(i==bA->nr,PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 897 } 898 CHKERRQ(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 CHKERRQ(VecDestroy(&R[j])); 902 } 903 CHKERRQ(PetscFree(R)); 904 } 905 906 if (left) { 907 /* allocate L */ 908 CHKERRQ(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 CHKERRQ(MatCreateVecs(bA->m[i][j],NULL,&L[i])); 914 break; 915 } 916 } 917 PetscCheckFalse(j==bA->nc,PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 918 } 919 920 CHKERRQ(VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left)); 921 for (i=0; i<bA->nr; i++) { 922 CHKERRQ(VecDestroy(&L[i])); 923 } 924 925 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 938 if (isascii) { 939 940 CHKERRQ(PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL)); 941 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Matrix object: \n")); 942 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 943 CHKERRQ(PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n",bA->nr,bA->nc)); 944 945 CHKERRQ(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 CHKERRQ(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL \n",i,j)); 955 continue; 956 } 957 CHKERRQ(MatGetSize(bA->m[i][j],&NR,&NC)); 958 CHKERRQ(MatGetType(bA->m[i][j], &type)); 959 if (((PetscObject)bA->m[i][j])->name) CHKERRQ(PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name)); 960 if (((PetscObject)bA->m[i][j])->prefix) CHKERRQ(PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix)); 961 CHKERRQ(PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest)); 962 963 CHKERRQ(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 CHKERRQ(PetscViewerASCIIPushTab(viewer)); /* push1 */ 967 CHKERRQ(MatView(bA->m[i][j],viewer)); 968 CHKERRQ(PetscViewerASCIIPopTab(viewer)); /* pop1 */ 969 } 970 } 971 } 972 CHKERRQ(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 CHKERRQ(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 PetscCheckFalse(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 CHKERRQ(MatCopy(bA->m[i][j],bB->m[i][j],str)); 1005 } else PetscCheckFalse(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 CHKERRQ(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 PetscCheckFalse(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 CHKERRQ(MatAXPY(bY->m[i][j],a,bX->m[i][j],str)); 1028 } else if (bX->m[i][j]) { 1029 Mat M; 1030 1031 PetscCheckFalse(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 CHKERRQ(MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M)); 1033 CHKERRQ(MatNestSetSubMat(Y,i,j,M)); 1034 CHKERRQ(MatDestroy(&M)); 1035 } 1036 if (bY->m[i][j]) CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatDuplicate(bA->m[i][j],op,&b[i*nc+j])); 1057 } else { 1058 b[i*nc+j] = NULL; 1059 } 1060 } 1061 } 1062 CHKERRQ(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 CHKERRQ(MatDestroy(&b[i])); 1066 } 1067 CHKERRQ(PetscFree(b)); 1068 1069 CHKERRQ(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 1070 CHKERRQ(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 PetscCheckFalse(idxm >= bA->nr,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,idxm,bA->nr-1); 1081 PetscCheckFalse(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 CHKERRQ(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 PetscCheckFalse(idxm >= bA->nr,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,idxm,bA->nr-1); 1118 PetscCheckFalse(jdxm >= bA->nc,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT,jdxm,bA->nc-1); 1119 CHKERRQ(MatGetLocalSize(mat,&m,&n)); 1120 CHKERRQ(MatGetSize(mat,&M,&N)); 1121 CHKERRQ(ISGetLocalSize(bA->isglobal.row[idxm],&mi)); 1122 CHKERRQ(ISGetSize(bA->isglobal.row[idxm],&Mi)); 1123 CHKERRQ(ISGetLocalSize(bA->isglobal.col[jdxm],&ni)); 1124 CHKERRQ(ISGetSize(bA->isglobal.col[jdxm],&Ni)); 1125 PetscCheckFalse(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 PetscCheckFalse(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 CHKERRQ(PetscObjectReference((PetscObject)mat)); 1132 CHKERRQ(MatDestroy(&bA->m[idxm][jdxm])); 1133 bA->m[idxm][jdxm] = mat; 1134 CHKERRQ(PetscObjectStateIncrease((PetscObject)A)); 1135 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatReset_Nest(A)); 1369 1370 s->nr = nr; 1371 s->nc = nc; 1372 1373 /* Create space for submatrices */ 1374 CHKERRQ(PetscMalloc1(nr,&s->m)); 1375 for (i=0; i<nr; i++) { 1376 CHKERRQ(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 CHKERRQ(PetscObjectReference((PetscObject)a[i*nc+j])); 1383 } 1384 } 1385 } 1386 CHKERRQ(MatGetVecType(A,&vtype)); 1387 CHKERRQ(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 CHKERRQ(MatGetVecType(a[i*nc+j],&vtype)); 1396 sametype = PETSC_TRUE; 1397 } else if (sametype) { 1398 CHKERRQ(MatGetVecType(a[i*nc+j],&type)); 1399 CHKERRQ(PetscStrcmp(vtype,type,&sametype)); 1400 } 1401 } 1402 } 1403 } 1404 if (sametype) { /* propagate vectype */ 1405 CHKERRQ(MatSetVecType(A,vtype)); 1406 } 1407 } 1408 1409 CHKERRQ(MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col)); 1410 1411 CHKERRQ(PetscMalloc1(nr,&s->row_len)); 1412 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j])); 1421 } 1422 } 1423 } 1424 1425 CHKERRQ(MatNestGetSizes_Private(A,&m,&n,&M,&N)); 1426 1427 CHKERRQ(PetscLayoutSetSize(A->rmap,M)); 1428 CHKERRQ(PetscLayoutSetLocalSize(A->rmap,m)); 1429 CHKERRQ(PetscLayoutSetSize(A->cmap,N)); 1430 CHKERRQ(PetscLayoutSetLocalSize(A->cmap,n)); 1431 1432 CHKERRQ(PetscLayoutSetUp(A->rmap)); 1433 CHKERRQ(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 CHKERRQ(MatHasCongruentLayouts(A,&cong)); 1438 if (cong && nr != nc) cong = PETSC_FALSE; 1439 if (cong) { 1440 for (i = 0; cong && i < nr; i++) { 1441 CHKERRQ(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 CHKERRQ(PetscCalloc2(nr,&s->left,nc,&s->right)); 1452 CHKERRQ(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 PetscCheckFalse(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 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(ISGetLocalSize(islocal[i],&mi)); 1507 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1508 } else { 1509 CHKERRQ(ISGetLocalSize(isglobal[i],&mi)); 1510 } 1511 m += mi; 1512 } 1513 if (!flg) PetscFunctionReturn(0); 1514 1515 CHKERRQ(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 CHKERRQ(MatNestFindNonzeroSubMatRow(A,i,&sub)); 1525 } else { 1526 CHKERRQ(MatNestFindNonzeroSubMatCol(A,i,&sub)); 1527 } 1528 if (sub) { 1529 if (!colflg) { 1530 CHKERRQ(MatGetLocalToGlobalMapping(sub,&smap,NULL)); 1531 } else { 1532 CHKERRQ(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 CHKERRQ(ISGetIndices(isglobal[i],&ix2)); 1540 if (islocal[i]) { 1541 PetscInt *ilocal,*iremote; 1542 PetscInt mil,nleaves; 1543 1544 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m)); 1548 1549 /* PetscSFSetGraphLayout does not like negative indices */ 1550 CHKERRQ(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 CHKERRQ(ISGetLocalSize(isglobal[i],&mil)); 1558 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)A),&sf)); 1559 CHKERRQ(PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map)); 1560 CHKERRQ(PetscLayoutSetLocalSize(map,mil)); 1561 CHKERRQ(PetscLayoutSetUp(map)); 1562 CHKERRQ(PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote)); 1563 CHKERRQ(PetscLayoutDestroy(&map)); 1564 CHKERRQ(PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE)); 1565 CHKERRQ(PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE)); 1566 CHKERRQ(PetscSFDestroy(&sf)); 1567 CHKERRQ(PetscFree2(ilocal,iremote)); 1568 } else { 1569 CHKERRQ(ISGetLocalSize(isglobal[i],&mi)); 1570 for (j=0; j<mi; j++) ix[m+j] = ix2[i]; 1571 } 1572 CHKERRQ(ISRestoreIndices(isglobal[i],&ix2)); 1573 m += mi; 1574 } 1575 CHKERRQ(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 CHKERRQ(PetscMalloc1(nr,&vs->isglobal.row)); 1607 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatNestFindNonzeroSubMatRow(A,i,&sub)); 1619 PetscCheck(sub,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %" PetscInt_FMT,i); 1620 CHKERRQ(MatGetLocalSize(sub,&n,NULL)); 1621 PetscCheckFalse(n < 0,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1622 nsum += n; 1623 } 1624 CHKERRMPI(MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A))); 1625 offset -= nsum; 1626 for (i=0; i<vs->nr; i++) { 1627 CHKERRQ(MatNestFindNonzeroSubMatRow(A,i,&sub)); 1628 CHKERRQ(MatGetLocalSize(sub,&n,NULL)); 1629 CHKERRQ(MatGetBlockSizes(sub,&bs,NULL)); 1630 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i])); 1631 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatNestFindNonzeroSubMatCol(A,j,&sub)); 1648 PetscCheck(sub,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %" PetscInt_FMT,i); 1649 CHKERRQ(MatGetLocalSize(sub,NULL,&n)); 1650 PetscCheckFalse(n < 0,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1651 nsum += n; 1652 } 1653 CHKERRMPI(MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A))); 1654 offset -= nsum; 1655 for (j=0; j<vs->nc; j++) { 1656 CHKERRQ(MatNestFindNonzeroSubMatCol(A,j,&sub)); 1657 CHKERRQ(MatGetLocalSize(sub,NULL,&n)); 1658 CHKERRQ(MatGetBlockSizes(sub,NULL,&bs)); 1659 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j])); 1660 CHKERRQ(ISSetBlockSize(vs->isglobal.col[j],bs)); 1661 offset += n; 1662 } 1663 } 1664 1665 /* Set up the local ISs */ 1666 CHKERRQ(PetscMalloc1(vs->nr,&vs->islocal.row)); 1667 CHKERRQ(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 CHKERRQ(MatNestFindNonzeroSubMatRow(A,i,&sub)); 1673 if (sub) CHKERRQ(MatGetLocalToGlobalMapping(sub,&rmap,NULL)); 1674 if (rmap) { 1675 CHKERRQ(MatGetBlockSizes(sub,&bs,NULL)); 1676 CHKERRQ(ISLocalToGlobalMappingGetSize(rmap,&nlocal)); 1677 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc)); 1678 CHKERRQ(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 CHKERRQ(MatNestFindNonzeroSubMatCol(A,i,&sub)); 1691 if (sub) CHKERRQ(MatGetLocalToGlobalMapping(sub,NULL,&cmap)); 1692 if (cmap) { 1693 CHKERRQ(MatGetBlockSizes(sub,NULL,&bs)); 1694 CHKERRQ(ISLocalToGlobalMappingGetSize(cmap,&nlocal)); 1695 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc)); 1696 CHKERRQ(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 CHKERRQ(MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap)); 1709 CHKERRQ(MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap)); 1710 if (rmap && cmap) CHKERRQ(MatSetLocalToGlobalMapping(A,rmap,cmap)); 1711 CHKERRQ(ISLocalToGlobalMappingDestroy(&rmap)); 1712 CHKERRQ(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 CHKERRQ(MatGetSize(B,&M,&N)); 1722 CHKERRQ(MatGetLocalSize(B,&m,&n)); 1723 CHKERRQ(ISGetSize(vs->isglobal.row[i],&Mi)); 1724 CHKERRQ(ISGetSize(vs->isglobal.col[j],&Ni)); 1725 CHKERRQ(ISGetLocalSize(vs->isglobal.row[i],&mi)); 1726 CHKERRQ(ISGetLocalSize(vs->isglobal.col[j],&ni)); 1727 PetscCheckFalse(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 PetscCheckFalse(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 CHKERRQ(MatCreate(comm,&A)); 1772 CHKERRQ(MatSetType(A,MATNEST)); 1773 A->preallocated = PETSC_TRUE; 1774 CHKERRQ(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 CHKERRQ(MatGetSize(A,&nr,&nc)); 1792 if (reuse == MAT_REUSE_MATRIX) { 1793 PetscInt rnr; 1794 1795 CHKERRQ(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 PetscCheckFalse(rnr != nr,PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows"); 1798 CHKERRQ(MatSeqAIJGetArray(*newmat,&vv)); 1799 } 1800 /* extract CSR for nested SeqAIJ matrices */ 1801 nnz = 0; 1802 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans)); 1812 if (istrans) { 1813 Mat Bt; 1814 1815 CHKERRQ(MatTransposeGetMat(B,&Bt)); 1816 CHKERRQ(MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j])); 1817 B = trans[i*nest->nc+j]; 1818 } 1819 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscMalloc1(nr+1,&ii)); 1832 CHKERRQ(PetscMalloc1(nnz,&jj)); 1833 CHKERRQ(PetscMalloc1(nnz,&vv)); 1834 } else { 1835 PetscCheckFalse(nnz != ii[nr],PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros"); 1836 } 1837 1838 /* new row pointer */ 1839 CHKERRQ(PetscArrayzero(ii,nr+1)); 1840 for (i=0; i<nest->nr; ++i) { 1841 PetscInt ncr,rst; 1842 1843 CHKERRQ(ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL)); 1844 CHKERRQ(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 CHKERRQ(PetscCalloc1(nr,&ci)); 1861 for (i=0; i<nest->nr; ++i) { 1862 PetscInt ncr,rst; 1863 1864 CHKERRQ(ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL)); 1865 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSeqAIJRestoreArray(B,&avv[k])); 1902 CHKERRQ(MatDestroy(&trans[k])); 1903 } 1904 } 1905 } 1906 CHKERRQ(PetscFree4(aii,ajj,avv,trans)); 1907 1908 /* finalize newmat */ 1909 if (reuse == MAT_INITIAL_MATRIX) { 1910 CHKERRQ(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat)); 1911 } else if (reuse == MAT_INPLACE_MATRIX) { 1912 Mat B; 1913 1914 CHKERRQ(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B)); 1915 CHKERRQ(MatHeaderReplace(A,&B)); 1916 } 1917 CHKERRQ(MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY)); 1918 CHKERRQ(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 CHKERRQ(ISAllGather(nest->isglobal.col[j], &bNis)); 1941 CHKERRQ(ISGetSize(bNis,&bN)); 1942 CHKERRQ(ISGetIndices(bNis,&bNindices)); 1943 for (i=0; i<nest->nr; ++i) { 1944 Mat B,D=NULL; 1945 PetscInt bm, br; 1946 const PetscInt *bmindices; 1947 B = nest->m[i][j]; 1948 if (!B) continue; 1949 CHKERRQ(PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&flg)); 1950 if (flg) { 1951 CHKERRQ(PetscTryMethod(B,"MatTransposeGetMat_C",(Mat,Mat*),(B,&D))); 1952 CHKERRQ(PetscTryMethod(B,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(B,&D))); 1953 CHKERRQ(MatConvert(B,((PetscObject)D)->type_name,MAT_INITIAL_MATRIX,&D)); 1954 B = D; 1955 } 1956 CHKERRQ(PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQSBAIJ,MATMPISBAIJ,"")); 1957 if (flg) { 1958 if (D) { 1959 CHKERRQ(MatConvert(D,MATBAIJ,MAT_INPLACE_MATRIX,&D)); 1960 } else { 1961 CHKERRQ(MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&D)); 1962 } 1963 B = D; 1964 } 1965 CHKERRQ(ISGetLocalSize(nest->isglobal.row[i],&bm)); 1966 CHKERRQ(ISGetIndices(nest->isglobal.row[i],&bmindices)); 1967 CHKERRQ(MatGetOwnershipRange(B,&rstart,NULL)); 1968 for (br = 0; br < bm; ++br) { 1969 PetscInt row = bmindices[br], brncols, *cols; 1970 const PetscInt *brcols; 1971 const PetscScalar *brcoldata; 1972 PetscScalar *vals = NULL; 1973 CHKERRQ(MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata)); 1974 CHKERRQ(PetscMalloc1(brncols,&cols)); 1975 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1976 /* 1977 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1978 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1979 */ 1980 if (a != 1.0) { 1981 CHKERRQ(PetscMalloc1(brncols,&vals)); 1982 for (k=0; k<brncols; k++) vals[k] = a * brcoldata[k]; 1983 CHKERRQ(MatSetValues(Y,1,&row,brncols,cols,vals,ADD_VALUES)); 1984 CHKERRQ(PetscFree(vals)); 1985 } else { 1986 CHKERRQ(MatSetValues(Y,1,&row,brncols,cols,brcoldata,ADD_VALUES)); 1987 } 1988 CHKERRQ(MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata)); 1989 CHKERRQ(PetscFree(cols)); 1990 } 1991 if (D) { 1992 CHKERRQ(MatDestroy(&D)); 1993 } 1994 CHKERRQ(ISRestoreIndices(nest->isglobal.row[i],&bmindices)); 1995 } 1996 CHKERRQ(ISRestoreIndices(bNis,&bNindices)); 1997 CHKERRQ(ISDestroy(&bNis)); 1998 } 1999 CHKERRQ(MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY)); 2000 CHKERRQ(MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY)); 2001 PetscFunctionReturn(0); 2002 } 2003 2004 PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 2005 { 2006 Mat_Nest *nest = (Mat_Nest*)A->data; 2007 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart,cstart,cend; 2008 PetscMPIInt size; 2009 Mat C; 2010 2011 PetscFunctionBegin; 2012 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A),&size)); 2013 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 2014 PetscInt nf; 2015 PetscBool fast; 2016 2017 CHKERRQ(PetscStrcmp(newtype,MATAIJ,&fast)); 2018 if (!fast) { 2019 CHKERRQ(PetscStrcmp(newtype,MATSEQAIJ,&fast)); 2020 } 2021 for (i=0; i<nest->nr && fast; ++i) { 2022 for (j=0; j<nest->nc && fast; ++j) { 2023 Mat B = nest->m[i][j]; 2024 if (B) { 2025 CHKERRQ(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast)); 2026 if (!fast) { 2027 PetscBool istrans; 2028 2029 CHKERRQ(PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans)); 2030 if (istrans) { 2031 Mat Bt; 2032 2033 CHKERRQ(MatTransposeGetMat(B,&Bt)); 2034 CHKERRQ(PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast)); 2035 } 2036 } 2037 } 2038 } 2039 } 2040 for (i=0, nf=0; i<nest->nr && fast; ++i) { 2041 CHKERRQ(PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast)); 2042 if (fast) { 2043 PetscInt f,s; 2044 2045 CHKERRQ(ISStrideGetInfo(nest->isglobal.row[i],&f,&s)); 2046 if (f != nf || s != 1) { fast = PETSC_FALSE; } 2047 else { 2048 CHKERRQ(ISGetSize(nest->isglobal.row[i],&f)); 2049 nf += f; 2050 } 2051 } 2052 } 2053 for (i=0, nf=0; i<nest->nc && fast; ++i) { 2054 CHKERRQ(PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast)); 2055 if (fast) { 2056 PetscInt f,s; 2057 2058 CHKERRQ(ISStrideGetInfo(nest->isglobal.col[i],&f,&s)); 2059 if (f != nf || s != 1) { fast = PETSC_FALSE; } 2060 else { 2061 CHKERRQ(ISGetSize(nest->isglobal.col[i],&f)); 2062 nf += f; 2063 } 2064 } 2065 } 2066 if (fast) { 2067 CHKERRQ(MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat)); 2068 PetscFunctionReturn(0); 2069 } 2070 } 2071 CHKERRQ(MatGetSize(A,&M,&N)); 2072 CHKERRQ(MatGetLocalSize(A,&m,&n)); 2073 CHKERRQ(MatGetOwnershipRangeColumn(A,&cstart,&cend)); 2074 if (reuse == MAT_REUSE_MATRIX) C = *newmat; 2075 else { 2076 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2077 CHKERRQ(MatSetType(C,newtype)); 2078 CHKERRQ(MatSetSizes(C,m,n,M,N)); 2079 } 2080 CHKERRQ(PetscMalloc1(2*m,&dnnz)); 2081 onnz = dnnz + m; 2082 for (k=0; k<m; k++) { 2083 dnnz[k] = 0; 2084 onnz[k] = 0; 2085 } 2086 for (j=0; j<nest->nc; ++j) { 2087 IS bNis; 2088 PetscInt bN; 2089 const PetscInt *bNindices; 2090 /* Using global column indices and ISAllGather() is not scalable. */ 2091 CHKERRQ(ISAllGather(nest->isglobal.col[j], &bNis)); 2092 CHKERRQ(ISGetSize(bNis, &bN)); 2093 CHKERRQ(ISGetIndices(bNis,&bNindices)); 2094 for (i=0; i<nest->nr; ++i) { 2095 PetscSF bmsf; 2096 PetscSFNode *iremote; 2097 Mat B; 2098 PetscInt bm, *sub_dnnz,*sub_onnz, br; 2099 const PetscInt *bmindices; 2100 B = nest->m[i][j]; 2101 if (!B) continue; 2102 CHKERRQ(ISGetLocalSize(nest->isglobal.row[i],&bm)); 2103 CHKERRQ(ISGetIndices(nest->isglobal.row[i],&bmindices)); 2104 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf)); 2105 CHKERRQ(PetscMalloc1(bm,&iremote)); 2106 CHKERRQ(PetscMalloc1(bm,&sub_dnnz)); 2107 CHKERRQ(PetscMalloc1(bm,&sub_onnz)); 2108 for (k = 0; k < bm; ++k) { 2109 sub_dnnz[k] = 0; 2110 sub_onnz[k] = 0; 2111 } 2112 /* 2113 Locate the owners for all of the locally-owned global row indices for this row block. 2114 These determine the roots of PetscSF used to communicate preallocation data to row owners. 2115 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 2116 */ 2117 CHKERRQ(MatGetOwnershipRange(B,&rstart,NULL)); 2118 for (br = 0; br < bm; ++br) { 2119 PetscInt row = bmindices[br], brncols, col; 2120 const PetscInt *brcols; 2121 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 2122 PetscMPIInt rowowner = 0; 2123 CHKERRQ(PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel)); 2124 /* how many roots */ 2125 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 2126 /* get nonzero pattern */ 2127 CHKERRQ(MatGetRow(B,br+rstart,&brncols,&brcols,NULL)); 2128 for (k=0; k<brncols; k++) { 2129 col = bNindices[brcols[k]]; 2130 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 2131 sub_dnnz[br]++; 2132 } else { 2133 sub_onnz[br]++; 2134 } 2135 } 2136 CHKERRQ(MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL)); 2137 } 2138 CHKERRQ(ISRestoreIndices(nest->isglobal.row[i],&bmindices)); 2139 /* bsf will have to take care of disposing of bedges. */ 2140 CHKERRQ(PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER)); 2141 CHKERRQ(PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM)); 2142 CHKERRQ(PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM)); 2143 CHKERRQ(PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM)); 2144 CHKERRQ(PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM)); 2145 CHKERRQ(PetscFree(sub_dnnz)); 2146 CHKERRQ(PetscFree(sub_onnz)); 2147 CHKERRQ(PetscSFDestroy(&bmsf)); 2148 } 2149 CHKERRQ(ISRestoreIndices(bNis,&bNindices)); 2150 CHKERRQ(ISDestroy(&bNis)); 2151 } 2152 /* Resize preallocation if overestimated */ 2153 for (i=0;i<m;i++) { 2154 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 2155 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 2156 } 2157 CHKERRQ(MatSeqAIJSetPreallocation(C,0,dnnz)); 2158 CHKERRQ(MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz)); 2159 CHKERRQ(PetscFree(dnnz)); 2160 CHKERRQ(MatAXPY_Dense_Nest(C,1.0,A)); 2161 if (reuse == MAT_INPLACE_MATRIX) { 2162 CHKERRQ(MatHeaderReplace(A,&C)); 2163 } else *newmat = C; 2164 PetscFunctionReturn(0); 2165 } 2166 2167 PetscErrorCode MatConvert_Nest_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 2168 { 2169 Mat B; 2170 PetscInt m,n,M,N; 2171 2172 PetscFunctionBegin; 2173 CHKERRQ(MatGetSize(A,&M,&N)); 2174 CHKERRQ(MatGetLocalSize(A,&m,&n)); 2175 if (reuse == MAT_REUSE_MATRIX) { 2176 B = *newmat; 2177 CHKERRQ(MatZeroEntries(B)); 2178 } else { 2179 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,&B)); 2180 } 2181 CHKERRQ(MatAXPY_Dense_Nest(B,1.0,A)); 2182 if (reuse == MAT_INPLACE_MATRIX) { 2183 CHKERRQ(MatHeaderReplace(A,&B)); 2184 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 2185 PetscFunctionReturn(0); 2186 } 2187 2188 PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has) 2189 { 2190 Mat_Nest *bA = (Mat_Nest*)mat->data; 2191 MatOperation opAdd; 2192 PetscInt i,j,nr = bA->nr,nc = bA->nc; 2193 PetscBool flg; 2194 PetscFunctionBegin; 2195 2196 *has = PETSC_FALSE; 2197 if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 2198 opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD); 2199 for (j=0; j<nc; j++) { 2200 for (i=0; i<nr; i++) { 2201 if (!bA->m[i][j]) continue; 2202 CHKERRQ(MatHasOperation(bA->m[i][j],opAdd,&flg)); 2203 if (!flg) PetscFunctionReturn(0); 2204 } 2205 } 2206 } 2207 if (((void**)mat->ops)[op]) *has = PETSC_TRUE; 2208 PetscFunctionReturn(0); 2209 } 2210 2211 /*MC 2212 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 2213 2214 Level: intermediate 2215 2216 Notes: 2217 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 2218 It allows the use of symmetric and block formats for parts of multi-physics simulations. 2219 It is usually used with DMComposite and DMCreateMatrix() 2220 2221 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 2222 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 2223 than the nest matrix. 2224 2225 .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(), 2226 VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(), 2227 MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats() 2228 M*/ 2229 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 2230 { 2231 Mat_Nest *s; 2232 2233 PetscFunctionBegin; 2234 CHKERRQ(PetscNewLog(A,&s)); 2235 A->data = (void*)s; 2236 2237 s->nr = -1; 2238 s->nc = -1; 2239 s->m = NULL; 2240 s->splitassembly = PETSC_FALSE; 2241 2242 CHKERRQ(PetscMemzero(A->ops,sizeof(*A->ops))); 2243 2244 A->ops->mult = MatMult_Nest; 2245 A->ops->multadd = MatMultAdd_Nest; 2246 A->ops->multtranspose = MatMultTranspose_Nest; 2247 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 2248 A->ops->transpose = MatTranspose_Nest; 2249 A->ops->assemblybegin = MatAssemblyBegin_Nest; 2250 A->ops->assemblyend = MatAssemblyEnd_Nest; 2251 A->ops->zeroentries = MatZeroEntries_Nest; 2252 A->ops->copy = MatCopy_Nest; 2253 A->ops->axpy = MatAXPY_Nest; 2254 A->ops->duplicate = MatDuplicate_Nest; 2255 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 2256 A->ops->destroy = MatDestroy_Nest; 2257 A->ops->view = MatView_Nest; 2258 A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 2259 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 2260 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 2261 A->ops->getdiagonal = MatGetDiagonal_Nest; 2262 A->ops->diagonalscale = MatDiagonalScale_Nest; 2263 A->ops->scale = MatScale_Nest; 2264 A->ops->shift = MatShift_Nest; 2265 A->ops->diagonalset = MatDiagonalSet_Nest; 2266 A->ops->setrandom = MatSetRandom_Nest; 2267 A->ops->hasoperation = MatHasOperation_Nest; 2268 A->ops->missingdiagonal = MatMissingDiagonal_Nest; 2269 2270 A->spptr = NULL; 2271 A->assembled = PETSC_FALSE; 2272 2273 /* expose Nest api's */ 2274 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest)); 2275 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest)); 2276 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest)); 2277 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest)); 2278 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest)); 2279 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest)); 2280 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest)); 2281 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest)); 2282 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ)); 2283 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ)); 2284 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ)); 2285 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS)); 2286 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",MatConvert_Nest_Dense)); 2287 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",MatConvert_Nest_Dense)); 2288 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense)); 2289 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense)); 2290 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense)); 2291 2292 CHKERRQ(PetscObjectChangeTypeName((PetscObject)A,MATNEST)); 2293 PetscFunctionReturn(0); 2294 } 2295