1 /* Program usage: mpiexec -n 1 rosenbrock2 [-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 TaoConvergedReason reason; 38 PetscInt its, recycled_its = 0, oneshot_its = 0; 39 40 /* Initialize TAO and PETSc */ 41 PetscFunctionBeginUser; 42 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 43 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 44 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors"); 45 46 /* Initialize problem parameters */ 47 user.n = 2; 48 user.alpha = 99.0; 49 user.chained = PETSC_FALSE; 50 /* Check for command line arguments to override defaults */ 51 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &user.n, &flg)); 52 PetscCall(PetscOptionsGetReal(NULL, NULL, "-alpha", &user.alpha, &flg)); 53 PetscCall(PetscOptionsGetBool(NULL, NULL, "-chained", &user.chained, &flg)); 54 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg)); 55 56 /* Allocate vectors for the solution and gradient */ 57 PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.n, &x)); 58 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, 2, user.n, user.n, 1, NULL, &H)); 59 60 /* The TAO code begins here */ 61 62 /* Create TAO solver with desired solution method */ 63 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao)); 64 PetscCall(TaoSetType(tao, TAOLMVM)); 65 66 /* Set solution vec and an initial guess */ 67 PetscCall(VecSet(x, zero)); 68 PetscCall(TaoSetSolution(tao, x)); 69 70 /* Set routines for function, gradient, hessian evaluation */ 71 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user)); 72 PetscCall(TaoSetHessian(tao, H, H, FormHessian, &user)); 73 74 /* Check for TAO command line options */ 75 PetscCall(TaoSetFromOptions(tao)); 76 77 /* Solve the problem */ 78 PetscCall(TaoSetTolerances(tao, 1.e-5, 0.0, 0.0)); 79 PetscCall(TaoSetMaximumIterations(tao, 5)); 80 PetscCall(TaoLMVMRecycle(tao, PETSC_TRUE)); 81 reason = TAO_CONTINUE_ITERATING; 82 while (reason != TAO_CONVERGED_GATOL) { 83 PetscCall(TaoSolve(tao)); 84 PetscCall(TaoGetConvergedReason(tao, &reason)); 85 PetscCall(TaoGetIterationNumber(tao, &its)); 86 recycled_its += its; 87 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n")); 88 } 89 90 /* Disable recycling and solve again! */ 91 PetscCall(TaoSetMaximumIterations(tao, 100)); 92 PetscCall(TaoLMVMRecycle(tao, PETSC_FALSE)); 93 PetscCall(VecSet(x, zero)); 94 PetscCall(TaoSolve(tao)); 95 PetscCall(TaoGetConvergedReason(tao, &reason)); 96 PetscCheck(reason == TAO_CONVERGED_GATOL, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Solution failed to converge!"); 97 PetscCall(TaoGetIterationNumber(tao, &oneshot_its)); 98 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n")); 99 PetscCall(PetscPrintf(PETSC_COMM_SELF, "recycled its: %" PetscInt_FMT " | oneshot its: %" PetscInt_FMT "\n", recycled_its, oneshot_its)); 100 PetscCheck(recycled_its == oneshot_its, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "LMVM recycling does not work!"); 101 102 PetscCall(TaoDestroy(&tao)); 103 PetscCall(VecDestroy(&x)); 104 PetscCall(MatDestroy(&H)); 105 106 PetscCall(PetscFinalize()); 107 return 0; 108 } 109 110 /* -------------------------------------------------------------------- */ 111 /* 112 FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X). 113 114 Input Parameters: 115 . tao - the Tao context 116 . X - input vector 117 . ptr - optional user-defined context, as set by TaoSetFunctionGradient() 118 119 Output Parameters: 120 . G - vector containing the newly evaluated gradient 121 . f - function value 122 123 Note: 124 Some optimization methods ask for the function and the gradient evaluation 125 at the same time. Evaluating both at once may be more efficient that 126 evaluating each separately. 127 */ 128 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr) 129 { 130 AppCtx *user = (AppCtx *)ptr; 131 PetscInt i, nn = user->n / 2; 132 PetscReal ff = 0, t1, t2, alpha = user->alpha; 133 PetscScalar *g; 134 const PetscScalar *x; 135 136 PetscFunctionBeginUser; 137 /* Get pointers to vector data */ 138 PetscCall(VecGetArrayRead(X, &x)); 139 PetscCall(VecGetArray(G, &g)); 140 141 /* Compute G(X) */ 142 if (user->chained) { 143 g[0] = 0; 144 for (i = 0; i < user->n - 1; i++) { 145 t1 = x[i + 1] - x[i] * x[i]; 146 ff += PetscSqr(1 - x[i]) + alpha * t1 * t1; 147 g[i] += -2 * (1 - x[i]) + 2 * alpha * t1 * (-2 * x[i]); 148 g[i + 1] = 2 * alpha * t1; 149 } 150 } else { 151 for (i = 0; i < nn; i++) { 152 t1 = x[2 * i + 1] - x[2 * i] * x[2 * i]; 153 t2 = 1 - x[2 * i]; 154 ff += alpha * t1 * t1 + t2 * t2; 155 g[2 * i] = -4 * alpha * t1 * x[2 * i] - 2.0 * t2; 156 g[2 * i + 1] = 2 * alpha * t1; 157 } 158 } 159 160 /* Restore vectors */ 161 PetscCall(VecRestoreArrayRead(X, &x)); 162 PetscCall(VecRestoreArray(G, &g)); 163 *f = ff; 164 165 PetscCall(PetscLogFlops(15.0 * nn)); 166 PetscFunctionReturn(PETSC_SUCCESS); 167 } 168 169 /* ------------------------------------------------------------------- */ 170 /* 171 FormHessian - Evaluates Hessian matrix. 172 173 Input Parameters: 174 . tao - the Tao context 175 . x - input vector 176 . ptr - optional user-defined context, as set by TaoSetHessian() 177 178 Output Parameters: 179 . H - Hessian matrix 180 181 Note: Providing the Hessian may not be necessary. Only some solvers 182 require this matrix. 183 */ 184 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 185 { 186 AppCtx *user = (AppCtx *)ptr; 187 PetscInt i, ind[2]; 188 PetscReal alpha = user->alpha; 189 PetscReal v[2][2]; 190 const PetscScalar *x; 191 PetscBool assembled; 192 193 PetscFunctionBeginUser; 194 /* Zero existing matrix entries */ 195 PetscCall(MatAssembled(H, &assembled)); 196 if (assembled) PetscCall(MatZeroEntries(H)); 197 198 /* Get a pointer to vector data */ 199 PetscCall(VecGetArrayRead(X, &x)); 200 201 /* Compute H(X) entries */ 202 if (user->chained) { 203 PetscCall(MatZeroEntries(H)); 204 for (i = 0; i < user->n - 1; i++) { 205 PetscScalar t1 = x[i + 1] - x[i] * x[i]; 206 v[0][0] = 2 + 2 * alpha * (t1 * (-2) - 2 * x[i]); 207 v[0][1] = 2 * alpha * (-2 * x[i]); 208 v[1][0] = 2 * alpha * (-2 * x[i]); 209 v[1][1] = 2 * alpha * t1; 210 ind[0] = i; 211 ind[1] = i + 1; 212 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], ADD_VALUES)); 213 } 214 } else { 215 for (i = 0; i < user->n / 2; i++) { 216 v[1][1] = 2 * alpha; 217 v[0][0] = -4 * alpha * (x[2 * i + 1] - 3 * x[2 * i] * x[2 * i]) + 2; 218 v[1][0] = v[0][1] = -4.0 * alpha * x[2 * i]; 219 ind[0] = 2 * i; 220 ind[1] = 2 * i + 1; 221 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], INSERT_VALUES)); 222 } 223 } 224 PetscCall(VecRestoreArrayRead(X, &x)); 225 226 /* Assemble matrix */ 227 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 228 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 229 PetscCall(PetscLogFlops(9.0 * user->n / 2.0)); 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 /*TEST 234 235 build: 236 requires: !complex 237 238 test: 239 args: -tao_type lmvm -tao_monitor 240 requires: !single 241 242 TEST*/ 243