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