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