Actual source code: admm.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1:  #include <../src/tao/constrained/impls/admm/admm.h>
  2:  #include <petsctao.h>
  3:  #include <petsc/private/petscimpl.h>

  5: /* Updates terminating criteria
  6:  *
  7:  * 1  ||r_k|| = ||Ax+Bz-c|| =< catol_admm* max{||Ax||,||Bz||,||c||}
  8:  *
  9:  * 2. Updates dual residual, d_k
 10:  *
 11:  * 3. ||d_k|| = ||mu*A^T*B(z_k-z_{k-1})|| =< gatol_admm * ||A^Ty||   */

 13: static PetscBool cited = PETSC_FALSE;
 14: static const char citation[] =
 15:   "@misc{xu2017adaptive,\n"
 16:   "   title={Adaptive Relaxed ADMM: Convergence Theory and Practical Implementation},\n"
 17:   "   author={Zheng Xu and Mario A. T. Figueiredo and Xiaoming Yuan and Christoph Studer and Tom Goldstein},\n"
 18:   "   year={2017},\n"
 19:   "   eprint={1704.02712},\n"
 20:   "   archivePrefix={arXiv},\n"
 21:   "   primaryClass={cs.CV}\n"
 22:   "}  \n";

 24: static PetscErrorCode TaoADMMToleranceUpdate(Tao tao)
 25: {
 26:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;
 28:   PetscReal      Axnorm,Bznorm,ATynorm,temp;
 29:   Vec            tempJR,tempL;
 30:   Tao            mis;

 33:   mis    = am->subsolverX;
 34:   tempJR = am->workJacobianRight;
 35:   tempL  = am->workLeft;
 36:   /* ATy */
 37:   TaoComputeJacobianEquality(mis, am->y, mis->jacobian_equality, mis->jacobian_equality_pre);
 38:   MatMultTranspose(mis->jacobian_equality,am->y,tempJR);
 39:   VecNorm(tempJR,NORM_2,&ATynorm);
 40:   /* dualres = mu * ||AT(Bz-Bzold)||_2 */
 41:   VecWAXPY(tempJR,-1.,am->Bzold,am->Bz);
 42:   MatMultTranspose(mis->jacobian_equality,tempJR,tempL);
 43:   VecNorm(tempL,NORM_2,&am->dualres);
 44:   am->dualres *= am->mu;

 46:   /* ||Ax||_2, ||Bz||_2 */
 47:   VecNorm(am->Ax,NORM_2,&Axnorm);
 48:   VecNorm(am->Bz,NORM_2,&Bznorm);

 50:   /* Set catol to be catol_admm *  max{||Ax||,||Bz||,||c||} *
 51:    * Set gatol to be gatol_admm *  ||A^Ty|| *
 52:    * while cnorm is ||r_k||_2, and gnorm is ||d_k||_2 */
 53:   temp = am->catol_admm * PetscMax(Axnorm, (!am->const_norm) ? Bznorm : PetscMax(Bznorm,am->const_norm));
 54:   TaoSetConstraintTolerances(tao,temp,PETSC_DEFAULT);
 55:   TaoSetTolerances(tao, am->gatol_admm*ATynorm, PETSC_DEFAULT,PETSC_DEFAULT);
 56:   return(0);
 57: }

 59: /* Penaly Update for Adaptive ADMM. */
 60: static PetscErrorCode AdaptiveADMMPenaltyUpdate(Tao tao)
 61: {
 62:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;
 64:   PetscReal      ydiff_norm, yhatdiff_norm, Axdiff_norm, Bzdiff_norm, Axyhat, Bzy, a_sd, a_mg, a_k, b_sd, b_mg, b_k;
 65:   PetscBool      hflag, gflag;
 66:   Vec            tempJR,tempJR2;

 69:   tempJR  = am->workJacobianRight;
 70:   tempJR2 = am->workJacobianRight2;
 71:   hflag   = PETSC_FALSE;
 72:   gflag   = PETSC_FALSE;
 73:   a_k     = -1;
 74:   b_k     = -1;

 76:   VecWAXPY(tempJR,-1.,am->Axold,am->Ax);
 77:   VecWAXPY(tempJR2,-1.,am->yhatold,am->yhat);
 78:   VecNorm(tempJR,NORM_2,&Axdiff_norm);
 79:   VecNorm(tempJR2,NORM_2,&yhatdiff_norm);
 80:   VecDot(tempJR,tempJR2,&Axyhat);

 82:   VecWAXPY(tempJR,-1.,am->Bz0,am->Bz);
 83:   VecWAXPY(tempJR2,-1.,am->y,am->y0);
 84:   VecNorm(tempJR,NORM_2,&Bzdiff_norm);
 85:   VecNorm(tempJR2,NORM_2,&ydiff_norm);
 86:   VecDot(tempJR,tempJR2,&Bzy);

 88:   if (Axyhat > am->orthval*Axdiff_norm*yhatdiff_norm + am->mueps) {
 89:     hflag = PETSC_TRUE;
 90:     a_sd  = PetscSqr(yhatdiff_norm)/Axyhat; /* alphaSD */
 91:     a_mg  = Axyhat/PetscSqr(Axdiff_norm);   /* alphaMG */
 92:     a_k   = (a_mg/a_sd) > 0.5 ? a_mg : a_sd - 0.5*a_mg;
 93:   }
 94:   if (Bzy > am->orthval*Bzdiff_norm*ydiff_norm + am->mueps) {
 95:     gflag = PETSC_TRUE;
 96:     b_sd  = PetscSqr(ydiff_norm)/Bzy;  /* betaSD */
 97:     b_mg  = Bzy/PetscSqr(Bzdiff_norm); /* betaMG */
 98:     b_k   = (b_mg/b_sd) > 0.5 ? b_mg : b_sd - 0.5*b_mg;
 99:   }
100:   am->muold = am->mu;
101:   if (gflag && hflag) {
102:     am->mu = PetscSqrtReal(a_k*b_k);
103:   } else if (hflag) {
104:     am->mu = a_k;
105:   } else if (gflag) {
106:     am->mu = b_k;
107:   }
108:   if (am->mu > am->muold) {
109:     am->mu = am->muold;
110:   }
111:   if (am->mu < am->mumin) {
112:     am->mu = am->mumin;
113:   }
114:   return(0);
115: }

117: static PetscErrorCode  TaoADMMSetRegularizerType_ADMM(Tao tao, TaoADMMRegularizerType type)
118: {
119:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

122:   am->regswitch = type;
123:   return(0);
124: }

126: static PetscErrorCode  TaoADMMGetRegularizerType_ADMM(Tao tao, TaoADMMRegularizerType *type)
127: {
128:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

131:   *type = am->regswitch;
132:   return(0);
133: }

135: static PetscErrorCode  TaoADMMSetUpdateType_ADMM(Tao tao, TaoADMMUpdateType type)
136: {
137:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

140:   am->update = type;
141:   return(0);
142: }

144: static PetscErrorCode  TaoADMMGetUpdateType_ADMM(Tao tao, TaoADMMUpdateType *type)
145: {
146:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

149:   *type = am->update;
150:   return(0);
151: }

153: /* This routine updates Jacobians with new x,z vectors,
154:  * and then updates Ax and Bz vectors, then computes updated residual vector*/
155: static PetscErrorCode ADMMUpdateConstraintResidualVector(Tao tao, Vec x, Vec z, Vec Ax, Vec Bz, Vec residual)
156: {
157:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;
158:   Tao            mis,reg;

162:   mis  = am->subsolverX;
163:   reg  = am->subsolverZ;
164:   TaoComputeJacobianEquality(mis, x, mis->jacobian_equality, mis->jacobian_equality_pre);
165:   MatMult(mis->jacobian_equality,x,Ax);
166:   TaoComputeJacobianEquality(reg, z, reg->jacobian_equality, reg->jacobian_equality_pre);
167:   MatMult(reg->jacobian_equality,z,Bz);

169:   VecWAXPY(residual,1.,Bz,Ax);
170:   if (am->constraint != NULL) {
171:     VecAXPY(residual,-1.,am->constraint);
172:   }
173:   return(0);
174: }

176: /* Updates Augmented Lagrangians to given routines *
177:  * For subsolverX, routine needs to be ComputeObjectiveAndGraidnet
178:  * Separate Objective and Gradient routines are not supported.  */
179: static PetscErrorCode SubObjGradUpdate(Tao tao, Vec x, PetscReal *f, Vec g, void *ptr)
180: {
181:   Tao            parent = (Tao)ptr;
182:   TAO_ADMM       *am    = (TAO_ADMM*)parent->data;
184:   PetscReal      temp,temp2;
185:   Vec            tempJR;

188:   tempJR = am->workJacobianRight;
189:   ADMMUpdateConstraintResidualVector(parent, x, am->subsolverZ->solution, am->Ax, am->Bz, am->residual);
190:   (*am->ops->misfitobjgrad)(am->subsolverX,x,f,g,am->misfitobjgradP);

192:   am->last_misfit_val = *f;
193:   /* Objective  Add + yT(Ax+Bz-c) + mu/2*||Ax+Bz-c||_2^2 */
194:   VecTDot(am->residual,am->y,&temp);
195:   VecTDot(am->residual,am->residual,&temp2);
196:   *f   += temp + (am->mu/2)*temp2;

198:   /* Gradient. Add + mu*AT(Ax+Bz-c) + yTA*/
199:   MatMultTranspose(tao->jacobian_equality,am->residual,tempJR);
200:   VecAXPY(g,am->mu,tempJR);
201:   MatMultTranspose(tao->jacobian_equality,am->y,tempJR);
202:   VecAXPY(g,1.,tempJR);
203:   return(0);
204: }

206: /* Updates Augmented Lagrangians to given routines
207:  * For subsolverZ, routine needs to be ComputeObjectiveAndGraidnet
208:  * Separate Objective and Gradient routines are not supported.  */
209: static PetscErrorCode RegObjGradUpdate(Tao tao, Vec z, PetscReal *f, Vec g, void *ptr)
210: {
211:   Tao            parent = (Tao)ptr;
212:   TAO_ADMM       *am    = (TAO_ADMM*)parent->data;
214:   PetscReal      temp,temp2;
215:   Vec            tempJR;

218:   tempJR = am->workJacobianRight;
219:   ADMMUpdateConstraintResidualVector(parent, am->subsolverX->solution, z, am->Ax, am->Bz, am->residual);
220:   (*am->ops->regobjgrad)(am->subsolverZ,z,f,g,am->regobjgradP);
221:   am->last_reg_val= *f;
222:   /* Objective  Add  + yT(Ax+Bz-c) + mu/2*||Ax+Bz-c||_2^2 */
223:   VecTDot(am->residual,am->y,&temp);
224:   VecTDot(am->residual,am->residual,&temp2);
225:   *f   += temp + (am->mu/2)*temp2;

227:   /* Gradient. Add + mu*BT(Ax+Bz-c) + yTB*/
228:   MatMultTranspose(am->subsolverZ->jacobian_equality,am->residual,tempJR);
229:   VecAXPY(g,am->mu,tempJR);
230:   MatMultTranspose(am->subsolverZ->jacobian_equality,am->y,tempJR);
231:   VecAXPY(g,1.,tempJR);
232:   return(0);
233: }

235: /* Computes epsilon padded L1 norm lambda*sum(sqrt(x^2+eps^2)-eps */
236: static PetscErrorCode ADMML1EpsilonNorm(Tao tao, Vec x, PetscReal eps, PetscReal *norm)
237: {
238:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;
239:   PetscInt       N;

243:   VecGetSize(am->workLeft,&N);
244:   VecPointwiseMult(am->workLeft,x,x);
245:   VecShift(am->workLeft,am->l1epsilon*am->l1epsilon);
246:   VecSqrtAbs(am->workLeft);
247:   VecSum(am->workLeft,norm);
248:   *norm += N*am->l1epsilon;
249:   *norm *= am->lambda;
250:   return(0);
251: }

253: static PetscErrorCode ADMMInternalHessianUpdate(Mat H, Mat Constraint, PetscBool Identity, void *ptr)
254: {
255:   TAO_ADMM       *am = (TAO_ADMM*)ptr;

259:   switch (am->update) {
260:   case (TAO_ADMM_UPDATE_BASIC):
261:     break;
262:   case (TAO_ADMM_UPDATE_ADAPTIVE):
263:   case (TAO_ADMM_UPDATE_ADAPTIVE_RELAXED):
264:     if (H && (am->muold != am->mu)) {
265:       if (!Identity) {
266:         MatAXPY(H,am->mu-am->muold,Constraint,DIFFERENT_NONZERO_PATTERN);
267:       } else {
268:         MatShift(H,am->mu-am->muold);
269:       }
270:     }
271:     break;
272:   }
273:   return(0);
274: }

276: /* Updates Hessian - adds second derivative of augmented Lagrangian
277:  * H \gets H + \rho*ATA
278:  * Here, \rho does not change in TAO_ADMM_UPDATE_BASIC - thus no-op
279:  * For ADAPTAIVE,ADAPTIVE_RELAXED,
280:  * H \gets H + (\rho-\rhoold)*ATA
281:  * Here, we assume that A is linear constraint i.e., doesnt change.
282:  * Thus, for both ADAPTIVE, and RELAXED, ATA matrix is pre-set (except for A=I (null case)) see TaoSetUp_ADMM */
283: static PetscErrorCode SubHessianUpdate(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
284: {
285:   Tao            parent = (Tao)ptr;
286:   TAO_ADMM       *am    = (TAO_ADMM*)parent->data;

290:   if (am->Hxchange) {
291:     /* Case where Hessian gets updated with respect to x vector input. */
292:     (*am->ops->misfithess)(am->subsolverX,x,H,Hpre,am->misfithessP);
293:     ADMMInternalHessianUpdate(am->subsolverX->hessian,am->ATA,am->xJI,am);
294:   } else if (am->Hxbool) {
295:     /* Hessian doesn't get updated. H(x) = c */
296:     /* Update Lagrangian only only per TAO call */
297:     ADMMInternalHessianUpdate(am->subsolverX->hessian,am->ATA,am->xJI,am);
298:     am->Hxbool = PETSC_FALSE;
299:   }
300:   return(0);
301: }

303: /* Same as SubHessianUpdate, except for B matrix instead of A matrix */
304: static PetscErrorCode RegHessianUpdate(Tao tao, Vec z, Mat H, Mat Hpre, void *ptr)
305: {
306:   Tao            parent = (Tao)ptr;
307:   TAO_ADMM       *am    = (TAO_ADMM*)parent->data;


312:   if (am->Hzchange) {
313:     /* Case where Hessian gets updated with respect to x vector input. */
314:     (*am->ops->reghess)(am->subsolverZ,z,H,Hpre,am->reghessP);
315:     ADMMInternalHessianUpdate(am->subsolverZ->hessian,am->BTB,am->zJI,am);
316:   } else if (am->Hzbool) {
317:     /* Hessian doesn't get updated. H(x) = c */
318:     /* Update Lagrangian only only per TAO call */
319:     ADMMInternalHessianUpdate(am->subsolverZ->hessian,am->BTB,am->zJI,am);
320:     am->Hzbool = PETSC_FALSE;
321:   }
322:   return(0);
323: }

325: /* Shell Matrix routine for A matrix.
326:  * This gets used when user puts NULL for
327:  * TaoSetJacobianEqualityRoutine(tao, NULL,NULL, ...)
328:  * Essentially sets A=I*/
329: static PetscErrorCode JacobianIdentity(Mat mat,Vec in,Vec out)
330: {

334:   VecCopy(in,out);
335:   return(0);
336: }

338: /* Shell Matrix routine for B matrix.
339:  * This gets used when user puts NULL for
340:  * TaoADMMSetRegularizerConstraintJacobian(tao, NULL,NULL, ...)
341:  * Sets B=-I */
342: static PetscErrorCode JacobianIdentityB(Mat mat,Vec in,Vec out)
343: {

347:   VecCopy(in,out);
348:   VecScale(out,-1.);
349:   return(0);
350: }

352: /* Solve f(x) + g(z) s.t. Ax + Bz = c */
353: static PetscErrorCode TaoSolve_ADMM(Tao tao)
354: {
355:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;
357:   PetscInt       N;
358:   PetscReal      reg_func;
359:   PetscBool      is_reg_shell;
360:   Vec            tempL;

363:   if (am->regswitch != TAO_ADMM_REGULARIZER_SOFT_THRESH) {
364:     if (!am->subsolverX->ops->computejacobianequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Must call TaoADMMSetMisfitConstraintJacobian() first");
365:     if (!am->subsolverZ->ops->computejacobianequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Must call TaoADMMSetRegularizerConstraintJacobian() first");
366:     if (am->constraint != NULL) {
367:       VecNorm(am->constraint,NORM_2,&am->const_norm);
368:     }
369:   }
370:   tempL = am->workLeft;
371:   VecGetSize(tempL,&N);

373:   if (am->Hx && am->ops->misfithess) {
374:     TaoSetHessianRoutine(am->subsolverX, am->Hx, am->Hx, SubHessianUpdate, tao);
375:   }

377:   if (!am->zJI) {
378:     /* Currently, B is assumed to be a linear system, i.e., not getting updated*/
379:     MatTransposeMatMult(am->JB,am->JB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(am->BTB));
380:   }
381:   if (!am->xJI) {
382:     /* Currently, A is assumed to be a linear system, i.e., not getting updated*/
383:     MatTransposeMatMult(am->subsolverX->jacobian_equality,am->subsolverX->jacobian_equality,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(am->ATA));
384:   }

386:   is_reg_shell = PETSC_FALSE;

388:   PetscObjectTypeCompare((PetscObject)am->subsolverZ, TAOSHELL, &is_reg_shell);

390:   if (!is_reg_shell) {
391:     switch (am->regswitch) {
392:     case (TAO_ADMM_REGULARIZER_USER):
393:       if (!am->ops->regobjgrad) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Must call TaoADMMSetRegularizerObjectiveAndGradientRoutine() first if one wishes to use TAO_ADMM_REGULARIZER_USER with non-TAOSHELL type");
394:       break;
395:     case (TAO_ADMM_REGULARIZER_SOFT_THRESH):
396:       /* Soft Threshold. */
397:       break;
398:     }
399:     if (am->ops->regobjgrad) {
400:       TaoSetObjectiveAndGradientRoutine(am->subsolverZ, RegObjGradUpdate, tao);
401:     }
402:     if (am->Hz && am->ops->reghess) {
403:       TaoSetHessianRoutine(am->subsolverZ, am->Hz, am->Hzpre, RegHessianUpdate, tao);
404:     }
405:   }

407:   switch (am->update) {
408:   case TAO_ADMM_UPDATE_BASIC:
409:     if (am->subsolverX->hessian) {
410:       /* In basic case, Hessian does not get updated w.r.t. to spectral penalty
411:        * Here, when A is set, i.e., am->xJI, add mu*ATA to Hessian*/
412:       if (!am->xJI) {
413:         MatAXPY(am->subsolverX->hessian,am->mu,am->ATA,DIFFERENT_NONZERO_PATTERN);
414:       } else {
415:         MatShift(am->subsolverX->hessian,am->mu);
416:       }
417:     }
418:     if (am->subsolverZ->hessian && am->regswitch == TAO_ADMM_REGULARIZER_USER) {
419:       if (am->regswitch == TAO_ADMM_REGULARIZER_USER && !am->zJI) {
420:         MatAXPY(am->subsolverZ->hessian,am->mu,am->BTB,DIFFERENT_NONZERO_PATTERN);
421:       } else {
422:         MatShift(am->subsolverZ->hessian,am->mu);
423:       }
424:     }
425:     break;
426:   case TAO_ADMM_UPDATE_ADAPTIVE:
427:   case TAO_ADMM_UPDATE_ADAPTIVE_RELAXED:
428:     break;
429:   }

431:   PetscCitationsRegister(citation,&cited);
432:   tao->reason = TAO_CONTINUE_ITERATING;

434:   while (tao->reason == TAO_CONTINUE_ITERATING) {
435:     if (tao->ops->update) {
436:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
437:     }
438:     VecCopy(am->Bz, am->Bzold);

440:     /* x update */
441:     TaoSolve(am->subsolverX);
442:     TaoComputeJacobianEquality(am->subsolverX, am->subsolverX->solution, am->subsolverX->jacobian_equality, am->subsolverX->jacobian_equality_pre);
443:     MatMult(am->subsolverX->jacobian_equality, am->subsolverX->solution,am->Ax);

445:     am->Hxbool = PETSC_TRUE;

447:     /* z update */
448:     switch (am->regswitch) {
449:     case TAO_ADMM_REGULARIZER_USER:
450:       TaoSolve(am->subsolverZ);
451:       break;
452:     case TAO_ADMM_REGULARIZER_SOFT_THRESH:
453:       /* L1 assumes A,B jacobians are identity nxn matrix */
454:       VecWAXPY(am->workJacobianRight,1/am->mu,am->y,am->Ax);
455:       TaoSoftThreshold(am->workJacobianRight,-am->lambda/am->mu,am->lambda/am->mu,am->subsolverZ->solution);
456:       break;
457:     }
458:     am->Hzbool = PETSC_TRUE;
459:     /* Returns Ax + Bz - c with updated Ax,Bz vectors */
460:     ADMMUpdateConstraintResidualVector(tao, am->subsolverX->solution, am->subsolverZ->solution, am->Ax, am->Bz, am->residual);
461:     /* Dual variable, y += y + mu*(Ax+Bz-c) */
462:     VecWAXPY(am->y, am->mu, am->residual, am->yold);

464:     /* stopping tolerance update */
465:     TaoADMMToleranceUpdate(tao);

467:     /* Updating Spectral Penalty */
468:     switch (am->update) {
469:     case TAO_ADMM_UPDATE_BASIC:
470:       am->muold = am->mu;
471:       break;
472:     case TAO_ADMM_UPDATE_ADAPTIVE:
473:     case TAO_ADMM_UPDATE_ADAPTIVE_RELAXED:
474:       if (tao->niter == 0) {
475:         VecCopy(am->y, am->y0);
476:         VecWAXPY(am->residual, 1., am->Ax, am->Bzold);
477:         if (am->constraint) {
478:           VecAXPY(am->residual, -1., am->constraint);
479:         }
480:         VecWAXPY(am->yhatold, -am->mu, am->residual, am->yold);
481:         VecCopy(am->Ax, am->Axold);
482:         VecCopy(am->Bz, am->Bz0);
483:         am->muold = am->mu;
484:       } else if (tao->niter % am->T == 1) {
485:         /* we have compute Bzold in a previous iteration, and we computed Ax above */
486:         VecWAXPY(am->residual, 1., am->Ax, am->Bzold);
487:         if (am->constraint) {
488:           VecAXPY(am->residual, -1., am->constraint);
489:         }
490:         VecWAXPY(am->yhat, -am->mu, am->residual, am->yold);
491:         AdaptiveADMMPenaltyUpdate(tao);
492:         VecCopy(am->Ax, am->Axold);
493:         VecCopy(am->Bz, am->Bz0);
494:         VecCopy(am->yhat, am->yhatold);
495:         VecCopy(am->y, am->y0);
496:       } else {
497:         am->muold = am->mu;
498:       }
499:       break;
500:     default:
501:       break;
502:     }
503:     tao->niter++;

505:     /* Calculate original function values. misfit part was done in TaoADMMToleranceUpdate*/
506:     switch (am->regswitch) {
507:     case TAO_ADMM_REGULARIZER_USER:
508:       if (is_reg_shell) {
509:         ADMML1EpsilonNorm(tao,am->subsolverZ->solution,am->l1epsilon,&reg_func);
510:       } else {
511:         (*am->ops->regobjgrad)(am->subsolverZ,am->subsolverX->solution,&reg_func,tempL,am->regobjgradP);
512:       }
513:       break;
514:     case TAO_ADMM_REGULARIZER_SOFT_THRESH:
515:       ADMML1EpsilonNorm(tao,am->subsolverZ->solution,am->l1epsilon,&reg_func);
516:       break;
517:     }
518:     VecCopy(am->y,am->yold);
519:     ADMMUpdateConstraintResidualVector(tao, am->subsolverX->solution, am->subsolverZ->solution, am->Ax, am->Bz, am->residual);
520:     VecNorm(am->residual,NORM_2,&am->resnorm);
521:     TaoLogConvergenceHistory(tao,am->last_misfit_val + reg_func,am->dualres,am->resnorm,tao->ksp_its);

523:     TaoMonitor(tao,tao->niter,am->last_misfit_val + reg_func,am->dualres,am->resnorm,1.0);
524:     (*tao->ops->convergencetest)(tao,tao->cnvP);
525:   }
526:   /* Update vectors */
527:   VecCopy(am->subsolverX->solution,tao->solution);
528:   VecCopy(am->subsolverX->gradient,tao->gradient);
529:   PetscObjectCompose((PetscObject)am->subsolverX,"TaoGetADMMParentTao_ADMM", NULL);
530:   PetscObjectCompose((PetscObject)am->subsolverZ,"TaoGetADMMParentTao_ADMM", NULL);
531:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetRegularizerType_C",NULL);
532:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetRegularizerType_C",NULL);
533:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetUpdateType_C",NULL);
534:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetUpdateType_C",NULL);
535:   return(0);
536: }

538: static PetscErrorCode TaoSetFromOptions_ADMM(PetscOptionItems *PetscOptionsObject,Tao tao)
539: {
540:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

544:   PetscOptionsHead(PetscOptionsObject,"ADMM problem that solves f(x) in a form of f(x) + g(z) subject to x - z = 0. Norm 1 and 2 are supported. Different subsolver routines can be selected. ");
545:   PetscOptionsReal("-tao_admm_regularizer_coefficient","regularizer constant","",am->lambda,&am->lambda,NULL);
546:   PetscOptionsReal("-tao_admm_spectral_penalty","Constant for Augmented Lagrangian term.","",am->mu,&am->mu,NULL);
547:   PetscOptionsReal("-tao_admm_relaxation_parameter","x relaxation parameter for Z update.","",am->gamma,&am->gamma,NULL);
548:   PetscOptionsReal("-tao_admm_tolerance_update_factor","ADMM dynamic tolerance update factor.","",am->tol,&am->tol,NULL);
549:   PetscOptionsReal("-tao_admm_spectral_penalty_update_factor","ADMM spectral penalty update curvature safeguard value.","",am->orthval,&am->orthval,NULL);
550:   PetscOptionsReal("-tao_admm_minimum_spectral_penalty","Set ADMM minimum spectral penalty.","",am->mumin,&am->mumin,NULL);
551:   PetscOptionsEnum("-tao_admm_dual_update","Lagrangian dual update policy","TaoADMMUpdateType",
552:                           TaoADMMUpdateTypes,(PetscEnum)am->update,(PetscEnum*)&am->update,NULL);
553:   PetscOptionsEnum("-tao_admm_regularizer_type","ADMM regularizer update rule","TaoADMMRegularizerType",
554:                           TaoADMMRegularizerTypes,(PetscEnum)am->regswitch,(PetscEnum*)&am->regswitch,NULL);
555:   PetscOptionsTail();
556:   TaoSetFromOptions(am->subsolverX);
557:   if (am->regswitch != TAO_ADMM_REGULARIZER_SOFT_THRESH) {
558:     TaoSetFromOptions(am->subsolverZ);
559:   }
560:   return(0);
561: }

563: static PetscErrorCode TaoView_ADMM(Tao tao,PetscViewer viewer)
564: {
565:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

569:   PetscViewerASCIIPushTab(viewer);
570:   TaoView(am->subsolverX,viewer);
571:   TaoView(am->subsolverZ,viewer);
572:   PetscViewerASCIIPopTab(viewer);
573:   return(0);
574: }

576: static PetscErrorCode TaoSetUp_ADMM(Tao tao)
577: {
578:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;
580:   PetscInt       n,N,M;

583:   VecGetLocalSize(tao->solution,&n);
584:   VecGetSize(tao->solution,&N);
585:   /* If Jacobian is given as NULL, it means Jacobian is identity matrix with size of solution vector */
586:   if (!am->JB) {
587:     am->zJI   = PETSC_TRUE;
588:     MatCreateShell(PETSC_COMM_WORLD,n,n,PETSC_DETERMINE,PETSC_DETERMINE,NULL,&am->JB);
589:     MatShellSetOperation(am->JB,MATOP_MULT,(void (*)(void))JacobianIdentityB);
590:     MatShellSetOperation(am->JB,MATOP_MULT_TRANSPOSE,(void (*)(void))JacobianIdentityB);
591:     MatShellSetOperation(am->JB,MATOP_TRANSPOSE_MAT_MULT,(void (*)(void))JacobianIdentityB);
592:     am->JBpre = am->JB;
593:   }
594:   if (!am->JA) {
595:     am->xJI   = PETSC_TRUE;
596:     MatCreateShell(PETSC_COMM_WORLD,n,n,PETSC_DETERMINE,PETSC_DETERMINE,NULL,&am->JA);
597:     MatShellSetOperation(am->JA,MATOP_MULT,(void (*)(void))JacobianIdentity);
598:     MatShellSetOperation(am->JA,MATOP_MULT_TRANSPOSE,(void (*)(void))JacobianIdentity);
599:     MatShellSetOperation(am->JA,MATOP_TRANSPOSE_MAT_MULT,(void (*)(void))JacobianIdentity);
600:     am->JApre = am->JA;
601:   }
602:   MatCreateVecs(am->JA,NULL,&am->Ax);
603:   if (!tao->gradient) {
604:     VecDuplicate(tao->solution,&tao->gradient);
605:   }
606:   TaoSetInitialVector(am->subsolverX, tao->solution);
607:   if (!am->z) {
608:     VecDuplicate(tao->solution,&am->z);
609:     VecSet(am->z,0.0);
610:   }
611:   TaoSetInitialVector(am->subsolverZ, am->z);
612:   if (!am->workLeft) {
613:     VecDuplicate(tao->solution,&am->workLeft);
614:   }
615:   if (!am->Axold) {
616:     VecDuplicate(am->Ax,&am->Axold);
617:   }
618:   if (!am->workJacobianRight) {
619:     VecDuplicate(am->Ax,&am->workJacobianRight);
620:   }
621:   if (!am->workJacobianRight2) {
622:     VecDuplicate(am->Ax,&am->workJacobianRight2);
623:   }
624:   if (!am->Bz) {
625:     VecDuplicate(am->Ax,&am->Bz);
626:   }
627:   if (!am->Bzold) {
628:     VecDuplicate(am->Ax,&am->Bzold);
629:   }
630:   if (!am->Bz0) {
631:     VecDuplicate(am->Ax,&am->Bz0);
632:   }
633:   if (!am->y) {
634:     VecDuplicate(am->Ax,&am->y);
635:     VecSet(am->y,0.0);
636:   }
637:   if (!am->yold) {
638:     VecDuplicate(am->Ax,&am->yold);
639:     VecSet(am->yold,0.0);
640:   }
641:   if (!am->y0) {
642:     VecDuplicate(am->Ax,&am->y0);
643:     VecSet(am->y0,0.0);
644:   }
645:   if (!am->yhat) {
646:     VecDuplicate(am->Ax,&am->yhat);
647:     VecSet(am->yhat,0.0);
648:   }
649:   if (!am->yhatold) {
650:     VecDuplicate(am->Ax,&am->yhatold);
651:     VecSet(am->yhatold,0.0);
652:   }
653:   if (!am->residual) {
654:     VecDuplicate(am->Ax,&am->residual);
655:     VecSet(am->residual,0.0);
656:   }
657:   if (!am->constraint) {
658:     am->constraint = NULL;
659:   } else {
660:     VecGetSize(am->constraint,&M);
661:     if (M != N) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"Solution vector and constraint vector must be of same size!");
662:   }

664:   /* Save changed tao tolerance for adaptive tolerance */
665:   if (tao->gatol_changed) {
666:     am->gatol_admm = tao->gatol;
667:   }
668:   if (tao->catol_changed) {
669:     am->catol_admm = tao->catol;
670:   }

672:   /*Update spectral and dual elements to X subsolver */
673:   TaoSetObjectiveAndGradientRoutine(am->subsolverX, SubObjGradUpdate, tao);
674:   TaoSetJacobianEqualityRoutine(am->subsolverX,am->JA,am->JApre, am->ops->misfitjac, am->misfitjacobianP);
675:   TaoSetJacobianEqualityRoutine(am->subsolverZ,am->JB,am->JBpre, am->ops->regjac, am->regjacobianP);
676:   return(0);
677: }

679: static PetscErrorCode TaoDestroy_ADMM(Tao tao)
680: {
681:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

685:   VecDestroy(&am->z);
686:   VecDestroy(&am->Ax);
687:   VecDestroy(&am->Axold);
688:   VecDestroy(&am->Bz);
689:   VecDestroy(&am->Bzold);
690:   VecDestroy(&am->Bz0);
691:   VecDestroy(&am->residual);
692:   VecDestroy(&am->y);
693:   VecDestroy(&am->yold);
694:   VecDestroy(&am->y0);
695:   VecDestroy(&am->yhat);
696:   VecDestroy(&am->yhatold);
697:   VecDestroy(&am->workLeft);
698:   VecDestroy(&am->workJacobianRight);
699:   VecDestroy(&am->workJacobianRight2);

701:   MatDestroy(&am->JA);
702:   MatDestroy(&am->JB);
703:   if (!am->xJI) {
704:     MatDestroy(&am->JApre);
705:   }
706:   if (!am->zJI) {
707:     MatDestroy(&am->JBpre);
708:   }
709:   if (am->Hx) {
710:     MatDestroy(&am->Hx);
711:     MatDestroy(&am->Hxpre);
712:   }
713:   if (am->Hz) {
714:     MatDestroy(&am->Hz);
715:     MatDestroy(&am->Hzpre);
716:   }
717:   MatDestroy(&am->ATA);
718:   MatDestroy(&am->BTB);
719:   TaoDestroy(&am->subsolverX);
720:   TaoDestroy(&am->subsolverZ);
721:   am->parent = NULL;
722:   PetscFree(tao->data);
723:   return(0);
724: }

726: /*MC

728:   TAOADMM - Alternating direction method of multipliers method fo solving linear problems with
729:             constraints. in a min_x f(x) + g(z)  s.t. Ax+Bz=c.
730:             This algorithm employs two sub Tao solvers, of which type can be specificed
731:             by the user. User need to provide ObjectiveAndGradient routine, and/or HessianRoutine for both subsolvers.
732:             Hessians can be given boolean flag determining whether they change with respect to a input vector. This can be set via
733:             TaoADMMSet{Misfit,Regularizer}HessianChangeStatus.
734:             Second subsolver does support TAOSHELL. It should be noted that L1-norm is used for objective value for TAOSHELL type.
735:             There is option to set regularizer option, and currently soft-threshold is implemented. For spectral penalty update,
736:             currently there are baisc option and adaptive option.
737:             Constraint is set at Ax+Bz=c, and A and B can be set with TaoADMMSet{Misfit,Regularizer}ConstraintJacobian.
738:             c can be set with TaoADMMSetConstraintVectorRHS.
739:             The user can also provide regularizer weight for second subsolver.

741:   References:
742: .    1. - Xu, Zheng and Figueiredo, Mario A. T. and Yuan, Xiaoming and Studer, Christoph and Goldstein, Tom
743:           "Adaptive Relaxed ADMM: Convergence Theory and Practical Implementation"
744:           The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), July, 2017.

746:   Options Database Keys:
747: + -tao_admm_regularizer_coefficient        - regularizer constant (default 1.e-6)
748: . -tao_admm_spectral_penalty               - Constant for Augmented Lagrangian term (default 1.)
749: . -tao_admm_relaxation_parameter           - relaxation parameter for Z update (default 1.)
750: . -tao_admm_tolerance_update_factor        - ADMM dynamic tolerance update factor (default 1.e-12)
751: . -tao_admm_spectral_penalty_update_factor - ADMM spectral penalty update curvature safeguard value (default 0.2)
752: . -tao_admm_minimum_spectral_penalty       - Set ADMM minimum spectral penalty (default 0)
753: . -tao_admm_dual_update                    - Lagrangian dual update policy ("basic","adaptive","adaptive-relaxed") (default "basic")
754: - -tao_admm_regularizer_type               - ADMM regularizer update rule ("user","soft-threshold") (default "soft-threshold")

756:   Level: beginner

758: .seealso: TaoADMMSetMisfitHessianChangeStatus(), TaoADMMSetRegHessianChangeStatus(), TaoADMMGetSpectralPenalty(), 
759:           TaoADMMGetMisfitSubsolver(), TaoADMMGetRegularizationSubsolver(), TaoADMMSetConstraintVectorRHS(),
760:           TaoADMMSetMinimumSpectralPenalty(), TaoADMMSetRegularizerCoefficient(),
761:           TaoADMMSetRegularizerConstraintJacobian(), TaoADMMSetMisfitConstraintJacobian(),
762:           TaoADMMSetMisfitObjectiveAndGradientRoutine(), TaoADMMSetMisfitHessianRoutine(),
763:           TaoADMMSetRegularizerObjectiveAndGradientRoutine(), TaoADMMSetRegularizerHessianRoutine(),
764:           TaoGetADMMParentTao(), TaoADMMGetDualVector(), TaoADMMSetRegularizerType(),
765:           TaoADMMGetRegularizerType(), TaoADMMSetUpdateType(), TaoADMMGetUpdateType()  
766: M*/

768: PETSC_EXTERN PetscErrorCode TaoCreate_ADMM(Tao tao)
769: {
770:   TAO_ADMM       *am;

774:   PetscNewLog(tao,&am);

776:   tao->ops->destroy        = TaoDestroy_ADMM;
777:   tao->ops->setup          = TaoSetUp_ADMM;
778:   tao->ops->setfromoptions = TaoSetFromOptions_ADMM;
779:   tao->ops->view           = TaoView_ADMM;
780:   tao->ops->solve          = TaoSolve_ADMM;

782:   tao->data           = (void*)am;
783:   am->l1epsilon       = 1e-6;
784:   am->lambda          = 1e-4;
785:   am->mu              = 1.;
786:   am->muold           = 0.;
787:   am->mueps           = PETSC_MACHINE_EPSILON;
788:   am->mumin           = 0.;
789:   am->orthval         = 0.2;
790:   am->T               = 2;
791:   am->parent          = tao;
792:   am->update          = TAO_ADMM_UPDATE_BASIC;
793:   am->regswitch       = TAO_ADMM_REGULARIZER_SOFT_THRESH;
794:   am->tol             = PETSC_SMALL;
795:   am->const_norm      = 0;
796:   am->resnorm         = 0;
797:   am->dualres         = 0;
798:   am->ops->regobjgrad = 0;
799:   am->ops->reghess    = 0;
800:   am->gamma           = 1;
801:   am->regobjgradP     = 0;
802:   am->reghessP        = 0;
803:   am->gatol_admm      = 1e-8;
804:   am->catol_admm      = 0;
805:   am->Hxchange        = PETSC_TRUE;
806:   am->Hzchange        = PETSC_TRUE;
807:   am->Hzbool          = PETSC_TRUE;
808:   am->Hxbool          = PETSC_TRUE;

810:   TaoCreate(PetscObjectComm((PetscObject)tao),&am->subsolverX);
811:   TaoSetOptionsPrefix(am->subsolverX,"misfit_");
812:   PetscObjectIncrementTabLevel((PetscObject)am->subsolverX,(PetscObject)tao,1);
813:   TaoCreate(PetscObjectComm((PetscObject)tao),&am->subsolverZ);
814:   TaoSetOptionsPrefix(am->subsolverZ,"reg_");
815:   PetscObjectIncrementTabLevel((PetscObject)am->subsolverZ,(PetscObject)tao,1);

817:   TaoSetType(am->subsolverX,TAONLS);
818:   TaoSetType(am->subsolverZ,TAONLS);
819:   PetscObjectCompose((PetscObject)am->subsolverX,"TaoGetADMMParentTao_ADMM", (PetscObject) tao);
820:   PetscObjectCompose((PetscObject)am->subsolverZ,"TaoGetADMMParentTao_ADMM", (PetscObject) tao);
821:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetRegularizerType_C",TaoADMMSetRegularizerType_ADMM);
822:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetRegularizerType_C",TaoADMMGetRegularizerType_ADMM);
823:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMSetUpdateType_C",TaoADMMSetUpdateType_ADMM);
824:   PetscObjectComposeFunction((PetscObject)tao,"TaoADMMGetUpdateType_C",TaoADMMGetUpdateType_ADMM);
825:   return(0);
826: }

828: /*@
829:   TaoADMMSetMisfitHessianChangeStatus - Set boolean that determines  whether Hessian matrix of misfit subsolver changes with respect to input vector.

831:   Collective on Tao

833:   Input Parameters:
834: +  tao - the Tao solver context. 
835: -  b - the Hessian matrix change status boolean, PETSC_FALSE  when the Hessian matrix does not change, TRUE otherwise.

837:   Level: advanced
838:   
839: .seealso: TAOADMM

841: @*/
842: PetscErrorCode TaoADMMSetMisfitHessianChangeStatus(Tao tao, PetscBool b)
843: {
844:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

847:   am->Hxchange = b;
848:   return(0);
849: }

851: /*@
852:   TaoADMMSetRegHessianChangeStatus - Set boolean that determines whether Hessian matrix of regularization subsolver changes with respect to input vector.
853:   

855:   Collective on Tao

857:   Input Parameters:
858: +  tao - the Tao solver context
859: -  b - the Hessian matrix change status boolean, PETSC_FALSE when the Hessian matrix does not change, TRUE otherwise.

861:   Level: advanced
862:   
863: .seealso: TAOADMM

865: @*/
866: PetscErrorCode TaoADMMSetRegHessianChangeStatus(Tao tao, PetscBool b)
867: {
868:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

871:   am->Hzchange = b;
872:   return(0);
873: }

875: /*@
876:   TaoADMMSetSpectralPenalty - Set the spectral penalty (mu) value

878:   Collective on Tao

880:   Input Parameters:
881: +  tao - the Tao solver context
882: -  mu - spectral penalty

884:   Level: advanced

886: .seealso: TaoADMMSetMinimumSpectralPenalty(), TAOADMM  
887: @*/
888: PetscErrorCode TaoADMMSetSpectralPenalty(Tao tao, PetscReal mu)
889: {
890:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

893:   am->mu = mu;
894:   return(0);
895: }

897: /*@
898:   TaoADMMGetSpectralPenalty - Get the spectral penalty (mu) value

900:   Collective on Tao

902:   Input Parameter:
903: .  tao - the Tao solver context

905:   Output Parameter:
906: .  mu - spectral penalty

908:   Level: advanced

910: .seealso: TaoADMMSetMinimumSpectralPenalty(), TaoADMMSetSpectralPenalty(), TAOADMM
911: @*/
912: PetscErrorCode TaoADMMGetSpectralPenalty(Tao tao, PetscReal *mu)
913: {
914:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

919:   *mu = am->mu;
920:   return(0);
921: }

923: /*@
924:   TaoADMMGetMisfitSubsolver - Get the pointer to the misfit subsolver inside ADMM

926:   Collective on Tao

928:   Input Parameter:
929: .  tao - the Tao solver context

931:    Output Parameter:
932: .  misfit - the Tao subsolver context

934:   Level: advanced
935:   
936: .seealso: TAOADMM

938: @*/
939: PetscErrorCode TaoADMMGetMisfitSubsolver(Tao tao, Tao *misfit)
940: {
941:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

944:   *misfit = am->subsolverX;
945:   return(0);
946: }

948: /*@
949:   TaoADMMGetRegularizationSubsolver - Get the pointer to the regularization subsolver inside ADMM

951:   Collective on Tao

953:   Input Parameter:
954: .  tao - the Tao solver context

956:   Output Parameter:
957: .  reg - the Tao subsolver context

959:   Level: advanced
960:   
961: .seealso: TAOADMM

963: @*/
964: PetscErrorCode TaoADMMGetRegularizationSubsolver(Tao tao, Tao *reg)
965: {
966:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

969:   *reg = am->subsolverZ;
970:   return(0);
971: }

973: /*@
974:   TaoADMMSetConstraintVectorRHS - Set the RHS constraint vector for ADMM

976:   Collective on Tao

978:   Input Parameters:
979: + tao - the Tao solver context
980: - c - RHS vector

982:   Level: advanced
983:   
984: .seealso: TAOADMM

986: @*/
987: PetscErrorCode TaoADMMSetConstraintVectorRHS(Tao tao, Vec c)
988: {
989:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

992:   am->constraint = c;
993:   return(0);
994: }

996: /*@
997:   TaoADMMSetMinimumSpectralPenalty - Set the minimum value for the spectral penalty

999:   Collective on Tao

1001:   Input Parameters:
1002: +  tao - the Tao solver context
1003: -  mu  - minimum spectral penalty value

1005:   Level: advanced

1007: .seealso: TaoADMMGetSpectralPenalty(), TAOADMM  
1008: @*/
1009: PetscErrorCode TaoADMMSetMinimumSpectralPenalty(Tao tao, PetscReal mu)
1010: {
1011:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

1014:   am->mumin= mu;
1015:   return(0);
1016: }

1018: /*@
1019:   TaoADMMSetRegularizerCoefficient - Set the regularization coefficient lambda for L1 norm regularization case

1021:   Collective on Tao

1023:   Input Parameters:
1024: +  tao - the Tao solver context
1025: -  lambda - L1-norm regularizer coefficient

1027:   Level: advanced

1029: .seealso: TaoADMMSetMisfitConstraintJacobian(), TaoADMMSetRegularizerConstraintJacobian(), TAOADMM

1031: @*/
1032: PetscErrorCode TaoADMMSetRegularizerCoefficient(Tao tao, PetscReal lambda)
1033: {
1034:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

1037:   am->lambda = lambda;
1038:   return(0);
1039: }

1041: /*@C
1042:   TaoADMMSetMisfitConstraintJacobian - Set the constraint matrix B for the ADMM algorithm. Matrix B constrains the z variable.

1044:   Collective on Tao

1046:   Input Parameters:
1047: + tao - the Tao solver context
1048: . J - user-created regularizer constraint Jacobian matrix
1049: . Jpre - user-created regularizer Jacobian constraint preconditioner matrix
1050: . func - function pointer for the regularizer constraint Jacobian update function
1051: - ctx - user context for the regularizer Hessian

1053:   Level: advanced

1055: .seealso: TaoADMMSetRegularizerCoefficient(), TaoADMMSetRegularizerConstraintJacobian(), TAOADMM

1057: @*/
1058: PetscErrorCode TaoADMMSetMisfitConstraintJacobian(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1059: {
1060:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

1065:   if (J) {
1068:   }
1069:   if (Jpre) {
1072:   }
1073:   if (ctx)  am->misfitjacobianP = ctx;
1074:   if (func) am->ops->misfitjac  = func;

1076:   if (J) {
1077:     PetscObjectReference((PetscObject)J);
1078:     MatDestroy(&am->JA);
1079:     am->JA = J;
1080:   }
1081:   if (Jpre) {
1082:     PetscObjectReference((PetscObject)Jpre);
1083:     MatDestroy(&am->JApre);
1084:     am->JApre = Jpre;
1085:   }
1086:   return(0);
1087: }

1089: /*@C
1090:   TaoADMMSetRegularizerConstraintJacobian - Set the constraint matrix B for ADMM algorithm. Matrix B constraints z variable.

1092:   Collective on Tao

1094:   Input Parameters:
1095: + tao - the Tao solver context
1096: . J - user-created regularizer constraint Jacobian matrix
1097: . Jpre - user-created regularizer Jacobian constraint preconditioner matrix
1098: . func - function pointer for the regularizer constraint Jacobian update function
1099: - ctx - user context for the regularizer Hessian

1101:   Level: advanced

1103: .seealso: TaoADMMSetRegularizerCoefficient(), TaoADMMSetMisfitConstraintJacobian(), TAOADMM

1105: @*/
1106: PetscErrorCode TaoADMMSetRegularizerConstraintJacobian(Tao tao, Mat J, Mat Jpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1107: {
1108:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

1113:   if (J) {
1116:   }
1117:   if (Jpre) {
1120:   }
1121:   if (ctx)  am->regjacobianP = ctx;
1122:   if (func) am->ops->regjac  = func;

1124:   if (J) {
1125:     PetscObjectReference((PetscObject)J);
1126:     MatDestroy(&am->JB);
1127:     am->JB = J;
1128:   }
1129:   if (Jpre) {
1130:     PetscObjectReference((PetscObject)Jpre);
1131:     MatDestroy(&am->JBpre);
1132:     am->JBpre = Jpre;
1133:   }
1134:   return(0);
1135: }

1137: /*@C
1138:    TaoADMMSetMisfitObjectiveAndGradientRoutine - Sets the user-defined misfit call-back function

1140:    Collective on tao

1142:    Input Parameters:
1143: +    tao - the Tao context
1144: .    func - function pointer for the misfit value and gradient evaluation
1145: -    ctx - user context for the misfit

1147:    Level: advanced

1149: .seealso: TAOADMM

1151: @*/
1152: PetscErrorCode TaoADMMSetMisfitObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao, Vec, PetscReal*, Vec, void*), void *ctx)
1153: {
1154:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

1158:   am->misfitobjgradP     = ctx;
1159:   am->ops->misfitobjgrad = func;
1160:   return(0);
1161: }

1163: /*@C
1164:    TaoADMMSetMisfitHessianRoutine - Sets the user-defined misfit Hessian call-back
1165:    function into the algorithm, to be used for subsolverX.

1167:    Collective on tao

1169:    Input Parameters:
1170:    + tao - the Tao context
1171:    . H - user-created matrix for the Hessian of the misfit term
1172:    . Hpre - user-created matrix for the preconditioner of Hessian of the misfit term
1173:    . func - function pointer for the misfit Hessian evaluation
1174:    - ctx - user context for the misfit Hessian

1176:    Level: advanced
1177:    
1178: .seealso: TAOADMM

1180: @*/
1181: PetscErrorCode TaoADMMSetMisfitHessianRoutine(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1182: {
1183:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

1188:   if (H) {
1191:   }
1192:   if (Hpre) {
1195:   }
1196:   if (ctx) {
1197:     am->misfithessP = ctx;
1198:   }
1199:   if (func) {
1200:     am->ops->misfithess = func;
1201:   }
1202:   if (H) {
1203:     PetscObjectReference((PetscObject)H);
1204:     MatDestroy(&am->Hx);
1205:     am->Hx = H;
1206:   }
1207:   if (Hpre) {
1208:     PetscObjectReference((PetscObject)Hpre);
1209:     MatDestroy(&am->Hxpre);
1210:     am->Hxpre = Hpre;
1211:   }
1212:   return(0);
1213: }

1215: /*@C
1216:    TaoADMMSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back function

1218:    Collective on tao

1220:    Input Parameters:
1221:    + tao - the Tao context
1222:    . func - function pointer for the regularizer value and gradient evaluation
1223:    - ctx - user context for the regularizer

1225:    Level: advanced
1226:    
1227: .seealso: TAOADMM

1229: @*/
1230: PetscErrorCode TaoADMMSetRegularizerObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao, Vec, PetscReal*, Vec, void*), void *ctx)
1231: {
1232:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

1236:   am->regobjgradP     = ctx;
1237:   am->ops->regobjgrad = func;
1238:   return(0);
1239: }

1241: /*@C
1242:    TaoADMMSetRegularizerHessianRoutine - Sets the user-defined regularizer Hessian call-back
1243:    function, to be used for subsolverZ.

1245:    Collective on tao

1247:    Input Parameters:
1248:    + tao - the Tao context
1249:    . H - user-created matrix for the Hessian of the regularization term
1250:    . Hpre - user-created matrix for the preconditioner of Hessian of the regularization term
1251:    . func - function pointer for the regularizer Hessian evaluation
1252:    - ctx - user context for the regularizer Hessian

1254:    Level: advanced
1255:    
1256: .seealso: TAOADMM

1258: @*/
1259: PetscErrorCode TaoADMMSetRegularizerHessianRoutine(Tao tao, Mat H, Mat Hpre, PetscErrorCode (*func)(Tao, Vec, Mat, Mat, void*), void *ctx)
1260: {
1261:   TAO_ADMM       *am = (TAO_ADMM*)tao->data;

1266:   if (H) {
1269:   }
1270:   if (Hpre) {
1273:   }
1274:   if (ctx) {
1275:     am->reghessP = ctx;
1276:   }
1277:   if (func) {
1278:     am->ops->reghess = func;
1279:   }
1280:   if (H) {
1281:     PetscObjectReference((PetscObject)H);
1282:     MatDestroy(&am->Hz);
1283:     am->Hz = H;
1284:   }
1285:   if (Hpre) {
1286:     PetscObjectReference((PetscObject)Hpre);
1287:     MatDestroy(&am->Hzpre);
1288:     am->Hzpre = Hpre;
1289:   }
1290:   return(0);
1291: }

1293: /*@
1294:    TaoGetADMMParentTao - Gets pointer to parent ADMM tao, used by inner subsolver.

1296:    Collective on tao

1298:    Input Parameter:
1299:    . tao - the Tao context

1301:    Output Parameter:
1302:    . admm_tao - the parent Tao context

1304:    Level: advanced
1305:    
1306: .seealso: TAOADMM

1308: @*/
1309: PetscErrorCode TaoGetADMMParentTao(Tao tao, Tao *admm_tao)
1310: {

1315:   PetscObjectQuery((PetscObject)tao,"TaoGetADMMParentTao_ADMM", (PetscObject*) admm_tao);
1316:   return(0);
1317: }

1319: /*@
1320:   TaoADMMGetDualVector - Returns the dual vector associated with the current TAOADMM state

1322:   Not Collective

1324:   Input Parameter:
1325:   . tao - the Tao context

1327:   Output Parameter:
1328:   . Y - the current solution

1330:   Level: intermediate
1331:   
1332: .seealso: TAOADMM

1334: @*/
1335: PetscErrorCode TaoADMMGetDualVector(Tao tao, Vec *Y)
1336: {
1337:   TAO_ADMM *am = (TAO_ADMM*)tao->data;

1341:   *Y = am->y;
1342:   return(0);
1343: }

1345: /*@
1346:   TaoADMMSetRegularizerType - Set regularizer type for ADMM routine

1348:   Not Collective

1350:   Input Parameter:
1351: + tao  - the Tao context
1352: - type - regularizer type

1354:   Options Database:
1355: .  -tao_admm_regularizer_type <admm_regularizer_user,admm_regularizer_soft_thresh>

1357:   Level: intermediate

1359: .seealso: TaoADMMGetRegularizerType(), TaoADMMRegularizerType, TAOADMM
1360: @*/
1361: PetscErrorCode TaoADMMSetRegularizerType(Tao tao, TaoADMMRegularizerType type)
1362: {

1368:   PetscTryMethod(tao,"TaoADMMSetRegularizerType_C",(Tao,TaoADMMRegularizerType),(tao,type));
1369:   return(0);
1370: }

1372: /*@
1373:    TaoADMMGetRegularizerType - Gets the type of regularizer routine for ADMM

1375:    Not Collective

1377:    Input Parameter:
1378: .  tao - the Tao context

1380:    Output Parameter:
1381: .  type - the type of regularizer

1383:   Options Database:
1384: .  -tao_admm_regularizer_type <admm_regularizer_user,admm_regularizer_soft_thresh>

1386:    Level: intermediate

1388: .seealso: TaoADMMSetRegularizerType(), TaoADMMRegularizerType, TAOADMM
1389: @*/
1390: PetscErrorCode TaoADMMGetRegularizerType(Tao tao, TaoADMMRegularizerType *type)
1391: {

1396:   PetscUseMethod(tao,"TaoADMMGetRegularizerType_C",(Tao,TaoADMMRegularizerType*),(tao,type));
1397:   return(0);
1398: }

1400: /*@
1401:   TaoADMMSetUpdateType - Set update routine for ADMM routine

1403:   Not Collective

1405:   Input Parameter:
1406: + tao  - the Tao context
1407: - type - spectral parameter update type

1409:   Level: intermediate

1411: .seealso: TaoADMMGetUpdateType(), TaoADMMUpdateType, TAOADMM
1412: @*/
1413: PetscErrorCode TaoADMMSetUpdateType(Tao tao, TaoADMMUpdateType type)
1414: {

1420:   PetscTryMethod(tao,"TaoADMMSetUpdateType_C",(Tao,TaoADMMUpdateType),(tao,type));
1421:   return(0);
1422: }

1424: /*@
1425:    TaoADMMGetUpdateType - Gets the type of spectral penalty update routine for ADMM

1427:    Not Collective

1429:    Input Parameter:
1430: .  tao - the Tao context

1432:    Output Parameter:
1433: .  type - the type of spectral penalty update routine

1435:    Level: intermediate

1437: .seealso: TaoADMMSetUpdateType(), TaoADMMUpdateType, TAOADMM
1438: @*/
1439: PetscErrorCode TaoADMMGetUpdateType(Tao tao, TaoADMMUpdateType *type)
1440: {

1445:   PetscUseMethod(tao,"TaoADMMGetUpdateType_C",(Tao,TaoADMMUpdateType*),(tao,type));
1446:   return(0);
1447: }