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 } PC_Redistribute; 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "PCView_Redistribute" 18 static PetscErrorCode PCView_Redistribute(PC pc,PetscViewer viewer) 19 { 20 PC_Redistribute *red = (PC_Redistribute*)pc->data; 21 PetscErrorCode ierr; 22 PetscTruth iascii,isstring; 23 24 PetscFunctionBegin; 25 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 26 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr); 27 if (iascii) { 28 ierr = PetscViewerASCIIPrintf(viewer," Redistribute preconditioner: \n");CHKERRQ(ierr); 29 ierr = KSPView(red->ksp,viewer);CHKERRQ(ierr); 30 } else if (isstring) { 31 ierr = PetscViewerStringSPrintf(viewer," Redistribute preconditioner");CHKERRQ(ierr); 32 ierr = KSPView(red->ksp,viewer);CHKERRQ(ierr); 33 } else { 34 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PC redistribute",((PetscObject)viewer)->type_name); 35 } 36 PetscFunctionReturn(0); 37 } 38 39 #include "private/matimpl.h" /*I "petscmat.h" I*/ 40 #undef __FUNCT__ 41 #define __FUNCT__ "PCSetUp_Redistribute" 42 static PetscErrorCode PCSetUp_Redistribute(PC pc) 43 { 44 PC_Redistribute *red = (PC_Redistribute*)pc->data; 45 PetscErrorCode ierr; 46 const char *prefix; 47 MPI_Comm comm; 48 PetscInt rstart,rend,i,nz,cnt,*rows,ncnt; 49 PetscMap *map,*nmap; 50 PetscMPIInt size,rank,imdex,tag,n; 51 PetscInt *source = PETSC_NULL; 52 PetscMPIInt *nprocs = PETSC_NULL,nrecvs; 53 PetscInt j,nsends; 54 PetscInt *owner = PETSC_NULL,*starts = PETSC_NULL,count,slen; 55 PetscInt *rvalues,*svalues,recvtotal; 56 PetscMPIInt *onodes1,*olengths1; 57 MPI_Request *send_waits = PETSC_NULL,*recv_waits = PETSC_NULL; 58 MPI_Status recv_status,*send_status; 59 60 PetscFunctionBegin; 61 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 62 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 63 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 64 65 if (!pc->setupcalled) { 66 if (!red->ksp) { 67 ierr = KSPCreate(comm,&red->ksp);CHKERRQ(ierr); 68 ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 69 ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr); 70 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 71 ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr); 72 ierr = KSPAppendOptionsPrefix(red->ksp,"redistribute_");CHKERRQ(ierr); 73 } 74 } else SETERRQ(PETSC_ERR_SUP,"Cannot yet re-setup a redistribute KSP"); 75 76 ierr = MatGetOwnershipRange(pc->mat,&rstart,&rend);CHKERRQ(ierr); 77 for (i=rstart; i<rend; i++) { 78 ierr = MatGetRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 79 if (nz > 1) cnt++; 80 } 81 ierr = PetscMalloc(cnt*sizeof(PetscInt),&rows);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) rows[cnt++] = i; 86 } 87 ierr = PetscMalloc(sizeof(PetscMap),&map);CHKERRQ(ierr); 88 ierr = PetscMapInitialize(comm,map);CHKERRQ(ierr); 89 ierr = PetscMapSetLocalSize(map,cnt);CHKERRQ(ierr); 90 ierr = PetscMapSetBlockSize(map,1);CHKERRQ(ierr); 91 ierr = PetscMapSetUp(map);CHKERRQ(ierr); 92 rstart = map->rstart; 93 rend = map->rend; 94 95 ierr = PetscMalloc(sizeof(PetscMap),&nmap);CHKERRQ(ierr); 96 ierr = PetscMapInitialize(comm,nmap);CHKERRQ(ierr); 97 ierr = MPI_Allreduce(&cnt,&ncnt,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 98 ierr = PetscMapSetSize(nmap,ncnt);CHKERRQ(ierr); 99 ierr = PetscMapSetBlockSize(nmap,1);CHKERRQ(ierr); 100 ierr = PetscMapSetUp(nmap);CHKERRQ(ierr); 101 102 /* this code is taken from VecScatterCreate_PtoS() */ 103 /* count number of contributors to each processor */ 104 ierr = PetscMalloc2(size,PetscMPIInt,&nprocs,cnt,PetscInt,&owner);CHKERRQ(ierr); 105 ierr = PetscMemzero(nprocs,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 106 j = 0; 107 nsends = 0; 108 for (i=rstart; i<rend; i++) { 109 if (i < nmap->range[j]) j = 0; 110 for (; j<size; j++) { 111 if (i < nmap->range[j+1]) { 112 if (!nprocs[j]++) nsends++; 113 owner[i] = j; 114 break; 115 } 116 } 117 } 118 /* inform other processors of number of messages and max length*/ 119 ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,nprocs,&nrecvs);CHKERRQ(ierr); 120 ierr = PetscGatherMessageLengths(comm,nsends,nrecvs,nprocs,&onodes1,&olengths1);CHKERRQ(ierr); 121 ierr = PetscSortMPIIntWithArray(nrecvs,onodes1,olengths1);CHKERRQ(ierr); 122 recvtotal = 0; for (i=0; i<nrecvs; i++) recvtotal += olengths1[i]; 123 124 /* post receives: rvalues - rows I will own; count - nu */ 125 ierr = PetscMalloc3(recvtotal,PetscInt,&rvalues,nrecvs,PetscInt,&source,nrecvs,MPI_Request,&recv_waits);CHKERRQ(ierr); 126 count = 0; 127 for (i=0; i<nrecvs; i++) { 128 ierr = MPI_Irecv((rvalues+count),olengths1[i],MPIU_INT,onodes1[i],tag,comm,recv_waits+i);CHKERRQ(ierr); 129 count += olengths1[i]; 130 } 131 132 /* do sends: 133 1) starts[i] gives the starting index in svalues for stuff going to 134 the ith processor 135 */ 136 ierr = PetscMalloc3(cnt,PetscInt,&svalues,nsends,MPI_Request,&send_waits,size+1,PetscInt,&starts);CHKERRQ(ierr); 137 starts[0] = 0; 138 for (i=1; i<size; i++) { starts[i] = starts[i-1] + nprocs[i-1];} 139 for (i=0; i<cnt; i++) { 140 if (owner[i] != rank) { 141 svalues[starts[owner[i]]++] = rows[i]; 142 } 143 } 144 145 starts[0] = 0; 146 for (i=1; i<size+1; i++) { starts[i] = starts[i-1] + nprocs[i-1];} 147 count = 0; 148 for (i=0; i<size; i++) { 149 if (nprocs[i]) { 150 ierr = MPI_Isend(svalues+starts[i],nprocs[i],MPIU_INT,i,tag,comm,send_waits+count++);CHKERRQ(ierr); 151 } 152 } 153 154 /* wait on receives */ 155 count = nrecvs; 156 slen = 0; 157 while (count) { 158 ierr = MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);CHKERRQ(ierr); 159 /* unpack receives into our local space */ 160 ierr = MPI_Get_count(&recv_status,MPIU_INT,&n);CHKERRQ(ierr); 161 slen += n; 162 count--; 163 } 164 if (slen != recvtotal) SETERRQ2(PETSC_ERR_PLIB,"Total message lengths %D not expected %D",slen,recvtotal); 165 166 ierr = PetscFree(olengths1);CHKERRQ(ierr); 167 ierr = PetscFree(onodes1);CHKERRQ(ierr); 168 ierr = PetscFree3(rvalues,source,recv_waits);CHKERRQ(ierr); 169 170 /* wait on sends */ 171 if (nsends) { 172 ierr = PetscMalloc(nsends*sizeof(MPI_Status),&send_status);CHKERRQ(ierr); 173 ierr = MPI_Waitall(nsends,send_waits,send_status);CHKERRQ(ierr); 174 ierr = PetscFree(send_status);CHKERRQ(ierr); 175 } 176 ierr = PetscFree3(svalues,send_waits,starts);CHKERRQ(ierr); 177 178 179 ierr = PetscFree(map);CHKERRQ(ierr); 180 ierr = PetscFree(nmap);CHKERRQ(ierr); 181 182 ierr = KSPSetUp(red->ksp);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "PCApply_Redistribute" 188 static PetscErrorCode PCApply_Redistribute(PC pc,Vec b,Vec x) 189 { 190 PC_Redistribute *red = (PC_Redistribute*)pc->data; 191 PetscErrorCode ierr; 192 193 PetscFunctionBegin; 194 ierr = VecScatterBegin(red->scatter,b,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 195 ierr = VecScatterEnd(red->scatter,b,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 196 ierr = KSPSolve(red->ksp,red->b,red->x);CHKERRQ(ierr); 197 ierr = VecScatterBegin(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 198 ierr = VecScatterEnd(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "PCDestroy_Redistribute" 204 static PetscErrorCode PCDestroy_Redistribute(PC pc) 205 { 206 PC_Redistribute *red = (PC_Redistribute*)pc->data; 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 if (red->scatter) {ierr = VecScatterDestroy(red->scatter);CHKERRQ(ierr);} 211 if (red->b) {ierr = VecDestroy(red->b);CHKERRQ(ierr);} 212 if (red->x) {ierr = VecDestroy(red->x);CHKERRQ(ierr);} 213 if (red->mat) {ierr = MatDestroy(red->mat);CHKERRQ(ierr);} 214 if (red->ksp) {ierr = KSPDestroy(red->ksp);CHKERRQ(ierr);} 215 ierr = PetscFree(red);CHKERRQ(ierr); 216 PetscFunctionReturn(0); 217 } 218 219 #undef __FUNCT__ 220 #define __FUNCT__ "PCSetFromOptions_Redistribute" 221 static PetscErrorCode PCSetFromOptions_Redistribute(PC pc) 222 { 223 PetscErrorCode ierr; 224 PC_Redistribute *red = (PC_Redistribute*)pc->data; 225 226 PetscFunctionBegin; 227 ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 /* -------------------------------------------------------------------------------------*/ 232 /*MC 233 PCREDISTRIBUTE - Redistributes a matrix for load balancing and then applys a KSP to that new matrix 234 235 Options for the redistribute preconditioners can be set with -redistribute_ksp_xxx <values> and -redistribute_pc_xxx <values> 236 237 Level: intermediate 238 239 .seealso: PCCreate(), PCSetType(), PCType (for list of available types) 240 M*/ 241 242 EXTERN_C_BEGIN 243 #undef __FUNCT__ 244 #define __FUNCT__ "PCCreate_Redistribute" 245 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Redistribute(PC pc) 246 { 247 PetscErrorCode ierr; 248 PC_Redistribute *red; 249 250 PetscFunctionBegin; 251 ierr = PetscNewLog(pc,PC_Redistribute,&red);CHKERRQ(ierr); 252 pc->data = (void*)red; 253 254 pc->ops->apply = PCApply_Redistribute; 255 pc->ops->applytranspose = 0; 256 pc->ops->setup = PCSetUp_Redistribute; 257 pc->ops->destroy = PCDestroy_Redistribute; 258 pc->ops->setfromoptions = PCSetFromOptions_Redistribute; 259 pc->ops->view = PCView_Redistribute; 260 PetscFunctionReturn(0); 261 } 262 EXTERN_C_END 263