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