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 PetscCheckFalse(!sell->garray,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 PetscCheckFalse(nonew == -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") 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 PetscCheckFalse(nonew == -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") 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 PetscCheck(im[i] < mat->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,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 PetscCheck(in[j] < mat->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,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 PetscCheckFalse(col < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") 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 PetscCheckFalse(mat->nooffprocentries,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Setting off process row %" PetscInt_FMT " 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; /* negative row */ 236 PetscCheck(idxm[i] < mat->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,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; /* negative column */ 241 PetscCheck(idxn[j] < mat->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,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 = PetscInfo(sell->A,"Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " 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 PetscCheckFalse(mat->was_assembled && !other_disassembled,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDisAssemble not implemented yet"); 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 PetscCheckFalse(nt != A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%" PetscInt_FMT ") and xx (%" PetscInt_FMT ")",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 PetscCheckFalse(A->rmap->N != A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block"); 473 PetscCheckFalse(A->rmap->rstart != A->cmap->rstart || A->rmap->rend != A->cmap->rend,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=%" PetscInt_FMT ", Cols=%" PetscInt_FMT,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 %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", 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 %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", 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 %" PetscInt_FMT " \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 %" PetscInt_FMT " \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 %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\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 = PetscInfo(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 SETERRQ(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 PetscCheckFalse(s1!=s3,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 PetscCheckFalse(s1!=s2,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 PetscCheckFalse(A->rmap->n != A->cmap->n,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 NULL, 1151 NULL, 1152 NULL 1153 }; 1154 1155 /* ----------------------------------------------------------------------------------------*/ 1156 1157 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1158 { 1159 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1160 PetscErrorCode ierr; 1161 1162 PetscFunctionBegin; 1163 ierr = MatStoreValues(sell->A);CHKERRQ(ierr); 1164 ierr = MatStoreValues(sell->B);CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167 1168 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1169 { 1170 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1171 PetscErrorCode ierr; 1172 1173 PetscFunctionBegin; 1174 ierr = MatRetrieveValues(sell->A);CHKERRQ(ierr); 1175 ierr = MatRetrieveValues(sell->B);CHKERRQ(ierr); 1176 PetscFunctionReturn(0); 1177 } 1178 1179 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[]) 1180 { 1181 Mat_MPISELL *b; 1182 PetscErrorCode ierr; 1183 1184 PetscFunctionBegin; 1185 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1186 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1187 b = (Mat_MPISELL*)B->data; 1188 1189 if (!B->preallocated) { 1190 /* Explicitly create 2 MATSEQSELL matrices. */ 1191 ierr = MatCreate(PETSC_COMM_SELF,&b->A);CHKERRQ(ierr); 1192 ierr = MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);CHKERRQ(ierr); 1193 ierr = MatSetBlockSizesFromMats(b->A,B,B);CHKERRQ(ierr); 1194 ierr = MatSetType(b->A,MATSEQSELL);CHKERRQ(ierr); 1195 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);CHKERRQ(ierr); 1196 ierr = MatCreate(PETSC_COMM_SELF,&b->B);CHKERRQ(ierr); 1197 ierr = MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);CHKERRQ(ierr); 1198 ierr = MatSetBlockSizesFromMats(b->B,B,B);CHKERRQ(ierr); 1199 ierr = MatSetType(b->B,MATSEQSELL);CHKERRQ(ierr); 1200 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);CHKERRQ(ierr); 1201 } 1202 1203 ierr = MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen);CHKERRQ(ierr); 1204 ierr = MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen);CHKERRQ(ierr); 1205 B->preallocated = PETSC_TRUE; 1206 B->was_assembled = PETSC_FALSE; 1207 /* 1208 critical for MatAssemblyEnd to work. 1209 MatAssemblyBegin checks it to set up was_assembled 1210 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1211 */ 1212 B->assembled = PETSC_FALSE; 1213 PetscFunctionReturn(0); 1214 } 1215 1216 PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 1217 { 1218 Mat mat; 1219 Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data; 1220 PetscErrorCode ierr; 1221 1222 PetscFunctionBegin; 1223 *newmat = NULL; 1224 ierr = MatCreate(PetscObjectComm((PetscObject)matin),&mat);CHKERRQ(ierr); 1225 ierr = MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);CHKERRQ(ierr); 1226 ierr = MatSetBlockSizesFromMats(mat,matin,matin);CHKERRQ(ierr); 1227 ierr = MatSetType(mat,((PetscObject)matin)->type_name);CHKERRQ(ierr); 1228 a = (Mat_MPISELL*)mat->data; 1229 1230 mat->factortype = matin->factortype; 1231 mat->assembled = PETSC_TRUE; 1232 mat->insertmode = NOT_SET_VALUES; 1233 mat->preallocated = PETSC_TRUE; 1234 1235 a->size = oldmat->size; 1236 a->rank = oldmat->rank; 1237 a->donotstash = oldmat->donotstash; 1238 a->roworiented = oldmat->roworiented; 1239 a->rowindices = NULL; 1240 a->rowvalues = NULL; 1241 a->getrowactive = PETSC_FALSE; 1242 1243 ierr = PetscLayoutReference(matin->rmap,&mat->rmap);CHKERRQ(ierr); 1244 ierr = PetscLayoutReference(matin->cmap,&mat->cmap);CHKERRQ(ierr); 1245 1246 if (oldmat->colmap) { 1247 #if defined(PETSC_USE_CTABLE) 1248 ierr = PetscTableCreateCopy(oldmat->colmap,&a->colmap);CHKERRQ(ierr); 1249 #else 1250 ierr = PetscMalloc1(mat->cmap->N,&a->colmap);CHKERRQ(ierr); 1251 ierr = PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt));CHKERRQ(ierr); 1252 ierr = PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N);CHKERRQ(ierr); 1253 #endif 1254 } else a->colmap = NULL; 1255 if (oldmat->garray) { 1256 PetscInt len; 1257 len = oldmat->B->cmap->n; 1258 ierr = PetscMalloc1(len+1,&a->garray);CHKERRQ(ierr); 1259 ierr = PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));CHKERRQ(ierr); 1260 if (len) { ierr = PetscArraycpy(a->garray,oldmat->garray,len);CHKERRQ(ierr); } 1261 } else a->garray = NULL; 1262 1263 ierr = VecDuplicate(oldmat->lvec,&a->lvec);CHKERRQ(ierr); 1264 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);CHKERRQ(ierr); 1265 ierr = VecScatterCopy(oldmat->Mvctx,&a->Mvctx);CHKERRQ(ierr); 1266 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);CHKERRQ(ierr); 1267 ierr = MatDuplicate(oldmat->A,cpvalues,&a->A);CHKERRQ(ierr); 1268 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);CHKERRQ(ierr); 1269 ierr = MatDuplicate(oldmat->B,cpvalues,&a->B);CHKERRQ(ierr); 1270 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);CHKERRQ(ierr); 1271 ierr = PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);CHKERRQ(ierr); 1272 *newmat = mat; 1273 PetscFunctionReturn(0); 1274 } 1275 1276 /*@C 1277 MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format. 1278 For good matrix assembly performance the user should preallocate the matrix storage by 1279 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1280 1281 Collective 1282 1283 Input Parameters: 1284 + B - the matrix 1285 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1286 (same value is used for all local rows) 1287 . d_nnz - array containing the number of nonzeros in the various rows of the 1288 DIAGONAL portion of the local submatrix (possibly different for each row) 1289 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1290 The size of this array is equal to the number of local rows, i.e 'm'. 1291 For matrices that will be factored, you must leave room for (and set) 1292 the diagonal entry even if it is zero. 1293 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1294 submatrix (same value is used for all local rows). 1295 - o_nnz - array containing the number of nonzeros in the various rows of the 1296 OFF-DIAGONAL portion of the local submatrix (possibly different for 1297 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1298 structure. The size of this array is equal to the number 1299 of local rows, i.e 'm'. 1300 1301 If the *_nnz parameter is given then the *_nz parameter is ignored 1302 1303 The stored row and column indices begin with zero. 1304 1305 The parallel matrix is partitioned such that the first m0 rows belong to 1306 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1307 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1308 1309 The DIAGONAL portion of the local submatrix of a processor can be defined 1310 as the submatrix which is obtained by extraction the part corresponding to 1311 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1312 first row that belongs to the processor, r2 is the last row belonging to 1313 the this processor, and c1-c2 is range of indices of the local part of a 1314 vector suitable for applying the matrix to. This is an mxn matrix. In the 1315 common case of a square matrix, the row and column ranges are the same and 1316 the DIAGONAL part is also square. The remaining portion of the local 1317 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1318 1319 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1320 1321 You can call MatGetInfo() to get information on how effective the preallocation was; 1322 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1323 You can also run with the option -info and look for messages with the string 1324 malloc in them to see if additional memory allocation was needed. 1325 1326 Example usage: 1327 1328 Consider the following 8x8 matrix with 34 non-zero values, that is 1329 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1330 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1331 as follows: 1332 1333 .vb 1334 1 2 0 | 0 3 0 | 0 4 1335 Proc0 0 5 6 | 7 0 0 | 8 0 1336 9 0 10 | 11 0 0 | 12 0 1337 ------------------------------------- 1338 13 0 14 | 15 16 17 | 0 0 1339 Proc1 0 18 0 | 19 20 21 | 0 0 1340 0 0 0 | 22 23 0 | 24 0 1341 ------------------------------------- 1342 Proc2 25 26 27 | 0 0 28 | 29 0 1343 30 0 0 | 31 32 33 | 0 34 1344 .ve 1345 1346 This can be represented as a collection of submatrices as: 1347 1348 .vb 1349 A B C 1350 D E F 1351 G H I 1352 .ve 1353 1354 Where the submatrices A,B,C are owned by proc0, D,E,F are 1355 owned by proc1, G,H,I are owned by proc2. 1356 1357 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1358 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1359 The 'M','N' parameters are 8,8, and have the same values on all procs. 1360 1361 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1362 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1363 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1364 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1365 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1366 matrix, ans [DF] as another SeqSELL matrix. 1367 1368 When d_nz, o_nz parameters are specified, d_nz storage elements are 1369 allocated for every row of the local diagonal submatrix, and o_nz 1370 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1371 One way to choose d_nz and o_nz is to use the max nonzerors per local 1372 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1373 In this case, the values of d_nz,o_nz are: 1374 .vb 1375 proc0 : dnz = 2, o_nz = 2 1376 proc1 : dnz = 3, o_nz = 2 1377 proc2 : dnz = 1, o_nz = 4 1378 .ve 1379 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1380 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1381 for proc3. i.e we are using 12+15+10=37 storage locations to store 1382 34 values. 1383 1384 When d_nnz, o_nnz parameters are specified, the storage is specified 1385 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1386 In the above case the values for d_nnz,o_nnz are: 1387 .vb 1388 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1389 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1390 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1391 .ve 1392 Here the space allocated is according to nz (or maximum values in the nnz 1393 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1394 1395 Level: intermediate 1396 1397 .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(), 1398 MATMPISELL, MatGetInfo(), PetscSplitOwnership() 1399 @*/ 1400 PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1401 { 1402 PetscErrorCode ierr; 1403 1404 PetscFunctionBegin; 1405 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1406 PetscValidType(B,1); 1407 ierr = PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 1408 PetscFunctionReturn(0); 1409 } 1410 1411 /*MC 1412 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1413 based on the sliced Ellpack format 1414 1415 Options Database Keys: 1416 . -mat_type sell - sets the matrix type to "seqsell" during a call to MatSetFromOptions() 1417 1418 Level: beginner 1419 1420 .seealso: MatCreateSell(), MATSEQSELL, MATSELL, MATSEQAIJ, MATAIJ, MATMPIAIJ 1421 M*/ 1422 1423 /*@C 1424 MatCreateSELL - Creates a sparse parallel matrix in SELL format. 1425 1426 Collective 1427 1428 Input Parameters: 1429 + comm - MPI communicator 1430 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1431 This value should be the same as the local size used in creating the 1432 y vector for the matrix-vector product y = Ax. 1433 . n - This value should be the same as the local size used in creating the 1434 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1435 calculated if N is given) For square matrices n is almost always m. 1436 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 1437 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 1438 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1439 (same value is used for all local rows) 1440 . d_rlen - array containing the number of nonzeros in the various rows of the 1441 DIAGONAL portion of the local submatrix (possibly different for each row) 1442 or NULL, if d_rlenmax is used to specify the nonzero structure. 1443 The size of this array is equal to the number of local rows, i.e 'm'. 1444 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1445 submatrix (same value is used for all local rows). 1446 - o_rlen - array containing the number of nonzeros in the various rows of the 1447 OFF-DIAGONAL portion of the local submatrix (possibly different for 1448 each row) or NULL, if o_rlenmax is used to specify the nonzero 1449 structure. The size of this array is equal to the number 1450 of local rows, i.e 'm'. 1451 1452 Output Parameter: 1453 . A - the matrix 1454 1455 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1456 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1457 [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation] 1458 1459 Notes: 1460 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1461 1462 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1463 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1464 storage requirements for this matrix. 1465 1466 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one 1467 processor than it must be used on all processors that share the object for 1468 that argument. 1469 1470 The user MUST specify either the local or global matrix dimensions 1471 (possibly both). 1472 1473 The parallel matrix is partitioned across processors such that the 1474 first m0 rows belong to process 0, the next m1 rows belong to 1475 process 1, the next m2 rows belong to process 2 etc.. where 1476 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1477 values corresponding to [m x N] submatrix. 1478 1479 The columns are logically partitioned with the n0 columns belonging 1480 to 0th partition, the next n1 columns belonging to the next 1481 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1482 1483 The DIAGONAL portion of the local submatrix on any given processor 1484 is the submatrix corresponding to the rows and columns m,n 1485 corresponding to the given processor. i.e diagonal matrix on 1486 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1487 etc. The remaining portion of the local submatrix [m x (N-n)] 1488 constitute the OFF-DIAGONAL portion. The example below better 1489 illustrates this concept. 1490 1491 For a square global matrix we define each processor's diagonal portion 1492 to be its local rows and the corresponding columns (a square submatrix); 1493 each processor's off-diagonal portion encompasses the remainder of the 1494 local matrix (a rectangular submatrix). 1495 1496 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1497 1498 When calling this routine with a single process communicator, a matrix of 1499 type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this 1500 type of communicator, use the construction mechanism: 1501 MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...); 1502 1503 Options Database Keys: 1504 - -mat_sell_oneindex - Internally use indexing starting at 1 1505 rather than 0. Note that when calling MatSetValues(), 1506 the user still MUST index entries starting at 0! 1507 1508 Example usage: 1509 1510 Consider the following 8x8 matrix with 34 non-zero values, that is 1511 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1512 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1513 as follows: 1514 1515 .vb 1516 1 2 0 | 0 3 0 | 0 4 1517 Proc0 0 5 6 | 7 0 0 | 8 0 1518 9 0 10 | 11 0 0 | 12 0 1519 ------------------------------------- 1520 13 0 14 | 15 16 17 | 0 0 1521 Proc1 0 18 0 | 19 20 21 | 0 0 1522 0 0 0 | 22 23 0 | 24 0 1523 ------------------------------------- 1524 Proc2 25 26 27 | 0 0 28 | 29 0 1525 30 0 0 | 31 32 33 | 0 34 1526 .ve 1527 1528 This can be represented as a collection of submatrices as: 1529 1530 .vb 1531 A B C 1532 D E F 1533 G H I 1534 .ve 1535 1536 Where the submatrices A,B,C are owned by proc0, D,E,F are 1537 owned by proc1, G,H,I are owned by proc2. 1538 1539 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1540 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1541 The 'M','N' parameters are 8,8, and have the same values on all procs. 1542 1543 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1544 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1545 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1546 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1547 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1548 matrix, ans [DF] as another SeqSELL matrix. 1549 1550 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1551 allocated for every row of the local diagonal submatrix, and o_rlenmax 1552 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1553 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1554 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1555 In this case, the values of d_rlenmax,o_rlenmax are: 1556 .vb 1557 proc0 : d_rlenmax = 2, o_rlenmax = 2 1558 proc1 : d_rlenmax = 3, o_rlenmax = 2 1559 proc2 : d_rlenmax = 1, o_rlenmax = 4 1560 .ve 1561 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1562 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1563 for proc3. i.e we are using 12+15+10=37 storage locations to store 1564 34 values. 1565 1566 When d_rlen, o_rlen parameters are specified, the storage is specified 1567 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1568 In the above case the values for d_nnz,o_nnz are: 1569 .vb 1570 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1571 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1572 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1573 .ve 1574 Here the space allocated is still 37 though there are 34 nonzeros because 1575 the allocation is always done according to rlenmax. 1576 1577 Level: intermediate 1578 1579 .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(), 1580 MATMPISELL, MatCreateMPISELLWithArrays() 1581 @*/ 1582 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) 1583 { 1584 PetscErrorCode ierr; 1585 PetscMPIInt size; 1586 1587 PetscFunctionBegin; 1588 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1589 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1590 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 1591 if (size > 1) { 1592 ierr = MatSetType(*A,MATMPISELL);CHKERRQ(ierr); 1593 ierr = MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen);CHKERRQ(ierr); 1594 } else { 1595 ierr = MatSetType(*A,MATSEQSELL);CHKERRQ(ierr); 1596 ierr = MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen);CHKERRQ(ierr); 1597 } 1598 PetscFunctionReturn(0); 1599 } 1600 1601 PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 1602 { 1603 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1604 PetscBool flg; 1605 PetscErrorCode ierr; 1606 1607 PetscFunctionBegin; 1608 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg);CHKERRQ(ierr); 1609 PetscCheckFalse(!flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input"); 1610 if (Ad) *Ad = a->A; 1611 if (Ao) *Ao = a->B; 1612 if (colmap) *colmap = a->garray; 1613 PetscFunctionReturn(0); 1614 } 1615 1616 /*@C 1617 MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns 1618 1619 Not Collective 1620 1621 Input Parameters: 1622 + A - the matrix 1623 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1624 - row, col - index sets of rows and columns to extract (or NULL) 1625 1626 Output Parameter: 1627 . A_loc - the local sequential matrix generated 1628 1629 Level: developer 1630 1631 .seealso: MatGetOwnershipRange(), MatMPISELLGetLocalMat() 1632 1633 @*/ 1634 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc) 1635 { 1636 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1637 PetscErrorCode ierr; 1638 PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx; 1639 IS isrowa,iscola; 1640 Mat *aloc; 1641 PetscBool match; 1642 1643 PetscFunctionBegin; 1644 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match);CHKERRQ(ierr); 1645 PetscCheckFalse(!match,PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input"); 1646 ierr = PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0);CHKERRQ(ierr); 1647 if (!row) { 1648 start = A->rmap->rstart; end = A->rmap->rend; 1649 ierr = ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa);CHKERRQ(ierr); 1650 } else { 1651 isrowa = *row; 1652 } 1653 if (!col) { 1654 start = A->cmap->rstart; 1655 cmap = a->garray; 1656 nzA = a->A->cmap->n; 1657 nzB = a->B->cmap->n; 1658 ierr = PetscMalloc1(nzA+nzB, &idx);CHKERRQ(ierr); 1659 ncols = 0; 1660 for (i=0; i<nzB; i++) { 1661 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1662 else break; 1663 } 1664 imark = i; 1665 for (i=0; i<nzA; i++) idx[ncols++] = start + i; 1666 for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; 1667 ierr = ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola);CHKERRQ(ierr); 1668 } else { 1669 iscola = *col; 1670 } 1671 if (scall != MAT_INITIAL_MATRIX) { 1672 ierr = PetscMalloc1(1,&aloc);CHKERRQ(ierr); 1673 aloc[0] = *A_loc; 1674 } 1675 ierr = MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc);CHKERRQ(ierr); 1676 *A_loc = aloc[0]; 1677 ierr = PetscFree(aloc);CHKERRQ(ierr); 1678 if (!row) { 1679 ierr = ISDestroy(&isrowa);CHKERRQ(ierr); 1680 } 1681 if (!col) { 1682 ierr = ISDestroy(&iscola);CHKERRQ(ierr); 1683 } 1684 ierr = PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0);CHKERRQ(ierr); 1685 PetscFunctionReturn(0); 1686 } 1687 1688 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1689 1690 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1691 { 1692 PetscErrorCode ierr; 1693 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1694 Mat B; 1695 Mat_MPIAIJ *b; 1696 1697 PetscFunctionBegin; 1698 PetscCheckFalse(!A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1699 1700 if (reuse == MAT_REUSE_MATRIX) { 1701 B = *newmat; 1702 } else { 1703 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1704 ierr = MatSetType(B,MATMPIAIJ);CHKERRQ(ierr); 1705 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1706 ierr = MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 1707 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1708 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1709 } 1710 b = (Mat_MPIAIJ*) B->data; 1711 1712 if (reuse == MAT_REUSE_MATRIX) { 1713 ierr = MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A);CHKERRQ(ierr); 1714 ierr = MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B);CHKERRQ(ierr); 1715 } else { 1716 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1717 ierr = MatDestroy(&b->B);CHKERRQ(ierr); 1718 ierr = MatDisAssemble_MPISELL(A);CHKERRQ(ierr); 1719 ierr = MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);CHKERRQ(ierr); 1720 ierr = MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);CHKERRQ(ierr); 1721 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1722 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1723 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1724 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1725 } 1726 1727 if (reuse == MAT_INPLACE_MATRIX) { 1728 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1729 } else { 1730 *newmat = B; 1731 } 1732 PetscFunctionReturn(0); 1733 } 1734 1735 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1736 { 1737 PetscErrorCode ierr; 1738 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 1739 Mat B; 1740 Mat_MPISELL *b; 1741 1742 PetscFunctionBegin; 1743 PetscCheckFalse(!A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1744 1745 if (reuse == MAT_REUSE_MATRIX) { 1746 B = *newmat; 1747 } else { 1748 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1749 ierr = MatSetType(B,MATMPISELL);CHKERRQ(ierr); 1750 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1751 ierr = MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 1752 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1753 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1754 } 1755 b = (Mat_MPISELL*) B->data; 1756 1757 if (reuse == MAT_REUSE_MATRIX) { 1758 ierr = MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A);CHKERRQ(ierr); 1759 ierr = MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B);CHKERRQ(ierr); 1760 } else { 1761 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1762 ierr = MatDestroy(&b->B);CHKERRQ(ierr); 1763 ierr = MatDisAssemble_MPIAIJ(A);CHKERRQ(ierr); 1764 ierr = MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A);CHKERRQ(ierr); 1765 ierr = MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B);CHKERRQ(ierr); 1766 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1767 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1768 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1769 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1770 } 1771 1772 if (reuse == MAT_INPLACE_MATRIX) { 1773 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1774 } else { 1775 *newmat = B; 1776 } 1777 PetscFunctionReturn(0); 1778 } 1779 1780 PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1781 { 1782 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 1783 PetscErrorCode ierr; 1784 Vec bb1=NULL; 1785 1786 PetscFunctionBegin; 1787 if (flag == SOR_APPLY_UPPER) { 1788 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1789 PetscFunctionReturn(0); 1790 } 1791 1792 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) { 1793 ierr = VecDuplicate(bb,&bb1);CHKERRQ(ierr); 1794 } 1795 1796 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1797 if (flag & SOR_ZERO_INITIAL_GUESS) { 1798 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1799 its--; 1800 } 1801 1802 while (its--) { 1803 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1804 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1805 1806 /* update rhs: bb1 = bb - B*x */ 1807 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 1808 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 1809 1810 /* local sweep */ 1811 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 1812 } 1813 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1814 if (flag & SOR_ZERO_INITIAL_GUESS) { 1815 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1816 its--; 1817 } 1818 while (its--) { 1819 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1820 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1821 1822 /* update rhs: bb1 = bb - B*x */ 1823 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 1824 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 1825 1826 /* local sweep */ 1827 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 1828 } 1829 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1830 if (flag & SOR_ZERO_INITIAL_GUESS) { 1831 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 1832 its--; 1833 } 1834 while (its--) { 1835 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1836 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1837 1838 /* update rhs: bb1 = bb - B*x */ 1839 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 1840 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 1841 1842 /* local sweep */ 1843 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 1844 } 1845 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported"); 1846 1847 ierr = VecDestroy(&bb1);CHKERRQ(ierr); 1848 1849 matin->factorerrortype = mat->A->factorerrortype; 1850 PetscFunctionReturn(0); 1851 } 1852 1853 /*MC 1854 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1855 1856 Options Database Keys: 1857 . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions() 1858 1859 Level: beginner 1860 1861 .seealso: MatCreateSELL() 1862 M*/ 1863 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1864 { 1865 Mat_MPISELL *b; 1866 PetscErrorCode ierr; 1867 PetscMPIInt size; 1868 1869 PetscFunctionBegin; 1870 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRMPI(ierr); 1871 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1872 B->data = (void*)b; 1873 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1874 B->assembled = PETSC_FALSE; 1875 B->insertmode = NOT_SET_VALUES; 1876 b->size = size; 1877 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);CHKERRMPI(ierr); 1878 /* build cache for off array entries formed */ 1879 ierr = MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);CHKERRQ(ierr); 1880 1881 b->donotstash = PETSC_FALSE; 1882 b->colmap = NULL; 1883 b->garray = NULL; 1884 b->roworiented = PETSC_TRUE; 1885 1886 /* stuff used for matrix vector multiply */ 1887 b->lvec = NULL; 1888 b->Mvctx = NULL; 1889 1890 /* stuff for MatGetRow() */ 1891 b->rowindices = NULL; 1892 b->rowvalues = NULL; 1893 b->getrowactive = PETSC_FALSE; 1894 1895 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL);CHKERRQ(ierr); 1896 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL);CHKERRQ(ierr); 1897 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL);CHKERRQ(ierr); 1898 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL);CHKERRQ(ierr); 1899 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ);CHKERRQ(ierr); 1900 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL);CHKERRQ(ierr); 1901 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMPISELL);CHKERRQ(ierr); 1902 PetscFunctionReturn(0); 1903 } 1904