Actual source code: redistribute.c

      2: /*
      3:   This file defines a "solve the problem redistributely on each subgroup of processor" preconditioner.
      4: */
      5: #include <petsc/private/pcimpl.h>
      6: #include <petscksp.h>

      8: typedef struct {
      9:   KSP         ksp;
     10:   Vec         x,b;
     11:   VecScatter  scatter;
     12:   IS          is;
     13:   PetscInt    dcnt,*drows;    /* these are the local rows that have only diagonal entry */
     14:   PetscScalar *diag;
     15:   Vec         work;
     16: } PC_Redistribute;

     18: static PetscErrorCode PCView_Redistribute(PC pc,PetscViewer viewer)
     19: {
     20:   PC_Redistribute *red = (PC_Redistribute*)pc->data;
     21:   PetscErrorCode  ierr;
     22:   PetscBool       iascii,isstring;
     23:   PetscInt        ncnt,N;

     26:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
     27:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
     28:   if (iascii) {
     29:     MPIU_Allreduce(&red->dcnt,&ncnt,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
     30:     MatGetSize(pc->pmat,&N,NULL);
     31:     PetscViewerASCIIPrintf(viewer,"    Number rows eliminated %D Percentage rows eliminated %g\n",ncnt,100.0*((PetscReal)ncnt)/((PetscReal)N));
     32:     PetscViewerASCIIPrintf(viewer,"  Redistribute preconditioner: \n");
     33:     KSPView(red->ksp,viewer);
     34:   } else if (isstring) {
     35:     PetscViewerStringSPrintf(viewer," Redistribute preconditioner");
     36:     KSPView(red->ksp,viewer);
     37:   }
     38:   return(0);
     39: }

     41: static PetscErrorCode PCSetUp_Redistribute(PC pc)
     42: {
     43:   PC_Redistribute          *red = (PC_Redistribute*)pc->data;
     44:   PetscErrorCode           ierr;
     45:   MPI_Comm                 comm;
     46:   PetscInt                 rstart,rend,i,nz,cnt,*rows,ncnt,dcnt,*drows;
     47:   PetscLayout              map,nmap;
     48:   PetscMPIInt              size,tag,n;
     49:   PETSC_UNUSED PetscMPIInt imdex;
     50:   PetscInt                 *source = NULL;
     51:   PetscMPIInt              *sizes = NULL,nrecvs;
     52:   PetscInt                 j,nsends;
     53:   PetscInt                 *owner = NULL,*starts = NULL,count,slen;
     54:   PetscInt                 *rvalues,*svalues,recvtotal;
     55:   PetscMPIInt              *onodes1,*olengths1;
     56:   MPI_Request              *send_waits = NULL,*recv_waits = NULL;
     57:   MPI_Status               recv_status,*send_status;
     58:   Vec                      tvec,diag;
     59:   Mat                      tmat;
     60:   const PetscScalar        *d,*values;
     61:   const PetscInt           *cols;

     64:   if (pc->setupcalled) {
     65:     KSPGetOperators(red->ksp,NULL,&tmat);
     66:     MatCreateSubMatrix(pc->pmat,red->is,red->is,MAT_REUSE_MATRIX,&tmat);
     67:     KSPSetOperators(red->ksp,tmat,tmat);
     68:   } else {
     69:     PetscInt NN;

     71:     PetscObjectGetComm((PetscObject)pc,&comm);
     72:     MPI_Comm_size(comm,&size);
     73:     PetscObjectGetNewTag((PetscObject)pc,&tag);

     75:     /* count non-diagonal rows on process */
     76:     MatGetOwnershipRange(pc->mat,&rstart,&rend);
     77:     cnt  = 0;
     78:     for (i=rstart; i<rend; i++) {
     79:       MatGetRow(pc->mat,i,&nz,&cols,&values);
     80:       for (PetscInt j=0; j<nz; j++) {
     81:         if (values[j] != 0 && cols[j] != i) {
     82:           cnt++;
     83:           break;
     84:         }
     85:       }
     86:       MatRestoreRow(pc->mat,i,&nz,&cols,&values);
     87:     }
     88:     PetscMalloc1(cnt,&rows);
     89:     PetscMalloc1(rend - rstart - cnt,&drows);

     91:     /* list non-diagonal rows on process */
     92:     cnt = 0; dcnt = 0;
     93:     for (i=rstart; i<rend; i++) {
     94:       PetscBool diagonly = PETSC_TRUE;
     95:       MatGetRow(pc->mat,i,&nz,&cols,&values);
     96:       for (PetscInt j=0; j<nz; j++) {
     97:         if (values[j] != 0 && cols[j] != i) {
     98:           diagonly = PETSC_FALSE;
     99:           break;
    100:         }
    101:       }
    102:       if (!diagonly) rows[cnt++] = i;
    103:       else drows[dcnt++] = i - rstart;
    104:       MatRestoreRow(pc->mat,i,&nz,&cols,&values);
    105:     }

    107:     /* create PetscLayout for non-diagonal rows on each process */
    108:     PetscLayoutCreate(comm,&map);
    109:     PetscLayoutSetLocalSize(map,cnt);
    110:     PetscLayoutSetBlockSize(map,1);
    111:     PetscLayoutSetUp(map);
    112:     rstart = map->rstart;
    113:     rend   = map->rend;

    115:     /* create PetscLayout for load-balanced non-diagonal rows on each process */
    116:     PetscLayoutCreate(comm,&nmap);
    117:     MPIU_Allreduce(&cnt,&ncnt,1,MPIU_INT,MPI_SUM,comm);
    118:     PetscLayoutSetSize(nmap,ncnt);
    119:     PetscLayoutSetBlockSize(nmap,1);
    120:     PetscLayoutSetUp(nmap);

    122:     MatGetSize(pc->pmat,&NN,NULL);
    123:     PetscInfo2(pc,"Number of diagonal rows eliminated %d, percentage eliminated %g\n",NN-ncnt,((PetscReal)(NN-ncnt))/((PetscReal)(NN)));

    125:     if (size > 1) {
    126:       /* the following block of code assumes MPI can send messages to self, which is not supported for MPI-uni hence we need to handle the size 1 case as a special case */
    127:       /*
    128:        this code is taken from VecScatterCreate_PtoS()
    129:        Determines what rows need to be moved where to
    130:        load balance the non-diagonal rows
    131:        */
    132:       /*  count number of contributors to each processor */
    133:       PetscMalloc2(size,&sizes,cnt,&owner);
    134:       PetscArrayzero(sizes,size);
    135:       j      = 0;
    136:       nsends = 0;
    137:       for (i=rstart; i<rend; i++) {
    138:         if (i < nmap->range[j]) j = 0;
    139:         for (; j<size; j++) {
    140:           if (i < nmap->range[j+1]) {
    141:             if (!sizes[j]++) nsends++;
    142:             owner[i-rstart] = j;
    143:             break;
    144:           }
    145:         }
    146:       }
    147:       /* inform other processors of number of messages and max length*/
    148:       PetscGatherNumberOfMessages(comm,NULL,sizes,&nrecvs);
    149:       PetscGatherMessageLengths(comm,nsends,nrecvs,sizes,&onodes1,&olengths1);
    150:       PetscSortMPIIntWithArray(nrecvs,onodes1,olengths1);
    151:       recvtotal = 0; for (i=0; i<nrecvs; i++) recvtotal += olengths1[i];

    153:       /* post receives:  rvalues - rows I will own; count - nu */
    154:       PetscMalloc3(recvtotal,&rvalues,nrecvs,&source,nrecvs,&recv_waits);
    155:       count = 0;
    156:       for (i=0; i<nrecvs; i++) {
    157:         MPI_Irecv((rvalues+count),olengths1[i],MPIU_INT,onodes1[i],tag,comm,recv_waits+i);
    158:         count += olengths1[i];
    159:       }

    161:       /* do sends:
    162:        1) starts[i] gives the starting index in svalues for stuff going to
    163:        the ith processor
    164:        */
    165:       PetscMalloc3(cnt,&svalues,nsends,&send_waits,size,&starts);
    166:       starts[0] = 0;
    167:       for (i=1; i<size; i++) starts[i] = starts[i-1] + sizes[i-1];
    168:       for (i=0; i<cnt; i++)  svalues[starts[owner[i]]++] = rows[i];
    169:       for (i=0; i<cnt; i++)  rows[i] = rows[i] - rstart;
    170:       red->drows = drows;
    171:       red->dcnt  = dcnt;
    172:       PetscFree(rows);

    174:       starts[0] = 0;
    175:       for (i=1; i<size; i++) starts[i] = starts[i-1] + sizes[i-1];
    176:       count = 0;
    177:       for (i=0; i<size; i++) {
    178:         if (sizes[i]) {
    179:           MPI_Isend(svalues+starts[i],sizes[i],MPIU_INT,i,tag,comm,send_waits+count++);
    180:         }
    181:       }

    183:       /*  wait on receives */
    184:       count = nrecvs;
    185:       slen  = 0;
    186:       while (count) {
    187:         MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);
    188:         /* unpack receives into our local space */
    189:         MPI_Get_count(&recv_status,MPIU_INT,&n);
    190:         slen += n;
    191:         count--;
    192:       }
    193:       if (slen != recvtotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Total message lengths %D not expected %D",slen,recvtotal);
    194:       ISCreateGeneral(comm,slen,rvalues,PETSC_COPY_VALUES,&red->is);

    196:       /* free all work space */
    197:       PetscFree(olengths1);
    198:       PetscFree(onodes1);
    199:       PetscFree3(rvalues,source,recv_waits);
    200:       PetscFree2(sizes,owner);
    201:       if (nsends) {   /* wait on sends */
    202:         PetscMalloc1(nsends,&send_status);
    203:         MPI_Waitall(nsends,send_waits,send_status);
    204:         PetscFree(send_status);
    205:       }
    206:       PetscFree3(svalues,send_waits,starts);
    207:     } else {
    208:       ISCreateGeneral(comm,cnt,rows,PETSC_OWN_POINTER,&red->is);
    209:       red->drows = drows;
    210:       red->dcnt  = dcnt;
    211:       slen = cnt;
    212:     }
    213:     PetscLayoutDestroy(&map);
    214:     PetscLayoutDestroy(&nmap);

    216:     VecCreateMPI(comm,slen,PETSC_DETERMINE,&red->b);
    217:     VecDuplicate(red->b,&red->x);
    218:     MatCreateVecs(pc->pmat,&tvec,NULL);
    219:     VecScatterCreate(tvec,red->is,red->b,NULL,&red->scatter);
    220:     VecDestroy(&tvec);
    221:     MatCreateSubMatrix(pc->pmat,red->is,red->is,MAT_INITIAL_MATRIX,&tmat);
    222:     KSPSetOperators(red->ksp,tmat,tmat);
    223:     MatDestroy(&tmat);
    224:   }

    226:   /* get diagonal portion of matrix */
    227:   PetscFree(red->diag);
    228:   PetscMalloc1(red->dcnt,&red->diag);
    229:   MatCreateVecs(pc->pmat,&diag,NULL);
    230:   MatGetDiagonal(pc->pmat,diag);
    231:   VecGetArrayRead(diag,&d);
    232:   for (i=0; i<red->dcnt; i++) red->diag[i] = 1.0/d[red->drows[i]];
    233:   VecRestoreArrayRead(diag,&d);
    234:   VecDestroy(&diag);
    235:   KSPSetUp(red->ksp);
    236:   return(0);
    237: }

    239: static PetscErrorCode PCApply_Redistribute(PC pc,Vec b,Vec x)
    240: {
    241:   PC_Redistribute   *red = (PC_Redistribute*)pc->data;
    242:   PetscErrorCode    ierr;
    243:   PetscInt          dcnt   = red->dcnt,i;
    244:   const PetscInt    *drows = red->drows;
    245:   PetscScalar       *xwork;
    246:   const PetscScalar *bwork,*diag = red->diag;

    249:   if (!red->work) {
    250:     VecDuplicate(b,&red->work);
    251:   }
    252:   /* compute the rows of solution that have diagonal entries only */
    253:   VecSet(x,0.0);         /* x = diag(A)^{-1} b */
    254:   VecGetArray(x,&xwork);
    255:   VecGetArrayRead(b,&bwork);
    256:   for (i=0; i<dcnt; i++) xwork[drows[i]] = diag[i]*bwork[drows[i]];
    257:   PetscLogFlops(dcnt);
    258:   VecRestoreArray(red->work,&xwork);
    259:   VecRestoreArrayRead(b,&bwork);
    260:   /* update the right hand side for the reduced system with diagonal rows (and corresponding columns) removed */
    261:   MatMult(pc->pmat,x,red->work);
    262:   VecAYPX(red->work,-1.0,b);   /* red->work = b - A x */

    264:   VecScatterBegin(red->scatter,red->work,red->b,INSERT_VALUES,SCATTER_FORWARD);
    265:   VecScatterEnd(red->scatter,red->work,red->b,INSERT_VALUES,SCATTER_FORWARD);
    266:   KSPSolve(red->ksp,red->b,red->x);
    267:   KSPCheckSolve(red->ksp,pc,red->x);
    268:   VecScatterBegin(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);
    269:   VecScatterEnd(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);
    270:   return(0);
    271: }

    273: static PetscErrorCode PCDestroy_Redistribute(PC pc)
    274: {
    275:   PC_Redistribute *red = (PC_Redistribute*)pc->data;
    276:   PetscErrorCode  ierr;

    279:   VecScatterDestroy(&red->scatter);
    280:   ISDestroy(&red->is);
    281:   VecDestroy(&red->b);
    282:   VecDestroy(&red->x);
    283:   KSPDestroy(&red->ksp);
    284:   VecDestroy(&red->work);
    285:   PetscFree(red->drows);
    286:   PetscFree(red->diag);
    287:   PetscFree(pc->data);
    288:   return(0);
    289: }

    291: static PetscErrorCode PCSetFromOptions_Redistribute(PetscOptionItems *PetscOptionsObject,PC pc)
    292: {
    293:   PetscErrorCode  ierr;
    294:   PC_Redistribute *red = (PC_Redistribute*)pc->data;

    297:   KSPSetFromOptions(red->ksp);
    298:   return(0);
    299: }

    301: /*@
    302:    PCRedistributeGetKSP - Gets the KSP created by the PCREDISTRIBUTE

    304:    Not Collective

    306:    Input Parameter:
    307: .  pc - the preconditioner context

    309:    Output Parameter:
    310: .  innerksp - the inner KSP

    312:    Level: advanced

    314: @*/
    315: PetscErrorCode  PCRedistributeGetKSP(PC pc,KSP *innerksp)
    316: {
    317:   PC_Redistribute *red = (PC_Redistribute*)pc->data;

    322:   *innerksp = red->ksp;
    323:   return(0);
    324: }

    326: /* -------------------------------------------------------------------------------------*/
    327: /*MC
    328:      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

    330:      Options for the redistribute preconditioners can be set with -redistribute_ksp_xxx <values> and -redistribute_pc_xxx <values>

    332:      Notes:
    333:     Usually run this with -ksp_type preonly

    335:      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
    336:      -pc_type redistribute -redistribute_ksp_type cg -redistribute_pc_type bjacobi -redistribute_sub_pc_type icc to take advantage of the symmetry.

    338:      This does NOT call a partitioner to reorder rows to lower communication; the ordering of the rows in the original matrix and redistributed matrix is the same.

    340:      Developer Notes:
    341:     Should add an option to this preconditioner to use a partitioner to redistribute the rows to lower communication.

    343:    Level: intermediate

    345: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCRedistributeGetKSP()
    346: M*/

    348: PETSC_EXTERN PetscErrorCode PCCreate_Redistribute(PC pc)
    349: {
    350:   PetscErrorCode  ierr;
    351:   PC_Redistribute *red;
    352:   const char      *prefix;

    355:   PetscNewLog(pc,&red);
    356:   pc->data = (void*)red;

    358:   pc->ops->apply          = PCApply_Redistribute;
    359:   pc->ops->applytranspose = NULL;
    360:   pc->ops->setup          = PCSetUp_Redistribute;
    361:   pc->ops->destroy        = PCDestroy_Redistribute;
    362:   pc->ops->setfromoptions = PCSetFromOptions_Redistribute;
    363:   pc->ops->view           = PCView_Redistribute;

    365:   KSPCreate(PetscObjectComm((PetscObject)pc),&red->ksp);
    366:   KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);
    367:   PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);
    368:   PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);
    369:   PCGetOptionsPrefix(pc,&prefix);
    370:   KSPSetOptionsPrefix(red->ksp,prefix);
    371:   KSPAppendOptionsPrefix(red->ksp,"redistribute_");
    372:   return(0);
    373: }