1 #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/ 2 #include <../src/mat/impls/sell/mpi/mpisell.h> /*I "petscmat.h" I*/ 3 #include <petsc/private/vecimpl.h> 4 #include <petsc/private/isimpl.h> 5 #include <petscblaslapack.h> 6 #include <petscsf.h> 7 8 /*MC 9 MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices. 10 11 This matrix type is identical to MATSEQSELL when constructed with a single process communicator, 12 and MATMPISELL otherwise. As a result, for single process communicators, 13 MatSeqSELLSetPreallocation is supported, and similarly MatMPISELLSetPreallocation is supported 14 for communicators controlling multiple processes. It is recommended that you call both of 15 the above preallocation routines for simplicity. 16 17 Options Database Keys: 18 . -mat_type sell - sets the matrix type to "sell" during a call to MatSetFromOptions() 19 20 Level: beginner 21 22 .seealso: MatCreateSELL(), MatCreateSeqSELL(), MATSEQSELL, MATMPISELL 23 M*/ 24 25 PetscErrorCode MatDiagonalSet_MPISELL(Mat Y,Vec D,InsertMode is) 26 { 27 PetscErrorCode ierr; 28 Mat_MPISELL *sell=(Mat_MPISELL*)Y->data; 29 30 PetscFunctionBegin; 31 if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) { 32 ierr = MatDiagonalSet(sell->A,D,is);CHKERRQ(ierr); 33 } else { 34 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 35 } 36 PetscFunctionReturn(0); 37 } 38 39 /* 40 Local utility routine that creates a mapping from the global column 41 number to the local number in the off-diagonal part of the local 42 storage of the matrix. When PETSC_USE_CTABLE is used this is scalable at 43 a slightly higher hash table cost; without it it is not scalable (each processor 44 has an order N integer array but is fast to acess. 45 */ 46 PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat) 47 { 48 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 49 PetscErrorCode ierr; 50 PetscInt n=sell->B->cmap->n,i; 51 52 PetscFunctionBegin; 53 if (!sell->garray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MPISELL Matrix was assembled but is missing garray"); 54 #if defined(PETSC_USE_CTABLE) 55 ierr = PetscTableCreate(n,mat->cmap->N+1,&sell->colmap);CHKERRQ(ierr); 56 for (i=0; i<n; i++) { 57 ierr = PetscTableAdd(sell->colmap,sell->garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr); 58 } 59 #else 60 ierr = PetscCalloc1(mat->cmap->N+1,&sell->colmap);CHKERRQ(ierr); 61 ierr = PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N+1)*sizeof(PetscInt));CHKERRQ(ierr); 62 for (i=0; i<n; i++) sell->colmap[sell->garray[i]] = i+1; 63 #endif 64 PetscFunctionReturn(0); 65 } 66 67 #define MatSetValues_SeqSELL_A_Private(row,col,value,addv,orow,ocol) \ 68 { \ 69 if (col <= lastcol1) low1 = 0; \ 70 else high1 = nrow1; \ 71 lastcol1 = col; \ 72 while (high1-low1 > 5) { \ 73 t = (low1+high1)/2; \ 74 if (*(cp1+8*t) > col) high1 = t; \ 75 else low1 = t; \ 76 } \ 77 for (_i=low1; _i<high1; _i++) { \ 78 if (*(cp1+8*_i) > col) break; \ 79 if (*(cp1+8*_i) == col) { \ 80 if (addv == ADD_VALUES) *(vp1+8*_i) += value; \ 81 else *(vp1+8*_i) = value; \ 82 goto a_noinsert; \ 83 } \ 84 } \ 85 if (value == 0.0 && ignorezeroentries) {low1 = 0; high1 = nrow1;goto a_noinsert;} \ 86 if (nonew == 1) {low1 = 0; high1 = nrow1; goto a_noinsert;} \ 87 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", orow, ocol); \ 88 MatSeqXSELLReallocateSELL(A,am,1,nrow1,a->sliidx,row/8,row,col,a->colidx,a->val,cp1,vp1,nonew,MatScalar); \ 89 /* shift up all the later entries in this row */ \ 90 for (ii=nrow1-1; ii>=_i; ii--) { \ 91 *(cp1+8*(ii+1)) = *(cp1+8*ii); \ 92 *(vp1+8*(ii+1)) = *(vp1+8*ii); \ 93 } \ 94 *(cp1+8*_i) = col; \ 95 *(vp1+8*_i) = value; \ 96 a->nz++; nrow1++; A->nonzerostate++; \ 97 a_noinsert: ; \ 98 a->rlen[row] = nrow1; \ 99 } 100 101 #define MatSetValues_SeqSELL_B_Private(row,col,value,addv,orow,ocol) \ 102 { \ 103 if (col <= lastcol2) low2 = 0; \ 104 else high2 = nrow2; \ 105 lastcol2 = col; \ 106 while (high2-low2 > 5) { \ 107 t = (low2+high2)/2; \ 108 if (*(cp2+8*t) > col) high2 = t; \ 109 else low2 = t; \ 110 } \ 111 for (_i=low2; _i<high2; _i++) { \ 112 if (*(cp2+8*_i) > col) break; \ 113 if (*(cp2+8*_i) == col) { \ 114 if (addv == ADD_VALUES) *(vp2+8*_i) += value; \ 115 else *(vp2+8*_i) = value; \ 116 goto b_noinsert; \ 117 } \ 118 } \ 119 if (value == 0.0 && ignorezeroentries) {low2 = 0; high2 = nrow2; goto b_noinsert;} \ 120 if (nonew == 1) {low2 = 0; high2 = nrow2; goto b_noinsert;} \ 121 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", orow, ocol); \ 122 MatSeqXSELLReallocateSELL(B,bm,1,nrow2,b->sliidx,row/8,row,col,b->colidx,b->val,cp2,vp2,nonew,MatScalar); \ 123 /* shift up all the later entries in this row */ \ 124 for (ii=nrow2-1; ii>=_i; ii--) { \ 125 *(cp2+8*(ii+1)) = *(cp2+8*ii); \ 126 *(vp2+8*(ii+1)) = *(vp2+8*ii); \ 127 } \ 128 *(cp2+8*_i) = col; \ 129 *(vp2+8*_i) = value; \ 130 b->nz++; nrow2++; B->nonzerostate++; \ 131 b_noinsert: ; \ 132 b->rlen[row] = nrow2; \ 133 } 134 135 PetscErrorCode MatSetValues_MPISELL(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv) 136 { 137 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 138 PetscScalar value; 139 PetscErrorCode ierr; 140 PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend,shift1,shift2; 141 PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col; 142 PetscBool roworiented=sell->roworiented; 143 144 /* Some Variables required in the macro */ 145 Mat A=sell->A; 146 Mat_SeqSELL *a=(Mat_SeqSELL*)A->data; 147 PetscBool ignorezeroentries=a->ignorezeroentries,found; 148 Mat B=sell->B; 149 Mat_SeqSELL *b=(Mat_SeqSELL*)B->data; 150 PetscInt *cp1,*cp2,ii,_i,nrow1,nrow2,low1,high1,low2,high2,t,lastcol1,lastcol2; 151 MatScalar *vp1,*vp2; 152 153 PetscFunctionBegin; 154 for (i=0; i<m; i++) { 155 if (im[i] < 0) continue; 156 if (PetscUnlikelyDebug(im[i] >= mat->rmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],mat->rmap->N-1); 157 if (im[i] >= rstart && im[i] < rend) { 158 row = im[i] - rstart; 159 lastcol1 = -1; 160 shift1 = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 161 cp1 = a->colidx+shift1; 162 vp1 = a->val+shift1; 163 nrow1 = a->rlen[row]; 164 low1 = 0; 165 high1 = nrow1; 166 lastcol2 = -1; 167 shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 168 cp2 = b->colidx+shift2; 169 vp2 = b->val+shift2; 170 nrow2 = b->rlen[row]; 171 low2 = 0; 172 high2 = nrow2; 173 174 for (j=0; j<n; j++) { 175 if (roworiented) value = v[i*n+j]; 176 else value = v[i+j*m]; 177 if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue; 178 if (in[j] >= cstart && in[j] < cend) { 179 col = in[j] - cstart; 180 MatSetValue_SeqSELL_Private(A,row,col,value,addv,im[i],in[j],cp1,vp1,lastcol1,low1,high1); /* set one value */ 181 } else if (in[j] < 0) continue; 182 else if (PetscUnlikelyDebug(in[j] >= mat->cmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[j],mat->cmap->N-1); 183 else { 184 if (mat->was_assembled) { 185 if (!sell->colmap) { 186 ierr = MatCreateColmap_MPISELL_Private(mat);CHKERRQ(ierr); 187 } 188 #if defined(PETSC_USE_CTABLE) 189 ierr = PetscTableFind(sell->colmap,in[j]+1,&col);CHKERRQ(ierr); 190 col--; 191 #else 192 col = sell->colmap[in[j]] - 1; 193 #endif 194 if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) { 195 ierr = MatDisAssemble_MPISELL(mat);CHKERRQ(ierr); 196 col = in[j]; 197 /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */ 198 B = sell->B; 199 b = (Mat_SeqSELL*)B->data; 200 shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 201 cp2 = b->colidx+shift2; 202 vp2 = b->val+shift2; 203 nrow2 = b->rlen[row]; 204 low2 = 0; 205 high2 = nrow2; 206 } else if (col < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", im[i], in[j]); 207 } else col = in[j]; 208 MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */ 209 } 210 } 211 } else { 212 if (mat->nooffprocentries) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Setting off process row %D even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set",im[i]); 213 if (!sell->donotstash) { 214 mat->assembled = PETSC_FALSE; 215 if (roworiented) { 216 ierr = MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));CHKERRQ(ierr); 217 } else { 218 ierr = MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)));CHKERRQ(ierr); 219 } 220 } 221 } 222 } 223 PetscFunctionReturn(0); 224 } 225 226 PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 227 { 228 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 229 PetscErrorCode ierr; 230 PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend; 231 PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col; 232 233 PetscFunctionBegin; 234 for (i=0; i<m; i++) { 235 if (idxm[i] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",idxm[i]);*/ 236 if (idxm[i] >= mat->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm[i],mat->rmap->N-1); 237 if (idxm[i] >= rstart && idxm[i] < rend) { 238 row = idxm[i] - rstart; 239 for (j=0; j<n; j++) { 240 if (idxn[j] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",idxn[j]); */ 241 if (idxn[j] >= mat->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",idxn[j],mat->cmap->N-1); 242 if (idxn[j] >= cstart && idxn[j] < cend) { 243 col = idxn[j] - cstart; 244 ierr = MatGetValues(sell->A,1,&row,1,&col,v+i*n+j);CHKERRQ(ierr); 245 } else { 246 if (!sell->colmap) { 247 ierr = MatCreateColmap_MPISELL_Private(mat);CHKERRQ(ierr); 248 } 249 #if defined(PETSC_USE_CTABLE) 250 ierr = PetscTableFind(sell->colmap,idxn[j]+1,&col);CHKERRQ(ierr); 251 col--; 252 #else 253 col = sell->colmap[idxn[j]] - 1; 254 #endif 255 if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0; 256 else { 257 ierr = MatGetValues(sell->B,1,&row,1,&col,v+i*n+j);CHKERRQ(ierr); 258 } 259 } 260 } 261 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported"); 262 } 263 PetscFunctionReturn(0); 264 } 265 266 extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat,Vec,Vec); 267 268 PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode) 269 { 270 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 271 PetscErrorCode ierr; 272 PetscInt nstash,reallocs; 273 274 PetscFunctionBegin; 275 if (sell->donotstash || mat->nooffprocentries) PetscFunctionReturn(0); 276 277 ierr = MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);CHKERRQ(ierr); 278 ierr = MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);CHKERRQ(ierr); 279 ierr = PetscInfo2(sell->A,"Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode) 284 { 285 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 286 PetscErrorCode ierr; 287 PetscMPIInt n; 288 PetscInt i,flg; 289 PetscInt *row,*col; 290 PetscScalar *val; 291 PetscBool other_disassembled; 292 /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */ 293 PetscFunctionBegin; 294 if (!sell->donotstash && !mat->nooffprocentries) { 295 while (1) { 296 ierr = MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);CHKERRQ(ierr); 297 if (!flg) break; 298 299 for (i=0; i<n; i++) { /* assemble one by one */ 300 ierr = MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode);CHKERRQ(ierr); 301 } 302 } 303 ierr = MatStashScatterEnd_Private(&mat->stash);CHKERRQ(ierr); 304 } 305 ierr = MatAssemblyBegin(sell->A,mode);CHKERRQ(ierr); 306 ierr = MatAssemblyEnd(sell->A,mode);CHKERRQ(ierr); 307 308 /* 309 determine if any processor has disassembled, if so we must 310 also disassemble ourselfs, in order that we may reassemble. 311 */ 312 /* 313 if nonzero structure of submatrix B cannot change then we know that 314 no processor disassembled thus we can skip this stuff 315 */ 316 if (!((Mat_SeqSELL*)sell->B->data)->nonew) { 317 ierr = MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 318 if (mat->was_assembled && !other_disassembled) 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));CHKERRQ(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) { 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) { 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));CHKERRQ(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));CHKERRQ(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 799 PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr) 800 { 801 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 802 Mat a=sell->A,b=sell->B; 803 PetscErrorCode ierr; 804 PetscInt s1,s2,s3; 805 806 PetscFunctionBegin; 807 ierr = MatGetLocalSize(mat,&s2,&s3);CHKERRQ(ierr); 808 if (rr) { 809 ierr = VecGetLocalSize(rr,&s1);CHKERRQ(ierr); 810 if (s1!=s3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size"); 811 /* Overlap communication with computation. */ 812 ierr = VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 813 } 814 if (ll) { 815 ierr = VecGetLocalSize(ll,&s1);CHKERRQ(ierr); 816 if (s1!=s2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size"); 817 ierr = (*b->ops->diagonalscale)(b,ll,NULL);CHKERRQ(ierr); 818 } 819 /* scale the diagonal block */ 820 ierr = (*a->ops->diagonalscale)(a,ll,rr);CHKERRQ(ierr); 821 822 if (rr) { 823 /* Do a scatter end and then right scale the off-diagonal block */ 824 ierr = VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 825 ierr = (*b->ops->diagonalscale)(b,NULL,sell->lvec);CHKERRQ(ierr); 826 } 827 PetscFunctionReturn(0); 828 } 829 830 PetscErrorCode MatSetUnfactored_MPISELL(Mat A) 831 { 832 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 833 PetscErrorCode ierr; 834 835 PetscFunctionBegin; 836 ierr = MatSetUnfactored(a->A);CHKERRQ(ierr); 837 PetscFunctionReturn(0); 838 } 839 840 PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool *flag) 841 { 842 Mat_MPISELL *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data; 843 Mat a,b,c,d; 844 PetscBool flg; 845 PetscErrorCode ierr; 846 847 PetscFunctionBegin; 848 a = matA->A; b = matA->B; 849 c = matB->A; d = matB->B; 850 851 ierr = MatEqual(a,c,&flg);CHKERRQ(ierr); 852 if (flg) { 853 ierr = MatEqual(b,d,&flg);CHKERRQ(ierr); 854 } 855 ierr = MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 856 PetscFunctionReturn(0); 857 } 858 859 PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str) 860 { 861 PetscErrorCode ierr; 862 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 863 Mat_MPISELL *b=(Mat_MPISELL*)B->data; 864 865 PetscFunctionBegin; 866 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 867 if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) { 868 /* because of the column compression in the off-processor part of the matrix a->B, 869 the number of columns in a->B and b->B may be different, hence we cannot call 870 the MatCopy() directly on the two parts. If need be, we can provide a more 871 efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices 872 then copying the submatrices */ 873 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 874 } else { 875 ierr = MatCopy(a->A,b->A,str);CHKERRQ(ierr); 876 ierr = MatCopy(a->B,b->B,str);CHKERRQ(ierr); 877 } 878 PetscFunctionReturn(0); 879 } 880 881 PetscErrorCode MatSetUp_MPISELL(Mat A) 882 { 883 PetscErrorCode ierr; 884 885 PetscFunctionBegin; 886 ierr = MatMPISELLSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);CHKERRQ(ierr); 887 PetscFunctionReturn(0); 888 } 889 890 891 extern PetscErrorCode MatConjugate_SeqSELL(Mat); 892 893 PetscErrorCode MatConjugate_MPISELL(Mat mat) 894 { 895 #if defined(PETSC_USE_COMPLEX) 896 PetscErrorCode ierr; 897 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 898 899 PetscFunctionBegin; 900 ierr = MatConjugate_SeqSELL(sell->A);CHKERRQ(ierr); 901 ierr = MatConjugate_SeqSELL(sell->B);CHKERRQ(ierr); 902 #else 903 PetscFunctionBegin; 904 #endif 905 PetscFunctionReturn(0); 906 } 907 908 PetscErrorCode MatRealPart_MPISELL(Mat A) 909 { 910 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 911 PetscErrorCode ierr; 912 913 PetscFunctionBegin; 914 ierr = MatRealPart(a->A);CHKERRQ(ierr); 915 ierr = MatRealPart(a->B);CHKERRQ(ierr); 916 PetscFunctionReturn(0); 917 } 918 919 PetscErrorCode MatImaginaryPart_MPISELL(Mat A) 920 { 921 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 922 PetscErrorCode ierr; 923 924 PetscFunctionBegin; 925 ierr = MatImaginaryPart(a->A);CHKERRQ(ierr); 926 ierr = MatImaginaryPart(a->B);CHKERRQ(ierr); 927 PetscFunctionReturn(0); 928 } 929 930 PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values) 931 { 932 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 933 PetscErrorCode ierr; 934 935 PetscFunctionBegin; 936 ierr = MatInvertBlockDiagonal(a->A,values);CHKERRQ(ierr); 937 A->factorerrortype = a->A->factorerrortype; 938 PetscFunctionReturn(0); 939 } 940 941 static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx) 942 { 943 PetscErrorCode ierr; 944 Mat_MPISELL *sell=(Mat_MPISELL*)x->data; 945 946 PetscFunctionBegin; 947 ierr = MatSetRandom(sell->A,rctx);CHKERRQ(ierr); 948 ierr = MatSetRandom(sell->B,rctx);CHKERRQ(ierr); 949 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 950 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 951 PetscFunctionReturn(0); 952 } 953 954 PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A) 955 { 956 PetscErrorCode ierr; 957 958 PetscFunctionBegin; 959 ierr = PetscOptionsHead(PetscOptionsObject,"MPISELL options");CHKERRQ(ierr); 960 ierr = PetscOptionsTail();CHKERRQ(ierr); 961 PetscFunctionReturn(0); 962 } 963 964 PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a) 965 { 966 PetscErrorCode ierr; 967 Mat_MPISELL *msell=(Mat_MPISELL*)Y->data; 968 Mat_SeqSELL *sell=(Mat_SeqSELL*)msell->A->data; 969 970 PetscFunctionBegin; 971 if (!Y->preallocated) { 972 ierr = MatMPISELLSetPreallocation(Y,1,NULL,0,NULL);CHKERRQ(ierr); 973 } else if (!sell->nz) { 974 PetscInt nonew = sell->nonew; 975 ierr = MatSeqSELLSetPreallocation(msell->A,1,NULL);CHKERRQ(ierr); 976 sell->nonew = nonew; 977 } 978 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 979 PetscFunctionReturn(0); 980 } 981 982 PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool *missing,PetscInt *d) 983 { 984 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 985 PetscErrorCode ierr; 986 987 PetscFunctionBegin; 988 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices"); 989 ierr = MatMissingDiagonal(a->A,missing,d);CHKERRQ(ierr); 990 if (d) { 991 PetscInt rstart; 992 ierr = MatGetOwnershipRange(A,&rstart,NULL);CHKERRQ(ierr); 993 *d += rstart; 994 995 } 996 PetscFunctionReturn(0); 997 } 998 999 PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a) 1000 { 1001 PetscFunctionBegin; 1002 *a = ((Mat_MPISELL*)A->data)->A; 1003 PetscFunctionReturn(0); 1004 } 1005 1006 /* -------------------------------------------------------------------*/ 1007 static struct _MatOps MatOps_Values = {MatSetValues_MPISELL, 1008 NULL, 1009 NULL, 1010 MatMult_MPISELL, 1011 /* 4*/ MatMultAdd_MPISELL, 1012 MatMultTranspose_MPISELL, 1013 MatMultTransposeAdd_MPISELL, 1014 NULL, 1015 NULL, 1016 NULL, 1017 /*10*/ NULL, 1018 NULL, 1019 NULL, 1020 MatSOR_MPISELL, 1021 NULL, 1022 /*15*/ MatGetInfo_MPISELL, 1023 MatEqual_MPISELL, 1024 MatGetDiagonal_MPISELL, 1025 MatDiagonalScale_MPISELL, 1026 NULL, 1027 /*20*/ MatAssemblyBegin_MPISELL, 1028 MatAssemblyEnd_MPISELL, 1029 MatSetOption_MPISELL, 1030 MatZeroEntries_MPISELL, 1031 /*24*/ NULL, 1032 NULL, 1033 NULL, 1034 NULL, 1035 NULL, 1036 /*29*/ MatSetUp_MPISELL, 1037 NULL, 1038 NULL, 1039 MatGetDiagonalBlock_MPISELL, 1040 NULL, 1041 /*34*/ MatDuplicate_MPISELL, 1042 NULL, 1043 NULL, 1044 NULL, 1045 NULL, 1046 /*39*/ NULL, 1047 NULL, 1048 NULL, 1049 MatGetValues_MPISELL, 1050 MatCopy_MPISELL, 1051 /*44*/ NULL, 1052 MatScale_MPISELL, 1053 MatShift_MPISELL, 1054 MatDiagonalSet_MPISELL, 1055 NULL, 1056 /*49*/ MatSetRandom_MPISELL, 1057 NULL, 1058 NULL, 1059 NULL, 1060 NULL, 1061 /*54*/ MatFDColoringCreate_MPIXAIJ, 1062 NULL, 1063 MatSetUnfactored_MPISELL, 1064 NULL, 1065 NULL, 1066 /*59*/ NULL, 1067 MatDestroy_MPISELL, 1068 MatView_MPISELL, 1069 NULL, 1070 NULL, 1071 /*64*/ NULL, 1072 NULL, 1073 NULL, 1074 NULL, 1075 NULL, 1076 /*69*/ NULL, 1077 NULL, 1078 NULL, 1079 NULL, 1080 NULL, 1081 NULL, 1082 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */ 1083 MatSetFromOptions_MPISELL, 1084 NULL, 1085 NULL, 1086 NULL, 1087 /*80*/ NULL, 1088 NULL, 1089 NULL, 1090 /*83*/ NULL, 1091 NULL, 1092 NULL, 1093 NULL, 1094 NULL, 1095 NULL, 1096 /*89*/ NULL, 1097 NULL, 1098 NULL, 1099 NULL, 1100 NULL, 1101 /*94*/ NULL, 1102 NULL, 1103 NULL, 1104 NULL, 1105 NULL, 1106 /*99*/ NULL, 1107 NULL, 1108 NULL, 1109 MatConjugate_MPISELL, 1110 NULL, 1111 /*104*/NULL, 1112 MatRealPart_MPISELL, 1113 MatImaginaryPart_MPISELL, 1114 NULL, 1115 NULL, 1116 /*109*/NULL, 1117 NULL, 1118 NULL, 1119 NULL, 1120 MatMissingDiagonal_MPISELL, 1121 /*114*/NULL, 1122 NULL, 1123 MatGetGhosts_MPISELL, 1124 NULL, 1125 NULL, 1126 /*119*/NULL, 1127 NULL, 1128 NULL, 1129 NULL, 1130 NULL, 1131 /*124*/NULL, 1132 NULL, 1133 MatInvertBlockDiagonal_MPISELL, 1134 NULL, 1135 NULL, 1136 /*129*/NULL, 1137 NULL, 1138 NULL, 1139 NULL, 1140 NULL, 1141 /*134*/NULL, 1142 NULL, 1143 NULL, 1144 NULL, 1145 NULL, 1146 /*139*/NULL, 1147 NULL, 1148 NULL, 1149 MatFDColoringSetUp_MPIXAIJ, 1150 NULL, 1151 /*144*/NULL 1152 }; 1153 1154 /* ----------------------------------------------------------------------------------------*/ 1155 1156 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1157 { 1158 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1159 PetscErrorCode ierr; 1160 1161 PetscFunctionBegin; 1162 ierr = MatStoreValues(sell->A);CHKERRQ(ierr); 1163 ierr = MatStoreValues(sell->B);CHKERRQ(ierr); 1164 PetscFunctionReturn(0); 1165 } 1166 1167 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1168 { 1169 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1170 PetscErrorCode ierr; 1171 1172 PetscFunctionBegin; 1173 ierr = MatRetrieveValues(sell->A);CHKERRQ(ierr); 1174 ierr = MatRetrieveValues(sell->B);CHKERRQ(ierr); 1175 PetscFunctionReturn(0); 1176 } 1177 1178 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[]) 1179 { 1180 Mat_MPISELL *b; 1181 PetscErrorCode ierr; 1182 1183 PetscFunctionBegin; 1184 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1185 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1186 b = (Mat_MPISELL*)B->data; 1187 1188 if (!B->preallocated) { 1189 /* Explicitly create 2 MATSEQSELL matrices. */ 1190 ierr = MatCreate(PETSC_COMM_SELF,&b->A);CHKERRQ(ierr); 1191 ierr = MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);CHKERRQ(ierr); 1192 ierr = MatSetBlockSizesFromMats(b->A,B,B);CHKERRQ(ierr); 1193 ierr = MatSetType(b->A,MATSEQSELL);CHKERRQ(ierr); 1194 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);CHKERRQ(ierr); 1195 ierr = MatCreate(PETSC_COMM_SELF,&b->B);CHKERRQ(ierr); 1196 ierr = MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);CHKERRQ(ierr); 1197 ierr = MatSetBlockSizesFromMats(b->B,B,B);CHKERRQ(ierr); 1198 ierr = MatSetType(b->B,MATSEQSELL);CHKERRQ(ierr); 1199 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);CHKERRQ(ierr); 1200 } 1201 1202 ierr = MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen);CHKERRQ(ierr); 1203 ierr = MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen);CHKERRQ(ierr); 1204 B->preallocated = PETSC_TRUE; 1205 B->was_assembled = PETSC_FALSE; 1206 /* 1207 critical for MatAssemblyEnd to work. 1208 MatAssemblyBegin checks it to set up was_assembled 1209 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1210 */ 1211 B->assembled = PETSC_FALSE; 1212 PetscFunctionReturn(0); 1213 } 1214 1215 PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 1216 { 1217 Mat mat; 1218 Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data; 1219 PetscErrorCode ierr; 1220 1221 PetscFunctionBegin; 1222 *newmat = NULL; 1223 ierr = MatCreate(PetscObjectComm((PetscObject)matin),&mat);CHKERRQ(ierr); 1224 ierr = MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);CHKERRQ(ierr); 1225 ierr = MatSetBlockSizesFromMats(mat,matin,matin);CHKERRQ(ierr); 1226 ierr = MatSetType(mat,((PetscObject)matin)->type_name);CHKERRQ(ierr); 1227 a = (Mat_MPISELL*)mat->data; 1228 1229 mat->factortype = matin->factortype; 1230 mat->assembled = PETSC_TRUE; 1231 mat->insertmode = NOT_SET_VALUES; 1232 mat->preallocated = PETSC_TRUE; 1233 1234 a->size = oldmat->size; 1235 a->rank = oldmat->rank; 1236 a->donotstash = oldmat->donotstash; 1237 a->roworiented = oldmat->roworiented; 1238 a->rowindices = NULL; 1239 a->rowvalues = NULL; 1240 a->getrowactive = PETSC_FALSE; 1241 1242 ierr = PetscLayoutReference(matin->rmap,&mat->rmap);CHKERRQ(ierr); 1243 ierr = PetscLayoutReference(matin->cmap,&mat->cmap);CHKERRQ(ierr); 1244 1245 if (oldmat->colmap) { 1246 #if defined(PETSC_USE_CTABLE) 1247 ierr = PetscTableCreateCopy(oldmat->colmap,&a->colmap);CHKERRQ(ierr); 1248 #else 1249 ierr = PetscMalloc1(mat->cmap->N,&a->colmap);CHKERRQ(ierr); 1250 ierr = PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt));CHKERRQ(ierr); 1251 ierr = PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N);CHKERRQ(ierr); 1252 #endif 1253 } else a->colmap = NULL; 1254 if (oldmat->garray) { 1255 PetscInt len; 1256 len = oldmat->B->cmap->n; 1257 ierr = PetscMalloc1(len+1,&a->garray);CHKERRQ(ierr); 1258 ierr = PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));CHKERRQ(ierr); 1259 if (len) { ierr = PetscArraycpy(a->garray,oldmat->garray,len);CHKERRQ(ierr); } 1260 } else a->garray = NULL; 1261 1262 ierr = VecDuplicate(oldmat->lvec,&a->lvec);CHKERRQ(ierr); 1263 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);CHKERRQ(ierr); 1264 ierr = VecScatterCopy(oldmat->Mvctx,&a->Mvctx);CHKERRQ(ierr); 1265 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);CHKERRQ(ierr); 1266 ierr = MatDuplicate(oldmat->A,cpvalues,&a->A);CHKERRQ(ierr); 1267 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);CHKERRQ(ierr); 1268 ierr = MatDuplicate(oldmat->B,cpvalues,&a->B);CHKERRQ(ierr); 1269 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);CHKERRQ(ierr); 1270 ierr = PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);CHKERRQ(ierr); 1271 *newmat = mat; 1272 PetscFunctionReturn(0); 1273 } 1274 1275 /*@C 1276 MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format. 1277 For good matrix assembly performance the user should preallocate the matrix storage by 1278 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1279 1280 Collective 1281 1282 Input Parameters: 1283 + B - the matrix 1284 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1285 (same value is used for all local rows) 1286 . d_nnz - array containing the number of nonzeros in the various rows of the 1287 DIAGONAL portion of the local submatrix (possibly different for each row) 1288 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1289 The size of this array is equal to the number of local rows, i.e 'm'. 1290 For matrices that will be factored, you must leave room for (and set) 1291 the diagonal entry even if it is zero. 1292 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1293 submatrix (same value is used for all local rows). 1294 - o_nnz - array containing the number of nonzeros in the various rows of the 1295 OFF-DIAGONAL portion of the local submatrix (possibly different for 1296 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1297 structure. The size of this array is equal to the number 1298 of local rows, i.e 'm'. 1299 1300 If the *_nnz parameter is given then the *_nz parameter is ignored 1301 1302 The stored row and column indices begin with zero. 1303 1304 The parallel matrix is partitioned such that the first m0 rows belong to 1305 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1306 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1307 1308 The DIAGONAL portion of the local submatrix of a processor can be defined 1309 as the submatrix which is obtained by extraction the part corresponding to 1310 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1311 first row that belongs to the processor, r2 is the last row belonging to 1312 the this processor, and c1-c2 is range of indices of the local part of a 1313 vector suitable for applying the matrix to. This is an mxn matrix. In the 1314 common case of a square matrix, the row and column ranges are the same and 1315 the DIAGONAL part is also square. The remaining portion of the local 1316 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1317 1318 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1319 1320 You can call MatGetInfo() to get information on how effective the preallocation was; 1321 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1322 You can also run with the option -info and look for messages with the string 1323 malloc in them to see if additional memory allocation was needed. 1324 1325 Example usage: 1326 1327 Consider the following 8x8 matrix with 34 non-zero values, that is 1328 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1329 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1330 as follows: 1331 1332 .vb 1333 1 2 0 | 0 3 0 | 0 4 1334 Proc0 0 5 6 | 7 0 0 | 8 0 1335 9 0 10 | 11 0 0 | 12 0 1336 ------------------------------------- 1337 13 0 14 | 15 16 17 | 0 0 1338 Proc1 0 18 0 | 19 20 21 | 0 0 1339 0 0 0 | 22 23 0 | 24 0 1340 ------------------------------------- 1341 Proc2 25 26 27 | 0 0 28 | 29 0 1342 30 0 0 | 31 32 33 | 0 34 1343 .ve 1344 1345 This can be represented as a collection of submatrices as: 1346 1347 .vb 1348 A B C 1349 D E F 1350 G H I 1351 .ve 1352 1353 Where the submatrices A,B,C are owned by proc0, D,E,F are 1354 owned by proc1, G,H,I are owned by proc2. 1355 1356 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1357 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1358 The 'M','N' parameters are 8,8, and have the same values on all procs. 1359 1360 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1361 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1362 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1363 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1364 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1365 matrix, ans [DF] as another SeqSELL matrix. 1366 1367 When d_nz, o_nz parameters are specified, d_nz storage elements are 1368 allocated for every row of the local diagonal submatrix, and o_nz 1369 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1370 One way to choose d_nz and o_nz is to use the max nonzerors per local 1371 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1372 In this case, the values of d_nz,o_nz are: 1373 .vb 1374 proc0 : dnz = 2, o_nz = 2 1375 proc1 : dnz = 3, o_nz = 2 1376 proc2 : dnz = 1, o_nz = 4 1377 .ve 1378 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1379 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1380 for proc3. i.e we are using 12+15+10=37 storage locations to store 1381 34 values. 1382 1383 When d_nnz, o_nnz parameters are specified, the storage is specified 1384 for every row, coresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1385 In the above case the values for d_nnz,o_nnz are: 1386 .vb 1387 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1388 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1389 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1390 .ve 1391 Here the space allocated is according to nz (or maximum values in the nnz 1392 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1393 1394 Level: intermediate 1395 1396 .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(), 1397 MATMPISELL, MatGetInfo(), PetscSplitOwnership() 1398 @*/ 1399 PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1400 { 1401 PetscErrorCode ierr; 1402 1403 PetscFunctionBegin; 1404 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1405 PetscValidType(B,1); 1406 ierr = PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 1407 PetscFunctionReturn(0); 1408 } 1409 1410 /*@C 1411 MatCreateSELL - Creates a sparse parallel matrix in SELL format. 1412 1413 Collective 1414 1415 Input Parameters: 1416 + comm - MPI communicator 1417 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1418 This value should be the same as the local size used in creating the 1419 y vector for the matrix-vector product y = Ax. 1420 . n - This value should be the same as the local size used in creating the 1421 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1422 calculated if N is given) For square matrices n is almost always m. 1423 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 1424 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 1425 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1426 (same value is used for all local rows) 1427 . d_rlen - array containing the number of nonzeros in the various rows of the 1428 DIAGONAL portion of the local submatrix (possibly different for each row) 1429 or NULL, if d_rlenmax is used to specify the nonzero structure. 1430 The size of this array is equal to the number of local rows, i.e 'm'. 1431 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1432 submatrix (same value is used for all local rows). 1433 - o_rlen - array containing the number of nonzeros in the various rows of the 1434 OFF-DIAGONAL portion of the local submatrix (possibly different for 1435 each row) or NULL, if o_rlenmax is used to specify the nonzero 1436 structure. The size of this array is equal to the number 1437 of local rows, i.e 'm'. 1438 1439 Output Parameter: 1440 . A - the matrix 1441 1442 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1443 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1444 [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation] 1445 1446 Notes: 1447 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1448 1449 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1450 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1451 storage requirements for this matrix. 1452 1453 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one 1454 processor than it must be used on all processors that share the object for 1455 that argument. 1456 1457 The user MUST specify either the local or global matrix dimensions 1458 (possibly both). 1459 1460 The parallel matrix is partitioned across processors such that the 1461 first m0 rows belong to process 0, the next m1 rows belong to 1462 process 1, the next m2 rows belong to process 2 etc.. where 1463 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1464 values corresponding to [m x N] submatrix. 1465 1466 The columns are logically partitioned with the n0 columns belonging 1467 to 0th partition, the next n1 columns belonging to the next 1468 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1469 1470 The DIAGONAL portion of the local submatrix on any given processor 1471 is the submatrix corresponding to the rows and columns m,n 1472 corresponding to the given processor. i.e diagonal matrix on 1473 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1474 etc. The remaining portion of the local submatrix [m x (N-n)] 1475 constitute the OFF-DIAGONAL portion. The example below better 1476 illustrates this concept. 1477 1478 For a square global matrix we define each processor's diagonal portion 1479 to be its local rows and the corresponding columns (a square submatrix); 1480 each processor's off-diagonal portion encompasses the remainder of the 1481 local matrix (a rectangular submatrix). 1482 1483 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1484 1485 When calling this routine with a single process communicator, a matrix of 1486 type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this 1487 type of communicator, use the construction mechanism: 1488 MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...); 1489 1490 Options Database Keys: 1491 - -mat_sell_oneindex - Internally use indexing starting at 1 1492 rather than 0. Note that when calling MatSetValues(), 1493 the user still MUST index entries starting at 0! 1494 1495 1496 Example usage: 1497 1498 Consider the following 8x8 matrix with 34 non-zero values, that is 1499 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1500 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1501 as follows: 1502 1503 .vb 1504 1 2 0 | 0 3 0 | 0 4 1505 Proc0 0 5 6 | 7 0 0 | 8 0 1506 9 0 10 | 11 0 0 | 12 0 1507 ------------------------------------- 1508 13 0 14 | 15 16 17 | 0 0 1509 Proc1 0 18 0 | 19 20 21 | 0 0 1510 0 0 0 | 22 23 0 | 24 0 1511 ------------------------------------- 1512 Proc2 25 26 27 | 0 0 28 | 29 0 1513 30 0 0 | 31 32 33 | 0 34 1514 .ve 1515 1516 This can be represented as a collection of submatrices as: 1517 1518 .vb 1519 A B C 1520 D E F 1521 G H I 1522 .ve 1523 1524 Where the submatrices A,B,C are owned by proc0, D,E,F are 1525 owned by proc1, G,H,I are owned by proc2. 1526 1527 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1528 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1529 The 'M','N' parameters are 8,8, and have the same values on all procs. 1530 1531 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1532 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1533 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1534 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1535 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1536 matrix, ans [DF] as another SeqSELL matrix. 1537 1538 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1539 allocated for every row of the local diagonal submatrix, and o_rlenmax 1540 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1541 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1542 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1543 In this case, the values of d_rlenmax,o_rlenmax are: 1544 .vb 1545 proc0 : d_rlenmax = 2, o_rlenmax = 2 1546 proc1 : d_rlenmax = 3, o_rlenmax = 2 1547 proc2 : d_rlenmax = 1, o_rlenmax = 4 1548 .ve 1549 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1550 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1551 for proc3. i.e we are using 12+15+10=37 storage locations to store 1552 34 values. 1553 1554 When d_rlen, o_rlen parameters are specified, the storage is specified 1555 for every row, coresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1556 In the above case the values for d_nnz,o_nnz are: 1557 .vb 1558 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1559 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1560 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1561 .ve 1562 Here the space allocated is still 37 though there are 34 nonzeros because 1563 the allocation is always done according to rlenmax. 1564 1565 Level: intermediate 1566 1567 .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(), 1568 MATMPISELL, MatCreateMPISELLWithArrays() 1569 @*/ 1570 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) 1571 { 1572 PetscErrorCode ierr; 1573 PetscMPIInt size; 1574 1575 PetscFunctionBegin; 1576 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1577 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1578 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 1579 if (size > 1) { 1580 ierr = MatSetType(*A,MATMPISELL);CHKERRQ(ierr); 1581 ierr = MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);CHKERRQ(ierr); 1582 } else { 1583 ierr = MatSetType(*A,MATSEQSELL);CHKERRQ(ierr); 1584 ierr = MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);CHKERRQ(ierr); 1585 } 1586 PetscFunctionReturn(0); 1587 } 1588 1589 PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 1590 { 1591 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1592 PetscBool flg; 1593 PetscErrorCode ierr; 1594 1595 PetscFunctionBegin; 1596 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);CHKERRQ(ierr); 1597 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input"); 1598 if (Ad) *Ad = a->A; 1599 if (Ao) *Ao = a->B; 1600 if (colmap) *colmap = a->garray; 1601 PetscFunctionReturn(0); 1602 } 1603 1604 /*@C 1605 MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns 1606 1607 Not Collective 1608 1609 Input Parameters: 1610 + A - the matrix 1611 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1612 - row, col - index sets of rows and columns to extract (or NULL) 1613 1614 Output Parameter: 1615 . A_loc - the local sequential matrix generated 1616 1617 Level: developer 1618 1619 .seealso: MatGetOwnershipRange(), MatMPISELLGetLocalMat() 1620 1621 @*/ 1622 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc) 1623 { 1624 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1625 PetscErrorCode ierr; 1626 PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx; 1627 IS isrowa,iscola; 1628 Mat *aloc; 1629 PetscBool match; 1630 1631 PetscFunctionBegin; 1632 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match);CHKERRQ(ierr); 1633 if (!match) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input"); 1634 ierr = PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0);CHKERRQ(ierr); 1635 if (!row) { 1636 start = A->rmap->rstart; end = A->rmap->rend; 1637 ierr = ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa);CHKERRQ(ierr); 1638 } else { 1639 isrowa = *row; 1640 } 1641 if (!col) { 1642 start = A->cmap->rstart; 1643 cmap = a->garray; 1644 nzA = a->A->cmap->n; 1645 nzB = a->B->cmap->n; 1646 ierr = PetscMalloc1(nzA+nzB, &idx);CHKERRQ(ierr); 1647 ncols = 0; 1648 for (i=0; i<nzB; i++) { 1649 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1650 else break; 1651 } 1652 imark = i; 1653 for (i=0; i<nzA; i++) idx[ncols++] = start + i; 1654 for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; 1655 ierr = ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola);CHKERRQ(ierr); 1656 } else { 1657 iscola = *col; 1658 } 1659 if (scall != MAT_INITIAL_MATRIX) { 1660 ierr = PetscMalloc1(1,&aloc);CHKERRQ(ierr); 1661 aloc[0] = *A_loc; 1662 } 1663 ierr = MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc);CHKERRQ(ierr); 1664 *A_loc = aloc[0]; 1665 ierr = PetscFree(aloc);CHKERRQ(ierr); 1666 if (!row) { 1667 ierr = ISDestroy(&isrowa);CHKERRQ(ierr); 1668 } 1669 if (!col) { 1670 ierr = ISDestroy(&iscola);CHKERRQ(ierr); 1671 } 1672 ierr = PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0);CHKERRQ(ierr); 1673 PetscFunctionReturn(0); 1674 } 1675 1676 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1677 1678 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1679 { 1680 PetscErrorCode ierr; 1681 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1682 Mat B; 1683 Mat_MPIAIJ *b; 1684 1685 PetscFunctionBegin; 1686 if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1687 1688 if (reuse == MAT_REUSE_MATRIX) { 1689 B = *newmat; 1690 } else { 1691 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1692 ierr = MatSetType(B,MATMPIAIJ);CHKERRQ(ierr); 1693 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1694 ierr = MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 1695 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1696 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1697 } 1698 b = (Mat_MPIAIJ*) B->data; 1699 1700 if (reuse == MAT_REUSE_MATRIX) { 1701 ierr = MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);CHKERRQ(ierr); 1702 ierr = MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);CHKERRQ(ierr); 1703 } else { 1704 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1705 ierr = MatDestroy(&b->B);CHKERRQ(ierr); 1706 ierr = MatDisAssemble_MPISELL(A);CHKERRQ(ierr); 1707 ierr = MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);CHKERRQ(ierr); 1708 ierr = MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);CHKERRQ(ierr); 1709 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1710 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1711 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1712 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1713 } 1714 1715 if (reuse == MAT_INPLACE_MATRIX) { 1716 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1717 } else { 1718 *newmat = B; 1719 } 1720 PetscFunctionReturn(0); 1721 } 1722 1723 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1724 { 1725 PetscErrorCode ierr; 1726 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 1727 Mat B; 1728 Mat_MPISELL *b; 1729 1730 PetscFunctionBegin; 1731 if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1732 1733 if (reuse == MAT_REUSE_MATRIX) { 1734 B = *newmat; 1735 } else { 1736 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1737 ierr = MatSetType(B,MATMPISELL);CHKERRQ(ierr); 1738 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1739 ierr = MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 1740 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1741 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1742 } 1743 b = (Mat_MPISELL*) B->data; 1744 1745 if (reuse == MAT_REUSE_MATRIX) { 1746 ierr = MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);CHKERRQ(ierr); 1747 ierr = MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);CHKERRQ(ierr); 1748 } else { 1749 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1750 ierr = MatDestroy(&b->B);CHKERRQ(ierr); 1751 ierr = MatDisAssemble_MPIAIJ(A);CHKERRQ(ierr); 1752 ierr = MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);CHKERRQ(ierr); 1753 ierr = MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);CHKERRQ(ierr); 1754 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1755 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1756 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1757 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1758 } 1759 1760 if (reuse == MAT_INPLACE_MATRIX) { 1761 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1762 } else { 1763 *newmat = B; 1764 } 1765 PetscFunctionReturn(0); 1766 } 1767 1768 PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1769 { 1770 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 1771 PetscErrorCode ierr; 1772 Vec bb1=NULL; 1773 1774 PetscFunctionBegin; 1775 if (flag == SOR_APPLY_UPPER) { 1776 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1777 PetscFunctionReturn(0); 1778 } 1779 1780 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) { 1781 ierr = VecDuplicate(bb,&bb1);CHKERRQ(ierr); 1782 } 1783 1784 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1785 if (flag & SOR_ZERO_INITIAL_GUESS) { 1786 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1787 its--; 1788 } 1789 1790 while (its--) { 1791 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1792 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1793 1794 /* update rhs: bb1 = bb - B*x */ 1795 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 1796 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 1797 1798 /* local sweep */ 1799 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 1800 } 1801 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1802 if (flag & SOR_ZERO_INITIAL_GUESS) { 1803 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1804 its--; 1805 } 1806 while (its--) { 1807 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1808 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1809 1810 /* update rhs: bb1 = bb - B*x */ 1811 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 1812 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 1813 1814 /* local sweep */ 1815 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 1816 } 1817 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1818 if (flag & SOR_ZERO_INITIAL_GUESS) { 1819 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1820 its--; 1821 } 1822 while (its--) { 1823 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1824 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1825 1826 /* update rhs: bb1 = bb - B*x */ 1827 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 1828 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 1829 1830 /* local sweep */ 1831 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 1832 } 1833 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported"); 1834 1835 ierr = VecDestroy(&bb1);CHKERRQ(ierr); 1836 1837 matin->factorerrortype = mat->A->factorerrortype; 1838 PetscFunctionReturn(0); 1839 } 1840 1841 /*MC 1842 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1843 1844 Options Database Keys: 1845 . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions() 1846 1847 Level: beginner 1848 1849 .seealso: MatCreateSELL() 1850 M*/ 1851 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1852 { 1853 Mat_MPISELL *b; 1854 PetscErrorCode ierr; 1855 PetscMPIInt size; 1856 1857 PetscFunctionBegin; 1858 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRMPI(ierr); 1859 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1860 B->data = (void*)b; 1861 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1862 B->assembled = PETSC_FALSE; 1863 B->insertmode = NOT_SET_VALUES; 1864 b->size = size; 1865 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);CHKERRMPI(ierr); 1866 /* build cache for off array entries formed */ 1867 ierr = MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);CHKERRQ(ierr); 1868 1869 b->donotstash = PETSC_FALSE; 1870 b->colmap = NULL; 1871 b->garray = NULL; 1872 b->roworiented = PETSC_TRUE; 1873 1874 /* stuff used for matrix vector multiply */ 1875 b->lvec = NULL; 1876 b->Mvctx = NULL; 1877 1878 /* stuff for MatGetRow() */ 1879 b->rowindices = NULL; 1880 b->rowvalues = NULL; 1881 b->getrowactive = PETSC_FALSE; 1882 1883 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);CHKERRQ(ierr); 1884 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);CHKERRQ(ierr); 1885 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);CHKERRQ(ierr); 1886 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);CHKERRQ(ierr); 1887 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);CHKERRQ(ierr); 1888 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);CHKERRQ(ierr); 1889 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);CHKERRQ(ierr); 1890 PetscFunctionReturn(0); 1891 } 1892