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,TAOBQNLS)); 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(TaoSetRecycleHistory(tao, PETSC_TRUE)); 79 reason = TAO_CONTINUE_ITERATING; 80 flg = PETSC_FALSE; 81 PetscCall(TaoGetRecycleHistory(tao, &flg)); 82 if (flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Recycle: enabled\n")); 83 while (reason != TAO_CONVERGED_GATOL) { 84 PetscCall(TaoSolve(tao)); 85 PetscCall(TaoGetConvergedReason(tao, &reason)); 86 PetscCall(TaoGetIterationNumber(tao, &its)); 87 recycled_its += its; 88 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n")); 89 } 90 91 /* Disable recycling and solve again! */ 92 PetscCall(TaoSetMaximumIterations(tao, 100)); 93 PetscCall(TaoSetRecycleHistory(tao, PETSC_FALSE)); 94 PetscCall(VecSet(x, zero)); 95 PetscCall(TaoGetRecycleHistory(tao, &flg)); 96 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Recycle: disabled\n")); 97 PetscCall(TaoSolve(tao)); 98 PetscCall(TaoGetConvergedReason(tao, &reason)); 99 PetscCheck(reason == TAO_CONVERGED_GATOL,PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Solution failed to converge!"); 100 PetscCall(TaoGetIterationNumber(tao, &oneshot_its)); 101 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n")); 102 PetscCall(PetscPrintf(PETSC_COMM_SELF, "recycled its: %" PetscInt_FMT " | oneshot its: %" PetscInt_FMT "\n", recycled_its, oneshot_its)); 103 PetscCheck(recycled_its == oneshot_its,PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Recycled solution does not match oneshot solution!"); 104 105 PetscCall(TaoDestroy(&tao)); 106 PetscCall(VecDestroy(&x)); 107 PetscCall(MatDestroy(&H)); 108 109 PetscCall(PetscFinalize()); 110 return 0; 111 } 112 113 /* -------------------------------------------------------------------- */ 114 /* 115 FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X). 116 117 Input Parameters: 118 . tao - the Tao context 119 . X - input vector 120 . ptr - optional user-defined context, as set by TaoSetFunctionGradient() 121 122 Output Parameters: 123 . G - vector containing the newly evaluated gradient 124 . f - function value 125 126 Note: 127 Some optimization methods ask for the function and the gradient evaluation 128 at the same time. Evaluating both at once may be more efficient than 129 evaluating each separately. 130 */ 131 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f, Vec G,void *ptr) 132 { 133 AppCtx *user = (AppCtx *) ptr; 134 PetscInt i,nn=user->n/2; 135 PetscReal ff=0,t1,t2,alpha=user->alpha; 136 PetscScalar *g; 137 const PetscScalar *x; 138 139 PetscFunctionBeginUser; 140 /* Get pointers to vector data */ 141 PetscCall(VecGetArrayRead(X,&x)); 142 PetscCall(VecGetArrayWrite(G,&g)); 143 144 /* Compute G(X) */ 145 if (user->chained) { 146 g[0] = 0; 147 for (i=0; i<user->n-1; i++) { 148 t1 = x[i+1] - x[i]*x[i]; 149 ff += PetscSqr(1 - x[i]) + alpha*t1*t1; 150 g[i] += -2*(1 - x[i]) + 2*alpha*t1*(-2*x[i]); 151 g[i+1] = 2*alpha*t1; 152 } 153 } else { 154 for (i=0; i<nn; i++) { 155 t1 = x[2*i+1]-x[2*i]*x[2*i]; t2= 1-x[2*i]; 156 ff += alpha*t1*t1 + t2*t2; 157 g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2; 158 g[2*i+1] = 2*alpha*t1; 159 } 160 } 161 162 /* Restore vectors */ 163 PetscCall(VecRestoreArrayRead(X,&x)); 164 PetscCall(VecRestoreArrayWrite(G,&g)); 165 *f = ff; 166 167 PetscCall(PetscLogFlops(15.0*nn)); 168 PetscFunctionReturn(0); 169 } 170 171 /* ------------------------------------------------------------------- */ 172 /* 173 FormHessian - Evaluates Hessian matrix. 174 175 Input Parameters: 176 . tao - the Tao context 177 . x - input vector 178 . ptr - optional user-defined context, as set by TaoSetHessian() 179 180 Output Parameters: 181 . H - Hessian matrix 182 183 Note: Providing the Hessian may not be necessary. Only some solvers 184 require this matrix. 185 */ 186 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 187 { 188 AppCtx *user = (AppCtx*)ptr; 189 PetscInt i, ind[2]; 190 PetscReal alpha=user->alpha; 191 PetscReal v[2][2]; 192 const PetscScalar *x; 193 PetscBool assembled; 194 195 PetscFunctionBeginUser; 196 /* Zero existing matrix entries */ 197 PetscCall(MatAssembled(H,&assembled)); 198 if (assembled || user->chained) PetscCall(MatZeroEntries(H)); 199 200 /* Get a pointer to vector data */ 201 PetscCall(VecGetArrayRead(X,&x)); 202 203 /* Compute H(X) entries */ 204 if (user->chained) { 205 for (i=0; i<user->n-1; i++) { 206 PetscScalar t1 = x[i+1] - x[i]*x[i]; 207 v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]); 208 v[0][1] = 2*alpha*(-2*x[i]); 209 v[1][0] = 2*alpha*(-2*x[i]); 210 v[1][1] = 2*alpha*t1; 211 ind[0] = i; ind[1] = i+1; 212 PetscCall(MatSetValues(H,2,ind,2,ind,v[0],ADD_VALUES)); 213 } 214 } else { 215 for (i=0; i<user->n/2; i++) { 216 v[1][1] = 2*alpha; 217 v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2; 218 v[1][0] = v[0][1] = -4.0*alpha*x[2*i]; 219 ind[0]=2*i; ind[1]=2*i+1; 220 PetscCall(MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES)); 221 } 222 } 223 PetscCall(VecRestoreArrayRead(X,&x)); 224 225 /* Assemble matrix */ 226 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 227 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 228 PetscCall(PetscLogFlops(9.0*user->n/2.0)); 229 PetscFunctionReturn(0); 230 } 231 232 /*TEST 233 234 build: 235 requires: !complex 236 237 test: 238 args: -tao_type bqnls -tao_monitor 239 requires: !single 240 241 TEST*/ 242