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. Use DIFFERENT_NONZERO_PATTERN",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 if (bY->m[i][j]) { 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 *ltog = NULL; 1414 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1415 if (islocal[i]) { 1416 ierr = ISGetLocalSize(islocal[i],&mi);CHKERRQ(ierr); 1417 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1418 } else { 1419 ierr = ISGetLocalSize(isglobal[i],&mi);CHKERRQ(ierr); 1420 } 1421 m += mi; 1422 } 1423 if (!flg) PetscFunctionReturn(0); 1424 1425 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1426 for (i=0,m=0; i<n; i++) { 1427 ISLocalToGlobalMapping smap = NULL; 1428 Mat sub = NULL; 1429 PetscSF sf; 1430 PetscLayout map; 1431 const PetscInt *ix2; 1432 1433 if (!colflg) { 1434 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1435 } else { 1436 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1437 } 1438 if (sub) { 1439 if (!colflg) { 1440 ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr); 1441 } else { 1442 ierr = MatGetLocalToGlobalMapping(sub,NULL,&smap);CHKERRQ(ierr); 1443 } 1444 } 1445 /* 1446 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1447 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1448 */ 1449 ierr = ISGetIndices(isglobal[i],&ix2);CHKERRQ(ierr); 1450 if (islocal[i]) { 1451 PetscInt *ilocal,*iremote; 1452 PetscInt mil,nleaves; 1453 1454 ierr = ISGetLocalSize(islocal[i],&mi);CHKERRQ(ierr); 1455 if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map"); 1456 for (j=0; j<mi; j++) ix[m+j] = j; 1457 ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr); 1458 1459 /* PetscSFSetGraphLayout does not like negative indices */ 1460 ierr = PetscMalloc2(mi,&ilocal,mi,&iremote);CHKERRQ(ierr); 1461 for (j=0, nleaves = 0; j<mi; j++) { 1462 if (ix[m+j] < 0) continue; 1463 ilocal[nleaves] = j; 1464 iremote[nleaves] = ix[m+j]; 1465 nleaves++; 1466 } 1467 ierr = ISGetLocalSize(isglobal[i],&mil);CHKERRQ(ierr); 1468 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);CHKERRQ(ierr); 1469 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);CHKERRQ(ierr); 1470 ierr = PetscLayoutSetLocalSize(map,mil);CHKERRQ(ierr); 1471 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 1472 ierr = PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);CHKERRQ(ierr); 1473 ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr); 1474 ierr = PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr); 1475 ierr = PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr); 1476 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1477 ierr = PetscFree2(ilocal,iremote);CHKERRQ(ierr); 1478 } else { 1479 ierr = ISGetLocalSize(isglobal[i],&mi);CHKERRQ(ierr); 1480 for (j=0; j<mi; j++) ix[m+j] = ix2[i]; 1481 } 1482 ierr = ISRestoreIndices(isglobal[i],&ix2);CHKERRQ(ierr); 1483 m += mi; 1484 } 1485 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1486 PetscFunctionReturn(0); 1487 } 1488 1489 1490 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1491 /* 1492 nprocessors = NP 1493 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1494 proc 0: => (g_0,h_0,) 1495 proc 1: => (g_1,h_1,) 1496 ... 1497 proc nprocs-1: => (g_NP-1,h_NP-1,) 1498 1499 proc 0: proc 1: proc nprocs-1: 1500 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1501 1502 proc 0: 1503 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1504 proc 1: 1505 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1506 1507 proc NP-1: 1508 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1509 */ 1510 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1511 { 1512 Mat_Nest *vs = (Mat_Nest*)A->data; 1513 PetscInt i,j,offset,n,nsum,bs; 1514 PetscErrorCode ierr; 1515 Mat sub = NULL; 1516 1517 PetscFunctionBegin; 1518 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr); 1519 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr); 1520 if (is_row) { /* valid IS is passed in */ 1521 /* refs on is[] are incremeneted */ 1522 for (i=0; i<vs->nr; i++) { 1523 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1524 1525 vs->isglobal.row[i] = is_row[i]; 1526 } 1527 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1528 nsum = 0; 1529 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1530 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1531 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1532 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1533 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1534 nsum += n; 1535 } 1536 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1537 offset -= nsum; 1538 for (i=0; i<vs->nr; i++) { 1539 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1540 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1541 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1542 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1543 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1544 offset += n; 1545 } 1546 } 1547 1548 if (is_col) { /* valid IS is passed in */ 1549 /* refs on is[] are incremeneted */ 1550 for (j=0; j<vs->nc; j++) { 1551 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1552 1553 vs->isglobal.col[j] = is_col[j]; 1554 } 1555 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1556 offset = A->cmap->rstart; 1557 nsum = 0; 1558 for (j=0; j<vs->nc; j++) { 1559 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1560 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1561 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1562 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1563 nsum += n; 1564 } 1565 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1566 offset -= nsum; 1567 for (j=0; j<vs->nc; j++) { 1568 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1569 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1570 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1571 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1572 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1573 offset += n; 1574 } 1575 } 1576 1577 /* Set up the local ISs */ 1578 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1579 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1580 for (i=0,offset=0; i<vs->nr; i++) { 1581 IS isloc; 1582 ISLocalToGlobalMapping rmap = NULL; 1583 PetscInt nlocal,bs; 1584 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1585 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1586 if (rmap) { 1587 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1588 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1589 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1590 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1591 } else { 1592 nlocal = 0; 1593 isloc = NULL; 1594 } 1595 vs->islocal.row[i] = isloc; 1596 offset += nlocal; 1597 } 1598 for (i=0,offset=0; i<vs->nc; i++) { 1599 IS isloc; 1600 ISLocalToGlobalMapping cmap = NULL; 1601 PetscInt nlocal,bs; 1602 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1603 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1604 if (cmap) { 1605 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1606 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1607 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1608 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1609 } else { 1610 nlocal = 0; 1611 isloc = NULL; 1612 } 1613 vs->islocal.col[i] = isloc; 1614 offset += nlocal; 1615 } 1616 1617 /* Set up the aggregate ISLocalToGlobalMapping */ 1618 { 1619 ISLocalToGlobalMapping rmap,cmap; 1620 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr); 1621 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr); 1622 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1623 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1624 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1625 } 1626 1627 #if defined(PETSC_USE_DEBUG) 1628 for (i=0; i<vs->nr; i++) { 1629 for (j=0; j<vs->nc; j++) { 1630 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1631 Mat B = vs->m[i][j]; 1632 if (!B) continue; 1633 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1634 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1635 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1636 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1637 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1638 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1639 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); 1640 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); 1641 } 1642 } 1643 #endif 1644 1645 /* Set A->assembled if all non-null blocks are currently assembled */ 1646 for (i=0; i<vs->nr; i++) { 1647 for (j=0; j<vs->nc; j++) { 1648 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1649 } 1650 } 1651 A->assembled = PETSC_TRUE; 1652 PetscFunctionReturn(0); 1653 } 1654 1655 /*@C 1656 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1657 1658 Collective on Mat 1659 1660 Input Parameter: 1661 + comm - Communicator for the new Mat 1662 . nr - number of nested row blocks 1663 . is_row - index sets for each nested row block, or NULL to make contiguous 1664 . nc - number of nested column blocks 1665 . is_col - index sets for each nested column block, or NULL to make contiguous 1666 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1667 1668 Output Parameter: 1669 . B - new matrix 1670 1671 Level: advanced 1672 1673 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(), 1674 MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(), 1675 MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats() 1676 @*/ 1677 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1678 { 1679 Mat A; 1680 PetscErrorCode ierr; 1681 1682 PetscFunctionBegin; 1683 *B = 0; 1684 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1685 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1686 A->preallocated = PETSC_TRUE; 1687 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1688 *B = A; 1689 PetscFunctionReturn(0); 1690 } 1691 1692 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1693 { 1694 Mat_Nest *nest = (Mat_Nest*)A->data; 1695 Mat *trans; 1696 PetscScalar **avv; 1697 PetscScalar *vv; 1698 PetscInt **aii,**ajj; 1699 PetscInt *ii,*jj,*ci; 1700 PetscInt nr,nc,nnz,i,j; 1701 PetscBool done; 1702 PetscErrorCode ierr; 1703 1704 PetscFunctionBegin; 1705 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 1706 if (reuse == MAT_REUSE_MATRIX) { 1707 PetscInt rnr; 1708 1709 ierr = MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);CHKERRQ(ierr); 1710 if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ"); 1711 if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows"); 1712 ierr = MatSeqAIJGetArray(*newmat,&vv);CHKERRQ(ierr); 1713 } 1714 /* extract CSR for nested SeqAIJ matrices */ 1715 nnz = 0; 1716 ierr = PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);CHKERRQ(ierr); 1717 for (i=0; i<nest->nr; ++i) { 1718 for (j=0; j<nest->nc; ++j) { 1719 Mat B = nest->m[i][j]; 1720 if (B) { 1721 PetscScalar *naa; 1722 PetscInt *nii,*njj,nnr; 1723 PetscBool istrans; 1724 1725 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1726 if (istrans) { 1727 Mat Bt; 1728 1729 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1730 ierr = MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);CHKERRQ(ierr); 1731 B = trans[i*nest->nc+j]; 1732 } 1733 ierr = MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);CHKERRQ(ierr); 1734 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ"); 1735 ierr = MatSeqAIJGetArray(B,&naa);CHKERRQ(ierr); 1736 nnz += nii[nnr]; 1737 1738 aii[i*nest->nc+j] = nii; 1739 ajj[i*nest->nc+j] = njj; 1740 avv[i*nest->nc+j] = naa; 1741 } 1742 } 1743 } 1744 if (reuse != MAT_REUSE_MATRIX) { 1745 ierr = PetscMalloc1(nr+1,&ii);CHKERRQ(ierr); 1746 ierr = PetscMalloc1(nnz,&jj);CHKERRQ(ierr); 1747 ierr = PetscMalloc1(nnz,&vv);CHKERRQ(ierr); 1748 } else { 1749 if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros"); 1750 } 1751 1752 /* new row pointer */ 1753 ierr = PetscArrayzero(ii,nr+1);CHKERRQ(ierr); 1754 for (i=0; i<nest->nr; ++i) { 1755 PetscInt ncr,rst; 1756 1757 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1758 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1759 for (j=0; j<nest->nc; ++j) { 1760 if (aii[i*nest->nc+j]) { 1761 PetscInt *nii = aii[i*nest->nc+j]; 1762 PetscInt ir; 1763 1764 for (ir=rst; ir<ncr+rst; ++ir) { 1765 ii[ir+1] += nii[1]-nii[0]; 1766 nii++; 1767 } 1768 } 1769 } 1770 } 1771 for (i=0; i<nr; i++) ii[i+1] += ii[i]; 1772 1773 /* construct CSR for the new matrix */ 1774 ierr = PetscCalloc1(nr,&ci);CHKERRQ(ierr); 1775 for (i=0; i<nest->nr; ++i) { 1776 PetscInt ncr,rst; 1777 1778 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1779 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1780 for (j=0; j<nest->nc; ++j) { 1781 if (aii[i*nest->nc+j]) { 1782 PetscScalar *nvv = avv[i*nest->nc+j]; 1783 PetscInt *nii = aii[i*nest->nc+j]; 1784 PetscInt *njj = ajj[i*nest->nc+j]; 1785 PetscInt ir,cst; 1786 1787 ierr = ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);CHKERRQ(ierr); 1788 for (ir=rst; ir<ncr+rst; ++ir) { 1789 PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir]; 1790 1791 for (ij=0;ij<rsize;ij++) { 1792 jj[ist+ij] = *njj+cst; 1793 vv[ist+ij] = *nvv; 1794 njj++; 1795 nvv++; 1796 } 1797 ci[ir] += rsize; 1798 nii++; 1799 } 1800 } 1801 } 1802 } 1803 ierr = PetscFree(ci);CHKERRQ(ierr); 1804 1805 /* restore info */ 1806 for (i=0; i<nest->nr; ++i) { 1807 for (j=0; j<nest->nc; ++j) { 1808 Mat B = nest->m[i][j]; 1809 if (B) { 1810 PetscInt nnr = 0, k = i*nest->nc+j; 1811 1812 B = (trans[k] ? trans[k] : B); 1813 ierr = MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);CHKERRQ(ierr); 1814 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ"); 1815 ierr = MatSeqAIJRestoreArray(B,&avv[k]);CHKERRQ(ierr); 1816 ierr = MatDestroy(&trans[k]);CHKERRQ(ierr); 1817 } 1818 } 1819 } 1820 ierr = PetscFree4(aii,ajj,avv,trans);CHKERRQ(ierr); 1821 1822 /* finalize newmat */ 1823 if (reuse == MAT_INITIAL_MATRIX) { 1824 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);CHKERRQ(ierr); 1825 } else if (reuse == MAT_INPLACE_MATRIX) { 1826 Mat B; 1827 1828 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);CHKERRQ(ierr); 1829 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1830 } 1831 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1832 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1833 { 1834 Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data); 1835 a->free_a = PETSC_TRUE; 1836 a->free_ij = PETSC_TRUE; 1837 } 1838 PetscFunctionReturn(0); 1839 } 1840 1841 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1842 { 1843 PetscErrorCode ierr; 1844 Mat_Nest *nest = (Mat_Nest*)A->data; 1845 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1846 PetscInt cstart,cend; 1847 PetscMPIInt size; 1848 Mat C; 1849 1850 PetscFunctionBegin; 1851 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1852 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1853 PetscInt nf; 1854 PetscBool fast; 1855 1856 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr); 1857 if (!fast) { 1858 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr); 1859 } 1860 for (i=0; i<nest->nr && fast; ++i) { 1861 for (j=0; j<nest->nc && fast; ++j) { 1862 Mat B = nest->m[i][j]; 1863 if (B) { 1864 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr); 1865 if (!fast) { 1866 PetscBool istrans; 1867 1868 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1869 if (istrans) { 1870 Mat Bt; 1871 1872 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1873 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);CHKERRQ(ierr); 1874 } 1875 } 1876 } 1877 } 1878 } 1879 for (i=0, nf=0; i<nest->nr && fast; ++i) { 1880 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1881 if (fast) { 1882 PetscInt f,s; 1883 1884 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr); 1885 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1886 else { 1887 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr); 1888 nf += f; 1889 } 1890 } 1891 } 1892 for (i=0, nf=0; i<nest->nc && fast; ++i) { 1893 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1894 if (fast) { 1895 PetscInt f,s; 1896 1897 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr); 1898 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1899 else { 1900 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr); 1901 nf += f; 1902 } 1903 } 1904 } 1905 if (fast) { 1906 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 } 1910 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1911 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1912 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 1913 switch (reuse) { 1914 case MAT_INITIAL_MATRIX: 1915 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1916 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1917 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1918 *newmat = C; 1919 break; 1920 case MAT_REUSE_MATRIX: 1921 C = *newmat; 1922 break; 1923 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1924 } 1925 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1926 onnz = dnnz + m; 1927 for (k=0; k<m; k++) { 1928 dnnz[k] = 0; 1929 onnz[k] = 0; 1930 } 1931 for (j=0; j<nest->nc; ++j) { 1932 IS bNis; 1933 PetscInt bN; 1934 const PetscInt *bNindices; 1935 /* Using global column indices and ISAllGather() is not scalable. */ 1936 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1937 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1938 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1939 for (i=0; i<nest->nr; ++i) { 1940 PetscSF bmsf; 1941 PetscSFNode *iremote; 1942 Mat B; 1943 PetscInt bm, *sub_dnnz,*sub_onnz, br; 1944 const PetscInt *bmindices; 1945 B = nest->m[i][j]; 1946 if (!B) continue; 1947 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1948 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1949 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1950 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 1951 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 1952 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 1953 for (k = 0; k < bm; ++k){ 1954 sub_dnnz[k] = 0; 1955 sub_onnz[k] = 0; 1956 } 1957 /* 1958 Locate the owners for all of the locally-owned global row indices for this row block. 1959 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1960 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1961 */ 1962 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1963 for (br = 0; br < bm; ++br) { 1964 PetscInt row = bmindices[br], brncols, col; 1965 const PetscInt *brcols; 1966 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1967 PetscMPIInt rowowner = 0; 1968 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1969 /* how many roots */ 1970 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1971 /* get nonzero pattern */ 1972 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1973 for (k=0; k<brncols; k++) { 1974 col = bNindices[brcols[k]]; 1975 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 1976 sub_dnnz[br]++; 1977 } else { 1978 sub_onnz[br]++; 1979 } 1980 } 1981 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1982 } 1983 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1984 /* bsf will have to take care of disposing of bedges. */ 1985 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1986 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1987 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1988 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1989 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1990 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 1991 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 1992 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1993 } 1994 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1995 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1996 } 1997 /* Resize preallocation if overestimated */ 1998 for (i=0;i<m;i++) { 1999 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 2000 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 2001 } 2002 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 2003 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 2004 ierr = PetscFree(dnnz);CHKERRQ(ierr); 2005 2006 /* Fill by row */ 2007 for (j=0; j<nest->nc; ++j) { 2008 /* Using global column indices and ISAllGather() is not scalable. */ 2009 IS bNis; 2010 PetscInt bN; 2011 const PetscInt *bNindices; 2012 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 2013 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 2014 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 2015 for (i=0; i<nest->nr; ++i) { 2016 Mat B; 2017 PetscInt bm, br; 2018 const PetscInt *bmindices; 2019 B = nest->m[i][j]; 2020 if (!B) continue; 2021 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 2022 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 2023 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 2024 for (br = 0; br < bm; ++br) { 2025 PetscInt row = bmindices[br], brncols, *cols; 2026 const PetscInt *brcols; 2027 const PetscScalar *brcoldata; 2028 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 2029 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 2030 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 2031 /* 2032 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 2033 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 2034 */ 2035 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 2036 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 2037 ierr = PetscFree(cols);CHKERRQ(ierr); 2038 } 2039 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 2040 } 2041 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 2042 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 2043 } 2044 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2045 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2046 PetscFunctionReturn(0); 2047 } 2048 2049 PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has) 2050 { 2051 Mat_Nest *bA = (Mat_Nest*)mat->data; 2052 PetscInt i,j,nr = bA->nr,nc = bA->nc; 2053 PetscBool flg; 2054 PetscErrorCode ierr; 2055 PetscFunctionBegin; 2056 2057 *has = PETSC_FALSE; 2058 if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MAT_MULT) { 2059 for (j=0; j<nc; j++) { 2060 for (i=0; i<nr; i++) { 2061 if (!bA->m[i][j]) continue; 2062 ierr = MatHasOperation(bA->m[i][j],op,&flg);CHKERRQ(ierr); 2063 if (!flg) PetscFunctionReturn(0); 2064 } 2065 } 2066 } 2067 if (((void**)mat->ops)[op] || (op == MATOP_MAT_MULT && flg)) *has = PETSC_TRUE; 2068 PetscFunctionReturn(0); 2069 } 2070 2071 /*MC 2072 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 2073 2074 Level: intermediate 2075 2076 Notes: 2077 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 2078 It allows the use of symmetric and block formats for parts of multi-physics simulations. 2079 It is usually used with DMComposite and DMCreateMatrix() 2080 2081 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 2082 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 2083 than the nest matrix. 2084 2085 .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(), 2086 VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(), 2087 MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats() 2088 M*/ 2089 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 2090 { 2091 Mat_Nest *s; 2092 PetscErrorCode ierr; 2093 2094 PetscFunctionBegin; 2095 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 2096 A->data = (void*)s; 2097 2098 s->nr = -1; 2099 s->nc = -1; 2100 s->m = NULL; 2101 s->splitassembly = PETSC_FALSE; 2102 2103 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 2104 2105 A->ops->mult = MatMult_Nest; 2106 A->ops->multadd = MatMultAdd_Nest; 2107 A->ops->multtranspose = MatMultTranspose_Nest; 2108 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 2109 A->ops->transpose = MatTranspose_Nest; 2110 A->ops->assemblybegin = MatAssemblyBegin_Nest; 2111 A->ops->assemblyend = MatAssemblyEnd_Nest; 2112 A->ops->zeroentries = MatZeroEntries_Nest; 2113 A->ops->copy = MatCopy_Nest; 2114 A->ops->axpy = MatAXPY_Nest; 2115 A->ops->duplicate = MatDuplicate_Nest; 2116 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 2117 A->ops->destroy = MatDestroy_Nest; 2118 A->ops->view = MatView_Nest; 2119 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 2120 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 2121 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 2122 A->ops->getdiagonal = MatGetDiagonal_Nest; 2123 A->ops->diagonalscale = MatDiagonalScale_Nest; 2124 A->ops->scale = MatScale_Nest; 2125 A->ops->shift = MatShift_Nest; 2126 A->ops->diagonalset = MatDiagonalSet_Nest; 2127 A->ops->setrandom = MatSetRandom_Nest; 2128 A->ops->hasoperation = MatHasOperation_Nest; 2129 2130 A->spptr = 0; 2131 A->assembled = PETSC_FALSE; 2132 2133 /* expose Nest api's */ 2134 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 2135 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 2136 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 2137 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 2138 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 2139 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 2140 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 2141 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 2142 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr); 2143 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr); 2144 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr); 2145 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr); 2146 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_mpidense_C",MatMatMult_Nest_Dense);CHKERRQ(ierr); 2147 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_seqdense_C",MatMatMult_Nest_Dense);CHKERRQ(ierr); 2148 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_dense_C", MatMatMult_Nest_Dense);CHKERRQ(ierr); 2149 2150 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 2151 PetscFunctionReturn(0); 2152 } 2153