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