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