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(); TaoSetObjectiveAndGradientRoutine(); 17 Routines: TaoSetHessianRoutine(); 18 Routines: TaoSetInitialVector(); 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 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 56 if (size >1) SETERRQ(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 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&user.n,&flg);CHKERRQ(ierr); 62 ierr = PetscOptionsGetReal(NULL,NULL,"-alpha",&user.alpha,&flg);CHKERRQ(ierr); 63 ierr = PetscOptionsGetBool(NULL,NULL,"-chained",&user.chained,&flg);CHKERRQ(ierr); 64 ierr = PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg);CHKERRQ(ierr); 65 66 /* Allocate vectors for the solution and gradient */ 67 ierr = VecCreateSeq(PETSC_COMM_SELF,user.n,&x);CHKERRQ(ierr); 68 ierr = MatCreateSeqBAIJ(PETSC_COMM_SELF,2,user.n,user.n,1,NULL,&H);CHKERRQ(ierr); 69 70 /* The TAO code begins here */ 71 72 /* Create TAO solver with desired solution method */ 73 ierr = TaoCreate(PETSC_COMM_SELF,&tao);CHKERRQ(ierr); 74 ierr = TaoSetType(tao,TAOBQNLS);CHKERRQ(ierr); 75 76 /* Set solution vec and an initial guess */ 77 ierr = VecSet(x, zero);CHKERRQ(ierr); 78 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr); 79 80 /* Set routines for function, gradient, hessian evaluation */ 81 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,&user);CHKERRQ(ierr); 82 ierr = TaoSetHessianRoutine(tao,H,H,FormHessian,&user);CHKERRQ(ierr); 83 84 /* Check for TAO command line options */ 85 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 86 87 /* Solve the problem */ 88 ierr = TaoSetTolerances(tao, 1.e-5, 0.0, 0.0);CHKERRQ(ierr); 89 ierr = TaoSetMaximumIterations(tao, 5);CHKERRQ(ierr); 90 ierr = TaoSetRecycleHistory(tao, PETSC_TRUE);CHKERRQ(ierr); 91 reason = TAO_CONTINUE_ITERATING; 92 flg = PETSC_FALSE; 93 ierr = TaoGetRecycleHistory(tao, &flg);CHKERRQ(ierr); 94 if (flg) {ierr = PetscPrintf(PETSC_COMM_SELF, "Recycle: enabled\n");CHKERRQ(ierr);} 95 while (reason != TAO_CONVERGED_GATOL) { 96 ierr = TaoSolve(tao);CHKERRQ(ierr); 97 ierr = TaoGetConvergedReason(tao, &reason);CHKERRQ(ierr); 98 ierr = TaoGetIterationNumber(tao, &its);CHKERRQ(ierr); 99 recycled_its += its; 100 ierr = PetscPrintf(PETSC_COMM_SELF, "-----------------------\n");CHKERRQ(ierr); 101 } 102 103 /* Disable recycling and solve again! */ 104 ierr = TaoSetMaximumIterations(tao, 100);CHKERRQ(ierr); 105 ierr = TaoSetRecycleHistory(tao, PETSC_FALSE);CHKERRQ(ierr); 106 ierr = VecSet(x, zero);CHKERRQ(ierr); 107 ierr = TaoGetRecycleHistory(tao, &flg);CHKERRQ(ierr); 108 if (!flg) {ierr = PetscPrintf(PETSC_COMM_SELF, "Recycle: disabled\n");CHKERRQ(ierr);} 109 ierr = TaoSolve(tao);CHKERRQ(ierr); 110 ierr = TaoGetConvergedReason(tao, &reason);CHKERRQ(ierr); 111 if (reason != TAO_CONVERGED_GATOL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Solution failed to converge!"); 112 ierr = TaoGetIterationNumber(tao, &oneshot_its);CHKERRQ(ierr); 113 ierr = PetscPrintf(PETSC_COMM_SELF, "-----------------------\n");CHKERRQ(ierr); 114 ierr = PetscPrintf(PETSC_COMM_SELF, "recycled its: %D | oneshot its: %D\n", recycled_its, oneshot_its);CHKERRQ(ierr); 115 if (recycled_its != oneshot_its) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Recycled solution does not match oneshot solution!"); 116 117 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 118 ierr = VecDestroy(&x);CHKERRQ(ierr); 119 ierr = MatDestroy(&H);CHKERRQ(ierr); 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 PetscErrorCode ierr; 148 PetscReal ff=0,t1,t2,alpha=user->alpha; 149 PetscScalar *g; 150 const PetscScalar *x; 151 152 PetscFunctionBeginUser; 153 /* Get pointers to vector data */ 154 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 155 ierr = VecGetArrayWrite(G,&g);CHKERRQ(ierr); 156 157 /* Compute G(X) */ 158 if (user->chained) { 159 g[0] = 0; 160 for (i=0; i<user->n-1; i++) { 161 t1 = x[i+1] - x[i]*x[i]; 162 ff += PetscSqr(1 - x[i]) + alpha*t1*t1; 163 g[i] += -2*(1 - x[i]) + 2*alpha*t1*(-2*x[i]); 164 g[i+1] = 2*alpha*t1; 165 } 166 } else { 167 for (i=0; i<nn; i++) { 168 t1 = x[2*i+1]-x[2*i]*x[2*i]; t2= 1-x[2*i]; 169 ff += alpha*t1*t1 + t2*t2; 170 g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2; 171 g[2*i+1] = 2*alpha*t1; 172 } 173 } 174 175 /* Restore vectors */ 176 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 177 ierr = VecRestoreArrayWrite(G,&g);CHKERRQ(ierr); 178 *f = ff; 179 180 ierr = PetscLogFlops(15.0*nn);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 /* ------------------------------------------------------------------- */ 185 /* 186 FormHessian - Evaluates Hessian matrix. 187 188 Input Parameters: 189 . tao - the Tao context 190 . x - input vector 191 . ptr - optional user-defined context, as set by TaoSetHessian() 192 193 Output Parameters: 194 . H - Hessian matrix 195 196 Note: Providing the Hessian may not be necessary. Only some solvers 197 require this matrix. 198 */ 199 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 200 { 201 AppCtx *user = (AppCtx*)ptr; 202 PetscErrorCode ierr; 203 PetscInt i, ind[2]; 204 PetscReal alpha=user->alpha; 205 PetscReal v[2][2]; 206 const PetscScalar *x; 207 PetscBool assembled; 208 209 PetscFunctionBeginUser; 210 /* Zero existing matrix entries */ 211 ierr = MatAssembled(H,&assembled);CHKERRQ(ierr); 212 if (assembled || user->chained) {ierr = MatZeroEntries(H);CHKERRQ(ierr);} 213 214 /* Get a pointer to vector data */ 215 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 216 217 /* Compute H(X) entries */ 218 if (user->chained) { 219 for (i=0; i<user->n-1; i++) { 220 PetscScalar t1 = x[i+1] - x[i]*x[i]; 221 v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]); 222 v[0][1] = 2*alpha*(-2*x[i]); 223 v[1][0] = 2*alpha*(-2*x[i]); 224 v[1][1] = 2*alpha*t1; 225 ind[0] = i; ind[1] = i+1; 226 ierr = MatSetValues(H,2,ind,2,ind,v[0],ADD_VALUES);CHKERRQ(ierr); 227 } 228 } else { 229 for (i=0; i<user->n/2; i++) { 230 v[1][1] = 2*alpha; 231 v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2; 232 v[1][0] = v[0][1] = -4.0*alpha*x[2*i]; 233 ind[0]=2*i; ind[1]=2*i+1; 234 ierr = MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES);CHKERRQ(ierr); 235 } 236 } 237 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 238 239 /* Assemble matrix */ 240 ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 241 ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 242 ierr = PetscLogFlops(9.0*user->n/2.0);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 /*TEST 247 248 build: 249 requires: !complex 250 251 test: 252 args: -tao_type bqnls -tao_monitor 253 requires: !single 254 255 TEST*/ 256