Actual source code: snescomposite.c

      2: /*
      3:       Defines a SNES that can consist of a collection of SNESes
      4: */
      5: #include <petsc/private/snesimpl.h>
      6: #include <petscblaslapack.h>

      8: const char *const        SNESCompositeTypes[]   = {"ADDITIVE","MULTIPLICATIVE","ADDITIVEOPTIMAL","SNESCompositeType","SNES_COMPOSITE",NULL};

     10: typedef struct _SNES_CompositeLink *SNES_CompositeLink;
     11: struct _SNES_CompositeLink {
     12:   SNES               snes;
     13:   PetscReal          dmp;
     14:   Vec                X;
     15:   SNES_CompositeLink next;
     16:   SNES_CompositeLink previous;
     17: };

     19: typedef struct {
     20:   SNES_CompositeLink head;
     21:   PetscInt           nsnes;
     22:   SNESCompositeType  type;
     23:   Vec                Xorig;
     24:   PetscInt           innerFailures; /* the number of inner failures we've seen */

     26:   /* context for ADDITIVEOPTIMAL */
     27:   Vec                *Xes,*Fes;      /* solution and residual vectors for the subsolvers */
     28:   PetscReal          *fnorms;        /* norms of the solutions */
     29:   PetscScalar        *h;             /* the matrix formed as q_ij = (rdot_i, rdot_j) */
     30:   PetscScalar        *g;             /* the dotproducts of the previous function with the candidate functions */
     31:   PetscBLASInt       n;              /* matrix dimension -- nsnes */
     32:   PetscBLASInt       nrhs;           /* the number of right hand sides */
     33:   PetscBLASInt       lda;            /* the padded matrix dimension */
     34:   PetscBLASInt       ldb;            /* the padded vector dimension */
     35:   PetscReal          *s;             /* the singular values */
     36:   PetscScalar        *beta;          /* the RHS and combination */
     37:   PetscReal          rcond;          /* the exit condition */
     38:   PetscBLASInt       rank;           /* the effective rank */
     39:   PetscScalar        *work;          /* the work vector */
     40:   PetscReal          *rwork;         /* the real work vector used for complex */
     41:   PetscBLASInt       lwork;          /* the size of the work vector */
     42:   PetscBLASInt       info;           /* the output condition */

     44:   PetscReal          rtol;           /* restart tolerance for accepting the combination */
     45:   PetscReal          stol;           /* restart tolerance for the combination */
     46: } SNES_Composite;

     48: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
     49: {
     50:   PetscErrorCode      ierr;
     51:   SNES_Composite      *jac = (SNES_Composite*)snes->data;
     52:   SNES_CompositeLink  next = jac->head;
     53:   Vec                 FSub;
     54:   SNESConvergedReason reason;

     57:   if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
     58:   if (snes->normschedule == SNES_NORM_ALWAYS) {
     59:     SNESSetInitialFunction(next->snes,F);
     60:   }
     61:   SNESSolve(next->snes,B,X);
     62:   SNESGetConvergedReason(next->snes,&reason);
     63:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
     64:     jac->innerFailures++;
     65:     if (jac->innerFailures >= snes->maxFailures) {
     66:       snes->reason = SNES_DIVERGED_INNER;
     67:       return(0);
     68:     }
     69:   }

     71:   while (next->next) {
     72:     /* only copy the function over in the case where the functions correspond */
     73:     if (next->snes->npcside== PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
     74:       SNESGetFunction(next->snes,&FSub,NULL,NULL);
     75:       next = next->next;
     76:       SNESSetInitialFunction(next->snes,FSub);
     77:     } else {
     78:       next = next->next;
     79:     }
     80:     SNESSolve(next->snes,B,X);
     81:     SNESGetConvergedReason(next->snes,&reason);
     82:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
     83:       jac->innerFailures++;
     84:       if (jac->innerFailures >= snes->maxFailures) {
     85:         snes->reason = SNES_DIVERGED_INNER;
     86:         return(0);
     87:       }
     88:     }
     89:   }
     90:   if (next->snes->npcside== PC_RIGHT) {
     91:     SNESGetFunction(next->snes,&FSub,NULL,NULL);
     92:     VecCopy(FSub,F);
     93:     if (fnorm) {
     94:       if (snes->xl && snes->xu) {
     95:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
     96:       } else {
     97:         VecNorm(F, NORM_2, fnorm);
     98:       }
     99:       SNESCheckFunctionNorm(snes,*fnorm);
    100:     }
    101:   } else if (snes->normschedule == SNES_NORM_ALWAYS) {
    102:     SNESComputeFunction(snes,X,F);
    103:     if (fnorm) {
    104:       if (snes->xl && snes->xu) {
    105:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
    106:       } else {
    107:         VecNorm(F, NORM_2, fnorm);
    108:       }
    109:       SNESCheckFunctionNorm(snes,*fnorm);
    110:     }
    111:   }
    112:   return(0);
    113: }

    115: static PetscErrorCode SNESCompositeApply_Additive(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
    116: {
    117:   PetscErrorCode      ierr;
    118:   SNES_Composite      *jac = (SNES_Composite*)snes->data;
    119:   SNES_CompositeLink  next = jac->head;
    120:   Vec                 Y,Xorig;
    121:   SNESConvergedReason reason;

    124:   Y = snes->vec_sol_update;
    125:   if (!jac->Xorig) {VecDuplicate(X,&jac->Xorig);}
    126:   Xorig = jac->Xorig;
    127:   VecCopy(X,Xorig);
    128:   if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
    129:   if (snes->normschedule == SNES_NORM_ALWAYS) {
    130:     SNESSetInitialFunction(next->snes,F);
    131:     while (next->next) {
    132:       next = next->next;
    133:       SNESSetInitialFunction(next->snes,F);
    134:     }
    135:   }
    136:   next = jac->head;
    137:   VecCopy(Xorig,Y);
    138:   SNESSolve(next->snes,B,Y);
    139:   SNESGetConvergedReason(next->snes,&reason);
    140:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    141:     jac->innerFailures++;
    142:     if (jac->innerFailures >= snes->maxFailures) {
    143:       snes->reason = SNES_DIVERGED_INNER;
    144:       return(0);
    145:     }
    146:   }
    147:   VecAXPY(Y,-1.0,Xorig);
    148:   VecAXPY(X,next->dmp,Y);
    149:   while (next->next) {
    150:     next = next->next;
    151:     VecCopy(Xorig,Y);
    152:     SNESSolve(next->snes,B,Y);
    153:     SNESGetConvergedReason(next->snes,&reason);
    154:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    155:       jac->innerFailures++;
    156:       if (jac->innerFailures >= snes->maxFailures) {
    157:         snes->reason = SNES_DIVERGED_INNER;
    158:         return(0);
    159:       }
    160:     }
    161:     VecAXPY(Y,-1.0,Xorig);
    162:     VecAXPY(X,next->dmp,Y);
    163:   }
    164:   if (snes->normschedule == SNES_NORM_ALWAYS) {
    165:     SNESComputeFunction(snes,X,F);
    166:     if (fnorm) {
    167:       if (snes->xl && snes->xu) {
    168:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
    169:       } else {
    170:         VecNorm(F, NORM_2, fnorm);
    171:       }
    172:       SNESCheckFunctionNorm(snes,*fnorm);
    173:     }
    174:   }
    175:   return(0);
    176: }

    178: /* approximately solve the overdetermined system:

    180:  2*F(x_i)\cdot F(\x_j)\alpha_i = 0
    181:  \alpha_i                      = 1

    183:  Which minimizes the L2 norm of the linearization of:
    184:  ||F(\sum_i \alpha_i*x_i)||^2

    186:  With the constraint that \sum_i\alpha_i = 1
    187:  Where x_i is the solution from the ith subsolver.
    188:  */
    189: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes,Vec X,Vec B,Vec F,PetscReal *fnorm)
    190: {
    191:   PetscErrorCode      ierr;
    192:   SNES_Composite      *jac = (SNES_Composite*)snes->data;
    193:   SNES_CompositeLink  next = jac->head;
    194:   Vec                 *Xes = jac->Xes,*Fes = jac->Fes;
    195:   PetscInt            i,j;
    196:   PetscScalar         tot,total,ftf;
    197:   PetscReal           min_fnorm;
    198:   PetscInt            min_i;
    199:   SNESConvergedReason reason;

    202:   if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");

    204:   if (snes->normschedule == SNES_NORM_ALWAYS) {
    205:     next = jac->head;
    206:     SNESSetInitialFunction(next->snes,F);
    207:     while (next->next) {
    208:       next = next->next;
    209:       SNESSetInitialFunction(next->snes,F);
    210:     }
    211:   }

    213:   next = jac->head;
    214:   i = 0;
    215:   VecCopy(X,Xes[i]);
    216:   SNESSolve(next->snes,B,Xes[i]);
    217:   SNESGetConvergedReason(next->snes,&reason);
    218:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    219:     jac->innerFailures++;
    220:     if (jac->innerFailures >= snes->maxFailures) {
    221:       snes->reason = SNES_DIVERGED_INNER;
    222:       return(0);
    223:     }
    224:   }
    225:   while (next->next) {
    226:     i++;
    227:     next = next->next;
    228:     VecCopy(X,Xes[i]);
    229:     SNESSolve(next->snes,B,Xes[i]);
    230:     SNESGetConvergedReason(next->snes,&reason);
    231:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    232:       jac->innerFailures++;
    233:       if (jac->innerFailures >= snes->maxFailures) {
    234:         snes->reason = SNES_DIVERGED_INNER;
    235:         return(0);
    236:       }
    237:     }
    238:   }

    240:   /* all the solutions are collected; combine optimally */
    241:   for (i=0;i<jac->n;i++) {
    242:     for (j=0;j<i+1;j++) {
    243:       VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
    244:     }
    245:     VecDotBegin(Fes[i],F,&jac->g[i]);
    246:   }

    248:   for (i=0;i<jac->n;i++) {
    249:     for (j=0;j<i+1;j++) {
    250:       VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
    251:       if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
    252:     }
    253:     VecDotEnd(Fes[i],F,&jac->g[i]);
    254:   }

    256:   ftf = (*fnorm)*(*fnorm);

    258:   for (i=0; i<jac->n; i++) {
    259:     for (j=i+1;j<jac->n;j++) {
    260:       jac->h[i + j*jac->n] = jac->h[j + i*jac->n];
    261:     }
    262:   }

    264:   for (i=0; i<jac->n; i++) {
    265:     for (j=0; j<jac->n; j++) {
    266:       jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
    267:     }
    268:     jac->beta[i] = ftf - jac->g[i];
    269:   }

    271:   jac->info  = 0;
    272:   jac->rcond = -1.;
    273:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
    274: #if defined(PETSC_USE_COMPLEX)
    275:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&jac->n,&jac->n,&jac->nrhs,jac->h,&jac->lda,jac->beta,&jac->lda,jac->s,&jac->rcond,&jac->rank,jac->work,&jac->lwork,jac->rwork,&jac->info));
    276: #else
    277:   PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&jac->n,&jac->n,&jac->nrhs,jac->h,&jac->lda,jac->beta,&jac->lda,jac->s,&jac->rcond,&jac->rank,jac->work,&jac->lwork,&jac->info));
    278: #endif
    279:   PetscFPTrapPop();
    280:   if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
    281:   if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
    282:   tot = 0.;
    283:   total = 0.;
    284:   for (i=0; i<jac->n; i++) {
    285:     if (snes->errorifnotconverged && PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
    286:     PetscInfo2(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
    287:     tot += jac->beta[i];
    288:     total += PetscAbsScalar(jac->beta[i]);
    289:   }
    290:   VecScale(X,(1. - tot));
    291:   VecMAXPY(X,jac->n,jac->beta,Xes);
    292:   SNESComputeFunction(snes,X,F);

    294:   if (snes->xl && snes->xu) {
    295:     SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
    296:   } else {
    297:     VecNorm(F, NORM_2, fnorm);
    298:   }

    300:   /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
    301:   min_fnorm = jac->fnorms[0];
    302:   min_i     = 0;
    303:   for (i=0; i<jac->n; i++) {
    304:     if (jac->fnorms[i] < min_fnorm) {
    305:       min_fnorm = jac->fnorms[i];
    306:       min_i     = i;
    307:     }
    308:   }

    310:   /* stagnation or divergence restart to the solution of the solver that failed the least */
    311:   if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
    312:     VecCopy(jac->Xes[min_i],X);
    313:     VecCopy(jac->Fes[min_i],F);
    314:     *fnorm = min_fnorm;
    315:   }
    316:   return(0);
    317: }

    319: static PetscErrorCode SNESSetUp_Composite(SNES snes)
    320: {
    321:   PetscErrorCode     ierr;
    322:   DM                 dm;
    323:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
    324:   SNES_CompositeLink next = jac->head;
    325:   PetscInt           n=0,i;
    326:   Vec                F;

    329:   SNESGetDM(snes,&dm);

    331:   if (snes->ops->computevariablebounds) {
    332:     /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
    333:     if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
    334:     if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
    335:     (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
    336:   }

    338:   while (next) {
    339:     n++;
    340:     SNESSetDM(next->snes,dm);
    341:     SNESSetApplicationContext(next->snes, snes->user);
    342:     if (snes->xl && snes->xu) {
    343:       if (snes->ops->computevariablebounds) {
    344:         SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
    345:       } else {
    346:         SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
    347:       }
    348:     }

    350:     next = next->next;
    351:   }
    352:   jac->nsnes = n;
    353:   SNESGetFunction(snes,&F,NULL,NULL);
    354:   if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
    355:     VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
    356:     PetscMalloc1(n,&jac->Fes);
    357:     PetscMalloc1(n,&jac->fnorms);
    358:     next = jac->head;
    359:     i = 0;
    360:     while (next) {
    361:       SNESGetFunction(next->snes,&F,NULL,NULL);
    362:       jac->Fes[i] = F;
    363:       PetscObjectReference((PetscObject)F);
    364:       next = next->next;
    365:       i++;
    366:     }
    367:     /* allocate the subspace direct solve area */
    368:     jac->nrhs  = 1;
    369:     jac->lda   = jac->nsnes;
    370:     jac->ldb   = jac->nsnes;
    371:     jac->n     = jac->nsnes;

    373:     PetscMalloc1(jac->n*jac->n,&jac->h);
    374:     PetscMalloc1(jac->n,&jac->beta);
    375:     PetscMalloc1(jac->n,&jac->s);
    376:     PetscMalloc1(jac->n,&jac->g);
    377:     jac->lwork = 12*jac->n;
    378: #if defined(PETSC_USE_COMPLEX)
    379:     PetscMalloc1(jac->lwork,&jac->rwork);
    380: #endif
    381:     PetscMalloc1(jac->lwork,&jac->work);
    382:   }

    384:   return(0);
    385: }

    387: static PetscErrorCode SNESReset_Composite(SNES snes)
    388: {
    389:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
    390:   PetscErrorCode   ierr;
    391:   SNES_CompositeLink next = jac->head;

    394:   while (next) {
    395:     SNESReset(next->snes);
    396:     next = next->next;
    397:   }
    398:   VecDestroy(&jac->Xorig);
    399:   if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
    400:   if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
    401:   PetscFree(jac->fnorms);
    402:   PetscFree(jac->h);
    403:   PetscFree(jac->s);
    404:   PetscFree(jac->g);
    405:   PetscFree(jac->beta);
    406:   PetscFree(jac->work);
    407:   PetscFree(jac->rwork);
    408:   return(0);
    409: }

    411: static PetscErrorCode SNESDestroy_Composite(SNES snes)
    412: {
    413:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
    414:   PetscErrorCode     ierr;
    415:   SNES_CompositeLink next = jac->head,next_tmp;

    418:   SNESReset_Composite(snes);
    419:   while (next) {
    420:     SNESDestroy(&next->snes);
    421:     next_tmp = next;
    422:     next     = next->next;
    423:     PetscFree(next_tmp);
    424:   }
    425:   PetscFree(snes->data);
    426:   return(0);
    427: }

    429: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
    430: {
    431:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
    432:   PetscErrorCode     ierr;
    433:   PetscInt           nmax = 8,i;
    434:   SNES_CompositeLink next;
    435:   char               *sneses[8];
    436:   PetscReal          dmps[8];
    437:   PetscBool          flg;

    440:   PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
    441:   PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
    442:   if (flg) {
    443:     SNESCompositeSetType(snes,jac->type);
    444:   }
    445:   PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
    446:   if (flg) {
    447:     for (i=0; i<nmax; i++) {
    448:       SNESCompositeAddSNES(snes,sneses[i]);
    449:       PetscFree(sneses[i]);   /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
    450:     }
    451:   }
    452:   PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
    453:   if (flg) {
    454:     for (i=0; i<nmax; i++) {
    455:       SNESCompositeSetDamping(snes,i,dmps[i]);
    456:     }
    457:   }
    458:   PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
    459:   PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
    460:   PetscOptionsTail();

    462:   next = jac->head;
    463:   while (next) {
    464:     SNESSetFromOptions(next->snes);
    465:     next = next->next;
    466:   }
    467:   return(0);
    468: }

    470: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
    471: {
    472:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
    473:   PetscErrorCode     ierr;
    474:   SNES_CompositeLink next = jac->head;
    475:   PetscBool          iascii;

    478:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    479:   if (iascii) {
    480:     PetscViewerASCIIPrintf(viewer,"  type - %s\n",SNESCompositeTypes[jac->type]);
    481:     PetscViewerASCIIPrintf(viewer,"  SNESes on composite preconditioner follow\n");
    482:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
    483:   }
    484:   if (iascii) {
    485:     PetscViewerASCIIPushTab(viewer);
    486:   }
    487:   while (next) {
    488:     SNESView(next->snes,viewer);
    489:     next = next->next;
    490:   }
    491:   if (iascii) {
    492:     PetscViewerASCIIPopTab(viewer);
    493:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
    494:   }
    495:   return(0);
    496: }

    498: /* ------------------------------------------------------------------------------*/

    500: static PetscErrorCode  SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
    501: {
    502:   SNES_Composite *jac = (SNES_Composite*)snes->data;

    505:   jac->type = type;
    506:   return(0);
    507: }

    509: static PetscErrorCode  SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
    510: {
    511:   SNES_Composite     *jac;
    512:   SNES_CompositeLink next,ilink;
    513:   PetscErrorCode     ierr;
    514:   PetscInt           cnt = 0;
    515:   const char         *prefix;
    516:   char               newprefix[20];
    517:   DM                 dm;

    520:   PetscNewLog(snes,&ilink);
    521:   ilink->next = NULL;
    522:   SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
    523:   PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
    524:   SNESGetDM(snes,&dm);
    525:   SNESSetDM(ilink->snes,dm);
    526:   SNESSetTolerances(ilink->snes,snes->abstol,snes->rtol,snes->stol,1,snes->max_funcs);
    527:   PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)ilink->snes);
    528:   jac  = (SNES_Composite*)snes->data;
    529:   next = jac->head;
    530:   if (!next) {
    531:     jac->head       = ilink;
    532:     ilink->previous = NULL;
    533:   } else {
    534:     cnt++;
    535:     while (next->next) {
    536:       next = next->next;
    537:       cnt++;
    538:     }
    539:     next->next      = ilink;
    540:     ilink->previous = next;
    541:   }
    542:   SNESGetOptionsPrefix(snes,&prefix);
    543:   SNESSetOptionsPrefix(ilink->snes,prefix);
    544:   PetscSNPrintf(newprefix,sizeof(newprefix),"sub_%d_",(int)cnt);
    545:   SNESAppendOptionsPrefix(ilink->snes,newprefix);
    546:   PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
    547:   SNESSetType(ilink->snes,type);
    548:   SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);

    550:   ilink->dmp = 1.0;
    551:   jac->nsnes++;
    552:   return(0);
    553: }

    555: static PetscErrorCode  SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
    556: {
    557:   SNES_Composite     *jac;
    558:   SNES_CompositeLink next;
    559:   PetscInt           i;

    562:   jac  = (SNES_Composite*)snes->data;
    563:   next = jac->head;
    564:   for (i=0; i<n; i++) {
    565:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
    566:     next = next->next;
    567:   }
    568:   *subsnes = next->snes;
    569:   return(0);
    570: }

    572: /* -------------------------------------------------------------------------------- */
    573: /*@C
    574:    SNESCompositeSetType - Sets the type of composite preconditioner.

    576:    Logically Collective on SNES

    578:    Input Parameter:
    579: +  snes - the preconditioner context
    580: -  type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE

    582:    Options Database Key:
    583: .  -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type

    585:    Level: Developer

    587: @*/
    588: PetscErrorCode  SNESCompositeSetType(SNES snes,SNESCompositeType type)
    589: {

    595:   PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
    596:   return(0);
    597: }

    599: /*@C
    600:    SNESCompositeAddSNES - Adds another SNES to the composite SNES.

    602:    Collective on SNES

    604:    Input Parameters:
    605: +  snes - the preconditioner context
    606: -  type - the type of the new preconditioner

    608:    Level: Developer

    610: @*/
    611: PetscErrorCode  SNESCompositeAddSNES(SNES snes,SNESType type)
    612: {

    617:   PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
    618:   return(0);
    619: }

    621: /*@
    622:    SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.

    624:    Not Collective

    626:    Input Parameter:
    627: +  snes - the preconditioner context
    628: -  n - the number of the snes requested

    630:    Output Parameters:
    631: .  subsnes - the SNES requested

    633:    Level: Developer

    635: .seealso: SNESCompositeAddSNES()
    636: @*/
    637: PetscErrorCode  SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
    638: {

    644:   PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
    645:   return(0);
    646: }

    648: /*@
    649:    SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.

    651:    Logically Collective on SNES

    653:    Input Parameter:
    654:    snes - the preconditioner context

    656:    Output Parameter:
    657:    n - the number of subsolvers

    659:    Level: Developer

    661: @*/
    662: PetscErrorCode  SNESCompositeGetNumber(SNES snes,PetscInt *n)
    663: {
    664:   SNES_Composite     *jac;
    665:   SNES_CompositeLink next;

    668:   jac  = (SNES_Composite*)snes->data;
    669:   next = jac->head;

    671:   *n = 0;
    672:   while (next) {
    673:     *n = *n + 1;
    674:     next = next->next;
    675:   }
    676:   return(0);
    677: }

    679: static PetscErrorCode  SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
    680: {
    681:   SNES_Composite     *jac;
    682:   SNES_CompositeLink next;
    683:   PetscInt           i;

    686:   jac  = (SNES_Composite*)snes->data;
    687:   next = jac->head;
    688:   for (i=0; i<n; i++) {
    689:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
    690:     next = next->next;
    691:   }
    692:   next->dmp = dmp;
    693:   return(0);
    694: }

    696: /*@
    697:    SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.

    699:    Not Collective

    701:    Input Parameter:
    702: +  snes - the preconditioner context
    703: .  n - the number of the snes requested
    704: -  dmp - the damping

    706:    Level: Developer

    708: .seealso: SNESCompositeAddSNES()
    709: @*/
    710: PetscErrorCode  SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
    711: {

    716:   PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
    717:   return(0);
    718: }

    720: static PetscErrorCode SNESSolve_Composite(SNES snes)
    721: {
    722:   Vec              F,X,B,Y;
    723:   PetscInt         i;
    724:   PetscReal        fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
    725:   PetscErrorCode   ierr;
    726:   SNESNormSchedule normtype;
    727:   SNES_Composite   *comp = (SNES_Composite*)snes->data;

    730:   X = snes->vec_sol;
    731:   F = snes->vec_func;
    732:   B = snes->vec_rhs;
    733:   Y = snes->vec_sol_update;

    735:   PetscObjectSAWsTakeAccess((PetscObject)snes);
    736:   snes->iter   = 0;
    737:   snes->norm   = 0.;
    738:   comp->innerFailures = 0;
    739:   PetscObjectSAWsGrantAccess((PetscObject)snes);
    740:   snes->reason = SNES_CONVERGED_ITERATING;
    741:   SNESGetNormSchedule(snes, &normtype);
    742:   if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
    743:     if (!snes->vec_func_init_set) {
    744:       SNESComputeFunction(snes,X,F);
    745:     } else snes->vec_func_init_set = PETSC_FALSE;

    747:     if (snes->xl && snes->xu) {
    748:       SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
    749:     } else {
    750:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
    751:     }
    752:     SNESCheckFunctionNorm(snes,fnorm);
    753:     PetscObjectSAWsTakeAccess((PetscObject)snes);
    754:     snes->iter = 0;
    755:     snes->norm = fnorm;
    756:     PetscObjectSAWsGrantAccess((PetscObject)snes);
    757:     SNESLogConvergenceHistory(snes,snes->norm,0);
    758:     SNESMonitor(snes,0,snes->norm);

    760:     /* test convergence */
    761:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
    762:     if (snes->reason) return(0);
    763:   } else {
    764:     PetscObjectSAWsGrantAccess((PetscObject)snes);
    765:     SNESLogConvergenceHistory(snes,snes->norm,0);
    766:     SNESMonitor(snes,0,snes->norm);
    767:   }

    769:   for (i = 0; i < snes->max_its; i++) {
    770:     /* Call general purpose update function */
    771:     if (snes->ops->update) {
    772:       (*snes->ops->update)(snes, snes->iter);
    773:     }

    775:     /* Copy the state before modification by application of the composite solver;
    776:        we will subtract the new state after application */
    777:     VecCopy(X, Y);

    779:     if (comp->type == SNES_COMPOSITE_ADDITIVE) {
    780:       SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
    781:     } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
    782:       SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
    783:     } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
    784:       SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
    785:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
    786:     if (snes->reason < 0) break;

    788:     /* Compute the solution update for convergence testing */
    789:     VecAYPX(Y, -1.0, X);

    791:     if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
    792:       SNESComputeFunction(snes,X,F);

    794:       if (snes->xl && snes->xu) {
    795:         VecNormBegin(X, NORM_2, &xnorm);
    796:         VecNormBegin(Y, NORM_2, &snorm);
    797:         SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
    798:         VecNormEnd(X, NORM_2, &xnorm);
    799:         VecNormEnd(Y, NORM_2, &snorm);
    800:       } else {
    801:         VecNormBegin(F, NORM_2, &fnorm);
    802:         VecNormBegin(X, NORM_2, &xnorm);
    803:         VecNormBegin(Y, NORM_2, &snorm);

    805:         VecNormEnd(F, NORM_2, &fnorm);
    806:         VecNormEnd(X, NORM_2, &xnorm);
    807:         VecNormEnd(Y, NORM_2, &snorm);
    808:       }
    809:       SNESCheckFunctionNorm(snes,fnorm);
    810:     } else if (normtype == SNES_NORM_ALWAYS) {
    811:       VecNormBegin(X, NORM_2, &xnorm);
    812:       VecNormBegin(Y, NORM_2, &snorm);
    813:       VecNormEnd(X, NORM_2, &xnorm);
    814:       VecNormEnd(Y, NORM_2, &snorm);
    815:     }
    816:     /* Monitor convergence */
    817:     PetscObjectSAWsTakeAccess((PetscObject)snes);
    818:     snes->iter = i+1;
    819:     snes->norm = fnorm;
    820:     snes->xnorm = xnorm;
    821:     snes->ynorm = snorm;
    822:     PetscObjectSAWsGrantAccess((PetscObject)snes);
    823:     SNESLogConvergenceHistory(snes,snes->norm,0);
    824:     SNESMonitor(snes,snes->iter,snes->norm);
    825:     /* Test for convergence */
    826:     if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
    827:     if (snes->reason) break;
    828:   }
    829:   if (normtype == SNES_NORM_ALWAYS) {
    830:     if (i == snes->max_its) {
    831:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
    832:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
    833:     }
    834:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
    835:   return(0);
    836: }

    838: /* -------------------------------------------------------------------------------------------*/

    840: /*MC
    841:      SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers

    843:    Options Database Keys:
    844: +  -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
    845: -  -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose

    847:    Level: intermediate

    849: .seealso:  SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
    850:            SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
    851:            SNESCompositeGetSNES()

    853:    References:
    854: .  1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
    855:    SIAM Review, 57(4), 2015

    857: M*/

    859: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
    860: {
    862:   SNES_Composite *jac;

    865:   PetscNewLog(snes,&jac);

    867:   snes->ops->solve           = SNESSolve_Composite;
    868:   snes->ops->setup           = SNESSetUp_Composite;
    869:   snes->ops->reset           = SNESReset_Composite;
    870:   snes->ops->destroy         = SNESDestroy_Composite;
    871:   snes->ops->setfromoptions  = SNESSetFromOptions_Composite;
    872:   snes->ops->view            = SNESView_Composite;

    874:   snes->usesksp        = PETSC_FALSE;

    876:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

    878:   snes->data = (void*)jac;
    879:   jac->type  = SNES_COMPOSITE_ADDITIVEOPTIMAL;
    880:   jac->Fes   = NULL;
    881:   jac->Xes   = NULL;
    882:   jac->fnorms = NULL;
    883:   jac->nsnes = 0;
    884:   jac->head  = NULL;
    885:   jac->stol  = 0.1;
    886:   jac->rtol  = 1.1;

    888:   jac->h     = NULL;
    889:   jac->s     = NULL;
    890:   jac->beta  = NULL;
    891:   jac->work  = NULL;
    892:   jac->rwork = NULL;

    894:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
    895:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
    896:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
    897:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
    898:   return(0);
    899: }