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