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