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