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