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