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