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 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; user.alpha = 99.0; 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 /* Test the LMVM matrix */ 75 if (test_lmvm) { 76 PetscCall(PetscOptionsSetValue(NULL, "-tao_type", "bqnktr")); 77 } 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 { 138 AppCtx *user = (AppCtx *) ptr; 139 PetscInt i,nn=user->n/2; 140 PetscReal ff=0,t1,t2,alpha=user->alpha; 141 PetscScalar *g; 142 const PetscScalar *x; 143 144 PetscFunctionBeginUser; 145 /* Get pointers to vector data */ 146 PetscCall(VecGetArrayRead(X,&x)); 147 PetscCall(VecGetArray(G,&g)); 148 149 /* Compute G(X) */ 150 if (user->chained) { 151 g[0] = 0; 152 for (i=0; i<user->n-1; i++) { 153 t1 = x[i+1] - x[i]*x[i]; 154 ff += PetscSqr(1 - x[i]) + alpha*t1*t1; 155 g[i] += -2*(1 - x[i]) + 2*alpha*t1*(-2*x[i]); 156 g[i+1] = 2*alpha*t1; 157 } 158 } else { 159 for (i=0; i<nn; i++) { 160 t1 = x[2*i+1]-x[2*i]*x[2*i]; 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 { 193 AppCtx *user = (AppCtx*)ptr; 194 PetscInt i, ind[2]; 195 PetscReal alpha=user->alpha; 196 PetscReal v[2][2]; 197 const PetscScalar *x; 198 PetscBool assembled; 199 200 PetscFunctionBeginUser; 201 /* Zero existing matrix entries */ 202 PetscCall(MatAssembled(H,&assembled)); 203 if (assembled) PetscCall(MatZeroEntries(H)); 204 205 /* Get a pointer to vector data */ 206 PetscCall(VecGetArrayRead(X,&x)); 207 208 /* Compute H(X) entries */ 209 if (user->chained) { 210 PetscCall(MatZeroEntries(H)); 211 for (i=0; i<user->n-1; i++) { 212 PetscScalar t1 = x[i+1] - x[i]*x[i]; 213 v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]); 214 v[0][1] = 2*alpha*(-2*x[i]); 215 v[1][0] = 2*alpha*(-2*x[i]); 216 v[1][1] = 2*alpha*t1; 217 ind[0] = i; 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; ind[1]=2*i+1; 226 PetscCall(MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES)); 227 } 228 } 229 PetscCall(VecRestoreArrayRead(X,&x)); 230 231 /* Assemble matrix */ 232 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 233 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 234 PetscCall(PetscLogFlops(9.0*user->n/2.0)); 235 PetscFunctionReturn(0); 236 } 237 238 /*TEST 239 240 build: 241 requires: !complex 242 243 test: 244 args: -tao_smonitor -tao_type nls -tao_gatol 1.e-4 245 requires: !single 246 247 test: 248 suffix: 2 249 args: -tao_smonitor -tao_type lmvm -tao_gatol 1.e-3 250 251 test: 252 suffix: 3 253 args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4 254 requires: !single 255 256 test: 257 suffix: 4 258 args: -tao_smonitor -tao_type ntr -tao_mf_hessian -tao_ntr_pc_type none -tao_gatol 1.e-4 259 260 test: 261 suffix: 5 262 args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 263 264 test: 265 suffix: 6 266 args: -tao_smonitor -tao_type bntl -tao_gatol 1.e-4 267 268 test: 269 suffix: 7 270 args: -tao_smonitor -tao_type bnls -tao_gatol 1.e-4 271 272 test: 273 suffix: 8 274 args: -tao_smonitor -tao_type bntr -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 275 276 test: 277 suffix: 9 278 args: -tao_smonitor -tao_type bntl -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 279 280 test: 281 suffix: 10 282 args: -tao_smonitor -tao_type bnls -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 283 284 test: 285 suffix: 11 286 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbroyden 287 288 test: 289 suffix: 12 290 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbadbroyden 291 292 test: 293 suffix: 13 294 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbroyden 295 296 test: 297 suffix: 14 298 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbfgs 299 300 test: 301 suffix: 15 302 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmdfp 303 304 test: 305 suffix: 16 306 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsr1 307 308 test: 309 suffix: 17 310 args: -tao_smonitor -tao_gatol 1e-4 -tao_type bqnls 311 312 test: 313 suffix: 18 314 args: -tao_smonitor -tao_gatol 1e-4 -tao_type blmvm 315 316 test: 317 suffix: 19 318 args: -tao_smonitor -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1 319 320 test: 321 suffix: 20 322 args: -tao_monitor -tao_gatol 1e-4 -tao_type blmvm -tao_ls_monitor 323 324 test: 325 suffix: 21 326 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbadbroyden 327 328 test: 329 suffix: 22 330 args: -tao_max_it 1 -tao_converged_reason 331 332 test: 333 suffix: 23 334 args: -tao_max_funcs 0 -tao_converged_reason 335 336 test: 337 suffix: 24 338 args: -tao_gatol 10 -tao_converged_reason 339 340 test: 341 suffix: 25 342 args: -tao_grtol 10 -tao_converged_reason 343 344 test: 345 suffix: 26 346 args: -tao_gttol 10 -tao_converged_reason 347 348 test: 349 suffix: 27 350 args: -tao_steptol 10 -tao_converged_reason 351 352 test: 353 suffix: 28 354 args: -tao_fmin 10 -tao_converged_reason 355 356 TEST*/ 357