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) continue; 179 else 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); 180 else { 181 if (mat->was_assembled) { 182 if (!sell->colmap) { 183 PetscCall(MatCreateColmap_MPISELL_Private(mat)); 184 } 185 #if defined(PETSC_USE_CTABLE) 186 PetscCall(PetscTableFind(sell->colmap,in[j]+1,&col)); 187 col--; 188 #else 189 col = sell->colmap[in[j]] - 1; 190 #endif 191 if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) { 192 PetscCall(MatDisAssemble_MPISELL(mat)); 193 col = in[j]; 194 /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */ 195 B = sell->B; 196 b = (Mat_SeqSELL*)B->data; 197 shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 198 cp2 = b->colidx+shift2; 199 vp2 = b->val+shift2; 200 nrow2 = b->rlen[row]; 201 low2 = 0; 202 high2 = nrow2; 203 } else 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]); 204 } else col = in[j]; 205 MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */ 206 } 207 } 208 } else { 209 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]); 210 if (!sell->donotstash) { 211 mat->assembled = PETSC_FALSE; 212 if (roworiented) { 213 PetscCall(MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 214 } else { 215 PetscCall(MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 216 } 217 } 218 } 219 } 220 PetscFunctionReturn(0); 221 } 222 223 PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 224 { 225 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 226 PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend; 227 PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col; 228 229 PetscFunctionBegin; 230 for (i=0; i<m; i++) { 231 if (idxm[i] < 0) continue; /* negative row */ 232 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); 233 if (idxm[i] >= rstart && idxm[i] < rend) { 234 row = idxm[i] - rstart; 235 for (j=0; j<n; j++) { 236 if (idxn[j] < 0) continue; /* negative column */ 237 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); 238 if (idxn[j] >= cstart && idxn[j] < cend) { 239 col = idxn[j] - cstart; 240 PetscCall(MatGetValues(sell->A,1,&row,1,&col,v+i*n+j)); 241 } else { 242 if (!sell->colmap) { 243 PetscCall(MatCreateColmap_MPISELL_Private(mat)); 244 } 245 #if defined(PETSC_USE_CTABLE) 246 PetscCall(PetscTableFind(sell->colmap,idxn[j]+1,&col)); 247 col--; 248 #else 249 col = sell->colmap[idxn[j]] - 1; 250 #endif 251 if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0; 252 else { 253 PetscCall(MatGetValues(sell->B,1,&row,1,&col,v+i*n+j)); 254 } 255 } 256 } 257 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported"); 258 } 259 PetscFunctionReturn(0); 260 } 261 262 extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat,Vec,Vec); 263 264 PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode) 265 { 266 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 267 PetscInt nstash,reallocs; 268 269 PetscFunctionBegin; 270 if (sell->donotstash || mat->nooffprocentries) PetscFunctionReturn(0); 271 272 PetscCall(MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range)); 273 PetscCall(MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs)); 274 PetscCall(PetscInfo(sell->A,"Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n",nstash,reallocs)); 275 PetscFunctionReturn(0); 276 } 277 278 PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode) 279 { 280 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 281 PetscMPIInt n; 282 PetscInt i,flg; 283 PetscInt *row,*col; 284 PetscScalar *val; 285 PetscBool other_disassembled; 286 /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */ 287 PetscFunctionBegin; 288 if (!sell->donotstash && !mat->nooffprocentries) { 289 while (1) { 290 PetscCall(MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg)); 291 if (!flg) break; 292 293 for (i=0; i<n; i++) { /* assemble one by one */ 294 PetscCall(MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode)); 295 } 296 } 297 PetscCall(MatStashScatterEnd_Private(&mat->stash)); 298 } 299 PetscCall(MatAssemblyBegin(sell->A,mode)); 300 PetscCall(MatAssemblyEnd(sell->A,mode)); 301 302 /* 303 determine if any processor has disassembled, if so we must 304 also disassemble ourselves, in order that we may reassemble. 305 */ 306 /* 307 if nonzero structure of submatrix B cannot change then we know that 308 no processor disassembled thus we can skip this stuff 309 */ 310 if (!((Mat_SeqSELL*)sell->B->data)->nonew) { 311 PetscCall(MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat))); 312 PetscCheck(!mat->was_assembled || other_disassembled,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDisAssemble not implemented yet"); 313 } 314 if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) { 315 PetscCall(MatSetUpMultiply_MPISELL(mat)); 316 } 317 /* 318 PetscCall(MatSetOption(sell->B,MAT_USE_INODES,PETSC_FALSE)); 319 */ 320 PetscCall(MatAssemblyBegin(sell->B,mode)); 321 PetscCall(MatAssemblyEnd(sell->B,mode)); 322 PetscCall(PetscFree2(sell->rowvalues,sell->rowindices)); 323 sell->rowvalues = NULL; 324 PetscCall(VecDestroy(&sell->diag)); 325 326 /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */ 327 if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) { 328 PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate; 329 PetscCall(MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat))); 330 } 331 PetscFunctionReturn(0); 332 } 333 334 PetscErrorCode MatZeroEntries_MPISELL(Mat A) 335 { 336 Mat_MPISELL *l=(Mat_MPISELL*)A->data; 337 338 PetscFunctionBegin; 339 PetscCall(MatZeroEntries(l->A)); 340 PetscCall(MatZeroEntries(l->B)); 341 PetscFunctionReturn(0); 342 } 343 344 PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy) 345 { 346 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 347 PetscInt nt; 348 349 PetscFunctionBegin; 350 PetscCall(VecGetLocalSize(xx,&nt)); 351 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); 352 PetscCall(VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 353 PetscCall((*a->A->ops->mult)(a->A,xx,yy)); 354 PetscCall(VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 355 PetscCall((*a->B->ops->multadd)(a->B,a->lvec,yy,yy)); 356 PetscFunctionReturn(0); 357 } 358 359 PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx) 360 { 361 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 362 363 PetscFunctionBegin; 364 PetscCall(MatMultDiagonalBlock(a->A,bb,xx)); 365 PetscFunctionReturn(0); 366 } 367 368 PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz) 369 { 370 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 371 372 PetscFunctionBegin; 373 PetscCall(VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 374 PetscCall((*a->A->ops->multadd)(a->A,xx,yy,zz)); 375 PetscCall(VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 376 PetscCall((*a->B->ops->multadd)(a->B,a->lvec,zz,zz)); 377 PetscFunctionReturn(0); 378 } 379 380 PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy) 381 { 382 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 383 384 PetscFunctionBegin; 385 /* do nondiagonal part */ 386 PetscCall((*a->B->ops->multtranspose)(a->B,xx,a->lvec)); 387 /* do local part */ 388 PetscCall((*a->A->ops->multtranspose)(a->A,xx,yy)); 389 /* add partial results together */ 390 PetscCall(VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE)); 391 PetscCall(VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE)); 392 PetscFunctionReturn(0); 393 } 394 395 PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f) 396 { 397 MPI_Comm comm; 398 Mat_MPISELL *Asell=(Mat_MPISELL*)Amat->data,*Bsell; 399 Mat Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs; 400 IS Me,Notme; 401 PetscInt M,N,first,last,*notme,i; 402 PetscMPIInt size; 403 404 PetscFunctionBegin; 405 /* Easy test: symmetric diagonal block */ 406 Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A; 407 PetscCall(MatIsTranspose(Adia,Bdia,tol,f)); 408 if (!*f) PetscFunctionReturn(0); 409 PetscCall(PetscObjectGetComm((PetscObject)Amat,&comm)); 410 PetscCallMPI(MPI_Comm_size(comm,&size)); 411 if (size == 1) PetscFunctionReturn(0); 412 413 /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */ 414 PetscCall(MatGetSize(Amat,&M,&N)); 415 PetscCall(MatGetOwnershipRange(Amat,&first,&last)); 416 PetscCall(PetscMalloc1(N-last+first,¬me)); 417 for (i=0; i<first; i++) notme[i] = i; 418 for (i=last; i<M; i++) notme[i-last+first] = i; 419 PetscCall(ISCreateGeneral(MPI_COMM_SELF,N-last+first,notme,PETSC_COPY_VALUES,&Notme)); 420 PetscCall(ISCreateStride(MPI_COMM_SELF,last-first,first,1,&Me)); 421 PetscCall(MatCreateSubMatrices(Amat,1,&Me,&Notme,MAT_INITIAL_MATRIX,&Aoffs)); 422 Aoff = Aoffs[0]; 423 PetscCall(MatCreateSubMatrices(Bmat,1,&Notme,&Me,MAT_INITIAL_MATRIX,&Boffs)); 424 Boff = Boffs[0]; 425 PetscCall(MatIsTranspose(Aoff,Boff,tol,f)); 426 PetscCall(MatDestroyMatrices(1,&Aoffs)); 427 PetscCall(MatDestroyMatrices(1,&Boffs)); 428 PetscCall(ISDestroy(&Me)); 429 PetscCall(ISDestroy(&Notme)); 430 PetscCall(PetscFree(notme)); 431 PetscFunctionReturn(0); 432 } 433 434 PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz) 435 { 436 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 437 438 PetscFunctionBegin; 439 /* do nondiagonal part */ 440 PetscCall((*a->B->ops->multtranspose)(a->B,xx,a->lvec)); 441 /* do local part */ 442 PetscCall((*a->A->ops->multtransposeadd)(a->A,xx,yy,zz)); 443 /* add partial results together */ 444 PetscCall(VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE)); 445 PetscCall(VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE)); 446 PetscFunctionReturn(0); 447 } 448 449 /* 450 This only works correctly for square matrices where the subblock A->A is the 451 diagonal block 452 */ 453 PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v) 454 { 455 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 456 457 PetscFunctionBegin; 458 PetscCheck(A->rmap->N == A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block"); 459 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"); 460 PetscCall(MatGetDiagonal(a->A,v)); 461 PetscFunctionReturn(0); 462 } 463 464 PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa) 465 { 466 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 467 468 PetscFunctionBegin; 469 PetscCall(MatScale(a->A,aa)); 470 PetscCall(MatScale(a->B,aa)); 471 PetscFunctionReturn(0); 472 } 473 474 PetscErrorCode MatDestroy_MPISELL(Mat mat) 475 { 476 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 477 478 PetscFunctionBegin; 479 #if defined(PETSC_USE_LOG) 480 PetscLogObjectState((PetscObject)mat,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT,mat->rmap->N,mat->cmap->N); 481 #endif 482 PetscCall(MatStashDestroy_Private(&mat->stash)); 483 PetscCall(VecDestroy(&sell->diag)); 484 PetscCall(MatDestroy(&sell->A)); 485 PetscCall(MatDestroy(&sell->B)); 486 #if defined(PETSC_USE_CTABLE) 487 PetscCall(PetscTableDestroy(&sell->colmap)); 488 #else 489 PetscCall(PetscFree(sell->colmap)); 490 #endif 491 PetscCall(PetscFree(sell->garray)); 492 PetscCall(VecDestroy(&sell->lvec)); 493 PetscCall(VecScatterDestroy(&sell->Mvctx)); 494 PetscCall(PetscFree2(sell->rowvalues,sell->rowindices)); 495 PetscCall(PetscFree(sell->ld)); 496 PetscCall(PetscFree(mat->data)); 497 498 PetscCall(PetscObjectChangeTypeName((PetscObject)mat,NULL)); 499 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL)); 500 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL)); 501 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL)); 502 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL)); 503 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL)); 504 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL)); 505 PetscFunctionReturn(0); 506 } 507 508 #include <petscdraw.h> 509 PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer) 510 { 511 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 512 PetscMPIInt rank=sell->rank,size=sell->size; 513 PetscBool isdraw,iascii,isbinary; 514 PetscViewer sviewer; 515 PetscViewerFormat format; 516 517 PetscFunctionBegin; 518 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 519 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 520 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 521 if (iascii) { 522 PetscCall(PetscViewerGetFormat(viewer,&format)); 523 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 524 MatInfo info; 525 PetscInt *inodes; 526 527 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank)); 528 PetscCall(MatGetInfo(mat,MAT_LOCAL,&info)); 529 PetscCall(MatInodeGetInodeSizes(sell->A,NULL,&inodes,NULL)); 530 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 531 if (!inodes) { 532 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", not using I-node routines\n", 533 rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory)); 534 } else { 535 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", using I-node routines\n", 536 rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory)); 537 } 538 PetscCall(MatGetInfo(sell->A,MAT_LOCAL,&info)); 539 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used)); 540 PetscCall(MatGetInfo(sell->B,MAT_LOCAL,&info)); 541 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used)); 542 PetscCall(PetscViewerFlush(viewer)); 543 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 544 PetscCall(PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n")); 545 PetscCall(VecScatterView(sell->Mvctx,viewer)); 546 PetscFunctionReturn(0); 547 } else if (format == PETSC_VIEWER_ASCII_INFO) { 548 PetscInt inodecount,inodelimit,*inodes; 549 PetscCall(MatInodeGetInodeSizes(sell->A,&inodecount,&inodes,&inodelimit)); 550 if (inodes) { 551 PetscCall(PetscViewerASCIIPrintf(viewer,"using I-node (on process 0) routines: found %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\n",inodecount,inodelimit)); 552 } else { 553 PetscCall(PetscViewerASCIIPrintf(viewer,"not using I-node (on process 0) routines\n")); 554 } 555 PetscFunctionReturn(0); 556 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 557 PetscFunctionReturn(0); 558 } 559 } else if (isbinary) { 560 if (size == 1) { 561 PetscCall(PetscObjectSetName((PetscObject)sell->A,((PetscObject)mat)->name)); 562 PetscCall(MatView(sell->A,viewer)); 563 } else { 564 /* PetscCall(MatView_MPISELL_Binary(mat,viewer)); */ 565 } 566 PetscFunctionReturn(0); 567 } else if (isdraw) { 568 PetscDraw draw; 569 PetscBool isnull; 570 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 571 PetscCall(PetscDrawIsNull(draw,&isnull)); 572 if (isnull) PetscFunctionReturn(0); 573 } 574 575 { 576 /* assemble the entire matrix onto first processor. */ 577 Mat A; 578 Mat_SeqSELL *Aloc; 579 PetscInt M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j; 580 MatScalar *aval; 581 PetscBool isnonzero; 582 583 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat),&A)); 584 if (rank == 0) { 585 PetscCall(MatSetSizes(A,M,N,M,N)); 586 } else { 587 PetscCall(MatSetSizes(A,0,0,M,N)); 588 } 589 /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */ 590 PetscCall(MatSetType(A,MATMPISELL)); 591 PetscCall(MatMPISELLSetPreallocation(A,0,NULL,0,NULL)); 592 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE)); 593 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)A)); 594 595 /* copy over the A part */ 596 Aloc = (Mat_SeqSELL*)sell->A->data; 597 acolidx = Aloc->colidx; aval = Aloc->val; 598 for (i=0; i<Aloc->totalslices; i++) { /* loop over slices */ 599 for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) { 600 isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]); 601 if (isnonzero) { /* check the mask bit */ 602 row = (i<<3)+(j&0x07) + mat->rmap->rstart; /* i<<3 is the starting row of this slice */ 603 col = *acolidx + mat->rmap->rstart; 604 PetscCall(MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES)); 605 } 606 aval++; acolidx++; 607 } 608 } 609 610 /* copy over the B part */ 611 Aloc = (Mat_SeqSELL*)sell->B->data; 612 acolidx = Aloc->colidx; aval = Aloc->val; 613 for (i=0; i<Aloc->totalslices; i++) { 614 for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) { 615 isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]); 616 if (isnonzero) { 617 row = (i<<3)+(j&0x07) + mat->rmap->rstart; 618 col = sell->garray[*acolidx]; 619 PetscCall(MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES)); 620 } 621 aval++; acolidx++; 622 } 623 } 624 625 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 626 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 627 /* 628 Everyone has to call to draw the matrix since the graphics waits are 629 synchronized across all processors that share the PetscDraw object 630 */ 631 PetscCall(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 632 if (rank == 0) { 633 PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISELL*)(A->data))->A,((PetscObject)mat)->name)); 634 PetscCall(MatView_SeqSELL(((Mat_MPISELL*)(A->data))->A,sviewer)); 635 } 636 PetscCall(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 637 PetscCall(PetscViewerFlush(viewer)); 638 PetscCall(MatDestroy(&A)); 639 } 640 PetscFunctionReturn(0); 641 } 642 643 PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer) 644 { 645 PetscBool iascii,isdraw,issocket,isbinary; 646 647 PetscFunctionBegin; 648 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 649 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 650 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 651 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket)); 652 if (iascii || isdraw || isbinary || issocket) { 653 PetscCall(MatView_MPISELL_ASCIIorDraworSocket(mat,viewer)); 654 } 655 PetscFunctionReturn(0); 656 } 657 658 PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[]) 659 { 660 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 661 662 PetscFunctionBegin; 663 PetscCall(MatGetSize(sell->B,NULL,nghosts)); 664 if (ghosts) *ghosts = sell->garray; 665 PetscFunctionReturn(0); 666 } 667 668 PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info) 669 { 670 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 671 Mat A=mat->A,B=mat->B; 672 PetscLogDouble isend[5],irecv[5]; 673 674 PetscFunctionBegin; 675 info->block_size = 1.0; 676 PetscCall(MatGetInfo(A,MAT_LOCAL,info)); 677 678 isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded; 679 isend[3] = info->memory; isend[4] = info->mallocs; 680 681 PetscCall(MatGetInfo(B,MAT_LOCAL,info)); 682 683 isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded; 684 isend[3] += info->memory; isend[4] += info->mallocs; 685 if (flag == MAT_LOCAL) { 686 info->nz_used = isend[0]; 687 info->nz_allocated = isend[1]; 688 info->nz_unneeded = isend[2]; 689 info->memory = isend[3]; 690 info->mallocs = isend[4]; 691 } else if (flag == MAT_GLOBAL_MAX) { 692 PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)matin))); 693 694 info->nz_used = irecv[0]; 695 info->nz_allocated = irecv[1]; 696 info->nz_unneeded = irecv[2]; 697 info->memory = irecv[3]; 698 info->mallocs = irecv[4]; 699 } else if (flag == MAT_GLOBAL_SUM) { 700 PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)matin))); 701 702 info->nz_used = irecv[0]; 703 info->nz_allocated = irecv[1]; 704 info->nz_unneeded = irecv[2]; 705 info->memory = irecv[3]; 706 info->mallocs = irecv[4]; 707 } 708 info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */ 709 info->fill_ratio_needed = 0; 710 info->factor_mallocs = 0; 711 PetscFunctionReturn(0); 712 } 713 714 PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg) 715 { 716 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 717 718 PetscFunctionBegin; 719 switch (op) { 720 case MAT_NEW_NONZERO_LOCATIONS: 721 case MAT_NEW_NONZERO_ALLOCATION_ERR: 722 case MAT_UNUSED_NONZERO_LOCATION_ERR: 723 case MAT_KEEP_NONZERO_PATTERN: 724 case MAT_NEW_NONZERO_LOCATION_ERR: 725 case MAT_USE_INODES: 726 case MAT_IGNORE_ZERO_ENTRIES: 727 MatCheckPreallocated(A,1); 728 PetscCall(MatSetOption(a->A,op,flg)); 729 PetscCall(MatSetOption(a->B,op,flg)); 730 break; 731 case MAT_ROW_ORIENTED: 732 MatCheckPreallocated(A,1); 733 a->roworiented = flg; 734 735 PetscCall(MatSetOption(a->A,op,flg)); 736 PetscCall(MatSetOption(a->B,op,flg)); 737 break; 738 case MAT_FORCE_DIAGONAL_ENTRIES: 739 case MAT_SORTED_FULL: 740 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 741 break; 742 case MAT_IGNORE_OFF_PROC_ENTRIES: 743 a->donotstash = flg; 744 break; 745 case MAT_SPD: 746 A->spd_set = PETSC_TRUE; 747 A->spd = flg; 748 if (flg) { 749 A->symmetric = PETSC_TRUE; 750 A->structurally_symmetric = PETSC_TRUE; 751 A->symmetric_set = PETSC_TRUE; 752 A->structurally_symmetric_set = PETSC_TRUE; 753 } 754 break; 755 case MAT_SYMMETRIC: 756 MatCheckPreallocated(A,1); 757 PetscCall(MatSetOption(a->A,op,flg)); 758 break; 759 case MAT_STRUCTURALLY_SYMMETRIC: 760 MatCheckPreallocated(A,1); 761 PetscCall(MatSetOption(a->A,op,flg)); 762 break; 763 case MAT_HERMITIAN: 764 MatCheckPreallocated(A,1); 765 PetscCall(MatSetOption(a->A,op,flg)); 766 break; 767 case MAT_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 }; 1116 1117 /* ----------------------------------------------------------------------------------------*/ 1118 1119 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1120 { 1121 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1122 1123 PetscFunctionBegin; 1124 PetscCall(MatStoreValues(sell->A)); 1125 PetscCall(MatStoreValues(sell->B)); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1130 { 1131 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1132 1133 PetscFunctionBegin; 1134 PetscCall(MatRetrieveValues(sell->A)); 1135 PetscCall(MatRetrieveValues(sell->B)); 1136 PetscFunctionReturn(0); 1137 } 1138 1139 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[]) 1140 { 1141 Mat_MPISELL *b; 1142 1143 PetscFunctionBegin; 1144 PetscCall(PetscLayoutSetUp(B->rmap)); 1145 PetscCall(PetscLayoutSetUp(B->cmap)); 1146 b = (Mat_MPISELL*)B->data; 1147 1148 if (!B->preallocated) { 1149 /* Explicitly create 2 MATSEQSELL matrices. */ 1150 PetscCall(MatCreate(PETSC_COMM_SELF,&b->A)); 1151 PetscCall(MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n)); 1152 PetscCall(MatSetBlockSizesFromMats(b->A,B,B)); 1153 PetscCall(MatSetType(b->A,MATSEQSELL)); 1154 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->A)); 1155 PetscCall(MatCreate(PETSC_COMM_SELF,&b->B)); 1156 PetscCall(MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N)); 1157 PetscCall(MatSetBlockSizesFromMats(b->B,B,B)); 1158 PetscCall(MatSetType(b->B,MATSEQSELL)); 1159 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->B)); 1160 } 1161 1162 PetscCall(MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen)); 1163 PetscCall(MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen)); 1164 B->preallocated = PETSC_TRUE; 1165 B->was_assembled = PETSC_FALSE; 1166 /* 1167 critical for MatAssemblyEnd to work. 1168 MatAssemblyBegin checks it to set up was_assembled 1169 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1170 */ 1171 B->assembled = PETSC_FALSE; 1172 PetscFunctionReturn(0); 1173 } 1174 1175 PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 1176 { 1177 Mat mat; 1178 Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data; 1179 1180 PetscFunctionBegin; 1181 *newmat = NULL; 1182 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin),&mat)); 1183 PetscCall(MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N)); 1184 PetscCall(MatSetBlockSizesFromMats(mat,matin,matin)); 1185 PetscCall(MatSetType(mat,((PetscObject)matin)->type_name)); 1186 a = (Mat_MPISELL*)mat->data; 1187 1188 mat->factortype = matin->factortype; 1189 mat->assembled = PETSC_TRUE; 1190 mat->insertmode = NOT_SET_VALUES; 1191 mat->preallocated = PETSC_TRUE; 1192 1193 a->size = oldmat->size; 1194 a->rank = oldmat->rank; 1195 a->donotstash = oldmat->donotstash; 1196 a->roworiented = oldmat->roworiented; 1197 a->rowindices = NULL; 1198 a->rowvalues = NULL; 1199 a->getrowactive = PETSC_FALSE; 1200 1201 PetscCall(PetscLayoutReference(matin->rmap,&mat->rmap)); 1202 PetscCall(PetscLayoutReference(matin->cmap,&mat->cmap)); 1203 1204 if (oldmat->colmap) { 1205 #if defined(PETSC_USE_CTABLE) 1206 PetscCall(PetscTableCreateCopy(oldmat->colmap,&a->colmap)); 1207 #else 1208 PetscCall(PetscMalloc1(mat->cmap->N,&a->colmap)); 1209 PetscCall(PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt))); 1210 PetscCall(PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N)); 1211 #endif 1212 } else a->colmap = NULL; 1213 if (oldmat->garray) { 1214 PetscInt len; 1215 len = oldmat->B->cmap->n; 1216 PetscCall(PetscMalloc1(len+1,&a->garray)); 1217 PetscCall(PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt))); 1218 if (len) PetscCall(PetscArraycpy(a->garray,oldmat->garray,len)); 1219 } else a->garray = NULL; 1220 1221 PetscCall(VecDuplicate(oldmat->lvec,&a->lvec)); 1222 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec)); 1223 PetscCall(VecScatterCopy(oldmat->Mvctx,&a->Mvctx)); 1224 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx)); 1225 PetscCall(MatDuplicate(oldmat->A,cpvalues,&a->A)); 1226 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A)); 1227 PetscCall(MatDuplicate(oldmat->B,cpvalues,&a->B)); 1228 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B)); 1229 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist)); 1230 *newmat = mat; 1231 PetscFunctionReturn(0); 1232 } 1233 1234 /*@C 1235 MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format. 1236 For good matrix assembly performance the user should preallocate the matrix storage by 1237 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1238 1239 Collective 1240 1241 Input Parameters: 1242 + B - the matrix 1243 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1244 (same value is used for all local rows) 1245 . d_nnz - array containing the number of nonzeros in the various rows of the 1246 DIAGONAL portion of the local submatrix (possibly different for each row) 1247 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1248 The size of this array is equal to the number of local rows, i.e 'm'. 1249 For matrices that will be factored, you must leave room for (and set) 1250 the diagonal entry even if it is zero. 1251 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1252 submatrix (same value is used for all local rows). 1253 - o_nnz - array containing the number of nonzeros in the various rows of the 1254 OFF-DIAGONAL portion of the local submatrix (possibly different for 1255 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1256 structure. The size of this array is equal to the number 1257 of local rows, i.e 'm'. 1258 1259 If the *_nnz parameter is given then the *_nz parameter is ignored 1260 1261 The stored row and column indices begin with zero. 1262 1263 The parallel matrix is partitioned such that the first m0 rows belong to 1264 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1265 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1266 1267 The DIAGONAL portion of the local submatrix of a processor can be defined 1268 as the submatrix which is obtained by extraction the part corresponding to 1269 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1270 first row that belongs to the processor, r2 is the last row belonging to 1271 the this processor, and c1-c2 is range of indices of the local part of a 1272 vector suitable for applying the matrix to. This is an mxn matrix. In the 1273 common case of a square matrix, the row and column ranges are the same and 1274 the DIAGONAL part is also square. The remaining portion of the local 1275 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1276 1277 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1278 1279 You can call MatGetInfo() to get information on how effective the preallocation was; 1280 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1281 You can also run with the option -info and look for messages with the string 1282 malloc in them to see if additional memory allocation was needed. 1283 1284 Example usage: 1285 1286 Consider the following 8x8 matrix with 34 non-zero values, that is 1287 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1288 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1289 as follows: 1290 1291 .vb 1292 1 2 0 | 0 3 0 | 0 4 1293 Proc0 0 5 6 | 7 0 0 | 8 0 1294 9 0 10 | 11 0 0 | 12 0 1295 ------------------------------------- 1296 13 0 14 | 15 16 17 | 0 0 1297 Proc1 0 18 0 | 19 20 21 | 0 0 1298 0 0 0 | 22 23 0 | 24 0 1299 ------------------------------------- 1300 Proc2 25 26 27 | 0 0 28 | 29 0 1301 30 0 0 | 31 32 33 | 0 34 1302 .ve 1303 1304 This can be represented as a collection of submatrices as: 1305 1306 .vb 1307 A B C 1308 D E F 1309 G H I 1310 .ve 1311 1312 Where the submatrices A,B,C are owned by proc0, D,E,F are 1313 owned by proc1, G,H,I are owned by proc2. 1314 1315 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1316 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1317 The 'M','N' parameters are 8,8, and have the same values on all procs. 1318 1319 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1320 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1321 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1322 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1323 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1324 matrix, ans [DF] as another SeqSELL matrix. 1325 1326 When d_nz, o_nz parameters are specified, d_nz storage elements are 1327 allocated for every row of the local diagonal submatrix, and o_nz 1328 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1329 One way to choose d_nz and o_nz is to use the max nonzerors per local 1330 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1331 In this case, the values of d_nz,o_nz are: 1332 .vb 1333 proc0 : dnz = 2, o_nz = 2 1334 proc1 : dnz = 3, o_nz = 2 1335 proc2 : dnz = 1, o_nz = 4 1336 .ve 1337 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1338 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1339 for proc3. i.e we are using 12+15+10=37 storage locations to store 1340 34 values. 1341 1342 When d_nnz, o_nnz parameters are specified, the storage is specified 1343 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1344 In the above case the values for d_nnz,o_nnz are: 1345 .vb 1346 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1347 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1348 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1349 .ve 1350 Here the space allocated is according to nz (or maximum values in the nnz 1351 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1352 1353 Level: intermediate 1354 1355 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreatesell()`, 1356 `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()` 1357 @*/ 1358 PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1359 { 1360 PetscFunctionBegin; 1361 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1362 PetscValidType(B,1); 1363 PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz)); 1364 PetscFunctionReturn(0); 1365 } 1366 1367 /*MC 1368 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1369 based on the sliced Ellpack format 1370 1371 Options Database Keys: 1372 . -mat_type sell - sets the matrix type to "seqsell" during a call to MatSetFromOptions() 1373 1374 Level: beginner 1375 1376 .seealso: `MatCreateSell()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ` 1377 M*/ 1378 1379 /*@C 1380 MatCreateSELL - Creates a sparse parallel matrix in SELL format. 1381 1382 Collective 1383 1384 Input Parameters: 1385 + comm - MPI communicator 1386 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1387 This value should be the same as the local size used in creating the 1388 y vector for the matrix-vector product y = Ax. 1389 . n - This value should be the same as the local size used in creating the 1390 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1391 calculated if N is given) For square matrices n is almost always m. 1392 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 1393 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 1394 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1395 (same value is used for all local rows) 1396 . d_rlen - array containing the number of nonzeros in the various rows of the 1397 DIAGONAL portion of the local submatrix (possibly different for each row) 1398 or NULL, if d_rlenmax is used to specify the nonzero structure. 1399 The size of this array is equal to the number of local rows, i.e 'm'. 1400 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1401 submatrix (same value is used for all local rows). 1402 - o_rlen - array containing the number of nonzeros in the various rows of the 1403 OFF-DIAGONAL portion of the local submatrix (possibly different for 1404 each row) or NULL, if o_rlenmax is used to specify the nonzero 1405 structure. The size of this array is equal to the number 1406 of local rows, i.e 'm'. 1407 1408 Output Parameter: 1409 . A - the matrix 1410 1411 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1412 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1413 [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation] 1414 1415 Notes: 1416 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1417 1418 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1419 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1420 storage requirements for this matrix. 1421 1422 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one 1423 processor than it must be used on all processors that share the object for 1424 that argument. 1425 1426 The user MUST specify either the local or global matrix dimensions 1427 (possibly both). 1428 1429 The parallel matrix is partitioned across processors such that the 1430 first m0 rows belong to process 0, the next m1 rows belong to 1431 process 1, the next m2 rows belong to process 2 etc.. where 1432 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1433 values corresponding to [m x N] submatrix. 1434 1435 The columns are logically partitioned with the n0 columns belonging 1436 to 0th partition, the next n1 columns belonging to the next 1437 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1438 1439 The DIAGONAL portion of the local submatrix on any given processor 1440 is the submatrix corresponding to the rows and columns m,n 1441 corresponding to the given processor. i.e diagonal matrix on 1442 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1443 etc. The remaining portion of the local submatrix [m x (N-n)] 1444 constitute the OFF-DIAGONAL portion. The example below better 1445 illustrates this concept. 1446 1447 For a square global matrix we define each processor's diagonal portion 1448 to be its local rows and the corresponding columns (a square submatrix); 1449 each processor's off-diagonal portion encompasses the remainder of the 1450 local matrix (a rectangular submatrix). 1451 1452 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1453 1454 When calling this routine with a single process communicator, a matrix of 1455 type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this 1456 type of communicator, use the construction mechanism: 1457 MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...); 1458 1459 Options Database Keys: 1460 - -mat_sell_oneindex - Internally use indexing starting at 1 1461 rather than 0. Note that when calling MatSetValues(), 1462 the user still MUST index entries starting at 0! 1463 1464 Example usage: 1465 1466 Consider the following 8x8 matrix with 34 non-zero values, that is 1467 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1468 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1469 as follows: 1470 1471 .vb 1472 1 2 0 | 0 3 0 | 0 4 1473 Proc0 0 5 6 | 7 0 0 | 8 0 1474 9 0 10 | 11 0 0 | 12 0 1475 ------------------------------------- 1476 13 0 14 | 15 16 17 | 0 0 1477 Proc1 0 18 0 | 19 20 21 | 0 0 1478 0 0 0 | 22 23 0 | 24 0 1479 ------------------------------------- 1480 Proc2 25 26 27 | 0 0 28 | 29 0 1481 30 0 0 | 31 32 33 | 0 34 1482 .ve 1483 1484 This can be represented as a collection of submatrices as: 1485 1486 .vb 1487 A B C 1488 D E F 1489 G H I 1490 .ve 1491 1492 Where the submatrices A,B,C are owned by proc0, D,E,F are 1493 owned by proc1, G,H,I are owned by proc2. 1494 1495 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1496 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1497 The 'M','N' parameters are 8,8, and have the same values on all procs. 1498 1499 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1500 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1501 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1502 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1503 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1504 matrix, ans [DF] as another SeqSELL matrix. 1505 1506 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1507 allocated for every row of the local diagonal submatrix, and o_rlenmax 1508 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1509 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1510 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1511 In this case, the values of d_rlenmax,o_rlenmax are: 1512 .vb 1513 proc0 : d_rlenmax = 2, o_rlenmax = 2 1514 proc1 : d_rlenmax = 3, o_rlenmax = 2 1515 proc2 : d_rlenmax = 1, o_rlenmax = 4 1516 .ve 1517 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1518 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1519 for proc3. i.e we are using 12+15+10=37 storage locations to store 1520 34 values. 1521 1522 When d_rlen, o_rlen parameters are specified, the storage is specified 1523 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1524 In the above case the values for d_nnz,o_nnz are: 1525 .vb 1526 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1527 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1528 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1529 .ve 1530 Here the space allocated is still 37 though there are 34 nonzeros because 1531 the allocation is always done according to rlenmax. 1532 1533 Level: intermediate 1534 1535 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MatMPISELLSetPreallocationSELL()`, 1536 `MATMPISELL`, `MatCreateMPISELLWithArrays()` 1537 @*/ 1538 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) 1539 { 1540 PetscMPIInt size; 1541 1542 PetscFunctionBegin; 1543 PetscCall(MatCreate(comm,A)); 1544 PetscCall(MatSetSizes(*A,m,n,M,N)); 1545 PetscCallMPI(MPI_Comm_size(comm,&size)); 1546 if (size > 1) { 1547 PetscCall(MatSetType(*A,MATMPISELL)); 1548 PetscCall(MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen)); 1549 } else { 1550 PetscCall(MatSetType(*A,MATSEQSELL)); 1551 PetscCall(MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen)); 1552 } 1553 PetscFunctionReturn(0); 1554 } 1555 1556 PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 1557 { 1558 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1559 PetscBool flg; 1560 1561 PetscFunctionBegin; 1562 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg)); 1563 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input"); 1564 if (Ad) *Ad = a->A; 1565 if (Ao) *Ao = a->B; 1566 if (colmap) *colmap = a->garray; 1567 PetscFunctionReturn(0); 1568 } 1569 1570 /*@C 1571 MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns 1572 1573 Not Collective 1574 1575 Input Parameters: 1576 + A - the matrix 1577 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1578 - row, col - index sets of rows and columns to extract (or NULL) 1579 1580 Output Parameter: 1581 . A_loc - the local sequential matrix generated 1582 1583 Level: developer 1584 1585 .seealso: `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()` 1586 1587 @*/ 1588 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc) 1589 { 1590 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1591 PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx; 1592 IS isrowa,iscola; 1593 Mat *aloc; 1594 PetscBool match; 1595 1596 PetscFunctionBegin; 1597 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match)); 1598 PetscCheck(match,PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input"); 1599 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0)); 1600 if (!row) { 1601 start = A->rmap->rstart; end = A->rmap->rend; 1602 PetscCall(ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa)); 1603 } else { 1604 isrowa = *row; 1605 } 1606 if (!col) { 1607 start = A->cmap->rstart; 1608 cmap = a->garray; 1609 nzA = a->A->cmap->n; 1610 nzB = a->B->cmap->n; 1611 PetscCall(PetscMalloc1(nzA+nzB, &idx)); 1612 ncols = 0; 1613 for (i=0; i<nzB; i++) { 1614 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1615 else break; 1616 } 1617 imark = i; 1618 for (i=0; i<nzA; i++) idx[ncols++] = start + i; 1619 for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; 1620 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola)); 1621 } else { 1622 iscola = *col; 1623 } 1624 if (scall != MAT_INITIAL_MATRIX) { 1625 PetscCall(PetscMalloc1(1,&aloc)); 1626 aloc[0] = *A_loc; 1627 } 1628 PetscCall(MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc)); 1629 *A_loc = aloc[0]; 1630 PetscCall(PetscFree(aloc)); 1631 if (!row) { 1632 PetscCall(ISDestroy(&isrowa)); 1633 } 1634 if (!col) { 1635 PetscCall(ISDestroy(&iscola)); 1636 } 1637 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0)); 1638 PetscFunctionReturn(0); 1639 } 1640 1641 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1642 1643 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1644 { 1645 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1646 Mat B; 1647 Mat_MPIAIJ *b; 1648 1649 PetscFunctionBegin; 1650 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1651 1652 if (reuse == MAT_REUSE_MATRIX) { 1653 B = *newmat; 1654 } else { 1655 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1656 PetscCall(MatSetType(B,MATMPIAIJ)); 1657 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1658 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1659 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1660 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1661 } 1662 b = (Mat_MPIAIJ*) B->data; 1663 1664 if (reuse == MAT_REUSE_MATRIX) { 1665 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1666 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1667 } else { 1668 PetscCall(MatDestroy(&b->A)); 1669 PetscCall(MatDestroy(&b->B)); 1670 PetscCall(MatDisAssemble_MPISELL(A)); 1671 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1672 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1673 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1674 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1675 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1676 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1677 } 1678 1679 if (reuse == MAT_INPLACE_MATRIX) { 1680 PetscCall(MatHeaderReplace(A,&B)); 1681 } else { 1682 *newmat = B; 1683 } 1684 PetscFunctionReturn(0); 1685 } 1686 1687 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1688 { 1689 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 1690 Mat B; 1691 Mat_MPISELL *b; 1692 1693 PetscFunctionBegin; 1694 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1695 1696 if (reuse == MAT_REUSE_MATRIX) { 1697 B = *newmat; 1698 } else { 1699 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1700 PetscCall(MatSetType(B,MATMPISELL)); 1701 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1702 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1703 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1704 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1705 } 1706 b = (Mat_MPISELL*) B->data; 1707 1708 if (reuse == MAT_REUSE_MATRIX) { 1709 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1710 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1711 } else { 1712 PetscCall(MatDestroy(&b->A)); 1713 PetscCall(MatDestroy(&b->B)); 1714 PetscCall(MatDisAssemble_MPIAIJ(A)); 1715 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1716 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1717 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1718 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1719 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1720 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1721 } 1722 1723 if (reuse == MAT_INPLACE_MATRIX) { 1724 PetscCall(MatHeaderReplace(A,&B)); 1725 } else { 1726 *newmat = B; 1727 } 1728 PetscFunctionReturn(0); 1729 } 1730 1731 PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1732 { 1733 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 1734 Vec bb1=NULL; 1735 1736 PetscFunctionBegin; 1737 if (flag == SOR_APPLY_UPPER) { 1738 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1739 PetscFunctionReturn(0); 1740 } 1741 1742 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) { 1743 PetscCall(VecDuplicate(bb,&bb1)); 1744 } 1745 1746 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1747 if (flag & SOR_ZERO_INITIAL_GUESS) { 1748 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1749 its--; 1750 } 1751 1752 while (its--) { 1753 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1754 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1755 1756 /* update rhs: bb1 = bb - B*x */ 1757 PetscCall(VecScale(mat->lvec,-1.0)); 1758 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1759 1760 /* local sweep */ 1761 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx)); 1762 } 1763 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1764 if (flag & SOR_ZERO_INITIAL_GUESS) { 1765 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1766 its--; 1767 } 1768 while (its--) { 1769 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1770 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1771 1772 /* update rhs: bb1 = bb - B*x */ 1773 PetscCall(VecScale(mat->lvec,-1.0)); 1774 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1775 1776 /* local sweep */ 1777 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx)); 1778 } 1779 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1780 if (flag & SOR_ZERO_INITIAL_GUESS) { 1781 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1782 its--; 1783 } 1784 while (its--) { 1785 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1786 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1787 1788 /* update rhs: bb1 = bb - B*x */ 1789 PetscCall(VecScale(mat->lvec,-1.0)); 1790 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1791 1792 /* local sweep */ 1793 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx)); 1794 } 1795 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported"); 1796 1797 PetscCall(VecDestroy(&bb1)); 1798 1799 matin->factorerrortype = mat->A->factorerrortype; 1800 PetscFunctionReturn(0); 1801 } 1802 1803 /*MC 1804 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1805 1806 Options Database Keys: 1807 . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions() 1808 1809 Level: beginner 1810 1811 .seealso: `MatCreateSELL()` 1812 M*/ 1813 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1814 { 1815 Mat_MPISELL *b; 1816 PetscMPIInt size; 1817 1818 PetscFunctionBegin; 1819 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1820 PetscCall(PetscNewLog(B,&b)); 1821 B->data = (void*)b; 1822 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1823 B->assembled = PETSC_FALSE; 1824 B->insertmode = NOT_SET_VALUES; 1825 b->size = size; 1826 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank)); 1827 /* build cache for off array entries formed */ 1828 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash)); 1829 1830 b->donotstash = PETSC_FALSE; 1831 b->colmap = NULL; 1832 b->garray = NULL; 1833 b->roworiented = PETSC_TRUE; 1834 1835 /* stuff used for matrix vector multiply */ 1836 b->lvec = NULL; 1837 b->Mvctx = NULL; 1838 1839 /* stuff for MatGetRow() */ 1840 b->rowindices = NULL; 1841 b->rowvalues = NULL; 1842 b->getrowactive = PETSC_FALSE; 1843 1844 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL)); 1845 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL)); 1846 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL)); 1847 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL)); 1848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ)); 1849 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL)); 1850 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATMPISELL)); 1851 PetscFunctionReturn(0); 1852 } 1853