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