1 /* Program usage: mpiexec -n 1 rosenbrock1 [-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 KSP ksp; 51 PC pc; 52 Mat M; 53 Vec in, out, out2; 54 PetscReal mult_solve_dist; 55 56 /* Initialize TAO and PETSc */ 57 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 58 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 59 if (size >1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors"); 60 61 /* Initialize problem parameters */ 62 user.n = 2; user.alpha = 99.0; user.chained = PETSC_FALSE; 63 /* Check for command line arguments to override defaults */ 64 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&user.n,&flg);CHKERRQ(ierr); 65 ierr = PetscOptionsGetReal(NULL,NULL,"-alpha",&user.alpha,&flg);CHKERRQ(ierr); 66 ierr = PetscOptionsGetBool(NULL,NULL,"-chained",&user.chained,&flg);CHKERRQ(ierr); 67 ierr = PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg);CHKERRQ(ierr); 68 69 /* Allocate vectors for the solution and gradient */ 70 ierr = VecCreateSeq(PETSC_COMM_SELF,user.n,&x);CHKERRQ(ierr); 71 ierr = MatCreateSeqBAIJ(PETSC_COMM_SELF,2,user.n,user.n,1,NULL,&H);CHKERRQ(ierr); 72 73 /* The TAO code begins here */ 74 75 /* Create TAO solver with desired solution method */ 76 ierr = TaoCreate(PETSC_COMM_SELF,&tao);CHKERRQ(ierr); 77 ierr = TaoSetType(tao,TAOLMVM);CHKERRQ(ierr); 78 79 /* Set solution vec and an initial guess */ 80 ierr = VecSet(x, zero);CHKERRQ(ierr); 81 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr); 82 83 /* Set routines for function, gradient, hessian evaluation */ 84 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,&user);CHKERRQ(ierr); 85 ierr = TaoSetHessianRoutine(tao,H,H,FormHessian,&user);CHKERRQ(ierr); 86 87 /* Test the LMVM matrix */ 88 if (test_lmvm) { 89 ierr = PetscOptionsSetValue(NULL, "-tao_type", "bqnktr");CHKERRQ(ierr); 90 } 91 92 /* Check for TAO command line options */ 93 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 94 95 /* SOLVE THE APPLICATION */ 96 ierr = TaoSolve(tao);CHKERRQ(ierr); 97 98 /* Test the LMVM matrix */ 99 if (test_lmvm) { 100 ierr = TaoGetKSP(tao, &ksp);CHKERRQ(ierr); 101 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 102 ierr = PCLMVMGetMatLMVM(pc, &M);CHKERRQ(ierr); 103 ierr = VecDuplicate(x, &in);CHKERRQ(ierr); 104 ierr = VecDuplicate(x, &out);CHKERRQ(ierr); 105 ierr = VecDuplicate(x, &out2);CHKERRQ(ierr); 106 ierr = VecSet(in, 1.0);CHKERRQ(ierr); 107 ierr = MatMult(M, in, out);CHKERRQ(ierr); 108 ierr = MatSolve(M, out, out2);CHKERRQ(ierr); 109 ierr = VecAXPY(out2, -1.0, in);CHKERRQ(ierr); 110 ierr = VecNorm(out2, NORM_2, &mult_solve_dist);CHKERRQ(ierr); 111 if (mult_solve_dist < 1.e-11) { 112 ierr = PetscPrintf(PetscObjectComm((PetscObject)tao), "error between LMVM MatMult and MatSolve: < 1.e-11\n");CHKERRQ(ierr); 113 } else if (mult_solve_dist < 1.e-6) { 114 ierr = PetscPrintf(PetscObjectComm((PetscObject)tao), "error between LMVM MatMult and MatSolve: < 1.e-6\n");CHKERRQ(ierr); 115 } else { 116 ierr = PetscPrintf(PetscObjectComm((PetscObject)tao), "error between LMVM MatMult and MatSolve: %e\n", (double)mult_solve_dist);CHKERRQ(ierr); 117 } 118 ierr = VecDestroy(&in);CHKERRQ(ierr); 119 ierr = VecDestroy(&out);CHKERRQ(ierr); 120 ierr = VecDestroy(&out2);CHKERRQ(ierr); 121 } 122 123 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 124 ierr = VecDestroy(&x);CHKERRQ(ierr); 125 ierr = MatDestroy(&H);CHKERRQ(ierr); 126 127 ierr = PetscFinalize(); 128 return ierr; 129 } 130 131 /* -------------------------------------------------------------------- */ 132 /* 133 FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X). 134 135 Input Parameters: 136 . tao - the Tao context 137 . X - input vector 138 . ptr - optional user-defined context, as set by TaoSetFunctionGradient() 139 140 Output Parameters: 141 . G - vector containing the newly evaluated gradient 142 . f - function value 143 144 Note: 145 Some optimization methods ask for the function and the gradient evaluation 146 at the same time. Evaluating both at once may be more efficient that 147 evaluating each separately. 148 */ 149 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f, Vec G,void *ptr) 150 { 151 AppCtx *user = (AppCtx *) ptr; 152 PetscInt i,nn=user->n/2; 153 PetscErrorCode ierr; 154 PetscReal ff=0,t1,t2,alpha=user->alpha; 155 PetscScalar *g; 156 const PetscScalar *x; 157 158 PetscFunctionBeginUser; 159 /* Get pointers to vector data */ 160 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 161 ierr = VecGetArray(G,&g);CHKERRQ(ierr); 162 163 /* Compute G(X) */ 164 if (user->chained) { 165 g[0] = 0; 166 for (i=0; i<user->n-1; i++) { 167 t1 = x[i+1] - x[i]*x[i]; 168 ff += PetscSqr(1 - x[i]) + alpha*t1*t1; 169 g[i] += -2*(1 - x[i]) + 2*alpha*t1*(-2*x[i]); 170 g[i+1] = 2*alpha*t1; 171 } 172 } else { 173 for (i=0; i<nn; i++) { 174 t1 = x[2*i+1]-x[2*i]*x[2*i]; t2= 1-x[2*i]; 175 ff += alpha*t1*t1 + t2*t2; 176 g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2; 177 g[2*i+1] = 2*alpha*t1; 178 } 179 } 180 181 /* Restore vectors */ 182 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 183 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr); 184 *f = ff; 185 186 ierr = PetscLogFlops(15.0*nn);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 /* ------------------------------------------------------------------- */ 191 /* 192 FormHessian - Evaluates Hessian matrix. 193 194 Input Parameters: 195 . tao - the Tao context 196 . x - input vector 197 . ptr - optional user-defined context, as set by TaoSetHessian() 198 199 Output Parameters: 200 . H - Hessian matrix 201 202 Note: Providing the Hessian may not be necessary. Only some solvers 203 require this matrix. 204 */ 205 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 206 { 207 AppCtx *user = (AppCtx*)ptr; 208 PetscErrorCode ierr; 209 PetscInt i, ind[2]; 210 PetscReal alpha=user->alpha; 211 PetscReal v[2][2]; 212 const PetscScalar *x; 213 PetscBool assembled; 214 215 PetscFunctionBeginUser; 216 /* Zero existing matrix entries */ 217 ierr = MatAssembled(H,&assembled);CHKERRQ(ierr); 218 if (assembled) {ierr = MatZeroEntries(H);CHKERRQ(ierr);} 219 220 /* Get a pointer to vector data */ 221 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 222 223 /* Compute H(X) entries */ 224 if (user->chained) { 225 ierr = MatZeroEntries(H);CHKERRQ(ierr); 226 for (i=0; i<user->n-1; i++) { 227 PetscScalar t1 = x[i+1] - x[i]*x[i]; 228 v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]); 229 v[0][1] = 2*alpha*(-2*x[i]); 230 v[1][0] = 2*alpha*(-2*x[i]); 231 v[1][1] = 2*alpha*t1; 232 ind[0] = i; ind[1] = i+1; 233 ierr = MatSetValues(H,2,ind,2,ind,v[0],ADD_VALUES);CHKERRQ(ierr); 234 } 235 } else { 236 for (i=0; i<user->n/2; i++) { 237 v[1][1] = 2*alpha; 238 v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2; 239 v[1][0] = v[0][1] = -4.0*alpha*x[2*i]; 240 ind[0]=2*i; ind[1]=2*i+1; 241 ierr = MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES);CHKERRQ(ierr); 242 } 243 } 244 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 245 246 /* Assemble matrix */ 247 ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 248 ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 249 ierr = PetscLogFlops(9.0*user->n/2.0);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 /*TEST 254 255 build: 256 requires: !complex 257 258 test: 259 args: -tao_smonitor -tao_type nls -tao_gatol 1.e-4 260 requires: !single 261 262 test: 263 suffix: 2 264 args: -tao_smonitor -tao_type lmvm -tao_gatol 1.e-3 265 266 test: 267 suffix: 3 268 args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4 269 requires: !single 270 271 test: 272 suffix: 4 273 args: -tao_smonitor -tao_type ntr -tao_mf_hessian -tao_ntr_pc_type none -tao_gatol 1.e-4 274 275 test: 276 suffix: 5 277 args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 278 279 test: 280 suffix: 6 281 args: -tao_smonitor -tao_type bntl -tao_gatol 1.e-4 282 283 test: 284 suffix: 7 285 args: -tao_smonitor -tao_type bnls -tao_gatol 1.e-4 286 287 test: 288 suffix: 8 289 args: -tao_smonitor -tao_type bntr -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 290 291 test: 292 suffix: 9 293 args: -tao_smonitor -tao_type bntl -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 294 295 test: 296 suffix: 10 297 args: -tao_smonitor -tao_type bnls -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 298 299 test: 300 suffix: 11 301 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbroyden 302 303 test: 304 suffix: 12 305 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbadbroyden 306 307 test: 308 suffix: 13 309 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbroyden 310 311 test: 312 suffix: 14 313 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbfgs 314 315 test: 316 suffix: 15 317 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmdfp 318 319 test: 320 suffix: 16 321 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsr1 322 323 test: 324 suffix: 17 325 args: -tao_smonitor -tao_gatol 1e-4 -tao_type bqnls 326 327 test: 328 suffix: 18 329 args: -tao_smonitor -tao_gatol 1e-4 -tao_type blmvm 330 331 test: 332 suffix: 19 333 args: -tao_smonitor -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1 334 335 test: 336 suffix: 20 337 args: -tao_monitor -tao_gatol 1e-4 -tao_type blmvm -tao_ls_monitor 338 339 test: 340 suffix: 21 341 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbadbroyden 342 343 test: 344 suffix: 22 345 args: -tao_max_it 1 -tao_converged_reason 346 347 test: 348 suffix: 23 349 args: -tao_max_funcs 0 -tao_converged_reason 350 351 test: 352 suffix: 24 353 args: -tao_gatol 10 -tao_converged_reason 354 355 test: 356 suffix: 25 357 args: -tao_grtol 10 -tao_converged_reason 358 359 test: 360 suffix: 26 361 args: -tao_gttol 10 -tao_converged_reason 362 363 test: 364 suffix: 27 365 args: -tao_steptol 10 -tao_converged_reason 366 367 test: 368 suffix: 28 369 args: -tao_fmin 10 -tao_converged_reason 370 371 TEST*/ 372