Actual source code: brgn.c

petsc-master 2019-12-13
Report Typos and Errors
  1:  #include <../src/tao/leastsquares/impls/brgn/brgn.h>

  3: #define BRGN_REGULARIZATION_USER    0
  4: #define BRGN_REGULARIZATION_L2PROX  1
  5: #define BRGN_REGULARIZATION_L2PURE  2
  6: #define BRGN_REGULARIZATION_L1DICT  3
  7: #define BRGN_REGULARIZATION_TYPES   4

  9: static const char *BRGN_REGULARIZATION_TABLE[64] = {"user","l2prox","l2pure","l1dict"};

 11: static PetscErrorCode GNHessianProd(Mat H,Vec in,Vec out)
 12: {
 13:   TAO_BRGN              *gn;
 14:   PetscErrorCode        ierr;
 15: 
 17:   MatShellGetContext(H,&gn);
 18:   MatMult(gn->subsolver->ls_jac,in,gn->r_work);
 19:   MatMultTranspose(gn->subsolver->ls_jac,gn->r_work,out);
 20:   switch (gn->reg_type) {
 21:   case BRGN_REGULARIZATION_USER:
 22:     MatMult(gn->Hreg,in,gn->x_work);
 23:     VecAXPY(out,gn->lambda,gn->x_work);
 24:     break;
 25:   case BRGN_REGULARIZATION_L2PURE:
 26:     VecAXPY(out,gn->lambda,in);
 27:     break;
 28:   case BRGN_REGULARIZATION_L2PROX:
 29:     VecAXPY(out,gn->lambda,in);
 30:     break;
 31:   case BRGN_REGULARIZATION_L1DICT:
 32:     /* out = out + lambda*D'*(diag.*(D*in)) */
 33:     if (gn->D) {
 34:       MatMult(gn->D,in,gn->y);/* y = D*in */
 35:     } else {
 36:       VecCopy(in,gn->y);
 37:     }
 38:     VecPointwiseMult(gn->y_work,gn->diag,gn->y);   /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
 39:     if (gn->D) {
 40:       MatMultTranspose(gn->D,gn->y_work,gn->x_work); /* x_work = D'*(diag.*(D*in)) */
 41:     } else {
 42:       VecCopy(gn->y_work,gn->x_work);
 43:     }
 44:     VecAXPY(out,gn->lambda,gn->x_work);
 45:     break;
 46:   }
 47:   return(0);
 48: }

 50: static PetscErrorCode GNObjectiveGradientEval(Tao tao,Vec X,PetscReal *fcn,Vec G,void *ptr)
 51: {
 52:   TAO_BRGN              *gn = (TAO_BRGN *)ptr;
 53:   PetscInt              K;                    /* dimension of D*X */
 54:   PetscScalar           yESum;
 55:   PetscErrorCode        ierr;
 56:   PetscReal             f_reg;
 57: 
 59:   /* compute objective *fcn*/
 60:   /* compute first term 0.5*||ls_res||_2^2 */
 61:   TaoComputeResidual(tao,X,tao->ls_res);
 62:   VecDot(tao->ls_res,tao->ls_res,fcn);
 63:   *fcn *= 0.5;
 64:   /* compute gradient G */
 65:   TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);
 66:   MatMultTranspose(tao->ls_jac,tao->ls_res,G);
 67:   /* add the regularization contribution */
 68:   switch (gn->reg_type) {
 69:   case BRGN_REGULARIZATION_USER:
 70:     (*gn->regularizerobjandgrad)(tao,X,&f_reg,gn->x_work,gn->reg_obj_ctx);
 71:     *fcn += gn->lambda*f_reg;
 72:     VecAXPY(G,gn->lambda,gn->x_work);
 73:     break;
 74:   case BRGN_REGULARIZATION_L2PURE:
 75:     /* compute f = f + lambda*0.5*xk'*xk */
 76:     VecDot(X,X,&f_reg);
 77:     *fcn += gn->lambda*0.5*f_reg;
 78:     /* compute G = G + lambda*xk */
 79:     VecAXPY(G,gn->lambda,X);
 80:     break;
 81:   case BRGN_REGULARIZATION_L2PROX:
 82:     /* compute f = f + lambda*0.5*(xk - xkm1)'*(xk - xkm1) */
 83:     VecAXPBYPCZ(gn->x_work,1.0,-1.0,0.0,X,gn->x_old);
 84:     VecDot(gn->x_work,gn->x_work,&f_reg);
 85:     *fcn += gn->lambda*0.5*f_reg;
 86:     /* compute G = G + lambda*(xk - xkm1) */
 87:     VecAXPBYPCZ(G,gn->lambda,-gn->lambda,1.0,X,gn->x_old);
 88:     break;
 89:   case BRGN_REGULARIZATION_L1DICT:
 90:     /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/
 91:     if (gn->D) {
 92:       MatMult(gn->D,X,gn->y);/* y = D*x */
 93:     } else {
 94:       VecCopy(X,gn->y);
 95:     }
 96:     VecPointwiseMult(gn->y_work,gn->y,gn->y);
 97:     VecShift(gn->y_work,gn->epsilon*gn->epsilon);
 98:     VecSqrtAbs(gn->y_work);  /* gn->y_work = sqrt(y.^2+epsilon^2) */
 99:     VecSum(gn->y_work,&yESum);
100:     VecGetSize(gn->y,&K);
101:     *fcn += gn->lambda*(yESum - K*gn->epsilon);
102:     /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */
103:     VecPointwiseDivide(gn->y_work,gn->y,gn->y_work); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
104:     if (gn->D) {
105:       MatMultTranspose(gn->D,gn->y_work,gn->x_work);
106:     } else {
107:       VecCopy(gn->y_work,gn->x_work);
108:     }
109:     VecAXPY(G,gn->lambda,gn->x_work);
110:     break;
111:   }
112:   return(0);
113: }

115: static PetscErrorCode GNComputeHessian(Tao tao,Vec X,Mat H,Mat Hpre,void *ptr)
116: {
117:   TAO_BRGN              *gn = (TAO_BRGN *)ptr;
119: 
121:   TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);

123:   switch (gn->reg_type) {
124:   case BRGN_REGULARIZATION_USER:
125:     (*gn->regularizerhessian)(tao,X,gn->Hreg,gn->reg_hess_ctx);
126:     break;
127:   case BRGN_REGULARIZATION_L2PURE:
128:     break;
129:   case BRGN_REGULARIZATION_L2PROX:
130:     break;
131:   case BRGN_REGULARIZATION_L1DICT:
132:     /* 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 */
133:     if (gn->D) {
134:       MatMult(gn->D,X,gn->y);/* y = D*x */
135:     } else {
136:       VecCopy(X,gn->y);
137:     }
138:     VecPointwiseMult(gn->y_work,gn->y,gn->y);
139:     VecShift(gn->y_work,gn->epsilon*gn->epsilon);
140:     VecCopy(gn->y_work,gn->diag);                  /* gn->diag = y.^2+epsilon^2 */
141:     VecSqrtAbs(gn->y_work);                        /* gn->y_work = sqrt(y.^2+epsilon^2) */
142:     VecPointwiseMult(gn->diag,gn->y_work,gn->diag);/* gn->diag = sqrt(y.^2+epsilon^2).^3 */
143:     VecReciprocal(gn->diag);
144:     VecScale(gn->diag,gn->epsilon*gn->epsilon);
145:     break;
146:   }
147:   return(0);
148: }

150: static PetscErrorCode GNHookFunction(Tao tao,PetscInt iter, void *ctx)
151: {
152:   TAO_BRGN              *gn = (TAO_BRGN *)ctx;
153:   PetscErrorCode        ierr;
154: 
156:   /* Update basic tao information from the subsolver */
157:   gn->parent->nfuncs = tao->nfuncs;
158:   gn->parent->ngrads = tao->ngrads;
159:   gn->parent->nfuncgrads = tao->nfuncgrads;
160:   gn->parent->nhess = tao->nhess;
161:   gn->parent->niter = tao->niter;
162:   gn->parent->ksp_its = tao->ksp_its;
163:   gn->parent->ksp_tot_its = tao->ksp_tot_its;
164:   TaoGetConvergedReason(tao,&gn->parent->reason);
165:   /* Update the solution vectors */
166:   if (iter == 0) {
167:     VecSet(gn->x_old,0.0);
168:   } else {
169:     VecCopy(tao->solution,gn->x_old);
170:     VecCopy(tao->solution,gn->parent->solution);
171:   }
172:   /* Update the gradient */
173:   VecCopy(tao->gradient,gn->parent->gradient);
174:   /* Call general purpose update function */
175:   if (gn->parent->ops->update) {
176:     (*gn->parent->ops->update)(gn->parent,gn->parent->niter,gn->parent->user_update);
177:   }
178:   return(0);
179: }

181: static PetscErrorCode TaoSolve_BRGN(Tao tao)
182: {
183:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
184:   PetscErrorCode        ierr;

187:   TaoSolve(gn->subsolver);
188:   /* Update basic tao information from the subsolver */
189:   tao->nfuncs = gn->subsolver->nfuncs;
190:   tao->ngrads = gn->subsolver->ngrads;
191:   tao->nfuncgrads = gn->subsolver->nfuncgrads;
192:   tao->nhess = gn->subsolver->nhess;
193:   tao->niter = gn->subsolver->niter;
194:   tao->ksp_its = gn->subsolver->ksp_its;
195:   tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
196:   TaoGetConvergedReason(gn->subsolver,&tao->reason);
197:   /* Update vectors */
198:   VecCopy(gn->subsolver->solution,tao->solution);
199:   VecCopy(gn->subsolver->gradient,tao->gradient);
200:   return(0);
201: }

203: static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao)
204: {
205:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
206:   PetscErrorCode        ierr;

209:   PetscOptionsHead(PetscOptionsObject,"least-squares problems with regularizer: ||f(x)||^2 + lambda*g(x), g(x) = ||xk-xkm1||^2 or ||Dx||_1 or user defined function.");
210:   PetscOptionsReal("-tao_brgn_regularizer_weight","regularizer weight (default 1e-4)","",gn->lambda,&gn->lambda,NULL);
211:   PetscOptionsReal("-tao_brgn_l1_smooth_epsilon","L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)","",gn->epsilon,&gn->epsilon,NULL);
212:   PetscOptionsEList("-tao_brgn_regularization_type","regularization type", "",BRGN_REGULARIZATION_TABLE,BRGN_REGULARIZATION_TYPES,BRGN_REGULARIZATION_TABLE[gn->reg_type],&gn->reg_type,NULL);
213:   PetscOptionsTail();
214:   TaoSetFromOptions(gn->subsolver);
215:   return(0);
216: }

218: static PetscErrorCode TaoView_BRGN(Tao tao,PetscViewer viewer)
219: {
220:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
221:   PetscErrorCode        ierr;

224:   PetscViewerASCIIPushTab(viewer);
225:   TaoView(gn->subsolver,viewer);
226:   PetscViewerASCIIPopTab(viewer);
227:   return(0);
228: }

230: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
231: {
232:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
233:   PetscErrorCode        ierr;
234:   PetscBool             is_bnls,is_bntr,is_bntl;
235:   PetscInt              i,n,N,K; /* dict has size K*N*/

238:   if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualRoutine() must be called before setup!");
239:   PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNLS,&is_bnls);
240:   PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTR,&is_bntr);
241:   PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTL,&is_bntl);
242:   if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualJacobianRoutine() must be called before setup!");
243:   if (!tao->gradient) {
244:     VecDuplicate(tao->solution,&tao->gradient);
245:   }
246:   if (!gn->x_work) {
247:     VecDuplicate(tao->solution,&gn->x_work);
248:   }
249:   if (!gn->r_work) {
250:     VecDuplicate(tao->ls_res,&gn->r_work);
251:   }
252:   if (!gn->x_old) {
253:     VecDuplicate(tao->solution,&gn->x_old);
254:     VecSet(gn->x_old,0.0);
255:   }
256: 
257:   if (BRGN_REGULARIZATION_L1DICT == gn->reg_type) {
258:     if (gn->D) {
259:       MatGetSize(gn->D,&K,&N); /* Shell matrices still must have sizes defined. K = N for identity matrix, K=N-1 or N for gradient matrix */
260:     } else {
261:       VecGetSize(tao->solution,&K); /* If user does not setup dict matrix, use identiy matrix, K=N */
262:     }
263:     if (!gn->y) {
264:       VecCreate(PETSC_COMM_SELF,&gn->y);
265:       VecSetSizes(gn->y,PETSC_DECIDE,K);
266:       VecSetFromOptions(gn->y);
267:       VecSet(gn->y,0.0);

269:     }
270:     if (!gn->y_work) {
271:       VecDuplicate(gn->y,&gn->y_work);
272:     }
273:     if (!gn->diag) {
274:       VecDuplicate(gn->y,&gn->diag);
275:       VecSet(gn->diag,0.0);
276:     }
277:   }

279:   if (!tao->setupcalled) {
280:     /* Hessian setup */
281:     VecGetLocalSize(tao->solution,&n);
282:     VecGetSize(tao->solution,&N);
283:     MatSetSizes(gn->H,n,n,N,N);
284:     MatSetType(gn->H,MATSHELL);
285:     MatSetUp(gn->H);
286:     MatShellSetOperation(gn->H,MATOP_MULT,(void (*)(void))GNHessianProd);
287:     MatShellSetContext(gn->H,(void*)gn);
288:     /* Subsolver setup,include initial vector and dicttionary D */
289:     TaoSetUpdate(gn->subsolver,GNHookFunction,(void*)gn);
290:     TaoSetInitialVector(gn->subsolver,tao->solution);
291:     if (tao->bounded) {
292:       TaoSetVariableBounds(gn->subsolver,tao->XL,tao->XU);
293:     }
294:     TaoSetResidualRoutine(gn->subsolver,tao->ls_res,tao->ops->computeresidual,tao->user_lsresP);
295:     TaoSetJacobianResidualRoutine(gn->subsolver,tao->ls_jac,tao->ls_jac,tao->ops->computeresidualjacobian,tao->user_lsjacP);
296:     TaoSetObjectiveAndGradientRoutine(gn->subsolver,GNObjectiveGradientEval,(void*)gn);
297:     TaoSetHessianRoutine(gn->subsolver,gn->H,gn->H,GNComputeHessian,(void*)gn);
298:     /* Propagate some options down */
299:     TaoSetTolerances(gn->subsolver,tao->gatol,tao->grtol,tao->gttol);
300:     TaoSetMaximumIterations(gn->subsolver,tao->max_it);
301:     TaoSetMaximumFunctionEvaluations(gn->subsolver,tao->max_funcs);
302:     for (i=0; i<tao->numbermonitors; ++i) {
303:       TaoSetMonitor(gn->subsolver,tao->monitor[i],tao->monitorcontext[i],tao->monitordestroy[i]);
304:       PetscObjectReference((PetscObject)(tao->monitorcontext[i]));
305:     }
306:     TaoSetUp(gn->subsolver);
307:   }
308:   return(0);
309: }

311: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
312: {
313:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
314:   PetscErrorCode        ierr;

317:   if (tao->setupcalled) {
318:     VecDestroy(&tao->gradient);
319:     VecDestroy(&gn->x_work);
320:     VecDestroy(&gn->r_work);
321:     VecDestroy(&gn->x_old);
322:     VecDestroy(&gn->diag);
323:     VecDestroy(&gn->y);
324:     VecDestroy(&gn->y_work);
325:   }
326:   MatDestroy(&gn->H);
327:   MatDestroy(&gn->D);
328:   MatDestroy(&gn->Hreg);
329:   TaoDestroy(&gn->subsolver);
330:   gn->parent = NULL;
331:   PetscFree(tao->data);
332:   return(0);
333: }

335: /*MC
336:   TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares 
337:             problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL 
338:             that constructs the Gauss-Newton problem with the user-provided least-squares 
339:             residual and Jacobian. The algorithm offers an L2-norm ("l2pure"), L2-norm proximal point ("l2prox") 
340:             regularizer, and L1-norm dictionary regularizer ("l1dict"), where we approximate the 
341:             L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
342:             The user can also provide own regularization function.

344:   Options Database Keys:
345: + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l2pure", "l1dict") (default "l2prox")
346: . -tao_brgn_regularizer_weight  - regularizer weight (default 1e-4)
347: - -tao_brgn_l1_smooth_epsilon   - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)

349:   Level: beginner
350: M*/
351: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
352: {
353:   TAO_BRGN       *gn;
355: 
357:   PetscNewLog(tao,&gn);
358: 
359:   tao->ops->destroy = TaoDestroy_BRGN;
360:   tao->ops->setup = TaoSetUp_BRGN;
361:   tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
362:   tao->ops->view = TaoView_BRGN;
363:   tao->ops->solve = TaoSolve_BRGN;
364: 
365:   tao->data = (void*)gn;
366:   gn->reg_type = BRGN_REGULARIZATION_L2PROX;
367:   gn->lambda = 1e-4;
368:   gn->epsilon = 1e-6;
369:   gn->parent = tao;
370: 
371:   MatCreate(PetscObjectComm((PetscObject)tao),&gn->H);
372:   MatSetOptionsPrefix(gn->H,"tao_brgn_hessian_");
373: 
374:   TaoCreate(PetscObjectComm((PetscObject)tao),&gn->subsolver);
375:   TaoSetType(gn->subsolver,TAOBNLS);
376:   TaoSetOptionsPrefix(gn->subsolver,"tao_brgn_subsolver_");
377:   return(0);
378: }

380: /*@
381:   TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN

383:   Collective on Tao

385:   Level: advanced
386:   
387:   Input Parameters:
388: +  tao - the Tao solver context
389: -  subsolver - the Tao sub-solver context
390: @*/
391: PetscErrorCode TaoBRGNGetSubsolver(Tao tao,Tao *subsolver)
392: {
393:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
394: 
396:   *subsolver = gn->subsolver;
397:   return(0);
398: }

400: /*@
401:   TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm

403:   Collective on Tao
404:   
405:   Input Parameters:
406: +  tao - the Tao solver context
407: -  lambda - L1-norm regularizer weight

409:   Level: beginner
410: @*/
411: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao,PetscReal lambda)
412: {
413:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
414: 
415:   /* Initialize lambda here */

418:   gn->lambda = lambda;
419:   return(0);
420: }

422: /*@
423:   TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm

425:   Collective on Tao
426:   
427:   Input Parameters:
428: +  tao - the Tao solver context
429: -  epsilon - L1-norm smooth approximation parameter

431:   Level: advanced
432: @*/
433: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao,PetscReal epsilon)
434: {
435:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
436: 
437:   /* Initialize epsilon here */

440:   gn->epsilon = epsilon;
441:   return(0);
442: }

444: /*@
445:    TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)

447:    Input Parameters:
448: +  tao  - the Tao context
449: -  dict - the user specified dictionary matrix.  We allow to set a null dictionary, which means identity matrix by default

451:     Level: advanced
452: @*/
453: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao,Mat dict)
454: {
455:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
459:   if (dict) {
462:     PetscObjectReference((PetscObject)dict);
463:   }
464:   MatDestroy(&gn->D);
465:   gn->D = dict;
466:   return(0);
467: }

469: /*@C
470:    TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back 
471:    function into the algorithm.

473:    Input Parameters:
474:    + tao - the Tao context
475:    . func - function pointer for the regularizer value and gradient evaluation
476:    - ctx - user context for the regularizer

478:    Level: advanced
479: @*/
480: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao,PetscErrorCode (*func)(Tao,Vec,PetscReal *,Vec,void*),void *ctx)
481: {
482:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;

486:   if (ctx) {
487:     gn->reg_obj_ctx = ctx;
488:   }
489:   if (func) {
490:     gn->regularizerobjandgrad = func;
491:   }
492:   return(0);
493: }

495: /*@C
496:    TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back 
497:    function into the algorithm.

499:    Input Parameters:
500:    + tao - the Tao context
501:    . Hreg - user-created matrix for the Hessian of the regularization term
502:    . func - function pointer for the regularizer Hessian evaluation
503:    - ctx - user context for the regularizer Hessian

505:    Level: advanced
506: @*/
507: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao,Mat Hreg,PetscErrorCode (*func)(Tao,Vec,Mat,void*),void *ctx)
508: {
509:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;

514:   if (Hreg) {
517:   } else SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONG,"NULL Hessian detected! User must provide valid Hessian for the regularizer.");
518:   if (ctx) {
519:     gn->reg_hess_ctx = ctx;
520:   }
521:   if (func) {
522:     gn->regularizerhessian = func;
523:   }
524:   if (Hreg) {
525:     PetscObjectReference((PetscObject)Hreg);
526:     MatDestroy(&gn->Hreg);
527:     gn->Hreg = Hreg;
528:   }
529:   return(0);
530: }