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 PetscLayout 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 Vec tvec; 60 Mat tmat; 61 62 PetscFunctionBegin; 63 if (pc->setupcalled) { 64 ierr = KSPGetOperators(red->ksp,PETSC_NULL,&tmat,PETSC_NULL);CHKERRQ(ierr); 65 ierr = MatGetSubMatrix(pc->pmat,red->is,red->is,MAT_REUSE_MATRIX,&tmat);CHKERRQ(ierr); 66 } else { 67 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 68 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 69 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 70 ierr = PetscObjectGetNewTag((PetscObject)pc,&tag);CHKERRQ(ierr); 71 72 /* count non-diagonal rows on process */ 73 ierr = MatGetOwnershipRange(pc->mat,&rstart,&rend);CHKERRQ(ierr); 74 cnt = 0; 75 for (i=rstart; i<rend; i++) { 76 ierr = MatGetRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 77 if (nz > 1) cnt++; 78 ierr = MatRestoreRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 79 } 80 ierr = PetscMalloc(cnt*sizeof(PetscInt),&rows);CHKERRQ(ierr); 81 82 /* list non-diagonal rows on process */ 83 cnt = 0; 84 for (i=rstart; i<rend; i++) { 85 ierr = MatGetRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 86 if (nz > 1) rows[cnt++] = i; 87 ierr = MatRestoreRow(pc->mat,i,&nz,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 88 } 89 90 /* create PetscLayout for non-diagonal rows on each process */ 91 ierr = PetscLayoutCreate(comm,&map);CHKERRQ(ierr); 92 ierr = PetscLayoutSetLocalSize(map,cnt);CHKERRQ(ierr); 93 ierr = PetscLayoutSetBlockSize(map,1);CHKERRQ(ierr); 94 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 95 rstart = map->rstart; 96 rend = map->rend; 97 98 /* create PetscLayout for load-balanced non-diagonal rows on each process */ 99 ierr = PetscLayoutCreate(comm,&nmap);CHKERRQ(ierr); 100 ierr = MPI_Allreduce(&cnt,&ncnt,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 101 ierr = PetscLayoutSetSize(nmap,ncnt);CHKERRQ(ierr); 102 ierr = PetscLayoutSetBlockSize(nmap,1);CHKERRQ(ierr); 103 ierr = PetscLayoutSetUp(nmap);CHKERRQ(ierr); 104 105 /* 106 this code is taken from VecScatterCreate_PtoS() 107 Determines what rows need to be moved where to 108 load balance the non-diagonal rows 109 */ 110 /* count number of contributors to each processor */ 111 ierr = PetscMalloc2(size,PetscMPIInt,&nprocs,cnt,PetscInt,&owner);CHKERRQ(ierr); 112 ierr = PetscMemzero(nprocs,size*sizeof(PetscMPIInt));CHKERRQ(ierr); 113 j = 0; 114 nsends = 0; 115 for (i=rstart; i<rend; i++) { 116 if (i < nmap->range[j]) j = 0; 117 for (; j<size; j++) { 118 if (i < nmap->range[j+1]) { 119 if (!nprocs[j]++) nsends++; 120 owner[i-rstart] = j; 121 break; 122 } 123 } 124 } 125 /* inform other processors of number of messages and max length*/ 126 ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,nprocs,&nrecvs);CHKERRQ(ierr); 127 ierr = PetscGatherMessageLengths(comm,nsends,nrecvs,nprocs,&onodes1,&olengths1);CHKERRQ(ierr); 128 ierr = PetscSortMPIIntWithArray(nrecvs,onodes1,olengths1);CHKERRQ(ierr); 129 recvtotal = 0; for (i=0; i<nrecvs; i++) recvtotal += olengths1[i]; 130 131 /* post receives: rvalues - rows I will own; count - nu */ 132 ierr = PetscMalloc3(recvtotal,PetscInt,&rvalues,nrecvs,PetscInt,&source,nrecvs,MPI_Request,&recv_waits);CHKERRQ(ierr); 133 count = 0; 134 for (i=0; i<nrecvs; i++) { 135 ierr = MPI_Irecv((rvalues+count),olengths1[i],MPIU_INT,onodes1[i],tag,comm,recv_waits+i);CHKERRQ(ierr); 136 count += olengths1[i]; 137 } 138 139 /* do sends: 140 1) starts[i] gives the starting index in svalues for stuff going to 141 the ith processor 142 */ 143 ierr = PetscMalloc3(cnt,PetscInt,&svalues,nsends,MPI_Request,&send_waits,size,PetscInt,&starts);CHKERRQ(ierr); 144 starts[0] = 0; 145 for (i=1; i<size; i++) { starts[i] = starts[i-1] + nprocs[i-1];} 146 for (i=0; i<cnt; i++) { 147 svalues[starts[owner[i]]++] = rows[i]; 148 } 149 ierr = PetscFree(rows);CHKERRQ(ierr); 150 starts[0] = 0; 151 for (i=1; i<size; i++) { starts[i] = starts[i-1] + nprocs[i-1];} 152 count = 0; 153 for (i=0; i<size; i++) { 154 if (nprocs[i]) { 155 ierr = MPI_Isend(svalues+starts[i],nprocs[i],MPIU_INT,i,tag,comm,send_waits+count++);CHKERRQ(ierr); 156 } 157 } 158 159 /* wait on receives */ 160 count = nrecvs; 161 slen = 0; 162 while (count) { 163 ierr = MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);CHKERRQ(ierr); 164 /* unpack receives into our local space */ 165 ierr = MPI_Get_count(&recv_status,MPIU_INT,&n);CHKERRQ(ierr); 166 slen += n; 167 count--; 168 } 169 if (slen != recvtotal) SETERRQ2(PETSC_ERR_PLIB,"Total message lengths %D not expected %D",slen,recvtotal); 170 171 ierr = ISCreateGeneral(comm,slen,rvalues,&red->is);CHKERRQ(ierr); 172 173 /* free up all work space */ 174 ierr = PetscFree(olengths1);CHKERRQ(ierr); 175 ierr = PetscFree(onodes1);CHKERRQ(ierr); 176 ierr = PetscFree3(rvalues,source,recv_waits);CHKERRQ(ierr); 177 ierr = PetscFree2(nprocs,owner);CHKERRQ(ierr); 178 if (nsends) { /* wait on sends */ 179 ierr = PetscMalloc(nsends*sizeof(MPI_Status),&send_status);CHKERRQ(ierr); 180 ierr = MPI_Waitall(nsends,send_waits,send_status);CHKERRQ(ierr); 181 ierr = PetscFree(send_status);CHKERRQ(ierr); 182 } 183 ierr = PetscFree3(svalues,send_waits,starts);CHKERRQ(ierr); 184 ierr = PetscLayoutDestroy(map);CHKERRQ(ierr); 185 ierr = PetscLayoutDestroy(nmap);CHKERRQ(ierr); 186 187 ierr = VecCreateMPI(comm,slen,PETSC_DETERMINE,&red->b);CHKERRQ(ierr); 188 ierr = VecDuplicate(red->b,&red->x);CHKERRQ(ierr); 189 ierr = MatGetVecs(pc->pmat,&tvec,PETSC_NULL);CHKERRQ(ierr); 190 ierr = VecScatterCreate(tvec,red->is,red->b,PETSC_NULL,&red->scatter);CHKERRQ(ierr); 191 ierr = VecDestroy(tvec);CHKERRQ(ierr); 192 ierr = MatGetSubMatrix(pc->pmat,red->is,red->is,MAT_INITIAL_MATRIX,&tmat);CHKERRQ(ierr); 193 } 194 ierr = KSPSetOperators(red->ksp,tmat,tmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 195 ierr = MatDestroy(tmat);CHKERRQ(ierr); 196 ierr = KSPSetUp(red->ksp);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "PCApply_Redistribute" 202 static PetscErrorCode PCApply_Redistribute(PC pc,Vec b,Vec x) 203 { 204 PC_Redistribute *red = (PC_Redistribute*)pc->data; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 ierr = VecScatterBegin(red->scatter,b,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 209 ierr = VecScatterEnd(red->scatter,b,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 210 ierr = KSPSolve(red->ksp,red->b,red->x);CHKERRQ(ierr); 211 ierr = VecScatterBegin(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 212 ierr = VecScatterEnd(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "PCDestroy_Redistribute" 218 static PetscErrorCode PCDestroy_Redistribute(PC pc) 219 { 220 PC_Redistribute *red = (PC_Redistribute*)pc->data; 221 PetscErrorCode ierr; 222 223 PetscFunctionBegin; 224 if (red->scatter) {ierr = VecScatterDestroy(red->scatter);CHKERRQ(ierr);} 225 if (red->is) {ierr = ISDestroy(red->is);CHKERRQ(ierr);} 226 if (red->b) {ierr = VecDestroy(red->b);CHKERRQ(ierr);} 227 if (red->x) {ierr = VecDestroy(red->x);CHKERRQ(ierr);} 228 if (red->mat) {ierr = MatDestroy(red->mat);CHKERRQ(ierr);} 229 if (red->ksp) {ierr = KSPDestroy(red->ksp);CHKERRQ(ierr);} 230 ierr = PetscFree(red);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "PCSetFromOptions_Redistribute" 236 static PetscErrorCode PCSetFromOptions_Redistribute(PC pc) 237 { 238 PetscErrorCode ierr; 239 PC_Redistribute *red = (PC_Redistribute*)pc->data; 240 241 PetscFunctionBegin; 242 ierr = KSPSetFromOptions(red->ksp);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 /* -------------------------------------------------------------------------------------*/ 247 /*MC 248 PCREDISTRIBUTE - Redistributes a matrix for load balancing and then applys a KSP to that new matrix 249 250 Options for the redistribute preconditioners can be set with -redistribute_ksp_xxx <values> and -redistribute_pc_xxx <values> 251 252 Level: intermediate 253 254 .seealso: PCCreate(), PCSetType(), PCType (for list of available types) 255 M*/ 256 257 EXTERN_C_BEGIN 258 #undef __FUNCT__ 259 #define __FUNCT__ "PCCreate_Redistribute" 260 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Redistribute(PC pc) 261 { 262 PetscErrorCode ierr; 263 PC_Redistribute *red; 264 const char *prefix; 265 266 PetscFunctionBegin; 267 ierr = PetscNewLog(pc,PC_Redistribute,&red);CHKERRQ(ierr); 268 pc->data = (void*)red; 269 270 pc->ops->apply = PCApply_Redistribute; 271 pc->ops->applytranspose = 0; 272 pc->ops->setup = PCSetUp_Redistribute; 273 pc->ops->destroy = PCDestroy_Redistribute; 274 pc->ops->setfromoptions = PCSetFromOptions_Redistribute; 275 pc->ops->view = PCView_Redistribute; 276 277 ierr = KSPCreate(((PetscObject)pc)->comm,&red->ksp);CHKERRQ(ierr); 278 ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 279 ierr = PetscLogObjectParent(pc,red->ksp);CHKERRQ(ierr); 280 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 281 ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr); 282 ierr = KSPAppendOptionsPrefix(red->ksp,"redistribute_");CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 EXTERN_C_END 286