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 NULL, 1119 NULL 1120 }; 1121 1122 /* ----------------------------------------------------------------------------------------*/ 1123 1124 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1125 { 1126 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1127 1128 PetscFunctionBegin; 1129 PetscCall(MatStoreValues(sell->A)); 1130 PetscCall(MatStoreValues(sell->B)); 1131 PetscFunctionReturn(0); 1132 } 1133 1134 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1135 { 1136 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1137 1138 PetscFunctionBegin; 1139 PetscCall(MatRetrieveValues(sell->A)); 1140 PetscCall(MatRetrieveValues(sell->B)); 1141 PetscFunctionReturn(0); 1142 } 1143 1144 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[]) 1145 { 1146 Mat_MPISELL *b; 1147 1148 PetscFunctionBegin; 1149 PetscCall(PetscLayoutSetUp(B->rmap)); 1150 PetscCall(PetscLayoutSetUp(B->cmap)); 1151 b = (Mat_MPISELL*)B->data; 1152 1153 if (!B->preallocated) { 1154 /* Explicitly create 2 MATSEQSELL matrices. */ 1155 PetscCall(MatCreate(PETSC_COMM_SELF,&b->A)); 1156 PetscCall(MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n)); 1157 PetscCall(MatSetBlockSizesFromMats(b->A,B,B)); 1158 PetscCall(MatSetType(b->A,MATSEQSELL)); 1159 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->A)); 1160 PetscCall(MatCreate(PETSC_COMM_SELF,&b->B)); 1161 PetscCall(MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N)); 1162 PetscCall(MatSetBlockSizesFromMats(b->B,B,B)); 1163 PetscCall(MatSetType(b->B,MATSEQSELL)); 1164 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->B)); 1165 } 1166 1167 PetscCall(MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen)); 1168 PetscCall(MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen)); 1169 B->preallocated = PETSC_TRUE; 1170 B->was_assembled = PETSC_FALSE; 1171 /* 1172 critical for MatAssemblyEnd to work. 1173 MatAssemblyBegin checks it to set up was_assembled 1174 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1175 */ 1176 B->assembled = PETSC_FALSE; 1177 PetscFunctionReturn(0); 1178 } 1179 1180 PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 1181 { 1182 Mat mat; 1183 Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data; 1184 1185 PetscFunctionBegin; 1186 *newmat = NULL; 1187 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin),&mat)); 1188 PetscCall(MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N)); 1189 PetscCall(MatSetBlockSizesFromMats(mat,matin,matin)); 1190 PetscCall(MatSetType(mat,((PetscObject)matin)->type_name)); 1191 a = (Mat_MPISELL*)mat->data; 1192 1193 mat->factortype = matin->factortype; 1194 mat->assembled = PETSC_TRUE; 1195 mat->insertmode = NOT_SET_VALUES; 1196 mat->preallocated = PETSC_TRUE; 1197 1198 a->size = oldmat->size; 1199 a->rank = oldmat->rank; 1200 a->donotstash = oldmat->donotstash; 1201 a->roworiented = oldmat->roworiented; 1202 a->rowindices = NULL; 1203 a->rowvalues = NULL; 1204 a->getrowactive = PETSC_FALSE; 1205 1206 PetscCall(PetscLayoutReference(matin->rmap,&mat->rmap)); 1207 PetscCall(PetscLayoutReference(matin->cmap,&mat->cmap)); 1208 1209 if (oldmat->colmap) { 1210 #if defined(PETSC_USE_CTABLE) 1211 PetscCall(PetscTableCreateCopy(oldmat->colmap,&a->colmap)); 1212 #else 1213 PetscCall(PetscMalloc1(mat->cmap->N,&a->colmap)); 1214 PetscCall(PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt))); 1215 PetscCall(PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N)); 1216 #endif 1217 } else a->colmap = NULL; 1218 if (oldmat->garray) { 1219 PetscInt len; 1220 len = oldmat->B->cmap->n; 1221 PetscCall(PetscMalloc1(len+1,&a->garray)); 1222 PetscCall(PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt))); 1223 if (len) PetscCall(PetscArraycpy(a->garray,oldmat->garray,len)); 1224 } else a->garray = NULL; 1225 1226 PetscCall(VecDuplicate(oldmat->lvec,&a->lvec)); 1227 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec)); 1228 PetscCall(VecScatterCopy(oldmat->Mvctx,&a->Mvctx)); 1229 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx)); 1230 PetscCall(MatDuplicate(oldmat->A,cpvalues,&a->A)); 1231 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A)); 1232 PetscCall(MatDuplicate(oldmat->B,cpvalues,&a->B)); 1233 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B)); 1234 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist)); 1235 *newmat = mat; 1236 PetscFunctionReturn(0); 1237 } 1238 1239 /*@C 1240 MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format. 1241 For good matrix assembly performance the user should preallocate the matrix storage by 1242 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1243 1244 Collective 1245 1246 Input Parameters: 1247 + B - the matrix 1248 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1249 (same value is used for all local rows) 1250 . d_nnz - array containing the number of nonzeros in the various rows of the 1251 DIAGONAL portion of the local submatrix (possibly different for each row) 1252 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1253 The size of this array is equal to the number of local rows, i.e 'm'. 1254 For matrices that will be factored, you must leave room for (and set) 1255 the diagonal entry even if it is zero. 1256 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1257 submatrix (same value is used for all local rows). 1258 - o_nnz - array containing the number of nonzeros in the various rows of the 1259 OFF-DIAGONAL portion of the local submatrix (possibly different for 1260 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1261 structure. The size of this array is equal to the number 1262 of local rows, i.e 'm'. 1263 1264 If the *_nnz parameter is given then the *_nz parameter is ignored 1265 1266 The stored row and column indices begin with zero. 1267 1268 The parallel matrix is partitioned such that the first m0 rows belong to 1269 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1270 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1271 1272 The DIAGONAL portion of the local submatrix of a processor can be defined 1273 as the submatrix which is obtained by extraction the part corresponding to 1274 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1275 first row that belongs to the processor, r2 is the last row belonging to 1276 the this processor, and c1-c2 is range of indices of the local part of a 1277 vector suitable for applying the matrix to. This is an mxn matrix. In the 1278 common case of a square matrix, the row and column ranges are the same and 1279 the DIAGONAL part is also square. The remaining portion of the local 1280 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1281 1282 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1283 1284 You can call MatGetInfo() to get information on how effective the preallocation was; 1285 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1286 You can also run with the option -info and look for messages with the string 1287 malloc in them to see if additional memory allocation was needed. 1288 1289 Example usage: 1290 1291 Consider the following 8x8 matrix with 34 non-zero values, that is 1292 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1293 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1294 as follows: 1295 1296 .vb 1297 1 2 0 | 0 3 0 | 0 4 1298 Proc0 0 5 6 | 7 0 0 | 8 0 1299 9 0 10 | 11 0 0 | 12 0 1300 ------------------------------------- 1301 13 0 14 | 15 16 17 | 0 0 1302 Proc1 0 18 0 | 19 20 21 | 0 0 1303 0 0 0 | 22 23 0 | 24 0 1304 ------------------------------------- 1305 Proc2 25 26 27 | 0 0 28 | 29 0 1306 30 0 0 | 31 32 33 | 0 34 1307 .ve 1308 1309 This can be represented as a collection of submatrices as: 1310 1311 .vb 1312 A B C 1313 D E F 1314 G H I 1315 .ve 1316 1317 Where the submatrices A,B,C are owned by proc0, D,E,F are 1318 owned by proc1, G,H,I are owned by proc2. 1319 1320 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1321 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1322 The 'M','N' parameters are 8,8, and have the same values on all procs. 1323 1324 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1325 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1326 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1327 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1328 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1329 matrix, ans [DF] as another SeqSELL matrix. 1330 1331 When d_nz, o_nz parameters are specified, d_nz storage elements are 1332 allocated for every row of the local diagonal submatrix, and o_nz 1333 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1334 One way to choose d_nz and o_nz is to use the max nonzerors per local 1335 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1336 In this case, the values of d_nz,o_nz are: 1337 .vb 1338 proc0 : dnz = 2, o_nz = 2 1339 proc1 : dnz = 3, o_nz = 2 1340 proc2 : dnz = 1, o_nz = 4 1341 .ve 1342 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1343 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1344 for proc3. i.e we are using 12+15+10=37 storage locations to store 1345 34 values. 1346 1347 When d_nnz, o_nnz parameters are specified, the storage is specified 1348 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1349 In the above case the values for d_nnz,o_nnz are: 1350 .vb 1351 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1352 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1353 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1354 .ve 1355 Here the space allocated is according to nz (or maximum values in the nnz 1356 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1357 1358 Level: intermediate 1359 1360 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreatesell()`, 1361 `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()` 1362 @*/ 1363 PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1364 { 1365 PetscFunctionBegin; 1366 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1367 PetscValidType(B,1); 1368 PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz)); 1369 PetscFunctionReturn(0); 1370 } 1371 1372 /*MC 1373 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1374 based on the sliced Ellpack format 1375 1376 Options Database Keys: 1377 . -mat_type sell - sets the matrix type to "seqsell" during a call to MatSetFromOptions() 1378 1379 Level: beginner 1380 1381 .seealso: `MatCreateSell()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ` 1382 M*/ 1383 1384 /*@C 1385 MatCreateSELL - Creates a sparse parallel matrix in SELL format. 1386 1387 Collective 1388 1389 Input Parameters: 1390 + comm - MPI communicator 1391 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1392 This value should be the same as the local size used in creating the 1393 y vector for the matrix-vector product y = Ax. 1394 . n - This value should be the same as the local size used in creating the 1395 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1396 calculated if N is given) For square matrices n is almost always m. 1397 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 1398 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 1399 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1400 (same value is used for all local rows) 1401 . d_rlen - array containing the number of nonzeros in the various rows of the 1402 DIAGONAL portion of the local submatrix (possibly different for each row) 1403 or NULL, if d_rlenmax is used to specify the nonzero structure. 1404 The size of this array is equal to the number of local rows, i.e 'm'. 1405 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1406 submatrix (same value is used for all local rows). 1407 - o_rlen - array containing the number of nonzeros in the various rows of the 1408 OFF-DIAGONAL portion of the local submatrix (possibly different for 1409 each row) or NULL, if o_rlenmax is used to specify the nonzero 1410 structure. The size of this array is equal to the number 1411 of local rows, i.e 'm'. 1412 1413 Output Parameter: 1414 . A - the matrix 1415 1416 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1417 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1418 [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation] 1419 1420 Notes: 1421 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1422 1423 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1424 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1425 storage requirements for this matrix. 1426 1427 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one 1428 processor than it must be used on all processors that share the object for 1429 that argument. 1430 1431 The user MUST specify either the local or global matrix dimensions 1432 (possibly both). 1433 1434 The parallel matrix is partitioned across processors such that the 1435 first m0 rows belong to process 0, the next m1 rows belong to 1436 process 1, the next m2 rows belong to process 2 etc.. where 1437 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1438 values corresponding to [m x N] submatrix. 1439 1440 The columns are logically partitioned with the n0 columns belonging 1441 to 0th partition, the next n1 columns belonging to the next 1442 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1443 1444 The DIAGONAL portion of the local submatrix on any given processor 1445 is the submatrix corresponding to the rows and columns m,n 1446 corresponding to the given processor. i.e diagonal matrix on 1447 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1448 etc. The remaining portion of the local submatrix [m x (N-n)] 1449 constitute the OFF-DIAGONAL portion. The example below better 1450 illustrates this concept. 1451 1452 For a square global matrix we define each processor's diagonal portion 1453 to be its local rows and the corresponding columns (a square submatrix); 1454 each processor's off-diagonal portion encompasses the remainder of the 1455 local matrix (a rectangular submatrix). 1456 1457 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1458 1459 When calling this routine with a single process communicator, a matrix of 1460 type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this 1461 type of communicator, use the construction mechanism: 1462 MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...); 1463 1464 Options Database Keys: 1465 - -mat_sell_oneindex - Internally use indexing starting at 1 1466 rather than 0. Note that when calling MatSetValues(), 1467 the user still MUST index entries starting at 0! 1468 1469 Example usage: 1470 1471 Consider the following 8x8 matrix with 34 non-zero values, that is 1472 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1473 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1474 as follows: 1475 1476 .vb 1477 1 2 0 | 0 3 0 | 0 4 1478 Proc0 0 5 6 | 7 0 0 | 8 0 1479 9 0 10 | 11 0 0 | 12 0 1480 ------------------------------------- 1481 13 0 14 | 15 16 17 | 0 0 1482 Proc1 0 18 0 | 19 20 21 | 0 0 1483 0 0 0 | 22 23 0 | 24 0 1484 ------------------------------------- 1485 Proc2 25 26 27 | 0 0 28 | 29 0 1486 30 0 0 | 31 32 33 | 0 34 1487 .ve 1488 1489 This can be represented as a collection of submatrices as: 1490 1491 .vb 1492 A B C 1493 D E F 1494 G H I 1495 .ve 1496 1497 Where the submatrices A,B,C are owned by proc0, D,E,F are 1498 owned by proc1, G,H,I are owned by proc2. 1499 1500 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1501 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1502 The 'M','N' parameters are 8,8, and have the same values on all procs. 1503 1504 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1505 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1506 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1507 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1508 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1509 matrix, ans [DF] as another SeqSELL matrix. 1510 1511 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1512 allocated for every row of the local diagonal submatrix, and o_rlenmax 1513 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1514 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1515 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1516 In this case, the values of d_rlenmax,o_rlenmax are: 1517 .vb 1518 proc0 : d_rlenmax = 2, o_rlenmax = 2 1519 proc1 : d_rlenmax = 3, o_rlenmax = 2 1520 proc2 : d_rlenmax = 1, o_rlenmax = 4 1521 .ve 1522 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1523 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1524 for proc3. i.e we are using 12+15+10=37 storage locations to store 1525 34 values. 1526 1527 When d_rlen, o_rlen parameters are specified, the storage is specified 1528 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1529 In the above case the values for d_nnz,o_nnz are: 1530 .vb 1531 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1532 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1533 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1534 .ve 1535 Here the space allocated is still 37 though there are 34 nonzeros because 1536 the allocation is always done according to rlenmax. 1537 1538 Level: intermediate 1539 1540 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MatMPISELLSetPreallocationSELL()`, 1541 `MATMPISELL`, `MatCreateMPISELLWithArrays()` 1542 @*/ 1543 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) 1544 { 1545 PetscMPIInt size; 1546 1547 PetscFunctionBegin; 1548 PetscCall(MatCreate(comm,A)); 1549 PetscCall(MatSetSizes(*A,m,n,M,N)); 1550 PetscCallMPI(MPI_Comm_size(comm,&size)); 1551 if (size > 1) { 1552 PetscCall(MatSetType(*A,MATMPISELL)); 1553 PetscCall(MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen)); 1554 } else { 1555 PetscCall(MatSetType(*A,MATSEQSELL)); 1556 PetscCall(MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen)); 1557 } 1558 PetscFunctionReturn(0); 1559 } 1560 1561 PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 1562 { 1563 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1564 PetscBool flg; 1565 1566 PetscFunctionBegin; 1567 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg)); 1568 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input"); 1569 if (Ad) *Ad = a->A; 1570 if (Ao) *Ao = a->B; 1571 if (colmap) *colmap = a->garray; 1572 PetscFunctionReturn(0); 1573 } 1574 1575 /*@C 1576 MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns 1577 1578 Not Collective 1579 1580 Input Parameters: 1581 + A - the matrix 1582 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1583 - row, col - index sets of rows and columns to extract (or NULL) 1584 1585 Output Parameter: 1586 . A_loc - the local sequential matrix generated 1587 1588 Level: developer 1589 1590 .seealso: `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()` 1591 1592 @*/ 1593 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc) 1594 { 1595 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1596 PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx; 1597 IS isrowa,iscola; 1598 Mat *aloc; 1599 PetscBool match; 1600 1601 PetscFunctionBegin; 1602 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match)); 1603 PetscCheck(match,PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input"); 1604 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0)); 1605 if (!row) { 1606 start = A->rmap->rstart; end = A->rmap->rend; 1607 PetscCall(ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa)); 1608 } else { 1609 isrowa = *row; 1610 } 1611 if (!col) { 1612 start = A->cmap->rstart; 1613 cmap = a->garray; 1614 nzA = a->A->cmap->n; 1615 nzB = a->B->cmap->n; 1616 PetscCall(PetscMalloc1(nzA+nzB, &idx)); 1617 ncols = 0; 1618 for (i=0; i<nzB; i++) { 1619 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1620 else break; 1621 } 1622 imark = i; 1623 for (i=0; i<nzA; i++) idx[ncols++] = start + i; 1624 for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; 1625 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola)); 1626 } else { 1627 iscola = *col; 1628 } 1629 if (scall != MAT_INITIAL_MATRIX) { 1630 PetscCall(PetscMalloc1(1,&aloc)); 1631 aloc[0] = *A_loc; 1632 } 1633 PetscCall(MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc)); 1634 *A_loc = aloc[0]; 1635 PetscCall(PetscFree(aloc)); 1636 if (!row) { 1637 PetscCall(ISDestroy(&isrowa)); 1638 } 1639 if (!col) { 1640 PetscCall(ISDestroy(&iscola)); 1641 } 1642 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0)); 1643 PetscFunctionReturn(0); 1644 } 1645 1646 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1647 1648 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1649 { 1650 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1651 Mat B; 1652 Mat_MPIAIJ *b; 1653 1654 PetscFunctionBegin; 1655 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1656 1657 if (reuse == MAT_REUSE_MATRIX) { 1658 B = *newmat; 1659 } else { 1660 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1661 PetscCall(MatSetType(B,MATMPIAIJ)); 1662 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1663 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1664 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1665 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1666 } 1667 b = (Mat_MPIAIJ*) B->data; 1668 1669 if (reuse == MAT_REUSE_MATRIX) { 1670 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1671 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1672 } else { 1673 PetscCall(MatDestroy(&b->A)); 1674 PetscCall(MatDestroy(&b->B)); 1675 PetscCall(MatDisAssemble_MPISELL(A)); 1676 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1677 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1678 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1679 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1680 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1681 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1682 } 1683 1684 if (reuse == MAT_INPLACE_MATRIX) { 1685 PetscCall(MatHeaderReplace(A,&B)); 1686 } else { 1687 *newmat = B; 1688 } 1689 PetscFunctionReturn(0); 1690 } 1691 1692 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1693 { 1694 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 1695 Mat B; 1696 Mat_MPISELL *b; 1697 1698 PetscFunctionBegin; 1699 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1700 1701 if (reuse == MAT_REUSE_MATRIX) { 1702 B = *newmat; 1703 } else { 1704 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1705 PetscCall(MatSetType(B,MATMPISELL)); 1706 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1707 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1708 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1709 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1710 } 1711 b = (Mat_MPISELL*) B->data; 1712 1713 if (reuse == MAT_REUSE_MATRIX) { 1714 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1715 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1716 } else { 1717 PetscCall(MatDestroy(&b->A)); 1718 PetscCall(MatDestroy(&b->B)); 1719 PetscCall(MatDisAssemble_MPIAIJ(A)); 1720 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1721 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1722 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1723 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1724 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1725 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1726 } 1727 1728 if (reuse == MAT_INPLACE_MATRIX) { 1729 PetscCall(MatHeaderReplace(A,&B)); 1730 } else { 1731 *newmat = B; 1732 } 1733 PetscFunctionReturn(0); 1734 } 1735 1736 PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1737 { 1738 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 1739 Vec bb1=NULL; 1740 1741 PetscFunctionBegin; 1742 if (flag == SOR_APPLY_UPPER) { 1743 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1744 PetscFunctionReturn(0); 1745 } 1746 1747 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) { 1748 PetscCall(VecDuplicate(bb,&bb1)); 1749 } 1750 1751 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1752 if (flag & SOR_ZERO_INITIAL_GUESS) { 1753 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1754 its--; 1755 } 1756 1757 while (its--) { 1758 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1759 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1760 1761 /* update rhs: bb1 = bb - B*x */ 1762 PetscCall(VecScale(mat->lvec,-1.0)); 1763 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1764 1765 /* local sweep */ 1766 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx)); 1767 } 1768 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1769 if (flag & SOR_ZERO_INITIAL_GUESS) { 1770 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1771 its--; 1772 } 1773 while (its--) { 1774 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1775 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1776 1777 /* update rhs: bb1 = bb - B*x */ 1778 PetscCall(VecScale(mat->lvec,-1.0)); 1779 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1780 1781 /* local sweep */ 1782 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx)); 1783 } 1784 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1785 if (flag & SOR_ZERO_INITIAL_GUESS) { 1786 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1787 its--; 1788 } 1789 while (its--) { 1790 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1791 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1792 1793 /* update rhs: bb1 = bb - B*x */ 1794 PetscCall(VecScale(mat->lvec,-1.0)); 1795 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1796 1797 /* local sweep */ 1798 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx)); 1799 } 1800 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported"); 1801 1802 PetscCall(VecDestroy(&bb1)); 1803 1804 matin->factorerrortype = mat->A->factorerrortype; 1805 PetscFunctionReturn(0); 1806 } 1807 1808 /*MC 1809 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1810 1811 Options Database Keys: 1812 . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions() 1813 1814 Level: beginner 1815 1816 .seealso: `MatCreateSELL()` 1817 M*/ 1818 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1819 { 1820 Mat_MPISELL *b; 1821 PetscMPIInt size; 1822 1823 PetscFunctionBegin; 1824 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1825 PetscCall(PetscNewLog(B,&b)); 1826 B->data = (void*)b; 1827 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1828 B->assembled = PETSC_FALSE; 1829 B->insertmode = NOT_SET_VALUES; 1830 b->size = size; 1831 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank)); 1832 /* build cache for off array entries formed */ 1833 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash)); 1834 1835 b->donotstash = PETSC_FALSE; 1836 b->colmap = NULL; 1837 b->garray = NULL; 1838 b->roworiented = PETSC_TRUE; 1839 1840 /* stuff used for matrix vector multiply */ 1841 b->lvec = NULL; 1842 b->Mvctx = NULL; 1843 1844 /* stuff for MatGetRow() */ 1845 b->rowindices = NULL; 1846 b->rowvalues = NULL; 1847 b->getrowactive = PETSC_FALSE; 1848 1849 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL)); 1850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL)); 1851 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL)); 1852 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL)); 1853 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ)); 1854 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL)); 1855 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATMPISELL)); 1856 PetscFunctionReturn(0); 1857 } 1858