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