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; user.alpha = 99.0; user.chained = PETSC_FALSE; 51 /* Check for command line arguments to override defaults */ 52 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&user.n,&flg)); 53 PetscCall(PetscOptionsGetReal(NULL,NULL,"-alpha",&user.alpha,&flg)); 54 PetscCall(PetscOptionsGetBool(NULL,NULL,"-chained",&user.chained,&flg)); 55 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg)); 56 57 /* Allocate vectors for the solution and gradient */ 58 PetscCall(VecCreateSeq(PETSC_COMM_SELF,user.n,&x)); 59 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF,2,user.n,user.n,1,NULL,&H)); 60 61 /* The TAO code begins here */ 62 63 /* Create TAO solver with desired solution method */ 64 PetscCall(TaoCreate(PETSC_COMM_SELF,&tao)); 65 PetscCall(TaoSetType(tao,TAOLMVM)); 66 67 /* Set solution vec and an initial guess */ 68 PetscCall(VecSet(x, zero)); 69 PetscCall(TaoSetSolution(tao,x)); 70 71 /* Set routines for function, gradient, hessian evaluation */ 72 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,&user)); 73 PetscCall(TaoSetHessian(tao,H,H,FormHessian,&user)); 74 75 /* Test the LMVM matrix */ 76 if (test_lmvm) { 77 PetscCall(PetscOptionsSetValue(NULL, "-tao_type", "bqnktr")); 78 } 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]; t2= 1-x[2*i]; 162 ff += alpha*t1*t1 + t2*t2; 163 g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2; 164 g[2*i+1] = 2*alpha*t1; 165 } 166 } 167 168 /* Restore vectors */ 169 PetscCall(VecRestoreArrayRead(X,&x)); 170 PetscCall(VecRestoreArray(G,&g)); 171 *f = ff; 172 173 PetscCall(PetscLogFlops(15.0*nn)); 174 PetscFunctionReturn(0); 175 } 176 177 /* ------------------------------------------------------------------- */ 178 /* 179 FormHessian - Evaluates Hessian matrix. 180 181 Input Parameters: 182 . tao - the Tao context 183 . x - input vector 184 . ptr - optional user-defined context, as set by TaoSetHessian() 185 186 Output Parameters: 187 . H - Hessian matrix 188 189 Note: Providing the Hessian may not be necessary. Only some solvers 190 require this matrix. 191 */ 192 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 193 { 194 AppCtx *user = (AppCtx*)ptr; 195 PetscInt i, ind[2]; 196 PetscReal alpha=user->alpha; 197 PetscReal v[2][2]; 198 const PetscScalar *x; 199 PetscBool assembled; 200 201 PetscFunctionBeginUser; 202 /* Zero existing matrix entries */ 203 PetscCall(MatAssembled(H,&assembled)); 204 if (assembled) PetscCall(MatZeroEntries(H)); 205 206 /* Get a pointer to vector data */ 207 PetscCall(VecGetArrayRead(X,&x)); 208 209 /* Compute H(X) entries */ 210 if (user->chained) { 211 PetscCall(MatZeroEntries(H)); 212 for (i=0; i<user->n-1; i++) { 213 PetscScalar t1 = x[i+1] - x[i]*x[i]; 214 v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]); 215 v[0][1] = 2*alpha*(-2*x[i]); 216 v[1][0] = 2*alpha*(-2*x[i]); 217 v[1][1] = 2*alpha*t1; 218 ind[0] = i; ind[1] = i+1; 219 PetscCall(MatSetValues(H,2,ind,2,ind,v[0],ADD_VALUES)); 220 } 221 } else { 222 for (i=0; i<user->n/2; i++) { 223 v[1][1] = 2*alpha; 224 v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2; 225 v[1][0] = v[0][1] = -4.0*alpha*x[2*i]; 226 ind[0]=2*i; 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