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