1 #include <../src/tao/leastsquares/impls/brgn/brgn.h> 2 3 /* Old code 4 static PetscErrorCode GNHessianProd(Mat H, Vec in, Vec out) 5 { 6 TAO_BRGN *gn; 7 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 ierr = MatShellGetContext(H, &gn);CHKERRQ(ierr); 11 ierr = MatMult(gn->subsolver->ls_jac, in, gn->r_work);CHKERRQ(ierr); 12 ierr = MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out);CHKERRQ(ierr); 13 ierr = VecAXPY(out, gn->lambda, in);CHKERRQ(ierr); 14 PetscFunctionReturn(0); 15 } 16 17 static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr) 18 { 19 TAO_BRGN *gn = (TAO_BRGN *)ptr; 20 PetscScalar xnorm2; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 ierr = TaoComputeResidual(tao, X, tao->ls_res);CHKERRQ(ierr); 25 ierr = VecDotBegin(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 26 ierr = VecAXPBYPCZ(gn->x_work, 1.0, -1.0, 0.0, X, gn->x_old);CHKERRQ(ierr); 27 ierr = VecDotBegin(gn->x_work, gn->x_work, &xnorm2);CHKERRQ(ierr); 28 ierr = VecDotEnd(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 29 ierr = VecDotEnd(gn->x_work, gn->x_work, &xnorm2);CHKERRQ(ierr); 30 *fcn = 0.5*(*fcn) + 0.5*gn->lambda*xnorm2; 31 32 ierr = TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);CHKERRQ(ierr); 33 ierr = MatMultTranspose(tao->ls_jac, tao->ls_res, G);CHKERRQ(ierr); 34 ierr = VecAXPBYPCZ(G, gn->lambda, -gn->lambda, 1.0, X, gn->x_old);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 */ 38 39 static PetscErrorCode GNHessianProd(Mat H, Vec in, Vec out) 40 { 41 TAO_BRGN *gn; 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 ierr = MatShellGetContext(H, &gn);CHKERRQ(ierr); 46 ierr = MatMult(gn->subsolver->ls_jac, in, gn->r_work);CHKERRQ(ierr); 47 ierr = MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out);CHKERRQ(ierr); 48 /* out = out + lambda*D'*(diag.*(D*in)) */ 49 ierr = MatMult(gn->D, in, gn->y);CHKERRQ(ierr); /* y = D*in */ 50 ierr = VecPointwiseMult(gn->y_work, gn->diag, gn->y);CHKERRQ(ierr); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */ 51 ierr = MatMultTranspose(gn->D, gn->y_work, gn->x_work);CHKERRQ(ierr); /* x_work = D'*(diag.*(D*in)) */ 52 ierr = VecAXPY(out, gn->lambda, gn->x_work);CHKERRQ(ierr); 53 54 PetscFunctionReturn(0); 55 } 56 57 static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr) 58 { 59 TAO_BRGN *gn = (TAO_BRGN *)ptr; 60 PetscInt Ny; /* dimension of D*X */ 61 PetscScalar yESum; 62 PetscErrorCode ierr; 63 64 PetscFunctionBegin; 65 /* compute objective */ 66 /* compute first term ||ls_res||^2 */ 67 ierr = TaoComputeResidual(tao, X, tao->ls_res);CHKERRQ(ierr); 68 ierr = VecDotBegin(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 69 ierr = VecDotEnd(tao->ls_res, tao->ls_res, fcn);CHKERRQ(ierr); 70 /* add the second term lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/ 71 ierr = MatMult(gn->D, X, gn->y);CHKERRQ(ierr); /* y = D*x */ 72 ierr = VecPointwiseMult(gn->y_work, gn->y, gn->y);CHKERRQ(ierr); 73 ierr = VecShift(gn->y_work, gn->epsilon*gn->epsilon);CHKERRQ(ierr); 74 ierr = VecSqrtAbs(gn->y_work);CHKERRQ(ierr); /* gn->y_work = sqrt(y.^2+epsilon^2) */ 75 ierr = VecSum(gn->y_work, &yESum);CHKERRQ(ierr);CHKERRQ(ierr); 76 ierr = VecGetSize(gn->y, &Ny);CHKERRQ(ierr); 77 *fcn = 0.5*(*fcn) + gn->lambda*(yESum - Ny*gn->epsilon); 78 79 /* compute gradient G */ 80 ierr = TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);CHKERRQ(ierr); 81 ierr = MatMultTranspose(tao->ls_jac, tao->ls_res, G);CHKERRQ(ierr); 82 /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)), where y = D*x */ 83 ierr = VecPointwiseDivide(gn->y_work, gn->y, gn->y_work);CHKERRQ(ierr); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */ 84 ierr = MatMultTranspose(gn->D, gn->y_work, gn->x_work);CHKERRQ(ierr); 85 ierr = VecAXPY(G, gn->lambda, gn->x_work);CHKERRQ(ierr); 86 87 PetscFunctionReturn(0); 88 } 89 90 91 static PetscErrorCode GNComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 92 { 93 TAO_BRGN *gn = (TAO_BRGN *)ptr; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);CHKERRQ(ierr); 98 99 /* calculate and store diagonal matrix as a vector: diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3* --> diag = epsilon^2 ./ sqrt(y.^2+epsilon^2).^3, where y = D*x */ 100 ierr = MatMult(gn->D, X, gn->y);CHKERRQ(ierr); /* y = D*x */ 101 ierr = VecPointwiseMult(gn->y_work, gn->y, gn->y);CHKERRQ(ierr); 102 ierr = VecShift(gn->y_work, gn->epsilon*gn->epsilon);CHKERRQ(ierr); 103 ierr = VecCopy(gn->y_work, gn->diag);CHKERRQ(ierr); /* gn->diag = y.^2+epsilon^2 */ 104 ierr = VecSqrtAbs(gn->y_work);CHKERRQ(ierr); /* gn->y_work = sqrt(y.^2+epsilon^2) */ 105 ierr = VecPointwiseMult(gn->diag, gn->y_work, gn->diag);CHKERRQ(ierr); /* gn->diag = sqrt(y.^2+epsilon^2).^3 */ 106 ierr = VecReciprocal(gn->diag);CHKERRQ(ierr); 107 ierr = VecScale(gn->diag, gn->epsilon*gn->epsilon);CHKERRQ(ierr); 108 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode GNHookFunction(Tao tao, PetscInt iter) 113 { 114 TAO_BRGN *gn = (TAO_BRGN *)tao->user_update; 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 /* Update basic tao information from the subsolver */ 119 gn->parent->nfuncs = tao->nfuncs; 120 gn->parent->ngrads = tao->ngrads; 121 gn->parent->nfuncgrads = tao->nfuncgrads; 122 gn->parent->nhess = tao->nhess; 123 gn->parent->niter = tao->niter; 124 gn->parent->ksp_its = tao->ksp_its; 125 gn->parent->ksp_tot_its = tao->ksp_tot_its; 126 ierr = TaoGetConvergedReason(tao, &gn->parent->reason);CHKERRQ(ierr); 127 /* Update the solution vectors */ 128 if (iter == 0) { 129 ierr = VecSet(gn->x_old, 0.0);CHKERRQ(ierr); 130 } else { 131 ierr = VecCopy(tao->solution, gn->x_old);CHKERRQ(ierr); 132 ierr = VecCopy(tao->solution, gn->parent->solution);CHKERRQ(ierr); 133 } 134 /* Update the gradient */ 135 ierr = VecCopy(tao->gradient, gn->parent->gradient);CHKERRQ(ierr); 136 /* Call general purpose update function */ 137 if (gn->parent->ops->update) { 138 ierr = (*gn->parent->ops->update)(gn->parent, gn->parent->niter);CHKERRQ(ierr); 139 } 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode TaoSolve_BRGN(Tao tao) 144 { 145 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = TaoSolve(gn->subsolver);CHKERRQ(ierr); 150 /* Update basic tao information from the subsolver */ 151 tao->nfuncs = gn->subsolver->nfuncs; 152 tao->ngrads = gn->subsolver->ngrads; 153 tao->nfuncgrads = gn->subsolver->nfuncgrads; 154 tao->nhess = gn->subsolver->nhess; 155 tao->niter = gn->subsolver->niter; 156 tao->ksp_its = gn->subsolver->ksp_its; 157 tao->ksp_tot_its = gn->subsolver->ksp_tot_its; 158 ierr = TaoGetConvergedReason(gn->subsolver, &tao->reason);CHKERRQ(ierr); 159 /* Update vectors */ 160 ierr = VecCopy(gn->subsolver->solution, tao->solution);CHKERRQ(ierr); 161 ierr = VecCopy(gn->subsolver->gradient, tao->gradient);CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 165 static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao) 166 { 167 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 /* old Tikhonov regularization code 172 ierr = PetscOptionsHead(PetscOptionsObject,"Gauss-Newton method for least-squares problems using Tikhonov regularization");CHKERRQ(ierr); 173 ierr = PetscOptionsReal("-tao_brgn_lambda", "Tikhonov regularization factor", "", gn->lambda, &gn->lambda, NULL);CHKERRQ(ierr); 174 */ 175 ierr = PetscOptionsHead(PetscOptionsObject,"least-squares problems with L1 regularizer: ||f(x)||^2 + lambda*||x||_1. Currently L1-norm is approximated with smooth form");CHKERRQ(ierr); 176 ierr = PetscOptionsReal("-tao_brgn_lambda", "L1-norm regularizer weight", "", gn->lambda, &gn->lambda, NULL);CHKERRQ(ierr); 177 ierr = PetscOptionsReal("-tao_brgn_epsilon", "L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon)", "", gn->epsilon, &gn->epsilon, NULL);CHKERRQ(ierr); 178 ierr = PetscOptionsTail();CHKERRQ(ierr); 179 ierr = TaoSetFromOptions(gn->subsolver);CHKERRQ(ierr); 180 PetscFunctionReturn(0); 181 } 182 183 static PetscErrorCode TaoView_BRGN(Tao tao, PetscViewer viewer) 184 { 185 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 190 ierr = TaoView(gn->subsolver, viewer);CHKERRQ(ierr); 191 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 static PetscErrorCode TaoSetUp_BRGN(Tao tao) 196 { 197 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 198 PetscErrorCode ierr; 199 PetscBool is_bnls, is_bntr, is_bntl; 200 PetscInt i, nx, Nx, Ny; /* XH: added Ny for size of y=D*x*/ 201 PetscScalar v = 1.0; /* XH: hack to set value of matrix */ 202 203 PetscFunctionBegin; 204 if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualRoutine() must be called before setup!"); 205 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls);CHKERRQ(ierr); 206 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr);CHKERRQ(ierr); 207 ierr = PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl);CHKERRQ(ierr); 208 if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualJacobianRoutine() must be called before setup!"); 209 if (!tao->gradient){ 210 ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); 211 } 212 if (!gn->x_work){ 213 ierr = VecDuplicate(tao->solution, &gn->x_work);CHKERRQ(ierr); 214 } 215 if (!gn->r_work){ 216 ierr = VecDuplicate(tao->ls_res, &gn->r_work);CHKERRQ(ierr); 217 } 218 if (!gn->x_old) { 219 ierr = VecDuplicate(tao->solution, &gn->x_old);CHKERRQ(ierr); 220 ierr = VecSet(gn->x_old, 0.0);CHKERRQ(ierr); 221 } 222 223 224 /* XH: hack: following works only when D is a squared matrix, to fix dimension of gn->diag/y/y_work, when D is not squared matrix */ 225 ierr = VecGetSize(tao->solution, &Nx);CHKERRQ(ierr); 226 /* Ny = Nx; */ /* for identity matrix */ 227 Ny = Nx - 1; /* for gradient matrix */ 228 if (!gn->y){ 229 ierr = VecCreate(PETSC_COMM_SELF,&gn->y);CHKERRQ(ierr); 230 ierr = VecSetSizes(gn->y,PETSC_DECIDE,Ny);CHKERRQ(ierr); 231 ierr = VecSetFromOptions(gn->y);CHKERRQ(ierr); 232 ierr = VecSet(gn->y,0.0);CHKERRQ(ierr); 233 234 } 235 if (!gn->y_work){ 236 ierr = VecDuplicate(gn->y, &gn->y_work);CHKERRQ(ierr); 237 } 238 if (!gn->diag){ 239 ierr = VecDuplicate(gn->y, &gn->diag);CHKERRQ(ierr); 240 ierr = VecSet(gn->diag, 0.0);CHKERRQ(ierr); 241 } 242 243 /* XH: hack: set gn->D as identity/gradient matrix here for text , how to set D matrix from user data?*/ 244 if (!gn->D){ 245 ierr = MatCreate(PETSC_COMM_SELF,&gn->D);CHKERRQ(ierr); 246 ierr = MatSetSizes(gn->D,PETSC_DECIDE,PETSC_DECIDE,Ny,Nx);CHKERRQ(ierr); 247 ierr = MatSetFromOptions(gn->D);CHKERRQ(ierr); 248 ierr = MatSetUp(gn->D);CHKERRQ(ierr); 249 /* identity matrix */ 250 /* 251 for (i=0; i<Ny; i++) { 252 ierr = MatSetValues(gn->D,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 253 } 254 */ 255 /* gradient matrix */ 256 /* [-1, 1, 0,...; 0, -1, 1, 0, ...] */ 257 for (i=0; i<Ny; i++) { 258 v = 1.0; 259 nx = i+1; 260 ierr = MatSetValues(gn->D,1,&i,1,&nx,&v,INSERT_VALUES);CHKERRQ(ierr); 261 v = -1.0; 262 ierr = MatSetValues(gn->D,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 263 } 264 ierr = MatAssemblyBegin(gn->D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 265 ierr = MatAssemblyEnd(gn->D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 266 } 267 268 /* XH: debug: check matrix */ 269 PetscPrintf(PETSC_COMM_SELF, "-------- Check D matrix. -------- \n"); 270 ierr = MatView(gn->D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 271 getchar(); 272 273 if (!tao->setupcalled) { 274 /* Hessian setup */ 275 ierr = VecGetLocalSize(tao->solution, &nx);CHKERRQ(ierr); 276 ierr = VecGetSize(tao->solution, &Nx);CHKERRQ(ierr); 277 ierr = MatSetSizes(gn->H, nx, nx, Nx, Nx);CHKERRQ(ierr); 278 ierr = MatSetType(gn->H, MATSHELL);CHKERRQ(ierr); 279 ierr = MatSetUp(gn->H);CHKERRQ(ierr); 280 ierr = MatShellSetOperation(gn->H, MATOP_MULT, (void (*)(void))GNHessianProd);CHKERRQ(ierr); 281 ierr = MatShellSetContext(gn->H, (void*)gn);CHKERRQ(ierr); 282 /* Subsolver setup */ 283 ierr = TaoSetUpdate(gn->subsolver, GNHookFunction, (void*)gn);CHKERRQ(ierr); 284 ierr = TaoSetInitialVector(gn->subsolver, tao->solution);CHKERRQ(ierr); 285 if (tao->bounded) { 286 ierr = TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU);CHKERRQ(ierr); 287 } 288 ierr = TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP);CHKERRQ(ierr); 289 ierr = TaoSetJacobianResidualRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP);CHKERRQ(ierr); 290 ierr = TaoSetObjectiveAndGradientRoutine(gn->subsolver, GNObjectiveGradientEval, (void*)gn);CHKERRQ(ierr); 291 ierr = TaoSetHessianRoutine(gn->subsolver, gn->H, gn->H, GNComputeHessian, (void*)gn);CHKERRQ(ierr); 292 /* Propagate some options down */ 293 ierr = TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol);CHKERRQ(ierr); 294 ierr = TaoSetMaximumIterations(gn->subsolver, tao->max_it);CHKERRQ(ierr); 295 ierr = TaoSetMaximumFunctionEvaluations(gn->subsolver, tao->max_funcs);CHKERRQ(ierr); 296 for (i=0; i<tao->numbermonitors; ++i) { 297 ierr = TaoSetMonitor(gn->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]);CHKERRQ(ierr); 298 ierr = PetscObjectReference((PetscObject)(tao->monitorcontext[i]));CHKERRQ(ierr); 299 } 300 ierr = TaoSetUp(gn->subsolver);CHKERRQ(ierr); 301 } 302 PetscFunctionReturn(0); 303 } 304 305 static PetscErrorCode TaoDestroy_BRGN(Tao tao) 306 { 307 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 if (tao->setupcalled) { 312 ierr = VecDestroy(&tao->gradient);CHKERRQ(ierr); 313 ierr = VecDestroy(&gn->x_work);CHKERRQ(ierr); 314 ierr = VecDestroy(&gn->r_work);CHKERRQ(ierr); 315 ierr = VecDestroy(&gn->x_old);CHKERRQ(ierr); 316 ierr = VecDestroy(&gn->diag);CHKERRQ(ierr); 317 ierr = VecDestroy(&gn->y);CHKERRQ(ierr); 318 ierr = VecDestroy(&gn->y_work);CHKERRQ(ierr); 319 } 320 ierr = MatDestroy(&gn->H);CHKERRQ(ierr); 321 ierr = MatDestroy(&gn->D);CHKERRQ(ierr); 322 ierr = TaoDestroy(&gn->subsolver);CHKERRQ(ierr); 323 gn->parent = NULL; 324 ierr = PetscFree(tao->data);CHKERRQ(ierr); 325 PetscFunctionReturn(0); 326 } 327 328 /*MC 329 TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares 330 problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL 331 that constructs the Guass-Newton problem with the user-provided least-squares 332 residual and Jacobian. The problem is regularized with an L2-norm proximal point 333 term. 334 335 Options Database Keys: 336 + -tao_bqnk_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop 337 . -tao_bqnk_init_type - trust radius initialization method ("constant", "direction", "interpolation") 338 . -tao_bqnk_update_type - trust radius update method ("step", "direction", "interpolation") 339 - -tao_bqnk_as_type - active-set estimation method ("none", "bertsekas") 340 341 Level: beginner 342 M*/ 343 PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao) 344 { 345 TAO_BRGN *gn; 346 PetscErrorCode ierr; 347 348 PetscFunctionBegin; 349 ierr = PetscNewLog(tao,&gn);CHKERRQ(ierr); 350 351 tao->ops->destroy = TaoDestroy_BRGN; 352 tao->ops->setup = TaoSetUp_BRGN; 353 tao->ops->setfromoptions = TaoSetFromOptions_BRGN; 354 tao->ops->view = TaoView_BRGN; 355 tao->ops->solve = TaoSolve_BRGN; 356 357 tao->data = (void*)gn; 358 gn->lambda = 1e-4; 359 gn->epsilon = 1e-6; 360 gn->parent = tao; 361 362 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &gn->H);CHKERRQ(ierr); 363 ierr = MatSetOptionsPrefix(gn->H, "tao_brgn_hessian_");CHKERRQ(ierr); 364 365 ierr = TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver);CHKERRQ(ierr); 366 ierr = TaoSetType(gn->subsolver, TAOBNLS);CHKERRQ(ierr); 367 ierr = TaoSetOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_");CHKERRQ(ierr); 368 PetscFunctionReturn(0); 369 } 370 371 /*@C 372 TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN 373 374 Collective on Tao 375 376 Level: developer 377 378 Input Parameters: 379 + tao - the Tao solver context 380 - subsolver - the Tao sub-solver context 381 @*/ 382 PetscErrorCode TaoBRGNGetSubsolver(Tao tao, Tao *subsolver) 383 { 384 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 385 386 PetscFunctionBegin; 387 *subsolver = gn->subsolver; 388 PetscFunctionReturn(0); 389 } 390 391 /*@C 392 TaoBRGNSetTikhonovLambda - Set the Tikhonov regularization factor for the Gauss-Newton least-squares algorithm 393 394 Collective on Tao 395 396 Level: developer 397 398 Input Parameters: 399 + tao - the Tao solver context 400 - lambda - Tikhonov regularization factor 401 @*/ 402 PetscErrorCode TaoBRGNSetTikhonovLambda(Tao tao, PetscReal lambda) 403 { 404 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 405 406 /* Initialize lambda here */ 407 408 PetscFunctionBegin; 409 gn->lambda = lambda; 410 PetscFunctionReturn(0); 411 } 412 413 /*@C 414 TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm 415 416 Collective on Tao 417 418 Level: developer 419 420 Input Parameters: 421 + tao - the Tao solver context 422 - epsilon - L1-norm smooth approximation parameter 423 @*/ 424 PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao, PetscReal epsilon) 425 { 426 TAO_BRGN *gn = (TAO_BRGN *)tao->data; 427 428 /* Initialize epsilon here */ 429 430 PetscFunctionBegin; 431 gn->epsilon = epsilon; 432 PetscFunctionReturn(0); 433 } 434 /* XH: done TaoBRGNSetL1SmoothEpsilon by copy TaoBRGNSetTikhonovLambda in peststao.h, brgn.c and zbrgnf.c. 435 maybe change the name of Tikhonov in TaoBRGNSetTikhonovLambda() etc, as lambda is no longer the Tikhonov regularizer weight but the L1 regularizer weight 436 Maybe change D*x to D(x), and A*x to A(x) as function handle 437 Maybe need to also keep y = D*x, to avoid duplicate frequent computation of D*x 438 */ 439