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