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