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