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