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, (char *)0, 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), "error between LMVM MatMult and MatSolve: < 1.e-11\n")); 101 } else if (mult_solve_dist < 1.e-6) { 102 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "error between LMVM MatMult and MatSolve: < 1.e-6\n")); 103 } else { 104 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "error between LMVM MatMult and MatSolve: %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 args: -tao_monitor_short -tao_type nls -tao_gatol 1.e-4 249 requires: !single 250 251 test: 252 suffix: 2 253 args: -tao_monitor_short -tao_type lmvm -tao_gatol 1.e-3 254 255 test: 256 suffix: 3 257 args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-4 258 requires: !single 259 260 test: 261 suffix: 4 262 args: -tao_monitor_short -tao_type ntr -tao_mf_hessian -tao_ntr_pc_type none -tao_gatol 1.e-4 263 264 test: 265 suffix: 5 266 args: -tao_monitor_short -tao_type bntr -tao_gatol 1.e-4 267 268 test: 269 suffix: 6 270 args: -tao_monitor_short -tao_type bntl -tao_gatol 1.e-4 271 272 test: 273 suffix: 7 274 args: -tao_monitor_short -tao_type bnls -tao_gatol 1.e-4 275 276 test: 277 suffix: 8 278 args: -tao_monitor_short -tao_type bntr -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 279 280 test: 281 suffix: 9 282 args: -tao_monitor_short -tao_type bntl -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 283 284 test: 285 suffix: 10 286 args: -tao_monitor_short -tao_type bnls -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 287 288 test: 289 suffix: 11 290 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbroyden 291 292 test: 293 suffix: 12 294 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbadbroyden 295 296 test: 297 suffix: 13 298 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbroyden 299 300 test: 301 suffix: 14 302 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbfgs 303 304 test: 305 suffix: 15 306 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmdfp 307 308 test: 309 suffix: 16 310 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsr1 311 312 test: 313 suffix: 17 314 args: -tao_monitor_short -tao_gatol 1e-4 -tao_type bqnls 315 316 test: 317 suffix: 18 318 args: -tao_monitor_short -tao_gatol 1e-4 -tao_type blmvm 319 320 test: 321 suffix: 19 322 args: -tao_monitor_short -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1 323 324 test: 325 suffix: 20 326 args: -tao_monitor -tao_gatol 1e-4 -tao_type blmvm -tao_ls_monitor 327 328 test: 329 suffix: 21 330 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbadbroyden 331 332 test: 333 suffix: 22 334 args: -tao_max_it 1 -tao_converged_reason 335 336 test: 337 suffix: 23 338 args: -tao_max_funcs 0 -tao_converged_reason 339 340 test: 341 suffix: 24 342 args: -tao_gatol 10 -tao_converged_reason 343 344 test: 345 suffix: 25 346 args: -tao_grtol 10 -tao_converged_reason 347 348 test: 349 suffix: 26 350 args: -tao_gttol 10 -tao_converged_reason 351 352 test: 353 suffix: 27 354 args: -tao_steptol 10 -tao_converged_reason 355 356 test: 357 suffix: 28 358 args: -tao_fmin 10 -tao_converged_reason 359 360 test: 361 suffix: snes 362 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 363 364 test: 365 suffix: snes_ls_armijo 366 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 367 368 test: 369 suffix: snes_tr_cgnegcurve_kmdc 370 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 371 requires: !single 372 373 test: 374 suffix: snes_ls_lmvm 375 args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtonls -snes_atol 1.e-4 -pc_type lmvm -tao_mf_hessian 376 377 TEST*/ 378