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