1 #include "ssls.h" 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "TaoSetUp_SSILS" 5 PetscErrorCode TaoSetUp_SSILS(TaoSolver tao) 6 { 7 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); 12 ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); 13 ierr = VecDuplicate(tao->solution,&ssls->ff);CHKERRQ(ierr); 14 ierr = VecDuplicate(tao->solution,&ssls->dpsi);CHKERRQ(ierr); 15 ierr = VecDuplicate(tao->solution,&ssls->da);CHKERRQ(ierr); 16 ierr = VecDuplicate(tao->solution,&ssls->db);CHKERRQ(ierr); 17 ierr = VecDuplicate(tao->solution,&ssls->t1);CHKERRQ(ierr); 18 ierr = VecDuplicate(tao->solution,&ssls->t2);CHKERRQ(ierr); 19 PetscFunctionReturn(0); 20 } 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "TaoDestroy_SSILS" 24 PetscErrorCode TaoDestroy_SSILS(TaoSolver tao) 25 { 26 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 ierr = VecDestroy(&ssls->ff);CHKERRQ(ierr); 31 ierr = VecDestroy(&ssls->dpsi);CHKERRQ(ierr); 32 ierr = VecDestroy(&ssls->da);CHKERRQ(ierr); 33 ierr = VecDestroy(&ssls->db);CHKERRQ(ierr); 34 ierr = VecDestroy(&ssls->t1);CHKERRQ(ierr); 35 ierr = VecDestroy(&ssls->t2);CHKERRQ(ierr); 36 ierr = PetscFree(tao->data);CHKERRQ(ierr); 37 tao->data = PETSC_NULL; 38 PetscFunctionReturn(0); 39 } 40 41 #undef __FUNCT__ 42 #define __FUNCT__ "TaoSolve_SSILS" 43 static PetscErrorCode TaoSolve_SSILS(TaoSolver tao) 44 { 45 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 46 PetscReal psi, ndpsi, normd, innerd, t=0; 47 PetscReal delta, rho; 48 PetscInt iter=0,kspits; 49 TaoSolverTerminationReason reason; 50 TaoLineSearchTerminationReason ls_reason; 51 PetscErrorCode ierr; 52 53 PetscFunctionBegin; 54 /* Assume that Setup has been called! 55 Set the structure for the Jacobian and create a linear solver. */ 56 delta = ssls->delta; 57 rho = ssls->rho; 58 59 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 60 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 61 ierr = TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);CHKERRQ(ierr); 62 ierr = TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);CHKERRQ(ierr); 63 64 /* Calculate the function value and fischer function value at the 65 current iterate */ 66 ierr = TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);CHKERRQ(ierr); 67 ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr); 68 69 while (1) { 70 ierr=PetscInfo3(tao, "iter: %D, merit: %g, ndpsi: %g\n",iter, (double)ssls->merit, (double)ndpsi);CHKERRQ(ierr); 71 /* Check the termination criteria */ 72 ierr = TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason);CHKERRQ(ierr); 73 if (reason!=TAO_CONTINUE_ITERATING) break; 74 75 /* Calculate direction. (Really negative of newton direction. Therefore, 76 rest of the code uses -d.) */ 77 ierr = KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre,ssls->matflag);CHKERRQ(ierr); 78 ierr = KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);CHKERRQ(ierr); 79 ierr = KSPGetIterationNumber(tao->ksp,&kspits);CHKERRQ(ierr); 80 tao->ksp_its+=kspits; 81 ierr = VecNorm(tao->stepdirection,NORM_2,&normd);CHKERRQ(ierr); 82 ierr = VecDot(tao->stepdirection,ssls->dpsi,&innerd);CHKERRQ(ierr); 83 84 /* Make sure that we have a descent direction */ 85 if (innerd <= delta*pow(normd, rho)) { 86 ierr = PetscInfo(tao, "newton direction not descent\n");CHKERRQ(ierr); 87 ierr = VecCopy(ssls->dpsi,tao->stepdirection);CHKERRQ(ierr); 88 ierr = VecDot(tao->stepdirection,ssls->dpsi,&innerd);CHKERRQ(ierr); 89 } 90 91 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 92 innerd = -innerd; 93 94 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); 95 ierr = TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);CHKERRQ(ierr); 96 ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr); 97 } 98 PetscFunctionReturn(0); 99 } 100 101 /* ---------------------------------------------------------- */ 102 EXTERN_C_BEGIN 103 #undef __FUNCT__ 104 #define __FUNCT__ "TaoCreate_SSILS" 105 PetscErrorCode TaoCreate_SSILS(TaoSolver tao) 106 { 107 TAO_SSLS *ssls; 108 PetscErrorCode ierr; 109 const char *armijo_type = TAOLINESEARCH_ARMIJO; 110 111 PetscFunctionBegin; 112 ierr = PetscNewLog(tao,&ssls);CHKERRQ(ierr); 113 tao->data = (void*)ssls; 114 tao->ops->solve=TaoSolve_SSILS; 115 tao->ops->setup=TaoSetUp_SSILS; 116 tao->ops->view=TaoView_SSLS; 117 tao->ops->setfromoptions=TaoSetFromOptions_SSLS; 118 tao->ops->destroy = TaoDestroy_SSILS; 119 120 ssls->delta = 1e-10; 121 ssls->rho = 2.1; 122 123 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr); 124 ierr = TaoLineSearchSetType(tao->linesearch,armijo_type);CHKERRQ(ierr); 125 ierr = TaoLineSearchSetFromOptions(tao->linesearch); 126 /* Note: linesearch objective and objectivegradient routines are set in solve routine */ 127 ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp);CHKERRQ(ierr); 128 129 tao->max_it = 2000; 130 tao->max_funcs = 4000; 131 tao->fatol = 0; tao->frtol = 0; tao->gttol=0; tao->grtol=0; 132 tao->gatol = 1.0e-16; 133 tao->fmin = 1.0e-8; 134 PetscFunctionReturn(0); 135 } 136 EXTERN_C_END 137 138