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