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