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 case MAT_SPD_ETERNAL: 750 break; 751 case MAT_SYMMETRIC: 752 MatCheckPreallocated(A,1); 753 PetscCall(MatSetOption(a->A,op,flg)); 754 break; 755 case MAT_STRUCTURALLY_SYMMETRIC: 756 MatCheckPreallocated(A,1); 757 PetscCall(MatSetOption(a->A,op,flg)); 758 break; 759 case MAT_HERMITIAN: 760 MatCheckPreallocated(A,1); 761 PetscCall(MatSetOption(a->A,op,flg)); 762 break; 763 case MAT_SYMMETRY_ETERNAL: 764 MatCheckPreallocated(A,1); 765 PetscCall(MatSetOption(a->A,op,flg)); 766 break; 767 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 768 MatCheckPreallocated(A,1); 769 PetscCall(MatSetOption(a->A,op,flg)); 770 break; 771 default: 772 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 773 } 774 PetscFunctionReturn(0); 775 } 776 777 PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr) 778 { 779 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 780 Mat a=sell->A,b=sell->B; 781 PetscInt s1,s2,s3; 782 783 PetscFunctionBegin; 784 PetscCall(MatGetLocalSize(mat,&s2,&s3)); 785 if (rr) { 786 PetscCall(VecGetLocalSize(rr,&s1)); 787 PetscCheck(s1==s3,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size"); 788 /* Overlap communication with computation. */ 789 PetscCall(VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD)); 790 } 791 if (ll) { 792 PetscCall(VecGetLocalSize(ll,&s1)); 793 PetscCheck(s1==s2,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size"); 794 PetscCall((*b->ops->diagonalscale)(b,ll,NULL)); 795 } 796 /* scale the diagonal block */ 797 PetscCall((*a->ops->diagonalscale)(a,ll,rr)); 798 799 if (rr) { 800 /* Do a scatter end and then right scale the off-diagonal block */ 801 PetscCall(VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD)); 802 PetscCall((*b->ops->diagonalscale)(b,NULL,sell->lvec)); 803 } 804 PetscFunctionReturn(0); 805 } 806 807 PetscErrorCode MatSetUnfactored_MPISELL(Mat A) 808 { 809 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 810 811 PetscFunctionBegin; 812 PetscCall(MatSetUnfactored(a->A)); 813 PetscFunctionReturn(0); 814 } 815 816 PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool *flag) 817 { 818 Mat_MPISELL *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data; 819 Mat a,b,c,d; 820 PetscBool flg; 821 822 PetscFunctionBegin; 823 a = matA->A; b = matA->B; 824 c = matB->A; d = matB->B; 825 826 PetscCall(MatEqual(a,c,&flg)); 827 if (flg) { 828 PetscCall(MatEqual(b,d,&flg)); 829 } 830 PetscCall(MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 831 PetscFunctionReturn(0); 832 } 833 834 PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str) 835 { 836 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 837 Mat_MPISELL *b=(Mat_MPISELL*)B->data; 838 839 PetscFunctionBegin; 840 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 841 if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) { 842 /* because of the column compression in the off-processor part of the matrix a->B, 843 the number of columns in a->B and b->B may be different, hence we cannot call 844 the MatCopy() directly on the two parts. If need be, we can provide a more 845 efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices 846 then copying the submatrices */ 847 PetscCall(MatCopy_Basic(A,B,str)); 848 } else { 849 PetscCall(MatCopy(a->A,b->A,str)); 850 PetscCall(MatCopy(a->B,b->B,str)); 851 } 852 PetscFunctionReturn(0); 853 } 854 855 PetscErrorCode MatSetUp_MPISELL(Mat A) 856 { 857 PetscFunctionBegin; 858 PetscCall(MatMPISELLSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL)); 859 PetscFunctionReturn(0); 860 } 861 862 extern PetscErrorCode MatConjugate_SeqSELL(Mat); 863 864 PetscErrorCode MatConjugate_MPISELL(Mat mat) 865 { 866 PetscFunctionBegin; 867 if (PetscDefined(USE_COMPLEX)) { 868 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 869 870 PetscCall(MatConjugate_SeqSELL(sell->A)); 871 PetscCall(MatConjugate_SeqSELL(sell->B)); 872 } 873 PetscFunctionReturn(0); 874 } 875 876 PetscErrorCode MatRealPart_MPISELL(Mat A) 877 { 878 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 879 880 PetscFunctionBegin; 881 PetscCall(MatRealPart(a->A)); 882 PetscCall(MatRealPart(a->B)); 883 PetscFunctionReturn(0); 884 } 885 886 PetscErrorCode MatImaginaryPart_MPISELL(Mat A) 887 { 888 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 889 890 PetscFunctionBegin; 891 PetscCall(MatImaginaryPart(a->A)); 892 PetscCall(MatImaginaryPart(a->B)); 893 PetscFunctionReturn(0); 894 } 895 896 PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values) 897 { 898 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 899 900 PetscFunctionBegin; 901 PetscCall(MatInvertBlockDiagonal(a->A,values)); 902 A->factorerrortype = a->A->factorerrortype; 903 PetscFunctionReturn(0); 904 } 905 906 static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx) 907 { 908 Mat_MPISELL *sell=(Mat_MPISELL*)x->data; 909 910 PetscFunctionBegin; 911 PetscCall(MatSetRandom(sell->A,rctx)); 912 PetscCall(MatSetRandom(sell->B,rctx)); 913 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 914 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 915 PetscFunctionReturn(0); 916 } 917 918 PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A) 919 { 920 PetscFunctionBegin; 921 PetscOptionsHeadBegin(PetscOptionsObject,"MPISELL options"); 922 PetscOptionsHeadEnd(); 923 PetscFunctionReturn(0); 924 } 925 926 PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a) 927 { 928 Mat_MPISELL *msell=(Mat_MPISELL*)Y->data; 929 Mat_SeqSELL *sell=(Mat_SeqSELL*)msell->A->data; 930 931 PetscFunctionBegin; 932 if (!Y->preallocated) { 933 PetscCall(MatMPISELLSetPreallocation(Y,1,NULL,0,NULL)); 934 } else if (!sell->nz) { 935 PetscInt nonew = sell->nonew; 936 PetscCall(MatSeqSELLSetPreallocation(msell->A,1,NULL)); 937 sell->nonew = nonew; 938 } 939 PetscCall(MatShift_Basic(Y,a)); 940 PetscFunctionReturn(0); 941 } 942 943 PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool *missing,PetscInt *d) 944 { 945 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 946 947 PetscFunctionBegin; 948 PetscCheck(A->rmap->n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices"); 949 PetscCall(MatMissingDiagonal(a->A,missing,d)); 950 if (d) { 951 PetscInt rstart; 952 PetscCall(MatGetOwnershipRange(A,&rstart,NULL)); 953 *d += rstart; 954 955 } 956 PetscFunctionReturn(0); 957 } 958 959 PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a) 960 { 961 PetscFunctionBegin; 962 *a = ((Mat_MPISELL*)A->data)->A; 963 PetscFunctionReturn(0); 964 } 965 966 /* -------------------------------------------------------------------*/ 967 static struct _MatOps MatOps_Values = {MatSetValues_MPISELL, 968 NULL, 969 NULL, 970 MatMult_MPISELL, 971 /* 4*/ MatMultAdd_MPISELL, 972 MatMultTranspose_MPISELL, 973 MatMultTransposeAdd_MPISELL, 974 NULL, 975 NULL, 976 NULL, 977 /*10*/ NULL, 978 NULL, 979 NULL, 980 MatSOR_MPISELL, 981 NULL, 982 /*15*/ MatGetInfo_MPISELL, 983 MatEqual_MPISELL, 984 MatGetDiagonal_MPISELL, 985 MatDiagonalScale_MPISELL, 986 NULL, 987 /*20*/ MatAssemblyBegin_MPISELL, 988 MatAssemblyEnd_MPISELL, 989 MatSetOption_MPISELL, 990 MatZeroEntries_MPISELL, 991 /*24*/ NULL, 992 NULL, 993 NULL, 994 NULL, 995 NULL, 996 /*29*/ MatSetUp_MPISELL, 997 NULL, 998 NULL, 999 MatGetDiagonalBlock_MPISELL, 1000 NULL, 1001 /*34*/ MatDuplicate_MPISELL, 1002 NULL, 1003 NULL, 1004 NULL, 1005 NULL, 1006 /*39*/ NULL, 1007 NULL, 1008 NULL, 1009 MatGetValues_MPISELL, 1010 MatCopy_MPISELL, 1011 /*44*/ NULL, 1012 MatScale_MPISELL, 1013 MatShift_MPISELL, 1014 MatDiagonalSet_MPISELL, 1015 NULL, 1016 /*49*/ MatSetRandom_MPISELL, 1017 NULL, 1018 NULL, 1019 NULL, 1020 NULL, 1021 /*54*/ MatFDColoringCreate_MPIXAIJ, 1022 NULL, 1023 MatSetUnfactored_MPISELL, 1024 NULL, 1025 NULL, 1026 /*59*/ NULL, 1027 MatDestroy_MPISELL, 1028 MatView_MPISELL, 1029 NULL, 1030 NULL, 1031 /*64*/ NULL, 1032 NULL, 1033 NULL, 1034 NULL, 1035 NULL, 1036 /*69*/ NULL, 1037 NULL, 1038 NULL, 1039 NULL, 1040 NULL, 1041 NULL, 1042 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */ 1043 MatSetFromOptions_MPISELL, 1044 NULL, 1045 NULL, 1046 NULL, 1047 /*80*/ NULL, 1048 NULL, 1049 NULL, 1050 /*83*/ NULL, 1051 NULL, 1052 NULL, 1053 NULL, 1054 NULL, 1055 NULL, 1056 /*89*/ NULL, 1057 NULL, 1058 NULL, 1059 NULL, 1060 NULL, 1061 /*94*/ NULL, 1062 NULL, 1063 NULL, 1064 NULL, 1065 NULL, 1066 /*99*/ NULL, 1067 NULL, 1068 NULL, 1069 MatConjugate_MPISELL, 1070 NULL, 1071 /*104*/NULL, 1072 MatRealPart_MPISELL, 1073 MatImaginaryPart_MPISELL, 1074 NULL, 1075 NULL, 1076 /*109*/NULL, 1077 NULL, 1078 NULL, 1079 NULL, 1080 MatMissingDiagonal_MPISELL, 1081 /*114*/NULL, 1082 NULL, 1083 MatGetGhosts_MPISELL, 1084 NULL, 1085 NULL, 1086 /*119*/NULL, 1087 NULL, 1088 NULL, 1089 NULL, 1090 NULL, 1091 /*124*/NULL, 1092 NULL, 1093 MatInvertBlockDiagonal_MPISELL, 1094 NULL, 1095 NULL, 1096 /*129*/NULL, 1097 NULL, 1098 NULL, 1099 NULL, 1100 NULL, 1101 /*134*/NULL, 1102 NULL, 1103 NULL, 1104 NULL, 1105 NULL, 1106 /*139*/NULL, 1107 NULL, 1108 NULL, 1109 MatFDColoringSetUp_MPIXAIJ, 1110 NULL, 1111 /*144*/NULL, 1112 NULL, 1113 NULL, 1114 NULL, 1115 NULL, 1116 NULL 1117 }; 1118 1119 /* ----------------------------------------------------------------------------------------*/ 1120 1121 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1122 { 1123 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1124 1125 PetscFunctionBegin; 1126 PetscCall(MatStoreValues(sell->A)); 1127 PetscCall(MatStoreValues(sell->B)); 1128 PetscFunctionReturn(0); 1129 } 1130 1131 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1132 { 1133 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1134 1135 PetscFunctionBegin; 1136 PetscCall(MatRetrieveValues(sell->A)); 1137 PetscCall(MatRetrieveValues(sell->B)); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[]) 1142 { 1143 Mat_MPISELL *b; 1144 1145 PetscFunctionBegin; 1146 PetscCall(PetscLayoutSetUp(B->rmap)); 1147 PetscCall(PetscLayoutSetUp(B->cmap)); 1148 b = (Mat_MPISELL*)B->data; 1149 1150 if (!B->preallocated) { 1151 /* Explicitly create 2 MATSEQSELL matrices. */ 1152 PetscCall(MatCreate(PETSC_COMM_SELF,&b->A)); 1153 PetscCall(MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n)); 1154 PetscCall(MatSetBlockSizesFromMats(b->A,B,B)); 1155 PetscCall(MatSetType(b->A,MATSEQSELL)); 1156 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->A)); 1157 PetscCall(MatCreate(PETSC_COMM_SELF,&b->B)); 1158 PetscCall(MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N)); 1159 PetscCall(MatSetBlockSizesFromMats(b->B,B,B)); 1160 PetscCall(MatSetType(b->B,MATSEQSELL)); 1161 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->B)); 1162 } 1163 1164 PetscCall(MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen)); 1165 PetscCall(MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen)); 1166 B->preallocated = PETSC_TRUE; 1167 B->was_assembled = PETSC_FALSE; 1168 /* 1169 critical for MatAssemblyEnd to work. 1170 MatAssemblyBegin checks it to set up was_assembled 1171 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1172 */ 1173 B->assembled = PETSC_FALSE; 1174 PetscFunctionReturn(0); 1175 } 1176 1177 PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 1178 { 1179 Mat mat; 1180 Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data; 1181 1182 PetscFunctionBegin; 1183 *newmat = NULL; 1184 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin),&mat)); 1185 PetscCall(MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N)); 1186 PetscCall(MatSetBlockSizesFromMats(mat,matin,matin)); 1187 PetscCall(MatSetType(mat,((PetscObject)matin)->type_name)); 1188 a = (Mat_MPISELL*)mat->data; 1189 1190 mat->factortype = matin->factortype; 1191 mat->assembled = PETSC_TRUE; 1192 mat->insertmode = NOT_SET_VALUES; 1193 mat->preallocated = PETSC_TRUE; 1194 1195 a->size = oldmat->size; 1196 a->rank = oldmat->rank; 1197 a->donotstash = oldmat->donotstash; 1198 a->roworiented = oldmat->roworiented; 1199 a->rowindices = NULL; 1200 a->rowvalues = NULL; 1201 a->getrowactive = PETSC_FALSE; 1202 1203 PetscCall(PetscLayoutReference(matin->rmap,&mat->rmap)); 1204 PetscCall(PetscLayoutReference(matin->cmap,&mat->cmap)); 1205 1206 if (oldmat->colmap) { 1207 #if defined(PETSC_USE_CTABLE) 1208 PetscCall(PetscTableCreateCopy(oldmat->colmap,&a->colmap)); 1209 #else 1210 PetscCall(PetscMalloc1(mat->cmap->N,&a->colmap)); 1211 PetscCall(PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt))); 1212 PetscCall(PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N)); 1213 #endif 1214 } else a->colmap = NULL; 1215 if (oldmat->garray) { 1216 PetscInt len; 1217 len = oldmat->B->cmap->n; 1218 PetscCall(PetscMalloc1(len+1,&a->garray)); 1219 PetscCall(PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt))); 1220 if (len) PetscCall(PetscArraycpy(a->garray,oldmat->garray,len)); 1221 } else a->garray = NULL; 1222 1223 PetscCall(VecDuplicate(oldmat->lvec,&a->lvec)); 1224 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec)); 1225 PetscCall(VecScatterCopy(oldmat->Mvctx,&a->Mvctx)); 1226 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx)); 1227 PetscCall(MatDuplicate(oldmat->A,cpvalues,&a->A)); 1228 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A)); 1229 PetscCall(MatDuplicate(oldmat->B,cpvalues,&a->B)); 1230 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B)); 1231 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist)); 1232 *newmat = mat; 1233 PetscFunctionReturn(0); 1234 } 1235 1236 /*@C 1237 MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format. 1238 For good matrix assembly performance the user should preallocate the matrix storage by 1239 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1240 1241 Collective 1242 1243 Input Parameters: 1244 + B - the matrix 1245 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1246 (same value is used for all local rows) 1247 . d_nnz - array containing the number of nonzeros in the various rows of the 1248 DIAGONAL portion of the local submatrix (possibly different for each row) 1249 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1250 The size of this array is equal to the number of local rows, i.e 'm'. 1251 For matrices that will be factored, you must leave room for (and set) 1252 the diagonal entry even if it is zero. 1253 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1254 submatrix (same value is used for all local rows). 1255 - o_nnz - array containing the number of nonzeros in the various rows of the 1256 OFF-DIAGONAL portion of the local submatrix (possibly different for 1257 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1258 structure. The size of this array is equal to the number 1259 of local rows, i.e 'm'. 1260 1261 If the *_nnz parameter is given then the *_nz parameter is ignored 1262 1263 The stored row and column indices begin with zero. 1264 1265 The parallel matrix is partitioned such that the first m0 rows belong to 1266 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1267 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1268 1269 The DIAGONAL portion of the local submatrix of a processor can be defined 1270 as the submatrix which is obtained by extraction the part corresponding to 1271 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1272 first row that belongs to the processor, r2 is the last row belonging to 1273 the this processor, and c1-c2 is range of indices of the local part of a 1274 vector suitable for applying the matrix to. This is an mxn matrix. In the 1275 common case of a square matrix, the row and column ranges are the same and 1276 the DIAGONAL part is also square. The remaining portion of the local 1277 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1278 1279 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1280 1281 You can call MatGetInfo() to get information on how effective the preallocation was; 1282 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1283 You can also run with the option -info and look for messages with the string 1284 malloc in them to see if additional memory allocation was needed. 1285 1286 Example usage: 1287 1288 Consider the following 8x8 matrix with 34 non-zero values, that is 1289 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1290 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1291 as follows: 1292 1293 .vb 1294 1 2 0 | 0 3 0 | 0 4 1295 Proc0 0 5 6 | 7 0 0 | 8 0 1296 9 0 10 | 11 0 0 | 12 0 1297 ------------------------------------- 1298 13 0 14 | 15 16 17 | 0 0 1299 Proc1 0 18 0 | 19 20 21 | 0 0 1300 0 0 0 | 22 23 0 | 24 0 1301 ------------------------------------- 1302 Proc2 25 26 27 | 0 0 28 | 29 0 1303 30 0 0 | 31 32 33 | 0 34 1304 .ve 1305 1306 This can be represented as a collection of submatrices as: 1307 1308 .vb 1309 A B C 1310 D E F 1311 G H I 1312 .ve 1313 1314 Where the submatrices A,B,C are owned by proc0, D,E,F are 1315 owned by proc1, G,H,I are owned by proc2. 1316 1317 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1318 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1319 The 'M','N' parameters are 8,8, and have the same values on all procs. 1320 1321 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1322 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1323 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1324 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1325 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1326 matrix, ans [DF] as another SeqSELL matrix. 1327 1328 When d_nz, o_nz parameters are specified, d_nz storage elements are 1329 allocated for every row of the local diagonal submatrix, and o_nz 1330 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1331 One way to choose d_nz and o_nz is to use the max nonzerors per local 1332 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1333 In this case, the values of d_nz,o_nz are: 1334 .vb 1335 proc0 : dnz = 2, o_nz = 2 1336 proc1 : dnz = 3, o_nz = 2 1337 proc2 : dnz = 1, o_nz = 4 1338 .ve 1339 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1340 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1341 for proc3. i.e we are using 12+15+10=37 storage locations to store 1342 34 values. 1343 1344 When d_nnz, o_nnz parameters are specified, the storage is specified 1345 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1346 In the above case the values for d_nnz,o_nnz are: 1347 .vb 1348 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1349 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1350 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1351 .ve 1352 Here the space allocated is according to nz (or maximum values in the nnz 1353 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1354 1355 Level: intermediate 1356 1357 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreatesell()`, 1358 `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()` 1359 @*/ 1360 PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1361 { 1362 PetscFunctionBegin; 1363 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1364 PetscValidType(B,1); 1365 PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz)); 1366 PetscFunctionReturn(0); 1367 } 1368 1369 /*MC 1370 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1371 based on the sliced Ellpack format 1372 1373 Options Database Keys: 1374 . -mat_type sell - sets the matrix type to "seqsell" during a call to MatSetFromOptions() 1375 1376 Level: beginner 1377 1378 .seealso: `MatCreateSell()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ` 1379 M*/ 1380 1381 /*@C 1382 MatCreateSELL - Creates a sparse parallel matrix in SELL format. 1383 1384 Collective 1385 1386 Input Parameters: 1387 + comm - MPI communicator 1388 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1389 This value should be the same as the local size used in creating the 1390 y vector for the matrix-vector product y = Ax. 1391 . n - This value should be the same as the local size used in creating the 1392 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1393 calculated if N is given) For square matrices n is almost always m. 1394 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 1395 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 1396 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1397 (same value is used for all local rows) 1398 . d_rlen - array containing the number of nonzeros in the various rows of the 1399 DIAGONAL portion of the local submatrix (possibly different for each row) 1400 or NULL, if d_rlenmax is used to specify the nonzero structure. 1401 The size of this array is equal to the number of local rows, i.e 'm'. 1402 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1403 submatrix (same value is used for all local rows). 1404 - o_rlen - array containing the number of nonzeros in the various rows of the 1405 OFF-DIAGONAL portion of the local submatrix (possibly different for 1406 each row) or NULL, if o_rlenmax is used to specify the nonzero 1407 structure. The size of this array is equal to the number 1408 of local rows, i.e 'm'. 1409 1410 Output Parameter: 1411 . A - the matrix 1412 1413 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1414 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1415 [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation] 1416 1417 Notes: 1418 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1419 1420 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1421 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1422 storage requirements for this matrix. 1423 1424 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one 1425 processor than it must be used on all processors that share the object for 1426 that argument. 1427 1428 The user MUST specify either the local or global matrix dimensions 1429 (possibly both). 1430 1431 The parallel matrix is partitioned across processors such that the 1432 first m0 rows belong to process 0, the next m1 rows belong to 1433 process 1, the next m2 rows belong to process 2 etc.. where 1434 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1435 values corresponding to [m x N] submatrix. 1436 1437 The columns are logically partitioned with the n0 columns belonging 1438 to 0th partition, the next n1 columns belonging to the next 1439 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1440 1441 The DIAGONAL portion of the local submatrix on any given processor 1442 is the submatrix corresponding to the rows and columns m,n 1443 corresponding to the given processor. i.e diagonal matrix on 1444 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1445 etc. The remaining portion of the local submatrix [m x (N-n)] 1446 constitute the OFF-DIAGONAL portion. The example below better 1447 illustrates this concept. 1448 1449 For a square global matrix we define each processor's diagonal portion 1450 to be its local rows and the corresponding columns (a square submatrix); 1451 each processor's off-diagonal portion encompasses the remainder of the 1452 local matrix (a rectangular submatrix). 1453 1454 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1455 1456 When calling this routine with a single process communicator, a matrix of 1457 type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this 1458 type of communicator, use the construction mechanism: 1459 MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...); 1460 1461 Options Database Keys: 1462 - -mat_sell_oneindex - Internally use indexing starting at 1 1463 rather than 0. Note that when calling MatSetValues(), 1464 the user still MUST index entries starting at 0! 1465 1466 Example usage: 1467 1468 Consider the following 8x8 matrix with 34 non-zero values, that is 1469 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1470 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1471 as follows: 1472 1473 .vb 1474 1 2 0 | 0 3 0 | 0 4 1475 Proc0 0 5 6 | 7 0 0 | 8 0 1476 9 0 10 | 11 0 0 | 12 0 1477 ------------------------------------- 1478 13 0 14 | 15 16 17 | 0 0 1479 Proc1 0 18 0 | 19 20 21 | 0 0 1480 0 0 0 | 22 23 0 | 24 0 1481 ------------------------------------- 1482 Proc2 25 26 27 | 0 0 28 | 29 0 1483 30 0 0 | 31 32 33 | 0 34 1484 .ve 1485 1486 This can be represented as a collection of submatrices as: 1487 1488 .vb 1489 A B C 1490 D E F 1491 G H I 1492 .ve 1493 1494 Where the submatrices A,B,C are owned by proc0, D,E,F are 1495 owned by proc1, G,H,I are owned by proc2. 1496 1497 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1498 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1499 The 'M','N' parameters are 8,8, and have the same values on all procs. 1500 1501 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1502 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1503 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1504 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1505 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1506 matrix, ans [DF] as another SeqSELL matrix. 1507 1508 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1509 allocated for every row of the local diagonal submatrix, and o_rlenmax 1510 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1511 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1512 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1513 In this case, the values of d_rlenmax,o_rlenmax are: 1514 .vb 1515 proc0 : d_rlenmax = 2, o_rlenmax = 2 1516 proc1 : d_rlenmax = 3, o_rlenmax = 2 1517 proc2 : d_rlenmax = 1, o_rlenmax = 4 1518 .ve 1519 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1520 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1521 for proc3. i.e we are using 12+15+10=37 storage locations to store 1522 34 values. 1523 1524 When d_rlen, o_rlen parameters are specified, the storage is specified 1525 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1526 In the above case the values for d_nnz,o_nnz are: 1527 .vb 1528 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1529 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1530 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1531 .ve 1532 Here the space allocated is still 37 though there are 34 nonzeros because 1533 the allocation is always done according to rlenmax. 1534 1535 Level: intermediate 1536 1537 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MatMPISELLSetPreallocationSELL()`, 1538 `MATMPISELL`, `MatCreateMPISELLWithArrays()` 1539 @*/ 1540 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) 1541 { 1542 PetscMPIInt size; 1543 1544 PetscFunctionBegin; 1545 PetscCall(MatCreate(comm,A)); 1546 PetscCall(MatSetSizes(*A,m,n,M,N)); 1547 PetscCallMPI(MPI_Comm_size(comm,&size)); 1548 if (size > 1) { 1549 PetscCall(MatSetType(*A,MATMPISELL)); 1550 PetscCall(MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen)); 1551 } else { 1552 PetscCall(MatSetType(*A,MATSEQSELL)); 1553 PetscCall(MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen)); 1554 } 1555 PetscFunctionReturn(0); 1556 } 1557 1558 PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 1559 { 1560 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1561 PetscBool flg; 1562 1563 PetscFunctionBegin; 1564 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg)); 1565 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input"); 1566 if (Ad) *Ad = a->A; 1567 if (Ao) *Ao = a->B; 1568 if (colmap) *colmap = a->garray; 1569 PetscFunctionReturn(0); 1570 } 1571 1572 /*@C 1573 MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns 1574 1575 Not Collective 1576 1577 Input Parameters: 1578 + A - the matrix 1579 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1580 - row, col - index sets of rows and columns to extract (or NULL) 1581 1582 Output Parameter: 1583 . A_loc - the local sequential matrix generated 1584 1585 Level: developer 1586 1587 .seealso: `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()` 1588 1589 @*/ 1590 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc) 1591 { 1592 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1593 PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx; 1594 IS isrowa,iscola; 1595 Mat *aloc; 1596 PetscBool match; 1597 1598 PetscFunctionBegin; 1599 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match)); 1600 PetscCheck(match,PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input"); 1601 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0)); 1602 if (!row) { 1603 start = A->rmap->rstart; end = A->rmap->rend; 1604 PetscCall(ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa)); 1605 } else { 1606 isrowa = *row; 1607 } 1608 if (!col) { 1609 start = A->cmap->rstart; 1610 cmap = a->garray; 1611 nzA = a->A->cmap->n; 1612 nzB = a->B->cmap->n; 1613 PetscCall(PetscMalloc1(nzA+nzB, &idx)); 1614 ncols = 0; 1615 for (i=0; i<nzB; i++) { 1616 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1617 else break; 1618 } 1619 imark = i; 1620 for (i=0; i<nzA; i++) idx[ncols++] = start + i; 1621 for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; 1622 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola)); 1623 } else { 1624 iscola = *col; 1625 } 1626 if (scall != MAT_INITIAL_MATRIX) { 1627 PetscCall(PetscMalloc1(1,&aloc)); 1628 aloc[0] = *A_loc; 1629 } 1630 PetscCall(MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc)); 1631 *A_loc = aloc[0]; 1632 PetscCall(PetscFree(aloc)); 1633 if (!row) { 1634 PetscCall(ISDestroy(&isrowa)); 1635 } 1636 if (!col) { 1637 PetscCall(ISDestroy(&iscola)); 1638 } 1639 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0)); 1640 PetscFunctionReturn(0); 1641 } 1642 1643 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1644 1645 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1646 { 1647 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1648 Mat B; 1649 Mat_MPIAIJ *b; 1650 1651 PetscFunctionBegin; 1652 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1653 1654 if (reuse == MAT_REUSE_MATRIX) { 1655 B = *newmat; 1656 } else { 1657 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1658 PetscCall(MatSetType(B,MATMPIAIJ)); 1659 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1660 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1661 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1662 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1663 } 1664 b = (Mat_MPIAIJ*) B->data; 1665 1666 if (reuse == MAT_REUSE_MATRIX) { 1667 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1668 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1669 } else { 1670 PetscCall(MatDestroy(&b->A)); 1671 PetscCall(MatDestroy(&b->B)); 1672 PetscCall(MatDisAssemble_MPISELL(A)); 1673 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1674 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1675 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1676 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1677 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1678 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1679 } 1680 1681 if (reuse == MAT_INPLACE_MATRIX) { 1682 PetscCall(MatHeaderReplace(A,&B)); 1683 } else { 1684 *newmat = B; 1685 } 1686 PetscFunctionReturn(0); 1687 } 1688 1689 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1690 { 1691 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 1692 Mat B; 1693 Mat_MPISELL *b; 1694 1695 PetscFunctionBegin; 1696 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1697 1698 if (reuse == MAT_REUSE_MATRIX) { 1699 B = *newmat; 1700 } else { 1701 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1702 PetscCall(MatSetType(B,MATMPISELL)); 1703 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1704 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1705 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1706 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1707 } 1708 b = (Mat_MPISELL*) B->data; 1709 1710 if (reuse == MAT_REUSE_MATRIX) { 1711 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1712 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1713 } else { 1714 PetscCall(MatDestroy(&b->A)); 1715 PetscCall(MatDestroy(&b->B)); 1716 PetscCall(MatDisAssemble_MPIAIJ(A)); 1717 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1718 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1719 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1720 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1721 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1722 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1723 } 1724 1725 if (reuse == MAT_INPLACE_MATRIX) { 1726 PetscCall(MatHeaderReplace(A,&B)); 1727 } else { 1728 *newmat = B; 1729 } 1730 PetscFunctionReturn(0); 1731 } 1732 1733 PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1734 { 1735 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 1736 Vec bb1=NULL; 1737 1738 PetscFunctionBegin; 1739 if (flag == SOR_APPLY_UPPER) { 1740 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1741 PetscFunctionReturn(0); 1742 } 1743 1744 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) { 1745 PetscCall(VecDuplicate(bb,&bb1)); 1746 } 1747 1748 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1749 if (flag & SOR_ZERO_INITIAL_GUESS) { 1750 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1751 its--; 1752 } 1753 1754 while (its--) { 1755 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1756 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1757 1758 /* update rhs: bb1 = bb - B*x */ 1759 PetscCall(VecScale(mat->lvec,-1.0)); 1760 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1761 1762 /* local sweep */ 1763 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx)); 1764 } 1765 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1766 if (flag & SOR_ZERO_INITIAL_GUESS) { 1767 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1768 its--; 1769 } 1770 while (its--) { 1771 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1772 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1773 1774 /* update rhs: bb1 = bb - B*x */ 1775 PetscCall(VecScale(mat->lvec,-1.0)); 1776 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1777 1778 /* local sweep */ 1779 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx)); 1780 } 1781 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1782 if (flag & SOR_ZERO_INITIAL_GUESS) { 1783 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1784 its--; 1785 } 1786 while (its--) { 1787 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1788 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1789 1790 /* update rhs: bb1 = bb - B*x */ 1791 PetscCall(VecScale(mat->lvec,-1.0)); 1792 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1793 1794 /* local sweep */ 1795 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx)); 1796 } 1797 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported"); 1798 1799 PetscCall(VecDestroy(&bb1)); 1800 1801 matin->factorerrortype = mat->A->factorerrortype; 1802 PetscFunctionReturn(0); 1803 } 1804 1805 /*MC 1806 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1807 1808 Options Database Keys: 1809 . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions() 1810 1811 Level: beginner 1812 1813 .seealso: `MatCreateSELL()` 1814 M*/ 1815 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1816 { 1817 Mat_MPISELL *b; 1818 PetscMPIInt size; 1819 1820 PetscFunctionBegin; 1821 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1822 PetscCall(PetscNewLog(B,&b)); 1823 B->data = (void*)b; 1824 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1825 B->assembled = PETSC_FALSE; 1826 B->insertmode = NOT_SET_VALUES; 1827 b->size = size; 1828 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank)); 1829 /* build cache for off array entries formed */ 1830 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash)); 1831 1832 b->donotstash = PETSC_FALSE; 1833 b->colmap = NULL; 1834 b->garray = NULL; 1835 b->roworiented = PETSC_TRUE; 1836 1837 /* stuff used for matrix vector multiply */ 1838 b->lvec = NULL; 1839 b->Mvctx = NULL; 1840 1841 /* stuff for MatGetRow() */ 1842 b->rowindices = NULL; 1843 b->rowvalues = NULL; 1844 b->getrowactive = PETSC_FALSE; 1845 1846 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL)); 1847 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL)); 1848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL)); 1849 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL)); 1850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ)); 1851 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL)); 1852 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATMPISELL)); 1853 PetscFunctionReturn(0); 1854 } 1855