1 #include <../src/tao/leastsquares/impls/brgn/brgn.h> 2 3 static PetscErrorCode GNHessianProd(Mat H, Vec in, Vec out) 4 { 5 TAO_BRGN *gn; 6 PetscErrorCode ierr; 7 8 PetscFunctionBegin; 9 ierr = MatShellGetContext(H, &gn);CHKERRQ(ierr); 10 ierr = MatMult(gn->subsolver->ls_jac, in, gn->r_work);CHKERRQ(ierr); 11 ierr = MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out);CHKERRQ(ierr); 12 ierr = VecAXPY(out, gn->lambda, in);CHKERRQ(ierr); 13 PetscFunctionReturn(0); 14 } 15 16 static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr) 17 { 18 TAO_BRGN *gn = (TAO_BRGN *)ptr; 19 PetscScalar xnorm2; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = TaoComputeResidual(tao, X, tao->ls_res);CHKERRQ(ierr); 24 ierr = VecDotBegin(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 25 ierr = VecAXPBY(gn->x_work, gn->lambda, 0.0, X);CHKERRQ(ierr); 26 ierr = VecDotBegin(gn->x_work, gn->x_work, &xnorm2);CHKERRQ(ierr); 27 ierr = VecDotEnd(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 28 ierr = VecDotEnd(gn->x_work, gn->x_work, &xnorm2);CHKERRQ(ierr); 29 *fcn += xnorm2; 30 31 ierr = TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);CHKERRQ(ierr); 32 ierr = MatMultTranspose(tao->ls_jac, tao->ls_res, G);CHKERRQ(ierr); 33 ierr = VecAXPY(G, gn->lambda, X);CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 static PetscErrorCode GNComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 38 { 39 TAO_BRGN *gn = (TAO_BRGN *)ptr; 40 PetscErrorCode ierr; 41 42 PetscFunctionBegin; 43 if (gn->explicit_H) { 44 ierr = MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_REUSE_MATRIX, PETSC_DEFAULT, &H);CHKERRQ(ierr); 45 ierr = MatShift(H, gn->lambda);CHKERRQ(ierr); 46 } 47 PetscFunctionReturn(0); 48 } 49 50 static PetscErrorCode TaoSolve_BRGN(Tao tao) 51 { 52 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 53 PetscErrorCode ierr; 54 55 PetscFunctionBegin; 56 ierr = TaoSolve(gn->subsolver);CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao) 61 { 62 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 ierr = PetscOptionsHead(PetscOptionsObject,"Gauss-Newton method for least-squares problems using Tikhonov regularization");CHKERRQ(ierr); 67 ierr = PetscOptionsReal("-tao_brgn_lambda", "Tikhonov regularization factor", "", gn->lambda, &gn->lambda, NULL);CHKERRQ(ierr); 68 ierr = PetscOptionsBool("-tao_brgn_explicit_hessian","(developer) explicitly perform MatTransposeMatMult for the least squares Jacobian","",gn->explicit_H,&gn->explicit_H,NULL);CHKERRQ(ierr); 69 ierr = PetscOptionsTail();CHKERRQ(ierr); 70 ierr = TaoSetFromOptions(gn->subsolver);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 static PetscErrorCode TaoView_BRGN(Tao tao, PetscViewer viewer) 75 { 76 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 ierr = TaoView(gn->subsolver, viewer);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 static PetscErrorCode TaoSetUp_BRGN(Tao tao) 85 { 86 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 87 PetscErrorCode ierr; 88 PetscBool is_bnls, is_bntr, is_bntl; 89 PetscInt i, nx, Nx; 90 91 PetscFunctionBegin; 92 if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualRoutine() must be called before setup!"); 93 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls);CHKERRQ(ierr); 94 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr);CHKERRQ(ierr); 95 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl);CHKERRQ(ierr); 96 if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualJacobianRoutine() must be called before setup!"); 97 if (!gn->x_work){ 98 ierr = VecDuplicate(tao->solution, &gn->x_work);CHKERRQ(ierr); 99 } 100 if (!gn->r_work){ 101 ierr = VecDuplicate(tao->ls_res, &gn->r_work);CHKERRQ(ierr); 102 } 103 /* Hessian setup */ 104 if (gn->explicit_H) { 105 ierr = MatDestroy(&gn->H);CHKERRQ(ierr); 106 ierr = MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &gn->H);CHKERRQ(ierr); 107 } else { 108 ierr = VecGetLocalSize(tao->solution, &nx);CHKERRQ(ierr); 109 ierr = VecGetSize(tao->solution, &Nx);CHKERRQ(ierr); 110 ierr = MatSetSizes(gn->H, nx, nx, Nx, Nx);CHKERRQ(ierr); 111 ierr = MatSetType(gn->H, MATSHELL);CHKERRQ(ierr); 112 ierr = MatSetUp(gn->H);CHKERRQ(ierr); 113 ierr = MatShellSetOperation(gn->H, MATOP_MULT, (void (*)(void))GNHessianProd);CHKERRQ(ierr); 114 ierr = MatShellSetContext(gn->H, (void*)gn);CHKERRQ(ierr); 115 } 116 /* Subsolver setup */ 117 ierr = TaoSetInitialVector(gn->subsolver, tao->solution);CHKERRQ(ierr); 118 ierr = TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol);CHKERRQ(ierr); 119 if (tao->bounded) { 120 ierr = TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU);CHKERRQ(ierr); 121 } 122 ierr = TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP);CHKERRQ(ierr); 123 ierr = TaoSetResidualJacobianRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP);CHKERRQ(ierr); 124 ierr = TaoSetObjectiveAndGradientRoutine(gn->subsolver, GNObjectiveGradientEval, (void*)gn);CHKERRQ(ierr); 125 ierr = TaoSetHessianRoutine(gn->subsolver, gn->H, gn->H, GNComputeHessian, (void*)gn);CHKERRQ(ierr); 126 for (i=0; i<tao->numbermonitors; ++i) { 127 ierr = TaoSetMonitor(gn->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]);CHKERRQ(ierr); 128 ierr = PetscObjectReference((PetscObject)(tao->monitorcontext[i]));CHKERRQ(ierr); 129 } 130 ierr = TaoSetUp(gn->subsolver);CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 134 static PetscErrorCode TaoDestroy_BRGN(Tao tao) 135 { 136 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 137 PetscErrorCode ierr; 138 139 PetscFunctionBegin; 140 if (tao->setupcalled) { 141 ierr = VecDestroy(&gn->x_work);CHKERRQ(ierr); 142 ierr = VecDestroy(&gn->r_work);CHKERRQ(ierr); 143 } 144 ierr = MatDestroy(&gn->H);CHKERRQ(ierr); 145 ierr = TaoDestroy(&gn->subsolver);CHKERRQ(ierr); 146 ierr = PetscFree(tao->data);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao) 151 { 152 TAO_BRGN *gn; 153 PetscErrorCode ierr; 154 155 PetscFunctionBegin; 156 ierr = PetscNewLog(tao,&gn);CHKERRQ(ierr); 157 158 tao->ops->destroy = TaoDestroy_BRGN; 159 tao->ops->setup = TaoSetUp_BRGN; 160 tao->ops->setfromoptions = TaoSetFromOptions_BRGN; 161 tao->ops->view = TaoView_BRGN; 162 tao->ops->solve = TaoSolve_BRGN; 163 164 tao->data = (void*)gn; 165 gn->lambda = 1.0; 166 gn->explicit_H = PETSC_FALSE; 167 gn->assembled_H = PETSC_FALSE; 168 169 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &gn->H);CHKERRQ(ierr); 170 ierr = MatSetOptionsPrefix(gn->H, "tao_brgn_hessian_");CHKERRQ(ierr); 171 ierr = PetscObjectIncrementTabLevel((PetscObject)gn->H, (PetscObject)tao, 1);CHKERRQ(ierr); 172 173 ierr = TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver);CHKERRQ(ierr); 174 ierr = TaoSetType(gn->subsolver, TAOBNLS);CHKERRQ(ierr); 175 ierr = TaoSetOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_");CHKERRQ(ierr); 176 ierr = PetscObjectIncrementTabLevel((PetscObject)gn->subsolver, (PetscObject)tao, 1);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 /*@C 181 TaoBRGNSetTikhonovLambda - Set the Tikhonov regularization factor for the Gauss-Newton least-squares algorithm 182 183 Collective on Tao 184 185 Level: developer 186 187 Input Parameters: 188 + tao - the Tao solver context 189 - lambda - Tikhonov regularization factor 190 @*/ 191 PetscErrorCode TaoBRGNSetTikhonovLambda(Tao tao, PetscReal lambda) 192 { 193 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 194 195 PetscFunctionBegin; 196 gn->lambda = lambda; 197 PetscFunctionReturn(0); 198 } 199