1 #define PETSCKSP_DLL 2 3 /* 4 This file defines a "solve the problem redistributely on each subgroup of processor" preconditioner. 5 */ 6 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 7 #include "petscksp.h" 8 9 typedef struct { 10 KSP ksp; 11 Vec x,b; 12 Mat mat; 13 VecScatter scatter; 14 IS is; 15 PetscInt dcnt,*drows; /* these are the local rows that have only diagonal entry */ 16 PetscScalar *diag; 17 Vec work; 18 } PC_Redistribute; 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "PCView_Redistribute" 22 static PetscErrorCode PCView_Redistribute(PC pc,PetscViewer viewer) 23 { 24 PC_Redistribute *red = (PC_Redistribute*)pc->data; 25 PetscErrorCode ierr; 26 PetscTruth iascii,isstring; 27 PetscInt ncnt,N; 28 29 PetscFunctionBegin; 30 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 31 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 32 if (iascii) { 33 ierr = MPI_Allreduce(&red->dcnt,&ncnt,1,MPIU_INT,MPI_SUM,((PetscObject)pc)->comm);CHKERRQ(ierr); 34 ierr = MatGetSize(pc->pmat,&N,PETSC_NULL);CHKERRQ(ierr); 35 ierr = PetscViewerASCIIPrintf(viewer," Number rows eliminated %D Percentage rows eliminated %g\n",ncnt,100.0*((PetscScalar)ncnt)/((PetscScalar)N));CHKERRQ(ierr); 36 ierr = PetscViewerASCIIPrintf(viewer," Redistribute preconditioner: \n");CHKERRQ(ierr); 37 ierr = KSPView(red->ksp,viewer);CHKERRQ(ierr); 38 } else if (isstring) { 39 ierr = PetscViewerStringSPrintf(viewer," Redistribute preconditioner");CHKERRQ(ierr); 40 ierr = KSPView(red->ksp,viewer);CHKERRQ(ierr); 41 } else { 42 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PC redistribute",((PetscObject)viewer)->type_name); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 #include "private/matimpl.h" /*I "petscmat.h" I*/ 48 #undef __FUNCT__ 49 #define __FUNCT__ "PCSetUp_Redistribute" 50 static PetscErrorCode PCSetUp_Redistribute(PC pc) 51 { 52 PC_Redistribute *red = (PC_Redistribute*)pc->data; 53 PetscErrorCode ierr; 54 MPI_Comm comm; 55 PetscInt rstart,rend,i,nz,cnt,*rows,ncnt,dcnt,*drows; 56 PetscLayout map,nmap; 57 PetscMPIInt size,rank,imdex,tag,n; 58 PetscInt *source = PETSC_NULL; 59 PetscMPIInt *nprocs = PETSC_NULL,nrecvs; 60 PetscInt j,nsends; 61 PetscInt *owner = PETSC_NULL,*starts = PETSC_NULL,count,slen; 62 PetscInt *rvalues,*svalues,recvtotal; 63 PetscMPIInt *onodes1,*olengths1; 64 MPI_Request *send_waits = PETSC_NULL,*recv_waits = PETSC_NULL; 65 MPI_Status recv_status,*send_status; 66 Vec tvec,diag; 67 Mat tmat; 68 const PetscScalar *d; 69 70 PetscFunctionBegin; 71 if (pc->setupcalled) { 72 ierr = KSPGetOperators(red->ksp,PETSC_NULL,&tmat,PETSC_NULL);CHKERRQ(ierr); 73 ierr = MatGetSubMatrix(pc->pmat,red->is,red->is,MAT_REUSE_MATRIX,&tmat);CHKERRQ(ierr); 74 } else { 75 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 76 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 77 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 78 ierr = PetscObjectGetNewTag((PetscObject)pc,&tag);CHKERRQ(ierr); 79 80 /* count non-diagonal rows on process */ 81 ierr = MatGetOwnershipRange(pc->mat,&rstart,&rend);CHKERRQ(ierr); 82 cnt = 0; 83 for (i=rstart; i<rend; i++) { 84 ierr = MatGetRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 85 if (nz > 1) cnt++; 86 ierr = MatRestoreRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 87 } 88 ierr = PetscMalloc(cnt*sizeof(PetscInt),&rows);CHKERRQ(ierr); 89 ierr = PetscMalloc((rend - rstart - cnt)*sizeof(PetscInt),&drows);CHKERRQ(ierr); 90 91 /* list non-diagonal rows on process */ 92 cnt = 0; dcnt = 0; 93 for (i=rstart; i<rend; i++) { 94 ierr = MatGetRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 95 if (nz > 1) rows[cnt++] = i; 96 else drows[dcnt++] = i - rstart; 97 ierr = MatRestoreRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 98 } 99 100 /* create PetscLayout for non-diagonal rows on each process */ 101 ierr = PetscLayoutCreate(comm,&map);CHKERRQ(ierr); 102 ierr = PetscLayoutSetLocalSize(map,cnt);CHKERRQ(ierr); 103 ierr = PetscLayoutSetBlockSize(map,1);CHKERRQ(ierr); 104 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 105 rstart = map->rstart; 106 rend = map->rend; 107 108 /* create PetscLayout for load-balanced non-diagonal rows on each process */ 109 ierr = PetscLayoutCreate(comm,&nmap);CHKERRQ(ierr); 110 ierr = MPI_Allreduce(&cnt,&ncnt,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 111 ierr = PetscLayoutSetSize(nmap,ncnt);CHKERRQ(ierr); 112 ierr = PetscLayoutSetBlockSize(nmap,1);CHKERRQ(ierr); 113 ierr = PetscLayoutSetUp(nmap);CHKERRQ(ierr); 114 115 ierr = PetscInfo2(pc,"Number of diagonal rows eliminated %d, percentage eliminated %g\n",pc->pmat->rmap->N-ncnt,((PetscScalar)(pc->pmat->rmap->N-ncnt))/((PetscScalar)(pc->pmat->rmap->N))); 116 /* 117 this code is taken from VecScatterCreate_PtoS() 118 Determines what rows need to be moved where to 119 load balance the non-diagonal rows 120 */ 121 /* count number of contributors to each processor */ 122 ierr = PetscMalloc2(size,PetscMPIInt,&nprocs,cnt,PetscInt,&owner);CHKERRQ(ierr); 123 ierr = PetscMemzero(nprocs,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 124 j = 0; 125 nsends = 0; 126 for (i=rstart; i<rend; i++) { 127 if (i < nmap->range[j]) j = 0; 128 for (; j<size; j++) { 129 if (i < nmap->range[j+1]) { 130 if (!nprocs[j]++) nsends++; 131 owner[i-rstart] = j; 132 break; 133 } 134 } 135 } 136 /* inform other processors of number of messages and max length*/ 137 ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,nprocs,&nrecvs);CHKERRQ(ierr); 138 ierr = PetscGatherMessageLengths(comm,nsends,nrecvs,nprocs,&onodes1,&olengths1);CHKERRQ(ierr); 139 ierr = PetscSortMPIIntWithArray(nrecvs,onodes1,olengths1);CHKERRQ(ierr); 140 recvtotal = 0; for (i=0; i<nrecvs; i++) recvtotal += olengths1[i]; 141 142 /* post receives: rvalues - rows I will own; count - nu */ 143 ierr = PetscMalloc3(recvtotal,PetscInt,&rvalues,nrecvs,PetscInt,&source,nrecvs,MPI_Request,&recv_waits);CHKERRQ(ierr); 144 count = 0; 145 for (i=0; i<nrecvs; i++) { 146 ierr = MPI_Irecv((rvalues+count),olengths1[i],MPIU_INT,onodes1[i],tag,comm,recv_waits+i);CHKERRQ(ierr); 147 count += olengths1[i]; 148 } 149 150 /* do sends: 151 1) starts[i] gives the starting index in svalues for stuff going to 152 the ith processor 153 */ 154 ierr = PetscMalloc3(cnt,PetscInt,&svalues,nsends,MPI_Request,&send_waits,size,PetscInt,&starts);CHKERRQ(ierr); 155 starts[0] = 0; 156 for (i=1; i<size; i++) { starts[i] = starts[i-1] + nprocs[i-1];} 157 for (i=0; i<cnt; i++) { 158 svalues[starts[owner[i]]++] = rows[i]; 159 } 160 for (i=0; i<cnt; i++) rows[i] = rows[i] - rstart; 161 red->drows = drows; 162 red->dcnt = dcnt; 163 ierr = PetscFree(rows);CHKERRQ(ierr); 164 165 starts[0] = 0; 166 for (i=1; i<size; i++) { starts[i] = starts[i-1] + nprocs[i-1];} 167 count = 0; 168 for (i=0; i<size; i++) { 169 if (nprocs[i]) { 170 ierr = MPI_Isend(svalues+starts[i],nprocs[i],MPIU_INT,i,tag,comm,send_waits+count++);CHKERRQ(ierr); 171 } 172 } 173 174 /* wait on receives */ 175 count = nrecvs; 176 slen = 0; 177 while (count) { 178 ierr = MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);CHKERRQ(ierr); 179 /* unpack receives into our local space */ 180 ierr = MPI_Get_count(&recv_status,MPIU_INT,&n);CHKERRQ(ierr); 181 slen += n; 182 count--; 183 } 184 if (slen != recvtotal) SETERRQ2(PETSC_ERR_PLIB,"Total message lengths %D not expected %D",slen,recvtotal); 185 186 ierr = ISCreateGeneral(comm,slen,rvalues,&red->is);CHKERRQ(ierr); 187 188 /* free up all work space */ 189 ierr = PetscFree(olengths1);CHKERRQ(ierr); 190 ierr = PetscFree(onodes1);CHKERRQ(ierr); 191 ierr = PetscFree3(rvalues,source,recv_waits);CHKERRQ(ierr); 192 ierr = PetscFree2(nprocs,owner);CHKERRQ(ierr); 193 if (nsends) { /* wait on sends */ 194 ierr = PetscMalloc(nsends*sizeof(MPI_Status),&send_status);CHKERRQ(ierr); 195 ierr = MPI_Waitall(nsends,send_waits,send_status);CHKERRQ(ierr); 196 ierr = PetscFree(send_status);CHKERRQ(ierr); 197 } 198 ierr = PetscFree3(svalues,send_waits,starts);CHKERRQ(ierr); 199 ierr = PetscLayoutDestroy(map);CHKERRQ(ierr); 200 ierr = PetscLayoutDestroy(nmap);CHKERRQ(ierr); 201 202 ierr = VecCreateMPI(comm,slen,PETSC_DETERMINE,&red->b);CHKERRQ(ierr); 203 ierr = VecDuplicate(red->b,&red->x);CHKERRQ(ierr); 204 ierr = MatGetVecs(pc->pmat,&tvec,PETSC_NULL);CHKERRQ(ierr); 205 ierr = VecScatterCreate(tvec,red->is,red->b,PETSC_NULL,&red->scatter);CHKERRQ(ierr); 206 ierr = VecDestroy(tvec);CHKERRQ(ierr); 207 ierr = MatGetSubMatrix(pc->pmat,red->is,red->is,MAT_INITIAL_MATRIX,&tmat);CHKERRQ(ierr); 208 } 209 ierr = KSPSetOperators(red->ksp,tmat,tmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 210 ierr = MatDestroy(tmat);CHKERRQ(ierr); 211 212 /* get diagonal portion of matrix */ 213 ierr = PetscMalloc(red->dcnt*sizeof(PetscScalar),&red->diag);CHKERRQ(ierr); 214 ierr = MatGetVecs(pc->pmat,&diag,PETSC_NULL);CHKERRQ(ierr); 215 ierr = MatGetDiagonal(pc->pmat,diag);CHKERRQ(ierr); 216 ierr = VecGetArray(diag,(PetscScalar**)&d);CHKERRQ(ierr); 217 for (i=0; i<red->dcnt; i++) { 218 red->diag[i] = 1.0/d[red->drows[i]]; 219 } 220 ierr = VecRestoreArray(diag,(PetscScalar**)&d);CHKERRQ(ierr); 221 ierr = VecDestroy(diag);CHKERRQ(ierr); 222 223 ierr = KSPSetUp(red->ksp);CHKERRQ(ierr); 224 PetscFunctionReturn(0); 225 } 226 227 #undef __FUNCT__ 228 #define __FUNCT__ "PCApply_Redistribute" 229 static PetscErrorCode PCApply_Redistribute(PC pc,Vec b,Vec x) 230 { 231 PC_Redistribute *red = (PC_Redistribute*)pc->data; 232 PetscErrorCode ierr; 233 PetscInt dcnt = red->dcnt,i; 234 const PetscInt *drows = red->drows; 235 PetscScalar *xwork; 236 const PetscScalar *bwork,*diag = red->diag; 237 238 PetscFunctionBegin; 239 if (!red->work) { 240 ierr = VecDuplicate(b,&red->work);CHKERRQ(ierr); 241 } 242 /* compute the rows of solution that have diagonal entries only */ 243 ierr = VecSet(x,0.0);CHKERRQ(ierr); /* x = diag(A)^{-1} b */ 244 ierr = VecGetArray(x,&xwork);CHKERRQ(ierr); 245 ierr = VecGetArray(b,(PetscScalar**)&bwork);CHKERRQ(ierr); 246 for (i=0; i<dcnt; i++) { 247 xwork[drows[i]] = diag[i]*bwork[drows[i]]; 248 } 249 ierr = PetscLogFlops(dcnt);CHKERRQ(ierr); 250 ierr = VecRestoreArray(red->work,&xwork);CHKERRQ(ierr); 251 ierr = VecRestoreArray(b,(PetscScalar**)&bwork);CHKERRQ(ierr); 252 /* update the right hand side for the reduced system with diagonal rows (and corresponding columns) removed */ 253 ierr = MatMult(pc->pmat,x,red->work);CHKERRQ(ierr); 254 ierr = VecAYPX(red->work,-1.0,b);CHKERRQ(ierr); /* red->work = b - A x */ 255 256 ierr = VecScatterBegin(red->scatter,red->work,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 257 ierr = VecScatterEnd(red->scatter,red->work,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 258 ierr = KSPSolve(red->ksp,red->b,red->x);CHKERRQ(ierr); 259 ierr = VecScatterBegin(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 260 ierr = VecScatterEnd(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "PCDestroy_Redistribute" 266 static PetscErrorCode PCDestroy_Redistribute(PC pc) 267 { 268 PC_Redistribute *red = (PC_Redistribute*)pc->data; 269 PetscErrorCode ierr; 270 271 PetscFunctionBegin; 272 if (red->scatter) {ierr = VecScatterDestroy(red->scatter);CHKERRQ(ierr);} 273 if (red->is) {ierr = ISDestroy(red->is);CHKERRQ(ierr);} 274 if (red->b) {ierr = VecDestroy(red->b);CHKERRQ(ierr);} 275 if (red->x) {ierr = VecDestroy(red->x);CHKERRQ(ierr);} 276 if (red->mat) {ierr = MatDestroy(red->mat);CHKERRQ(ierr);} 277 if (red->ksp) {ierr = KSPDestroy(red->ksp);CHKERRQ(ierr);} 278 if (red->work) {ierr = VecDestroy(red->work);CHKERRQ(ierr);} 279 ierr = PetscFree(red->drows); 280 ierr = PetscFree(red->diag); 281 ierr = PetscFree(red);CHKERRQ(ierr); 282 PetscFunctionReturn(0); 283 } 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "PCSetFromOptions_Redistribute" 287 static PetscErrorCode PCSetFromOptions_Redistribute(PC pc) 288 { 289 PetscErrorCode ierr; 290 PC_Redistribute *red = (PC_Redistribute*)pc->data; 291 292 PetscFunctionBegin; 293 ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 297 /* -------------------------------------------------------------------------------------*/ 298 /*MC 299 PCREDISTRIBUTE - Redistributes a matrix for load balancing, removing the rows that only have a diagonal entry and then applys a KSP to that new matrix 300 301 Options for the redistribute preconditioners can be set with -redistribute_ksp_xxx <values> and -redistribute_pc_xxx <values> 302 303 Notes: Usually run this with -ksp_type preonly 304 305 If you have used MatZeroRows() to eliminate (for example, Dirichlet) boundary conditions for a symmetric problem then you can use, for example, -ksp_type preonly 306 -pc_type redistribute -redistribute_ksp_type cg -redistribute_pc_type bjacobi -redistribute_sub_pc_type icc to take advantage of the symmetry. 307 308 Level: intermediate 309 310 .seealso: PCCreate(), PCSetType(), PCType (for list of available types) 311 M*/ 312 313 EXTERN_C_BEGIN 314 #undef __FUNCT__ 315 #define __FUNCT__ "PCCreate_Redistribute" 316 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Redistribute(PC pc) 317 { 318 PetscErrorCode ierr; 319 PC_Redistribute *red; 320 const char *prefix; 321 322 PetscFunctionBegin; 323 ierr = PetscNewLog(pc,PC_Redistribute,&red);CHKERRQ(ierr); 324 pc->data = (void*)red; 325 326 pc->ops->apply = PCApply_Redistribute; 327 pc->ops->applytranspose = 0; 328 pc->ops->setup = PCSetUp_Redistribute; 329 pc->ops->destroy = PCDestroy_Redistribute; 330 pc->ops->setfromoptions = PCSetFromOptions_Redistribute; 331 pc->ops->view = PCView_Redistribute; 332 333 ierr = KSPCreate(((PetscObject)pc)->comm,&red->ksp);CHKERRQ(ierr); 334 ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 335 ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr); 336 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 337 ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr); 338 ierr = KSPAppendOptionsPrefix(red->ksp,"redistribute_");CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 EXTERN_C_END 342