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