1 #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/ 2 #include <../src/mat/impls/sell/mpi/mpisell.h> /*I "petscmat.h" I*/ 3 #include <petsc/private/vecimpl.h> 4 #include <petsc/private/isimpl.h> 5 #include <petscblaslapack.h> 6 #include <petscsf.h> 7 8 /*MC 9 MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices. 10 11 This matrix type is identical to MATSEQSELL when constructed with a single process communicator, 12 and MATMPISELL otherwise. As a result, for single process communicators, 13 MatSeqSELLSetPreallocation is supported, and similarly MatMPISELLSetPreallocation is supported 14 for communicators controlling multiple processes. It is recommended that you call both of 15 the above preallocation routines for simplicity. 16 17 Options Database Keys: 18 . -mat_type sell - sets the matrix type to "sell" during a call to MatSetFromOptions() 19 20 Level: beginner 21 22 .seealso: `MatCreateSELL()`, `MatCreateSeqSELL()`, `MATSEQSELL`, `MATMPISELL` 23 M*/ 24 25 PetscErrorCode MatDiagonalSet_MPISELL(Mat Y,Vec D,InsertMode is) 26 { 27 Mat_MPISELL *sell=(Mat_MPISELL*)Y->data; 28 29 PetscFunctionBegin; 30 if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) { 31 PetscCall(MatDiagonalSet(sell->A,D,is)); 32 } else { 33 PetscCall(MatDiagonalSet_Default(Y,D,is)); 34 } 35 PetscFunctionReturn(0); 36 } 37 38 /* 39 Local utility routine that creates a mapping from the global column 40 number to the local number in the off-diagonal part of the local 41 storage of the matrix. When PETSC_USE_CTABLE is used this is scalable at 42 a slightly higher hash table cost; without it it is not scalable (each processor 43 has an order N integer array but is fast to acess. 44 */ 45 PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat) 46 { 47 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 48 PetscInt n=sell->B->cmap->n,i; 49 50 PetscFunctionBegin; 51 PetscCheck(sell->garray,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MPISELL Matrix was assembled but is missing garray"); 52 #if defined(PETSC_USE_CTABLE) 53 PetscCall(PetscTableCreate(n,mat->cmap->N+1,&sell->colmap)); 54 for (i=0; i<n; i++) { 55 PetscCall(PetscTableAdd(sell->colmap,sell->garray[i]+1,i+1,INSERT_VALUES)); 56 } 57 #else 58 PetscCall(PetscCalloc1(mat->cmap->N+1,&sell->colmap)); 59 PetscCall(PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N+1)*sizeof(PetscInt))); 60 for (i=0; i<n; i++) sell->colmap[sell->garray[i]] = i+1; 61 #endif 62 PetscFunctionReturn(0); 63 } 64 65 #define MatSetValues_SeqSELL_A_Private(row,col,value,addv,orow,ocol) \ 66 { \ 67 if (col <= lastcol1) low1 = 0; \ 68 else high1 = nrow1; \ 69 lastcol1 = col; \ 70 while (high1-low1 > 5) { \ 71 t = (low1+high1)/2; \ 72 if (*(cp1+8*t) > col) high1 = t; \ 73 else low1 = t; \ 74 } \ 75 for (_i=low1; _i<high1; _i++) { \ 76 if (*(cp1+8*_i) > col) break; \ 77 if (*(cp1+8*_i) == col) { \ 78 if (addv == ADD_VALUES) *(vp1+8*_i) += value; \ 79 else *(vp1+8*_i) = value; \ 80 goto a_noinsert; \ 81 } \ 82 } \ 83 if (value == 0.0 && ignorezeroentries) {low1 = 0; high1 = nrow1;goto a_noinsert;} \ 84 if (nonew == 1) {low1 = 0; high1 = nrow1; goto a_noinsert;} \ 85 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \ 86 MatSeqXSELLReallocateSELL(A,am,1,nrow1,a->sliidx,row/8,row,col,a->colidx,a->val,cp1,vp1,nonew,MatScalar); \ 87 /* shift up all the later entries in this row */ \ 88 for (ii=nrow1-1; ii>=_i; ii--) { \ 89 *(cp1+8*(ii+1)) = *(cp1+8*ii); \ 90 *(vp1+8*(ii+1)) = *(vp1+8*ii); \ 91 } \ 92 *(cp1+8*_i) = col; \ 93 *(vp1+8*_i) = value; \ 94 a->nz++; nrow1++; A->nonzerostate++; \ 95 a_noinsert: ; \ 96 a->rlen[row] = nrow1; \ 97 } 98 99 #define MatSetValues_SeqSELL_B_Private(row,col,value,addv,orow,ocol) \ 100 { \ 101 if (col <= lastcol2) low2 = 0; \ 102 else high2 = nrow2; \ 103 lastcol2 = col; \ 104 while (high2-low2 > 5) { \ 105 t = (low2+high2)/2; \ 106 if (*(cp2+8*t) > col) high2 = t; \ 107 else low2 = t; \ 108 } \ 109 for (_i=low2; _i<high2; _i++) { \ 110 if (*(cp2+8*_i) > col) break; \ 111 if (*(cp2+8*_i) == col) { \ 112 if (addv == ADD_VALUES) *(vp2+8*_i) += value; \ 113 else *(vp2+8*_i) = value; \ 114 goto b_noinsert; \ 115 } \ 116 } \ 117 if (value == 0.0 && ignorezeroentries) {low2 = 0; high2 = nrow2; goto b_noinsert;} \ 118 if (nonew == 1) {low2 = 0; high2 = nrow2; goto b_noinsert;} \ 119 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \ 120 MatSeqXSELLReallocateSELL(B,bm,1,nrow2,b->sliidx,row/8,row,col,b->colidx,b->val,cp2,vp2,nonew,MatScalar); \ 121 /* shift up all the later entries in this row */ \ 122 for (ii=nrow2-1; ii>=_i; ii--) { \ 123 *(cp2+8*(ii+1)) = *(cp2+8*ii); \ 124 *(vp2+8*(ii+1)) = *(vp2+8*ii); \ 125 } \ 126 *(cp2+8*_i) = col; \ 127 *(vp2+8*_i) = value; \ 128 b->nz++; nrow2++; B->nonzerostate++; \ 129 b_noinsert: ; \ 130 b->rlen[row] = nrow2; \ 131 } 132 133 PetscErrorCode MatSetValues_MPISELL(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv) 134 { 135 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 136 PetscScalar value; 137 PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend,shift1,shift2; 138 PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col; 139 PetscBool roworiented=sell->roworiented; 140 141 /* Some Variables required in the macro */ 142 Mat A=sell->A; 143 Mat_SeqSELL *a=(Mat_SeqSELL*)A->data; 144 PetscBool ignorezeroentries=a->ignorezeroentries,found; 145 Mat B=sell->B; 146 Mat_SeqSELL *b=(Mat_SeqSELL*)B->data; 147 PetscInt *cp1,*cp2,ii,_i,nrow1,nrow2,low1,high1,low2,high2,t,lastcol1,lastcol2; 148 MatScalar *vp1,*vp2; 149 150 PetscFunctionBegin; 151 for (i=0; i<m; i++) { 152 if (im[i] < 0) continue; 153 PetscCheck(im[i] < mat->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,im[i],mat->rmap->N-1); 154 if (im[i] >= rstart && im[i] < rend) { 155 row = im[i] - rstart; 156 lastcol1 = -1; 157 shift1 = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 158 cp1 = a->colidx+shift1; 159 vp1 = a->val+shift1; 160 nrow1 = a->rlen[row]; 161 low1 = 0; 162 high1 = nrow1; 163 lastcol2 = -1; 164 shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 165 cp2 = b->colidx+shift2; 166 vp2 = b->val+shift2; 167 nrow2 = b->rlen[row]; 168 low2 = 0; 169 high2 = nrow2; 170 171 for (j=0; j<n; j++) { 172 if (roworiented) value = v[i*n+j]; 173 else value = v[i+j*m]; 174 if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue; 175 if (in[j] >= cstart && in[j] < cend) { 176 col = in[j] - cstart; 177 MatSetValue_SeqSELL_Private(A,row,col,value,addv,im[i],in[j],cp1,vp1,lastcol1,low1,high1); /* set one value */ 178 } else if (in[j] < 0) { 179 continue; 180 } else { 181 PetscCheck(in[j] < mat->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[j],mat->cmap->N-1); 182 if (mat->was_assembled) { 183 if (!sell->colmap) { 184 PetscCall(MatCreateColmap_MPISELL_Private(mat)); 185 } 186 #if defined(PETSC_USE_CTABLE) 187 PetscCall(PetscTableFind(sell->colmap,in[j]+1,&col)); 188 col--; 189 #else 190 col = sell->colmap[in[j]] - 1; 191 #endif 192 if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) { 193 PetscCall(MatDisAssemble_MPISELL(mat)); 194 col = in[j]; 195 /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */ 196 B = sell->B; 197 b = (Mat_SeqSELL*)B->data; 198 shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 199 cp2 = b->colidx+shift2; 200 vp2 = b->val+shift2; 201 nrow2 = b->rlen[row]; 202 low2 = 0; 203 high2 = nrow2; 204 } else { 205 PetscCheck(col >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]); 206 } 207 } else col = in[j]; 208 MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */ 209 } 210 } 211 } else { 212 PetscCheck(!mat->nooffprocentries,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set",im[i]); 213 if (!sell->donotstash) { 214 mat->assembled = PETSC_FALSE; 215 if (roworiented) { 216 PetscCall(MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 217 } else { 218 PetscCall(MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 219 } 220 } 221 } 222 } 223 PetscFunctionReturn(0); 224 } 225 226 PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 227 { 228 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 229 PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend; 230 PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col; 231 232 PetscFunctionBegin; 233 for (i=0; i<m; i++) { 234 if (idxm[i] < 0) continue; /* negative row */ 235 PetscCheck(idxm[i] < mat->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,idxm[i],mat->rmap->N-1); 236 if (idxm[i] >= rstart && idxm[i] < rend) { 237 row = idxm[i] - rstart; 238 for (j=0; j<n; j++) { 239 if (idxn[j] < 0) continue; /* negative column */ 240 PetscCheck(idxn[j] < mat->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,idxn[j],mat->cmap->N-1); 241 if (idxn[j] >= cstart && idxn[j] < cend) { 242 col = idxn[j] - cstart; 243 PetscCall(MatGetValues(sell->A,1,&row,1,&col,v+i*n+j)); 244 } else { 245 if (!sell->colmap) { 246 PetscCall(MatCreateColmap_MPISELL_Private(mat)); 247 } 248 #if defined(PETSC_USE_CTABLE) 249 PetscCall(PetscTableFind(sell->colmap,idxn[j]+1,&col)); 250 col--; 251 #else 252 col = sell->colmap[idxn[j]] - 1; 253 #endif 254 if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0; 255 else { 256 PetscCall(MatGetValues(sell->B,1,&row,1,&col,v+i*n+j)); 257 } 258 } 259 } 260 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported"); 261 } 262 PetscFunctionReturn(0); 263 } 264 265 extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat,Vec,Vec); 266 267 PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode) 268 { 269 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 270 PetscInt nstash,reallocs; 271 272 PetscFunctionBegin; 273 if (sell->donotstash || mat->nooffprocentries) PetscFunctionReturn(0); 274 275 PetscCall(MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range)); 276 PetscCall(MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs)); 277 PetscCall(PetscInfo(sell->A,"Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n",nstash,reallocs)); 278 PetscFunctionReturn(0); 279 } 280 281 PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode) 282 { 283 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 284 PetscMPIInt n; 285 PetscInt i,flg; 286 PetscInt *row,*col; 287 PetscScalar *val; 288 PetscBool other_disassembled; 289 /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */ 290 PetscFunctionBegin; 291 if (!sell->donotstash && !mat->nooffprocentries) { 292 while (1) { 293 PetscCall(MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg)); 294 if (!flg) break; 295 296 for (i=0; i<n; i++) { /* assemble one by one */ 297 PetscCall(MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode)); 298 } 299 } 300 PetscCall(MatStashScatterEnd_Private(&mat->stash)); 301 } 302 PetscCall(MatAssemblyBegin(sell->A,mode)); 303 PetscCall(MatAssemblyEnd(sell->A,mode)); 304 305 /* 306 determine if any processor has disassembled, if so we must 307 also disassemble ourselves, in order that we may reassemble. 308 */ 309 /* 310 if nonzero structure of submatrix B cannot change then we know that 311 no processor disassembled thus we can skip this stuff 312 */ 313 if (!((Mat_SeqSELL*)sell->B->data)->nonew) { 314 PetscCall(MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat))); 315 PetscCheck(!mat->was_assembled || other_disassembled,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDisAssemble not implemented yet"); 316 } 317 if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) { 318 PetscCall(MatSetUpMultiply_MPISELL(mat)); 319 } 320 /* 321 PetscCall(MatSetOption(sell->B,MAT_USE_INODES,PETSC_FALSE)); 322 */ 323 PetscCall(MatAssemblyBegin(sell->B,mode)); 324 PetscCall(MatAssemblyEnd(sell->B,mode)); 325 PetscCall(PetscFree2(sell->rowvalues,sell->rowindices)); 326 sell->rowvalues = NULL; 327 PetscCall(VecDestroy(&sell->diag)); 328 329 /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */ 330 if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) { 331 PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate; 332 PetscCall(MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat))); 333 } 334 PetscFunctionReturn(0); 335 } 336 337 PetscErrorCode MatZeroEntries_MPISELL(Mat A) 338 { 339 Mat_MPISELL *l=(Mat_MPISELL*)A->data; 340 341 PetscFunctionBegin; 342 PetscCall(MatZeroEntries(l->A)); 343 PetscCall(MatZeroEntries(l->B)); 344 PetscFunctionReturn(0); 345 } 346 347 PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy) 348 { 349 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 350 PetscInt nt; 351 352 PetscFunctionBegin; 353 PetscCall(VecGetLocalSize(xx,&nt)); 354 PetscCheck(nt == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%" PetscInt_FMT ") and xx (%" PetscInt_FMT ")",A->cmap->n,nt); 355 PetscCall(VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 356 PetscCall((*a->A->ops->mult)(a->A,xx,yy)); 357 PetscCall(VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 358 PetscCall((*a->B->ops->multadd)(a->B,a->lvec,yy,yy)); 359 PetscFunctionReturn(0); 360 } 361 362 PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx) 363 { 364 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 365 366 PetscFunctionBegin; 367 PetscCall(MatMultDiagonalBlock(a->A,bb,xx)); 368 PetscFunctionReturn(0); 369 } 370 371 PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz) 372 { 373 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 374 375 PetscFunctionBegin; 376 PetscCall(VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 377 PetscCall((*a->A->ops->multadd)(a->A,xx,yy,zz)); 378 PetscCall(VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 379 PetscCall((*a->B->ops->multadd)(a->B,a->lvec,zz,zz)); 380 PetscFunctionReturn(0); 381 } 382 383 PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy) 384 { 385 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 386 387 PetscFunctionBegin; 388 /* do nondiagonal part */ 389 PetscCall((*a->B->ops->multtranspose)(a->B,xx,a->lvec)); 390 /* do local part */ 391 PetscCall((*a->A->ops->multtranspose)(a->A,xx,yy)); 392 /* add partial results together */ 393 PetscCall(VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE)); 394 PetscCall(VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE)); 395 PetscFunctionReturn(0); 396 } 397 398 PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f) 399 { 400 MPI_Comm comm; 401 Mat_MPISELL *Asell=(Mat_MPISELL*)Amat->data,*Bsell; 402 Mat Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs; 403 IS Me,Notme; 404 PetscInt M,N,first,last,*notme,i; 405 PetscMPIInt size; 406 407 PetscFunctionBegin; 408 /* Easy test: symmetric diagonal block */ 409 Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A; 410 PetscCall(MatIsTranspose(Adia,Bdia,tol,f)); 411 if (!*f) PetscFunctionReturn(0); 412 PetscCall(PetscObjectGetComm((PetscObject)Amat,&comm)); 413 PetscCallMPI(MPI_Comm_size(comm,&size)); 414 if (size == 1) PetscFunctionReturn(0); 415 416 /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */ 417 PetscCall(MatGetSize(Amat,&M,&N)); 418 PetscCall(MatGetOwnershipRange(Amat,&first,&last)); 419 PetscCall(PetscMalloc1(N-last+first,¬me)); 420 for (i=0; i<first; i++) notme[i] = i; 421 for (i=last; i<M; i++) notme[i-last+first] = i; 422 PetscCall(ISCreateGeneral(MPI_COMM_SELF,N-last+first,notme,PETSC_COPY_VALUES,&Notme)); 423 PetscCall(ISCreateStride(MPI_COMM_SELF,last-first,first,1,&Me)); 424 PetscCall(MatCreateSubMatrices(Amat,1,&Me,&Notme,MAT_INITIAL_MATRIX,&Aoffs)); 425 Aoff = Aoffs[0]; 426 PetscCall(MatCreateSubMatrices(Bmat,1,&Notme,&Me,MAT_INITIAL_MATRIX,&Boffs)); 427 Boff = Boffs[0]; 428 PetscCall(MatIsTranspose(Aoff,Boff,tol,f)); 429 PetscCall(MatDestroyMatrices(1,&Aoffs)); 430 PetscCall(MatDestroyMatrices(1,&Boffs)); 431 PetscCall(ISDestroy(&Me)); 432 PetscCall(ISDestroy(&Notme)); 433 PetscCall(PetscFree(notme)); 434 PetscFunctionReturn(0); 435 } 436 437 PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz) 438 { 439 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 440 441 PetscFunctionBegin; 442 /* do nondiagonal part */ 443 PetscCall((*a->B->ops->multtranspose)(a->B,xx,a->lvec)); 444 /* do local part */ 445 PetscCall((*a->A->ops->multtransposeadd)(a->A,xx,yy,zz)); 446 /* add partial results together */ 447 PetscCall(VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE)); 448 PetscCall(VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE)); 449 PetscFunctionReturn(0); 450 } 451 452 /* 453 This only works correctly for square matrices where the subblock A->A is the 454 diagonal block 455 */ 456 PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v) 457 { 458 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 459 460 PetscFunctionBegin; 461 PetscCheck(A->rmap->N == A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block"); 462 PetscCheck(A->rmap->rstart == A->cmap->rstart && A->rmap->rend == A->cmap->rend,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"row partition must equal col partition"); 463 PetscCall(MatGetDiagonal(a->A,v)); 464 PetscFunctionReturn(0); 465 } 466 467 PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa) 468 { 469 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 470 471 PetscFunctionBegin; 472 PetscCall(MatScale(a->A,aa)); 473 PetscCall(MatScale(a->B,aa)); 474 PetscFunctionReturn(0); 475 } 476 477 PetscErrorCode MatDestroy_MPISELL(Mat mat) 478 { 479 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 480 481 PetscFunctionBegin; 482 #if defined(PETSC_USE_LOG) 483 PetscLogObjectState((PetscObject)mat,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT,mat->rmap->N,mat->cmap->N); 484 #endif 485 PetscCall(MatStashDestroy_Private(&mat->stash)); 486 PetscCall(VecDestroy(&sell->diag)); 487 PetscCall(MatDestroy(&sell->A)); 488 PetscCall(MatDestroy(&sell->B)); 489 #if defined(PETSC_USE_CTABLE) 490 PetscCall(PetscTableDestroy(&sell->colmap)); 491 #else 492 PetscCall(PetscFree(sell->colmap)); 493 #endif 494 PetscCall(PetscFree(sell->garray)); 495 PetscCall(VecDestroy(&sell->lvec)); 496 PetscCall(VecScatterDestroy(&sell->Mvctx)); 497 PetscCall(PetscFree2(sell->rowvalues,sell->rowindices)); 498 PetscCall(PetscFree(sell->ld)); 499 PetscCall(PetscFree(mat->data)); 500 501 PetscCall(PetscObjectChangeTypeName((PetscObject)mat,NULL)); 502 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL)); 503 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL)); 504 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL)); 505 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL)); 506 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL)); 507 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL)); 508 PetscFunctionReturn(0); 509 } 510 511 #include <petscdraw.h> 512 PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer) 513 { 514 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 515 PetscMPIInt rank=sell->rank,size=sell->size; 516 PetscBool isdraw,iascii,isbinary; 517 PetscViewer sviewer; 518 PetscViewerFormat format; 519 520 PetscFunctionBegin; 521 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 522 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 523 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 524 if (iascii) { 525 PetscCall(PetscViewerGetFormat(viewer,&format)); 526 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 527 MatInfo info; 528 PetscInt *inodes; 529 530 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank)); 531 PetscCall(MatGetInfo(mat,MAT_LOCAL,&info)); 532 PetscCall(MatInodeGetInodeSizes(sell->A,NULL,&inodes,NULL)); 533 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 534 if (!inodes) { 535 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", not using I-node routines\n", 536 rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory)); 537 } else { 538 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", using I-node routines\n", 539 rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory)); 540 } 541 PetscCall(MatGetInfo(sell->A,MAT_LOCAL,&info)); 542 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used)); 543 PetscCall(MatGetInfo(sell->B,MAT_LOCAL,&info)); 544 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used)); 545 PetscCall(PetscViewerFlush(viewer)); 546 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 547 PetscCall(PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n")); 548 PetscCall(VecScatterView(sell->Mvctx,viewer)); 549 PetscFunctionReturn(0); 550 } else if (format == PETSC_VIEWER_ASCII_INFO) { 551 PetscInt inodecount,inodelimit,*inodes; 552 PetscCall(MatInodeGetInodeSizes(sell->A,&inodecount,&inodes,&inodelimit)); 553 if (inodes) { 554 PetscCall(PetscViewerASCIIPrintf(viewer,"using I-node (on process 0) routines: found %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\n",inodecount,inodelimit)); 555 } else { 556 PetscCall(PetscViewerASCIIPrintf(viewer,"not using I-node (on process 0) routines\n")); 557 } 558 PetscFunctionReturn(0); 559 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 560 PetscFunctionReturn(0); 561 } 562 } else if (isbinary) { 563 if (size == 1) { 564 PetscCall(PetscObjectSetName((PetscObject)sell->A,((PetscObject)mat)->name)); 565 PetscCall(MatView(sell->A,viewer)); 566 } else { 567 /* PetscCall(MatView_MPISELL_Binary(mat,viewer)); */ 568 } 569 PetscFunctionReturn(0); 570 } else if (isdraw) { 571 PetscDraw draw; 572 PetscBool isnull; 573 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 574 PetscCall(PetscDrawIsNull(draw,&isnull)); 575 if (isnull) PetscFunctionReturn(0); 576 } 577 578 { 579 /* assemble the entire matrix onto first processor. */ 580 Mat A; 581 Mat_SeqSELL *Aloc; 582 PetscInt M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j; 583 MatScalar *aval; 584 PetscBool isnonzero; 585 586 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat),&A)); 587 if (rank == 0) { 588 PetscCall(MatSetSizes(A,M,N,M,N)); 589 } else { 590 PetscCall(MatSetSizes(A,0,0,M,N)); 591 } 592 /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */ 593 PetscCall(MatSetType(A,MATMPISELL)); 594 PetscCall(MatMPISELLSetPreallocation(A,0,NULL,0,NULL)); 595 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE)); 596 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)A)); 597 598 /* copy over the A part */ 599 Aloc = (Mat_SeqSELL*)sell->A->data; 600 acolidx = Aloc->colidx; aval = Aloc->val; 601 for (i=0; i<Aloc->totalslices; i++) { /* loop over slices */ 602 for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) { 603 isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]); 604 if (isnonzero) { /* check the mask bit */ 605 row = (i<<3)+(j&0x07) + mat->rmap->rstart; /* i<<3 is the starting row of this slice */ 606 col = *acolidx + mat->rmap->rstart; 607 PetscCall(MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES)); 608 } 609 aval++; acolidx++; 610 } 611 } 612 613 /* copy over the B part */ 614 Aloc = (Mat_SeqSELL*)sell->B->data; 615 acolidx = Aloc->colidx; aval = Aloc->val; 616 for (i=0; i<Aloc->totalslices; i++) { 617 for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) { 618 isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]); 619 if (isnonzero) { 620 row = (i<<3)+(j&0x07) + mat->rmap->rstart; 621 col = sell->garray[*acolidx]; 622 PetscCall(MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES)); 623 } 624 aval++; acolidx++; 625 } 626 } 627 628 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 629 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 630 /* 631 Everyone has to call to draw the matrix since the graphics waits are 632 synchronized across all processors that share the PetscDraw object 633 */ 634 PetscCall(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 635 if (rank == 0) { 636 PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISELL*)(A->data))->A,((PetscObject)mat)->name)); 637 PetscCall(MatView_SeqSELL(((Mat_MPISELL*)(A->data))->A,sviewer)); 638 } 639 PetscCall(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 640 PetscCall(PetscViewerFlush(viewer)); 641 PetscCall(MatDestroy(&A)); 642 } 643 PetscFunctionReturn(0); 644 } 645 646 PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer) 647 { 648 PetscBool iascii,isdraw,issocket,isbinary; 649 650 PetscFunctionBegin; 651 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 652 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 653 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 654 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket)); 655 if (iascii || isdraw || isbinary || issocket) { 656 PetscCall(MatView_MPISELL_ASCIIorDraworSocket(mat,viewer)); 657 } 658 PetscFunctionReturn(0); 659 } 660 661 PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[]) 662 { 663 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 664 665 PetscFunctionBegin; 666 PetscCall(MatGetSize(sell->B,NULL,nghosts)); 667 if (ghosts) *ghosts = sell->garray; 668 PetscFunctionReturn(0); 669 } 670 671 PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info) 672 { 673 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 674 Mat A=mat->A,B=mat->B; 675 PetscLogDouble isend[5],irecv[5]; 676 677 PetscFunctionBegin; 678 info->block_size = 1.0; 679 PetscCall(MatGetInfo(A,MAT_LOCAL,info)); 680 681 isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded; 682 isend[3] = info->memory; isend[4] = info->mallocs; 683 684 PetscCall(MatGetInfo(B,MAT_LOCAL,info)); 685 686 isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded; 687 isend[3] += info->memory; isend[4] += info->mallocs; 688 if (flag == MAT_LOCAL) { 689 info->nz_used = isend[0]; 690 info->nz_allocated = isend[1]; 691 info->nz_unneeded = isend[2]; 692 info->memory = isend[3]; 693 info->mallocs = isend[4]; 694 } else if (flag == MAT_GLOBAL_MAX) { 695 PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)matin))); 696 697 info->nz_used = irecv[0]; 698 info->nz_allocated = irecv[1]; 699 info->nz_unneeded = irecv[2]; 700 info->memory = irecv[3]; 701 info->mallocs = irecv[4]; 702 } else if (flag == MAT_GLOBAL_SUM) { 703 PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)matin))); 704 705 info->nz_used = irecv[0]; 706 info->nz_allocated = irecv[1]; 707 info->nz_unneeded = irecv[2]; 708 info->memory = irecv[3]; 709 info->mallocs = irecv[4]; 710 } 711 info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */ 712 info->fill_ratio_needed = 0; 713 info->factor_mallocs = 0; 714 PetscFunctionReturn(0); 715 } 716 717 PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg) 718 { 719 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 720 721 PetscFunctionBegin; 722 switch (op) { 723 case MAT_NEW_NONZERO_LOCATIONS: 724 case MAT_NEW_NONZERO_ALLOCATION_ERR: 725 case MAT_UNUSED_NONZERO_LOCATION_ERR: 726 case MAT_KEEP_NONZERO_PATTERN: 727 case MAT_NEW_NONZERO_LOCATION_ERR: 728 case MAT_USE_INODES: 729 case MAT_IGNORE_ZERO_ENTRIES: 730 MatCheckPreallocated(A,1); 731 PetscCall(MatSetOption(a->A,op,flg)); 732 PetscCall(MatSetOption(a->B,op,flg)); 733 break; 734 case MAT_ROW_ORIENTED: 735 MatCheckPreallocated(A,1); 736 a->roworiented = flg; 737 738 PetscCall(MatSetOption(a->A,op,flg)); 739 PetscCall(MatSetOption(a->B,op,flg)); 740 break; 741 case MAT_FORCE_DIAGONAL_ENTRIES: 742 case MAT_SORTED_FULL: 743 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 744 break; 745 case MAT_IGNORE_OFF_PROC_ENTRIES: 746 a->donotstash = flg; 747 break; 748 case MAT_SPD: 749 A->spd_set = PETSC_TRUE; 750 A->spd = flg; 751 if (flg) { 752 A->symmetric = PETSC_TRUE; 753 A->structurally_symmetric = PETSC_TRUE; 754 A->symmetric_set = PETSC_TRUE; 755 A->structurally_symmetric_set = PETSC_TRUE; 756 } 757 break; 758 case MAT_SYMMETRIC: 759 MatCheckPreallocated(A,1); 760 PetscCall(MatSetOption(a->A,op,flg)); 761 break; 762 case MAT_STRUCTURALLY_SYMMETRIC: 763 MatCheckPreallocated(A,1); 764 PetscCall(MatSetOption(a->A,op,flg)); 765 break; 766 case MAT_HERMITIAN: 767 MatCheckPreallocated(A,1); 768 PetscCall(MatSetOption(a->A,op,flg)); 769 break; 770 case MAT_SYMMETRY_ETERNAL: 771 MatCheckPreallocated(A,1); 772 PetscCall(MatSetOption(a->A,op,flg)); 773 break; 774 default: 775 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 776 } 777 PetscFunctionReturn(0); 778 } 779 780 PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr) 781 { 782 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 783 Mat a=sell->A,b=sell->B; 784 PetscInt s1,s2,s3; 785 786 PetscFunctionBegin; 787 PetscCall(MatGetLocalSize(mat,&s2,&s3)); 788 if (rr) { 789 PetscCall(VecGetLocalSize(rr,&s1)); 790 PetscCheck(s1==s3,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size"); 791 /* Overlap communication with computation. */ 792 PetscCall(VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD)); 793 } 794 if (ll) { 795 PetscCall(VecGetLocalSize(ll,&s1)); 796 PetscCheck(s1==s2,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size"); 797 PetscCall((*b->ops->diagonalscale)(b,ll,NULL)); 798 } 799 /* scale the diagonal block */ 800 PetscCall((*a->ops->diagonalscale)(a,ll,rr)); 801 802 if (rr) { 803 /* Do a scatter end and then right scale the off-diagonal block */ 804 PetscCall(VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD)); 805 PetscCall((*b->ops->diagonalscale)(b,NULL,sell->lvec)); 806 } 807 PetscFunctionReturn(0); 808 } 809 810 PetscErrorCode MatSetUnfactored_MPISELL(Mat A) 811 { 812 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 813 814 PetscFunctionBegin; 815 PetscCall(MatSetUnfactored(a->A)); 816 PetscFunctionReturn(0); 817 } 818 819 PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool *flag) 820 { 821 Mat_MPISELL *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data; 822 Mat a,b,c,d; 823 PetscBool flg; 824 825 PetscFunctionBegin; 826 a = matA->A; b = matA->B; 827 c = matB->A; d = matB->B; 828 829 PetscCall(MatEqual(a,c,&flg)); 830 if (flg) { 831 PetscCall(MatEqual(b,d,&flg)); 832 } 833 PetscCall(MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 834 PetscFunctionReturn(0); 835 } 836 837 PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str) 838 { 839 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 840 Mat_MPISELL *b=(Mat_MPISELL*)B->data; 841 842 PetscFunctionBegin; 843 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 844 if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) { 845 /* because of the column compression in the off-processor part of the matrix a->B, 846 the number of columns in a->B and b->B may be different, hence we cannot call 847 the MatCopy() directly on the two parts. If need be, we can provide a more 848 efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices 849 then copying the submatrices */ 850 PetscCall(MatCopy_Basic(A,B,str)); 851 } else { 852 PetscCall(MatCopy(a->A,b->A,str)); 853 PetscCall(MatCopy(a->B,b->B,str)); 854 } 855 PetscFunctionReturn(0); 856 } 857 858 PetscErrorCode MatSetUp_MPISELL(Mat A) 859 { 860 PetscFunctionBegin; 861 PetscCall(MatMPISELLSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL)); 862 PetscFunctionReturn(0); 863 } 864 865 extern PetscErrorCode MatConjugate_SeqSELL(Mat); 866 867 PetscErrorCode MatConjugate_MPISELL(Mat mat) 868 { 869 PetscFunctionBegin; 870 if (PetscDefined(USE_COMPLEX)) { 871 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 872 873 PetscCall(MatConjugate_SeqSELL(sell->A)); 874 PetscCall(MatConjugate_SeqSELL(sell->B)); 875 } 876 PetscFunctionReturn(0); 877 } 878 879 PetscErrorCode MatRealPart_MPISELL(Mat A) 880 { 881 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 882 883 PetscFunctionBegin; 884 PetscCall(MatRealPart(a->A)); 885 PetscCall(MatRealPart(a->B)); 886 PetscFunctionReturn(0); 887 } 888 889 PetscErrorCode MatImaginaryPart_MPISELL(Mat A) 890 { 891 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 892 893 PetscFunctionBegin; 894 PetscCall(MatImaginaryPart(a->A)); 895 PetscCall(MatImaginaryPart(a->B)); 896 PetscFunctionReturn(0); 897 } 898 899 PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values) 900 { 901 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 902 903 PetscFunctionBegin; 904 PetscCall(MatInvertBlockDiagonal(a->A,values)); 905 A->factorerrortype = a->A->factorerrortype; 906 PetscFunctionReturn(0); 907 } 908 909 static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx) 910 { 911 Mat_MPISELL *sell=(Mat_MPISELL*)x->data; 912 913 PetscFunctionBegin; 914 PetscCall(MatSetRandom(sell->A,rctx)); 915 PetscCall(MatSetRandom(sell->B,rctx)); 916 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 917 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 918 PetscFunctionReturn(0); 919 } 920 921 PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A) 922 { 923 PetscFunctionBegin; 924 PetscOptionsHeadBegin(PetscOptionsObject,"MPISELL options"); 925 PetscOptionsHeadEnd(); 926 PetscFunctionReturn(0); 927 } 928 929 PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a) 930 { 931 Mat_MPISELL *msell=(Mat_MPISELL*)Y->data; 932 Mat_SeqSELL *sell=(Mat_SeqSELL*)msell->A->data; 933 934 PetscFunctionBegin; 935 if (!Y->preallocated) { 936 PetscCall(MatMPISELLSetPreallocation(Y,1,NULL,0,NULL)); 937 } else if (!sell->nz) { 938 PetscInt nonew = sell->nonew; 939 PetscCall(MatSeqSELLSetPreallocation(msell->A,1,NULL)); 940 sell->nonew = nonew; 941 } 942 PetscCall(MatShift_Basic(Y,a)); 943 PetscFunctionReturn(0); 944 } 945 946 PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool *missing,PetscInt *d) 947 { 948 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 949 950 PetscFunctionBegin; 951 PetscCheck(A->rmap->n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices"); 952 PetscCall(MatMissingDiagonal(a->A,missing,d)); 953 if (d) { 954 PetscInt rstart; 955 PetscCall(MatGetOwnershipRange(A,&rstart,NULL)); 956 *d += rstart; 957 958 } 959 PetscFunctionReturn(0); 960 } 961 962 PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a) 963 { 964 PetscFunctionBegin; 965 *a = ((Mat_MPISELL*)A->data)->A; 966 PetscFunctionReturn(0); 967 } 968 969 /* -------------------------------------------------------------------*/ 970 static struct _MatOps MatOps_Values = {MatSetValues_MPISELL, 971 NULL, 972 NULL, 973 MatMult_MPISELL, 974 /* 4*/ MatMultAdd_MPISELL, 975 MatMultTranspose_MPISELL, 976 MatMultTransposeAdd_MPISELL, 977 NULL, 978 NULL, 979 NULL, 980 /*10*/ NULL, 981 NULL, 982 NULL, 983 MatSOR_MPISELL, 984 NULL, 985 /*15*/ MatGetInfo_MPISELL, 986 MatEqual_MPISELL, 987 MatGetDiagonal_MPISELL, 988 MatDiagonalScale_MPISELL, 989 NULL, 990 /*20*/ MatAssemblyBegin_MPISELL, 991 MatAssemblyEnd_MPISELL, 992 MatSetOption_MPISELL, 993 MatZeroEntries_MPISELL, 994 /*24*/ NULL, 995 NULL, 996 NULL, 997 NULL, 998 NULL, 999 /*29*/ MatSetUp_MPISELL, 1000 NULL, 1001 NULL, 1002 MatGetDiagonalBlock_MPISELL, 1003 NULL, 1004 /*34*/ MatDuplicate_MPISELL, 1005 NULL, 1006 NULL, 1007 NULL, 1008 NULL, 1009 /*39*/ NULL, 1010 NULL, 1011 NULL, 1012 MatGetValues_MPISELL, 1013 MatCopy_MPISELL, 1014 /*44*/ NULL, 1015 MatScale_MPISELL, 1016 MatShift_MPISELL, 1017 MatDiagonalSet_MPISELL, 1018 NULL, 1019 /*49*/ MatSetRandom_MPISELL, 1020 NULL, 1021 NULL, 1022 NULL, 1023 NULL, 1024 /*54*/ MatFDColoringCreate_MPIXAIJ, 1025 NULL, 1026 MatSetUnfactored_MPISELL, 1027 NULL, 1028 NULL, 1029 /*59*/ NULL, 1030 MatDestroy_MPISELL, 1031 MatView_MPISELL, 1032 NULL, 1033 NULL, 1034 /*64*/ NULL, 1035 NULL, 1036 NULL, 1037 NULL, 1038 NULL, 1039 /*69*/ NULL, 1040 NULL, 1041 NULL, 1042 NULL, 1043 NULL, 1044 NULL, 1045 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */ 1046 MatSetFromOptions_MPISELL, 1047 NULL, 1048 NULL, 1049 NULL, 1050 /*80*/ NULL, 1051 NULL, 1052 NULL, 1053 /*83*/ NULL, 1054 NULL, 1055 NULL, 1056 NULL, 1057 NULL, 1058 NULL, 1059 /*89*/ NULL, 1060 NULL, 1061 NULL, 1062 NULL, 1063 NULL, 1064 /*94*/ NULL, 1065 NULL, 1066 NULL, 1067 NULL, 1068 NULL, 1069 /*99*/ NULL, 1070 NULL, 1071 NULL, 1072 MatConjugate_MPISELL, 1073 NULL, 1074 /*104*/NULL, 1075 MatRealPart_MPISELL, 1076 MatImaginaryPart_MPISELL, 1077 NULL, 1078 NULL, 1079 /*109*/NULL, 1080 NULL, 1081 NULL, 1082 NULL, 1083 MatMissingDiagonal_MPISELL, 1084 /*114*/NULL, 1085 NULL, 1086 MatGetGhosts_MPISELL, 1087 NULL, 1088 NULL, 1089 /*119*/NULL, 1090 NULL, 1091 NULL, 1092 NULL, 1093 NULL, 1094 /*124*/NULL, 1095 NULL, 1096 MatInvertBlockDiagonal_MPISELL, 1097 NULL, 1098 NULL, 1099 /*129*/NULL, 1100 NULL, 1101 NULL, 1102 NULL, 1103 NULL, 1104 /*134*/NULL, 1105 NULL, 1106 NULL, 1107 NULL, 1108 NULL, 1109 /*139*/NULL, 1110 NULL, 1111 NULL, 1112 MatFDColoringSetUp_MPIXAIJ, 1113 NULL, 1114 /*144*/NULL, 1115 NULL, 1116 NULL, 1117 NULL 1118 }; 1119 1120 /* ----------------------------------------------------------------------------------------*/ 1121 1122 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1123 { 1124 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1125 1126 PetscFunctionBegin; 1127 PetscCall(MatStoreValues(sell->A)); 1128 PetscCall(MatStoreValues(sell->B)); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1133 { 1134 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1135 1136 PetscFunctionBegin; 1137 PetscCall(MatRetrieveValues(sell->A)); 1138 PetscCall(MatRetrieveValues(sell->B)); 1139 PetscFunctionReturn(0); 1140 } 1141 1142 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[]) 1143 { 1144 Mat_MPISELL *b; 1145 1146 PetscFunctionBegin; 1147 PetscCall(PetscLayoutSetUp(B->rmap)); 1148 PetscCall(PetscLayoutSetUp(B->cmap)); 1149 b = (Mat_MPISELL*)B->data; 1150 1151 if (!B->preallocated) { 1152 /* Explicitly create 2 MATSEQSELL matrices. */ 1153 PetscCall(MatCreate(PETSC_COMM_SELF,&b->A)); 1154 PetscCall(MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n)); 1155 PetscCall(MatSetBlockSizesFromMats(b->A,B,B)); 1156 PetscCall(MatSetType(b->A,MATSEQSELL)); 1157 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->A)); 1158 PetscCall(MatCreate(PETSC_COMM_SELF,&b->B)); 1159 PetscCall(MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N)); 1160 PetscCall(MatSetBlockSizesFromMats(b->B,B,B)); 1161 PetscCall(MatSetType(b->B,MATSEQSELL)); 1162 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->B)); 1163 } 1164 1165 PetscCall(MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen)); 1166 PetscCall(MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen)); 1167 B->preallocated = PETSC_TRUE; 1168 B->was_assembled = PETSC_FALSE; 1169 /* 1170 critical for MatAssemblyEnd to work. 1171 MatAssemblyBegin checks it to set up was_assembled 1172 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1173 */ 1174 B->assembled = PETSC_FALSE; 1175 PetscFunctionReturn(0); 1176 } 1177 1178 PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 1179 { 1180 Mat mat; 1181 Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data; 1182 1183 PetscFunctionBegin; 1184 *newmat = NULL; 1185 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin),&mat)); 1186 PetscCall(MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N)); 1187 PetscCall(MatSetBlockSizesFromMats(mat,matin,matin)); 1188 PetscCall(MatSetType(mat,((PetscObject)matin)->type_name)); 1189 a = (Mat_MPISELL*)mat->data; 1190 1191 mat->factortype = matin->factortype; 1192 mat->assembled = PETSC_TRUE; 1193 mat->insertmode = NOT_SET_VALUES; 1194 mat->preallocated = PETSC_TRUE; 1195 1196 a->size = oldmat->size; 1197 a->rank = oldmat->rank; 1198 a->donotstash = oldmat->donotstash; 1199 a->roworiented = oldmat->roworiented; 1200 a->rowindices = NULL; 1201 a->rowvalues = NULL; 1202 a->getrowactive = PETSC_FALSE; 1203 1204 PetscCall(PetscLayoutReference(matin->rmap,&mat->rmap)); 1205 PetscCall(PetscLayoutReference(matin->cmap,&mat->cmap)); 1206 1207 if (oldmat->colmap) { 1208 #if defined(PETSC_USE_CTABLE) 1209 PetscCall(PetscTableCreateCopy(oldmat->colmap,&a->colmap)); 1210 #else 1211 PetscCall(PetscMalloc1(mat->cmap->N,&a->colmap)); 1212 PetscCall(PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt))); 1213 PetscCall(PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N)); 1214 #endif 1215 } else a->colmap = NULL; 1216 if (oldmat->garray) { 1217 PetscInt len; 1218 len = oldmat->B->cmap->n; 1219 PetscCall(PetscMalloc1(len+1,&a->garray)); 1220 PetscCall(PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt))); 1221 if (len) PetscCall(PetscArraycpy(a->garray,oldmat->garray,len)); 1222 } else a->garray = NULL; 1223 1224 PetscCall(VecDuplicate(oldmat->lvec,&a->lvec)); 1225 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec)); 1226 PetscCall(VecScatterCopy(oldmat->Mvctx,&a->Mvctx)); 1227 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx)); 1228 PetscCall(MatDuplicate(oldmat->A,cpvalues,&a->A)); 1229 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A)); 1230 PetscCall(MatDuplicate(oldmat->B,cpvalues,&a->B)); 1231 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B)); 1232 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist)); 1233 *newmat = mat; 1234 PetscFunctionReturn(0); 1235 } 1236 1237 /*@C 1238 MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format. 1239 For good matrix assembly performance the user should preallocate the matrix storage by 1240 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1241 1242 Collective 1243 1244 Input Parameters: 1245 + B - the matrix 1246 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1247 (same value is used for all local rows) 1248 . d_nnz - array containing the number of nonzeros in the various rows of the 1249 DIAGONAL portion of the local submatrix (possibly different for each row) 1250 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1251 The size of this array is equal to the number of local rows, i.e 'm'. 1252 For matrices that will be factored, you must leave room for (and set) 1253 the diagonal entry even if it is zero. 1254 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1255 submatrix (same value is used for all local rows). 1256 - o_nnz - array containing the number of nonzeros in the various rows of the 1257 OFF-DIAGONAL portion of the local submatrix (possibly different for 1258 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1259 structure. The size of this array is equal to the number 1260 of local rows, i.e 'm'. 1261 1262 If the *_nnz parameter is given then the *_nz parameter is ignored 1263 1264 The stored row and column indices begin with zero. 1265 1266 The parallel matrix is partitioned such that the first m0 rows belong to 1267 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1268 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1269 1270 The DIAGONAL portion of the local submatrix of a processor can be defined 1271 as the submatrix which is obtained by extraction the part corresponding to 1272 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1273 first row that belongs to the processor, r2 is the last row belonging to 1274 the this processor, and c1-c2 is range of indices of the local part of a 1275 vector suitable for applying the matrix to. This is an mxn matrix. In the 1276 common case of a square matrix, the row and column ranges are the same and 1277 the DIAGONAL part is also square. The remaining portion of the local 1278 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1279 1280 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1281 1282 You can call MatGetInfo() to get information on how effective the preallocation was; 1283 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1284 You can also run with the option -info and look for messages with the string 1285 malloc in them to see if additional memory allocation was needed. 1286 1287 Example usage: 1288 1289 Consider the following 8x8 matrix with 34 non-zero values, that is 1290 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1291 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1292 as follows: 1293 1294 .vb 1295 1 2 0 | 0 3 0 | 0 4 1296 Proc0 0 5 6 | 7 0 0 | 8 0 1297 9 0 10 | 11 0 0 | 12 0 1298 ------------------------------------- 1299 13 0 14 | 15 16 17 | 0 0 1300 Proc1 0 18 0 | 19 20 21 | 0 0 1301 0 0 0 | 22 23 0 | 24 0 1302 ------------------------------------- 1303 Proc2 25 26 27 | 0 0 28 | 29 0 1304 30 0 0 | 31 32 33 | 0 34 1305 .ve 1306 1307 This can be represented as a collection of submatrices as: 1308 1309 .vb 1310 A B C 1311 D E F 1312 G H I 1313 .ve 1314 1315 Where the submatrices A,B,C are owned by proc0, D,E,F are 1316 owned by proc1, G,H,I are owned by proc2. 1317 1318 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1319 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1320 The 'M','N' parameters are 8,8, and have the same values on all procs. 1321 1322 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1323 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1324 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1325 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1326 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1327 matrix, ans [DF] as another SeqSELL matrix. 1328 1329 When d_nz, o_nz parameters are specified, d_nz storage elements are 1330 allocated for every row of the local diagonal submatrix, and o_nz 1331 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1332 One way to choose d_nz and o_nz is to use the max nonzerors per local 1333 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1334 In this case, the values of d_nz,o_nz are: 1335 .vb 1336 proc0 : dnz = 2, o_nz = 2 1337 proc1 : dnz = 3, o_nz = 2 1338 proc2 : dnz = 1, o_nz = 4 1339 .ve 1340 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1341 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1342 for proc3. i.e we are using 12+15+10=37 storage locations to store 1343 34 values. 1344 1345 When d_nnz, o_nnz parameters are specified, the storage is specified 1346 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1347 In the above case the values for d_nnz,o_nnz are: 1348 .vb 1349 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1350 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1351 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1352 .ve 1353 Here the space allocated is according to nz (or maximum values in the nnz 1354 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1355 1356 Level: intermediate 1357 1358 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreatesell()`, 1359 `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()` 1360 @*/ 1361 PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1362 { 1363 PetscFunctionBegin; 1364 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1365 PetscValidType(B,1); 1366 PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz)); 1367 PetscFunctionReturn(0); 1368 } 1369 1370 /*MC 1371 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1372 based on the sliced Ellpack format 1373 1374 Options Database Keys: 1375 . -mat_type sell - sets the matrix type to "seqsell" during a call to MatSetFromOptions() 1376 1377 Level: beginner 1378 1379 .seealso: `MatCreateSell()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ` 1380 M*/ 1381 1382 /*@C 1383 MatCreateSELL - Creates a sparse parallel matrix in SELL format. 1384 1385 Collective 1386 1387 Input Parameters: 1388 + comm - MPI communicator 1389 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1390 This value should be the same as the local size used in creating the 1391 y vector for the matrix-vector product y = Ax. 1392 . n - This value should be the same as the local size used in creating the 1393 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1394 calculated if N is given) For square matrices n is almost always m. 1395 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 1396 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 1397 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1398 (same value is used for all local rows) 1399 . d_rlen - array containing the number of nonzeros in the various rows of the 1400 DIAGONAL portion of the local submatrix (possibly different for each row) 1401 or NULL, if d_rlenmax is used to specify the nonzero structure. 1402 The size of this array is equal to the number of local rows, i.e 'm'. 1403 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1404 submatrix (same value is used for all local rows). 1405 - o_rlen - array containing the number of nonzeros in the various rows of the 1406 OFF-DIAGONAL portion of the local submatrix (possibly different for 1407 each row) or NULL, if o_rlenmax is used to specify the nonzero 1408 structure. The size of this array is equal to the number 1409 of local rows, i.e 'm'. 1410 1411 Output Parameter: 1412 . A - the matrix 1413 1414 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1415 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1416 [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation] 1417 1418 Notes: 1419 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1420 1421 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1422 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1423 storage requirements for this matrix. 1424 1425 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one 1426 processor than it must be used on all processors that share the object for 1427 that argument. 1428 1429 The user MUST specify either the local or global matrix dimensions 1430 (possibly both). 1431 1432 The parallel matrix is partitioned across processors such that the 1433 first m0 rows belong to process 0, the next m1 rows belong to 1434 process 1, the next m2 rows belong to process 2 etc.. where 1435 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1436 values corresponding to [m x N] submatrix. 1437 1438 The columns are logically partitioned with the n0 columns belonging 1439 to 0th partition, the next n1 columns belonging to the next 1440 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1441 1442 The DIAGONAL portion of the local submatrix on any given processor 1443 is the submatrix corresponding to the rows and columns m,n 1444 corresponding to the given processor. i.e diagonal matrix on 1445 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1446 etc. The remaining portion of the local submatrix [m x (N-n)] 1447 constitute the OFF-DIAGONAL portion. The example below better 1448 illustrates this concept. 1449 1450 For a square global matrix we define each processor's diagonal portion 1451 to be its local rows and the corresponding columns (a square submatrix); 1452 each processor's off-diagonal portion encompasses the remainder of the 1453 local matrix (a rectangular submatrix). 1454 1455 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1456 1457 When calling this routine with a single process communicator, a matrix of 1458 type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this 1459 type of communicator, use the construction mechanism: 1460 MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...); 1461 1462 Options Database Keys: 1463 - -mat_sell_oneindex - Internally use indexing starting at 1 1464 rather than 0. Note that when calling MatSetValues(), 1465 the user still MUST index entries starting at 0! 1466 1467 Example usage: 1468 1469 Consider the following 8x8 matrix with 34 non-zero values, that is 1470 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1471 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1472 as follows: 1473 1474 .vb 1475 1 2 0 | 0 3 0 | 0 4 1476 Proc0 0 5 6 | 7 0 0 | 8 0 1477 9 0 10 | 11 0 0 | 12 0 1478 ------------------------------------- 1479 13 0 14 | 15 16 17 | 0 0 1480 Proc1 0 18 0 | 19 20 21 | 0 0 1481 0 0 0 | 22 23 0 | 24 0 1482 ------------------------------------- 1483 Proc2 25 26 27 | 0 0 28 | 29 0 1484 30 0 0 | 31 32 33 | 0 34 1485 .ve 1486 1487 This can be represented as a collection of submatrices as: 1488 1489 .vb 1490 A B C 1491 D E F 1492 G H I 1493 .ve 1494 1495 Where the submatrices A,B,C are owned by proc0, D,E,F are 1496 owned by proc1, G,H,I are owned by proc2. 1497 1498 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1499 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1500 The 'M','N' parameters are 8,8, and have the same values on all procs. 1501 1502 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1503 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1504 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1505 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1506 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1507 matrix, ans [DF] as another SeqSELL matrix. 1508 1509 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1510 allocated for every row of the local diagonal submatrix, and o_rlenmax 1511 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1512 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1513 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1514 In this case, the values of d_rlenmax,o_rlenmax are: 1515 .vb 1516 proc0 : d_rlenmax = 2, o_rlenmax = 2 1517 proc1 : d_rlenmax = 3, o_rlenmax = 2 1518 proc2 : d_rlenmax = 1, o_rlenmax = 4 1519 .ve 1520 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1521 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1522 for proc3. i.e we are using 12+15+10=37 storage locations to store 1523 34 values. 1524 1525 When d_rlen, o_rlen parameters are specified, the storage is specified 1526 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1527 In the above case the values for d_nnz,o_nnz are: 1528 .vb 1529 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1530 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1531 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1532 .ve 1533 Here the space allocated is still 37 though there are 34 nonzeros because 1534 the allocation is always done according to rlenmax. 1535 1536 Level: intermediate 1537 1538 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MatMPISELLSetPreallocationSELL()`, 1539 `MATMPISELL`, `MatCreateMPISELLWithArrays()` 1540 @*/ 1541 PetscErrorCode MatCreateSELL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[],Mat *A) 1542 { 1543 PetscMPIInt size; 1544 1545 PetscFunctionBegin; 1546 PetscCall(MatCreate(comm,A)); 1547 PetscCall(MatSetSizes(*A,m,n,M,N)); 1548 PetscCallMPI(MPI_Comm_size(comm,&size)); 1549 if (size > 1) { 1550 PetscCall(MatSetType(*A,MATMPISELL)); 1551 PetscCall(MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen)); 1552 } else { 1553 PetscCall(MatSetType(*A,MATSEQSELL)); 1554 PetscCall(MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen)); 1555 } 1556 PetscFunctionReturn(0); 1557 } 1558 1559 PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 1560 { 1561 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1562 PetscBool flg; 1563 1564 PetscFunctionBegin; 1565 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg)); 1566 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input"); 1567 if (Ad) *Ad = a->A; 1568 if (Ao) *Ao = a->B; 1569 if (colmap) *colmap = a->garray; 1570 PetscFunctionReturn(0); 1571 } 1572 1573 /*@C 1574 MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns 1575 1576 Not Collective 1577 1578 Input Parameters: 1579 + A - the matrix 1580 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1581 - row, col - index sets of rows and columns to extract (or NULL) 1582 1583 Output Parameter: 1584 . A_loc - the local sequential matrix generated 1585 1586 Level: developer 1587 1588 .seealso: `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()` 1589 1590 @*/ 1591 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc) 1592 { 1593 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1594 PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx; 1595 IS isrowa,iscola; 1596 Mat *aloc; 1597 PetscBool match; 1598 1599 PetscFunctionBegin; 1600 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match)); 1601 PetscCheck(match,PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input"); 1602 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0)); 1603 if (!row) { 1604 start = A->rmap->rstart; end = A->rmap->rend; 1605 PetscCall(ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa)); 1606 } else { 1607 isrowa = *row; 1608 } 1609 if (!col) { 1610 start = A->cmap->rstart; 1611 cmap = a->garray; 1612 nzA = a->A->cmap->n; 1613 nzB = a->B->cmap->n; 1614 PetscCall(PetscMalloc1(nzA+nzB, &idx)); 1615 ncols = 0; 1616 for (i=0; i<nzB; i++) { 1617 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1618 else break; 1619 } 1620 imark = i; 1621 for (i=0; i<nzA; i++) idx[ncols++] = start + i; 1622 for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; 1623 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola)); 1624 } else { 1625 iscola = *col; 1626 } 1627 if (scall != MAT_INITIAL_MATRIX) { 1628 PetscCall(PetscMalloc1(1,&aloc)); 1629 aloc[0] = *A_loc; 1630 } 1631 PetscCall(MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc)); 1632 *A_loc = aloc[0]; 1633 PetscCall(PetscFree(aloc)); 1634 if (!row) { 1635 PetscCall(ISDestroy(&isrowa)); 1636 } 1637 if (!col) { 1638 PetscCall(ISDestroy(&iscola)); 1639 } 1640 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0)); 1641 PetscFunctionReturn(0); 1642 } 1643 1644 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1645 1646 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1647 { 1648 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1649 Mat B; 1650 Mat_MPIAIJ *b; 1651 1652 PetscFunctionBegin; 1653 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1654 1655 if (reuse == MAT_REUSE_MATRIX) { 1656 B = *newmat; 1657 } else { 1658 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1659 PetscCall(MatSetType(B,MATMPIAIJ)); 1660 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1661 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1662 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1663 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1664 } 1665 b = (Mat_MPIAIJ*) B->data; 1666 1667 if (reuse == MAT_REUSE_MATRIX) { 1668 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1669 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1670 } else { 1671 PetscCall(MatDestroy(&b->A)); 1672 PetscCall(MatDestroy(&b->B)); 1673 PetscCall(MatDisAssemble_MPISELL(A)); 1674 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1675 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1676 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1677 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1678 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1679 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1680 } 1681 1682 if (reuse == MAT_INPLACE_MATRIX) { 1683 PetscCall(MatHeaderReplace(A,&B)); 1684 } else { 1685 *newmat = B; 1686 } 1687 PetscFunctionReturn(0); 1688 } 1689 1690 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1691 { 1692 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 1693 Mat B; 1694 Mat_MPISELL *b; 1695 1696 PetscFunctionBegin; 1697 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1698 1699 if (reuse == MAT_REUSE_MATRIX) { 1700 B = *newmat; 1701 } else { 1702 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1703 PetscCall(MatSetType(B,MATMPISELL)); 1704 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1705 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1706 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1707 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1708 } 1709 b = (Mat_MPISELL*) B->data; 1710 1711 if (reuse == MAT_REUSE_MATRIX) { 1712 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1713 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1714 } else { 1715 PetscCall(MatDestroy(&b->A)); 1716 PetscCall(MatDestroy(&b->B)); 1717 PetscCall(MatDisAssemble_MPIAIJ(A)); 1718 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1719 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1720 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1721 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1722 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1723 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1724 } 1725 1726 if (reuse == MAT_INPLACE_MATRIX) { 1727 PetscCall(MatHeaderReplace(A,&B)); 1728 } else { 1729 *newmat = B; 1730 } 1731 PetscFunctionReturn(0); 1732 } 1733 1734 PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1735 { 1736 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 1737 Vec bb1=NULL; 1738 1739 PetscFunctionBegin; 1740 if (flag == SOR_APPLY_UPPER) { 1741 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1742 PetscFunctionReturn(0); 1743 } 1744 1745 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) { 1746 PetscCall(VecDuplicate(bb,&bb1)); 1747 } 1748 1749 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1750 if (flag & SOR_ZERO_INITIAL_GUESS) { 1751 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1752 its--; 1753 } 1754 1755 while (its--) { 1756 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1757 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1758 1759 /* update rhs: bb1 = bb - B*x */ 1760 PetscCall(VecScale(mat->lvec,-1.0)); 1761 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1762 1763 /* local sweep */ 1764 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx)); 1765 } 1766 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1767 if (flag & SOR_ZERO_INITIAL_GUESS) { 1768 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1769 its--; 1770 } 1771 while (its--) { 1772 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1773 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1774 1775 /* update rhs: bb1 = bb - B*x */ 1776 PetscCall(VecScale(mat->lvec,-1.0)); 1777 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1778 1779 /* local sweep */ 1780 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx)); 1781 } 1782 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1783 if (flag & SOR_ZERO_INITIAL_GUESS) { 1784 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1785 its--; 1786 } 1787 while (its--) { 1788 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1789 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1790 1791 /* update rhs: bb1 = bb - B*x */ 1792 PetscCall(VecScale(mat->lvec,-1.0)); 1793 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1794 1795 /* local sweep */ 1796 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx)); 1797 } 1798 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported"); 1799 1800 PetscCall(VecDestroy(&bb1)); 1801 1802 matin->factorerrortype = mat->A->factorerrortype; 1803 PetscFunctionReturn(0); 1804 } 1805 1806 /*MC 1807 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1808 1809 Options Database Keys: 1810 . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions() 1811 1812 Level: beginner 1813 1814 .seealso: `MatCreateSELL()` 1815 M*/ 1816 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1817 { 1818 Mat_MPISELL *b; 1819 PetscMPIInt size; 1820 1821 PetscFunctionBegin; 1822 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1823 PetscCall(PetscNewLog(B,&b)); 1824 B->data = (void*)b; 1825 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1826 B->assembled = PETSC_FALSE; 1827 B->insertmode = NOT_SET_VALUES; 1828 b->size = size; 1829 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank)); 1830 /* build cache for off array entries formed */ 1831 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash)); 1832 1833 b->donotstash = PETSC_FALSE; 1834 b->colmap = NULL; 1835 b->garray = NULL; 1836 b->roworiented = PETSC_TRUE; 1837 1838 /* stuff used for matrix vector multiply */ 1839 b->lvec = NULL; 1840 b->Mvctx = NULL; 1841 1842 /* stuff for MatGetRow() */ 1843 b->rowindices = NULL; 1844 b->rowvalues = NULL; 1845 b->getrowactive = PETSC_FALSE; 1846 1847 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL)); 1848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL)); 1849 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL)); 1850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL)); 1851 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ)); 1852 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL)); 1853 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATMPISELL)); 1854 PetscFunctionReturn(0); 1855 } 1856