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