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