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