1 /* Program usage: mpiexec -n 1 rosenbrock2 [-help] [all TAO options] */ 2 3 /* Include "petsctao.h" so we can use TAO solvers. */ 4 #include <petsctao.h> 5 6 static char help[] = "This example demonstrates use of the TAO package to \n\ 7 solve an unconstrained minimization problem on a single processor. We \n\ 8 minimize the extended Rosenbrock function: \n\ 9 sum_{i=0}^{n/2-1} (alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2) \n\ 10 or the chained Rosenbrock function:\n\ 11 sum_{i=0}^{n-1} alpha*(x_{i+1} - x_i^2)^2 + (1 - x_i)^2\n"; 12 13 /*T 14 Concepts: TAO^Solving an unconstrained minimization problem 15 Routines: TaoCreate(); 16 Routines: TaoSetType(); TaoSetObjectiveAndGradient(); 17 Routines: TaoSetHessian(); 18 Routines: TaoSetSolution(); 19 Routines: TaoSetFromOptions(); 20 Routines: TaoSolve(); 21 Routines: TaoDestroy(); 22 Processors: 1 23 T*/ 24 25 /* 26 User-defined application context - contains data needed by the 27 application-provided call-back routines that evaluate the function, 28 gradient, and hessian. 29 */ 30 typedef struct { 31 PetscInt n; /* dimension */ 32 PetscReal alpha; /* condition parameter */ 33 PetscBool chained; 34 } AppCtx; 35 36 /* -------------- User-defined routines ---------- */ 37 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 38 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*); 39 40 int main(int argc,char **argv) 41 { 42 PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 43 PetscReal zero=0.0; 44 Vec x; /* solution vector */ 45 Mat H; 46 Tao tao; /* Tao solver context */ 47 PetscBool flg, test_lmvm = PETSC_FALSE; 48 PetscMPIInt size; /* number of processes running */ 49 AppCtx user; /* user-defined application context */ 50 TaoConvergedReason reason; 51 PetscInt its, recycled_its=0, oneshot_its=0; 52 53 /* Initialize TAO and PETSc */ 54 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 55 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 56 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors"); 57 58 /* Initialize problem parameters */ 59 user.n = 2; user.alpha = 99.0; user.chained = PETSC_FALSE; 60 /* Check for command line arguments to override defaults */ 61 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&user.n,&flg)); 62 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-alpha",&user.alpha,&flg)); 63 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-chained",&user.chained,&flg)); 64 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg)); 65 66 /* Allocate vectors for the solution and gradient */ 67 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,user.n,&x)); 68 CHKERRQ(MatCreateSeqBAIJ(PETSC_COMM_SELF,2,user.n,user.n,1,NULL,&H)); 69 70 /* The TAO code begins here */ 71 72 /* Create TAO solver with desired solution method */ 73 CHKERRQ(TaoCreate(PETSC_COMM_SELF,&tao)); 74 CHKERRQ(TaoSetType(tao,TAOLMVM)); 75 76 /* Set solution vec and an initial guess */ 77 CHKERRQ(VecSet(x, zero)); 78 CHKERRQ(TaoSetSolution(tao,x)); 79 80 /* Set routines for function, gradient, hessian evaluation */ 81 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,&user)); 82 CHKERRQ(TaoSetHessian(tao,H,H,FormHessian,&user)); 83 84 /* Check for TAO command line options */ 85 CHKERRQ(TaoSetFromOptions(tao)); 86 87 /* Solve the problem */ 88 CHKERRQ(TaoSetTolerances(tao, 1.e-5, 0.0, 0.0)); 89 CHKERRQ(TaoSetMaximumIterations(tao, 5)); 90 CHKERRQ(TaoLMVMRecycle(tao, PETSC_TRUE)); 91 reason = TAO_CONTINUE_ITERATING; 92 while (reason != TAO_CONVERGED_GATOL) { 93 CHKERRQ(TaoSolve(tao)); 94 CHKERRQ(TaoGetConvergedReason(tao, &reason)); 95 CHKERRQ(TaoGetIterationNumber(tao, &its)); 96 recycled_its += its; 97 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n")); 98 } 99 100 /* Disable recycling and solve again! */ 101 CHKERRQ(TaoSetMaximumIterations(tao, 100)); 102 CHKERRQ(TaoLMVMRecycle(tao, PETSC_FALSE)); 103 CHKERRQ(VecSet(x, zero)); 104 CHKERRQ(TaoSolve(tao)); 105 CHKERRQ(TaoGetConvergedReason(tao, &reason)); 106 PetscCheck(reason == TAO_CONVERGED_GATOL,PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Solution failed to converge!"); 107 CHKERRQ(TaoGetIterationNumber(tao, &oneshot_its)); 108 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n")); 109 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "recycled its: %D | oneshot its: %D\n", recycled_its, oneshot_its)); 110 PetscCheck(recycled_its == oneshot_its,PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "LMVM recycling does not work!"); 111 112 CHKERRQ(TaoDestroy(&tao)); 113 CHKERRQ(VecDestroy(&x)); 114 CHKERRQ(MatDestroy(&H)); 115 116 ierr = PetscFinalize(); 117 return ierr; 118 } 119 120 /* -------------------------------------------------------------------- */ 121 /* 122 FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X). 123 124 Input Parameters: 125 . tao - the Tao context 126 . X - input vector 127 . ptr - optional user-defined context, as set by TaoSetFunctionGradient() 128 129 Output Parameters: 130 . G - vector containing the newly evaluated gradient 131 . f - function value 132 133 Note: 134 Some optimization methods ask for the function and the gradient evaluation 135 at the same time. Evaluating both at once may be more efficient that 136 evaluating each separately. 137 */ 138 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f, Vec G,void *ptr) 139 { 140 AppCtx *user = (AppCtx *) ptr; 141 PetscInt i,nn=user->n/2; 142 PetscReal ff=0,t1,t2,alpha=user->alpha; 143 PetscScalar *g; 144 const PetscScalar *x; 145 146 PetscFunctionBeginUser; 147 /* Get pointers to vector data */ 148 CHKERRQ(VecGetArrayRead(X,&x)); 149 CHKERRQ(VecGetArray(G,&g)); 150 151 /* Compute G(X) */ 152 if (user->chained) { 153 g[0] = 0; 154 for (i=0; i<user->n-1; i++) { 155 t1 = x[i+1] - x[i]*x[i]; 156 ff += PetscSqr(1 - x[i]) + alpha*t1*t1; 157 g[i] += -2*(1 - x[i]) + 2*alpha*t1*(-2*x[i]); 158 g[i+1] = 2*alpha*t1; 159 } 160 } else { 161 for (i=0; i<nn; i++) { 162 t1 = x[2*i+1]-x[2*i]*x[2*i]; t2= 1-x[2*i]; 163 ff += alpha*t1*t1 + t2*t2; 164 g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2; 165 g[2*i+1] = 2*alpha*t1; 166 } 167 } 168 169 /* Restore vectors */ 170 CHKERRQ(VecRestoreArrayRead(X,&x)); 171 CHKERRQ(VecRestoreArray(G,&g)); 172 *f = ff; 173 174 CHKERRQ(PetscLogFlops(15.0*nn)); 175 PetscFunctionReturn(0); 176 } 177 178 /* ------------------------------------------------------------------- */ 179 /* 180 FormHessian - Evaluates Hessian matrix. 181 182 Input Parameters: 183 . tao - the Tao context 184 . x - input vector 185 . ptr - optional user-defined context, as set by TaoSetHessian() 186 187 Output Parameters: 188 . H - Hessian matrix 189 190 Note: Providing the Hessian may not be necessary. Only some solvers 191 require this matrix. 192 */ 193 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 194 { 195 AppCtx *user = (AppCtx*)ptr; 196 PetscInt i, ind[2]; 197 PetscReal alpha=user->alpha; 198 PetscReal v[2][2]; 199 const PetscScalar *x; 200 PetscBool assembled; 201 202 PetscFunctionBeginUser; 203 /* Zero existing matrix entries */ 204 CHKERRQ(MatAssembled(H,&assembled)); 205 if (assembled) CHKERRQ(MatZeroEntries(H)); 206 207 /* Get a pointer to vector data */ 208 CHKERRQ(VecGetArrayRead(X,&x)); 209 210 /* Compute H(X) entries */ 211 if (user->chained) { 212 CHKERRQ(MatZeroEntries(H)); 213 for (i=0; i<user->n-1; i++) { 214 PetscScalar t1 = x[i+1] - x[i]*x[i]; 215 v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]); 216 v[0][1] = 2*alpha*(-2*x[i]); 217 v[1][0] = 2*alpha*(-2*x[i]); 218 v[1][1] = 2*alpha*t1; 219 ind[0] = i; ind[1] = i+1; 220 CHKERRQ(MatSetValues(H,2,ind,2,ind,v[0],ADD_VALUES)); 221 } 222 } else { 223 for (i=0; i<user->n/2; i++) { 224 v[1][1] = 2*alpha; 225 v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2; 226 v[1][0] = v[0][1] = -4.0*alpha*x[2*i]; 227 ind[0]=2*i; ind[1]=2*i+1; 228 CHKERRQ(MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES)); 229 } 230 } 231 CHKERRQ(VecRestoreArrayRead(X,&x)); 232 233 /* Assemble matrix */ 234 CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 235 CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 236 CHKERRQ(PetscLogFlops(9.0*user->n/2.0)); 237 PetscFunctionReturn(0); 238 } 239 240 /*TEST 241 242 build: 243 requires: !complex 244 245 test: 246 args: -tao_type lmvm -tao_monitor 247 requires: !single 248 249 TEST*/ 250