1! Program usage: mpiexec -n 1 chwirut1f [-help] [all TAO options] 2! 3! Description: This example demonstrates use of the TAO package to solve a 4! nonlinear least-squares problem on a single processor. We minimize the 5! Chwirut function: 6! sum_{i=0}^{n/2-1} ( alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2) 7! 8! The C version of this code is test_chwirut1.c 9! 10!!/*T 11! Concepts: TAO^Solving an unconstrained minimization problem 12! Routines: TaoCreate(); 13! Routines: TaoSetType(); 14! Routines: TaoSetInitialVector(); 15! Routines: TaoSetResidualRoutine(); 16! Routines: TaoSetFromOptions(); 17! Routines: TaoSolve(); 18! Routines: TaoDestroy(); 19! Processors: 1 20!T*/ 21 22 23! 24! ---------------------------------------------------------------------- 25! 26#include "chwirut1f.h" 27 28! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 29! Variable declarations 30! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 31! 32! See additional variable declarations in the file chwirut1f.h 33 34 PetscErrorCode ierr ! used to check for functions returning nonzeros 35 Vec x ! solution vector 36 Vec f ! vector of functions 37 Tao tao ! Tao context 38 PetscInt nhist 39 PetscMPIInt size,rank ! number of processes running 40 PetscReal zero 41 PetscReal hist(100) ! objective value history 42 PetscReal resid(100)! residual history 43 PetscReal cnorm(100)! cnorm history 44 PetscInt lits(100) ! #ksp history 45 PetscInt oh 46 TaoConvergedReason reason 47 48! Note: Any user-defined Fortran routines (such as FormGradient) 49! MUST be declared as external. 50 51 external FormFunction 52 53 zero = 0.0 54 55! Initialize TAO and PETSc 56 call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 57 if (ierr .ne. 0) then 58 print*,'Unable to initialize PETSc' 59 stop 60 endif 61 62 call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) 63 call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) 64 if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only '); endif 65 66! Initialize problem parameters 67 m = 214 68 n = 3 69 70 71! Allocate vectors for the solution and gradient 72 call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr) 73 call VecCreateSeq(PETSC_COMM_SELF,m,f,ierr) 74 75 76! The TAO code begins here 77 78! Create TAO solver 79 call TaoCreate(PETSC_COMM_SELF,tao,ierr);CHKERRA(ierr) 80 call TaoSetType(tao,TAOPOUNDERS,ierr);CHKERRA(ierr) 81! Set routines for function, gradient, and hessian evaluation 82 83 call TaoSetResidualRoutine(tao,f, & 84 & FormFunction,0,ierr) 85 CHKERRA(ierr) 86 87! Optional: Set initial guess 88 call InitializeData() 89 call FormStartingPoint(x) 90 call TaoSetInitialVector(tao, x, ierr) 91 CHKERRA(ierr) 92 93 94! Check for TAO command line options 95 call TaoSetFromOptions(tao,ierr) 96 CHKERRA(ierr) 97 oh = 100 98 call TaoSetConvergenceHistory(tao,hist,resid,cnorm,lits, & 99 & oh,PETSC_TRUE,ierr) 100 CHKERRA(ierr) 101! SOLVE THE APPLICATION 102 call TaoSolve(tao,ierr) 103 CHKERRA(ierr) 104 call TaoGetConvergenceHistory(tao,nhist,ierr) 105 CHKERRA(ierr) 106 call TaoGetConvergedReason(tao, reason, ierr) 107 if (reason .le. 0) then 108 print *,'Tao failed.' 109 print *,'Try a different TAO method, adjust some parameters,' 110 print *,'or check the function evaluation routines.' 111 endif 112 113! Free TAO data structures 114 call TaoDestroy(tao,ierr) 115 116! Free PETSc data structures 117 call VecDestroy(x,ierr) 118 call VecDestroy(f,ierr) 119 120 call PetscFinalize(ierr) 121 122 end 123 124 125! -------------------------------------------------------------------- 126! FormFunction - Evaluates the function f(X) and gradient G(X) 127! 128! Input Parameters: 129! tao - the Tao context 130! X - input vector 131! dummy - not used 132! 133! Output Parameters: 134! f - function vector 135 136 subroutine FormFunction(tao, x, f, dummy, ierr) 137#include "chwirut1f.h" 138 139 Tao tao 140 Vec x,f 141 PetscErrorCode ierr 142 PetscInt dummy 143 144 PetscInt i 145 PetscScalar, pointer, dimension(:) :: x_v,f_v 146 147 ierr = 0 148 149! Get pointers to vector data 150 call VecGetArrayF90(x,x_v,ierr);CHKERRQ(ierr) 151 call VecGetArrayF90(f,f_v,ierr);CHKERRQ(ierr) 152 153 154! Compute F(X) 155 do i=0,m-1 156 f_v(i+1) = y(i) - exp(-x_v(1)*t(i))/(x_v(2) + x_v(3)*t(i)) 157 enddo 158 159 160! Restore vectors 161 call VecRestoreArrayF90(X,x_v,ierr);CHKERRQ(ierr) 162 call VecRestoreArrayF90(F,f_v,ierr);CHKERRQ(ierr) 163 164 165 return 166 end 167 168 subroutine FormStartingPoint(x) 169#include "chwirut1f.h" 170 171 Vec x 172 PetscScalar, pointer, dimension(:) :: x_v 173 PetscErrorCode ierr 174 175 call VecGetArrayF90(x,x_v,ierr) 176 x_v(1) = 0.15 177 x_v(2) = 0.008 178 x_v(3) = 0.01 179 call VecRestoreArrayF90(x,x_v,ierr) 180 return 181 end 182 183 subroutine InitializeData() 184#include "chwirut1f.h" 185 186 integer i 187 i=0 188 y(i) = 92.9000; t(i) = 0.5000; i=i+1 189 y(i) = 78.7000; t(i) = 0.6250; i=i+1 190 y(i) = 64.2000; t(i) = 0.7500; i=i+1 191 y(i) = 64.9000; t(i) = 0.8750; i=i+1 192 y(i) = 57.1000; t(i) = 1.0000; i=i+1 193 y(i) = 43.3000; t(i) = 1.2500; i=i+1 194 y(i) = 31.1000; t(i) = 1.7500; i=i+1 195 y(i) = 23.6000; t(i) = 2.2500; i=i+1 196 y(i) = 31.0500; t(i) = 1.7500; i=i+1 197 y(i) = 23.7750; t(i) = 2.2500; i=i+1 198 y(i) = 17.7375; t(i) = 2.7500; i=i+1 199 y(i) = 13.8000; t(i) = 3.2500; i=i+1 200 y(i) = 11.5875; t(i) = 3.7500; i=i+1 201 y(i) = 9.4125; t(i) = 4.2500; i=i+1 202 y(i) = 7.7250; t(i) = 4.7500; i=i+1 203 y(i) = 7.3500; t(i) = 5.2500; i=i+1 204 y(i) = 8.0250; t(i) = 5.7500; i=i+1 205 y(i) = 90.6000; t(i) = 0.5000; i=i+1 206 y(i) = 76.9000; t(i) = 0.6250; i=i+1 207 y(i) = 71.6000; t(i) = 0.7500; i=i+1 208 y(i) = 63.6000; t(i) = 0.8750; i=i+1 209 y(i) = 54.0000; t(i) = 1.0000; i=i+1 210 y(i) = 39.2000; t(i) = 1.2500; i=i+1 211 y(i) = 29.3000; t(i) = 1.7500; i=i+1 212 y(i) = 21.4000; t(i) = 2.2500; i=i+1 213 y(i) = 29.1750; t(i) = 1.7500; i=i+1 214 y(i) = 22.1250; t(i) = 2.2500; i=i+1 215 y(i) = 17.5125; t(i) = 2.7500; i=i+1 216 y(i) = 14.2500; t(i) = 3.2500; i=i+1 217 y(i) = 9.4500; t(i) = 3.7500; i=i+1 218 y(i) = 9.1500; t(i) = 4.2500; i=i+1 219 y(i) = 7.9125; t(i) = 4.7500; i=i+1 220 y(i) = 8.4750; t(i) = 5.2500; i=i+1 221 y(i) = 6.1125; t(i) = 5.7500; i=i+1 222 y(i) = 80.0000; t(i) = 0.5000; i=i+1 223 y(i) = 79.0000; t(i) = 0.6250; i=i+1 224 y(i) = 63.8000; t(i) = 0.7500; i=i+1 225 y(i) = 57.2000; t(i) = 0.8750; i=i+1 226 y(i) = 53.2000; t(i) = 1.0000; i=i+1 227 y(i) = 42.5000; t(i) = 1.2500; i=i+1 228 y(i) = 26.8000; t(i) = 1.7500; i=i+1 229 y(i) = 20.4000; t(i) = 2.2500; i=i+1 230 y(i) = 26.8500; t(i) = 1.7500; i=i+1 231 y(i) = 21.0000; t(i) = 2.2500; i=i+1 232 y(i) = 16.4625; t(i) = 2.7500; i=i+1 233 y(i) = 12.5250; t(i) = 3.2500; i=i+1 234 y(i) = 10.5375; t(i) = 3.7500; i=i+1 235 y(i) = 8.5875; t(i) = 4.2500; i=i+1 236 y(i) = 7.1250; t(i) = 4.7500; i=i+1 237 y(i) = 6.1125; t(i) = 5.2500; i=i+1 238 y(i) = 5.9625; t(i) = 5.7500; i=i+1 239 y(i) = 74.1000; t(i) = 0.5000; i=i+1 240 y(i) = 67.3000; t(i) = 0.6250; i=i+1 241 y(i) = 60.8000; t(i) = 0.7500; i=i+1 242 y(i) = 55.5000; t(i) = 0.8750; i=i+1 243 y(i) = 50.3000; t(i) = 1.0000; i=i+1 244 y(i) = 41.0000; t(i) = 1.2500; i=i+1 245 y(i) = 29.4000; t(i) = 1.7500; i=i+1 246 y(i) = 20.4000; t(i) = 2.2500; i=i+1 247 y(i) = 29.3625; t(i) = 1.7500; i=i+1 248 y(i) = 21.1500; t(i) = 2.2500; i=i+1 249 y(i) = 16.7625; t(i) = 2.7500; i=i+1 250 y(i) = 13.2000; t(i) = 3.2500; i=i+1 251 y(i) = 10.8750; t(i) = 3.7500; i=i+1 252 y(i) = 8.1750; t(i) = 4.2500; i=i+1 253 y(i) = 7.3500; t(i) = 4.7500; i=i+1 254 y(i) = 5.9625; t(i) = 5.2500; i=i+1 255 y(i) = 5.6250; t(i) = 5.7500; i=i+1 256 y(i) = 81.5000; t(i) = .5000; i=i+1 257 y(i) = 62.4000; t(i) = .7500; i=i+1 258 y(i) = 32.5000; t(i) = 1.5000; i=i+1 259 y(i) = 12.4100; t(i) = 3.0000; i=i+1 260 y(i) = 13.1200; t(i) = 3.0000; i=i+1 261 y(i) = 15.5600; t(i) = 3.0000; i=i+1 262 y(i) = 5.6300; t(i) = 6.0000; i=i+1 263 y(i) = 78.0000; t(i) = .5000; i=i+1 264 y(i) = 59.9000; t(i) = .7500; i=i+1 265 y(i) = 33.2000; t(i) = 1.5000; i=i+1 266 y(i) = 13.8400; t(i) = 3.0000; i=i+1 267 y(i) = 12.7500; t(i) = 3.0000; i=i+1 268 y(i) = 14.6200; t(i) = 3.0000; i=i+1 269 y(i) = 3.9400; t(i) = 6.0000; i=i+1 270 y(i) = 76.8000; t(i) = .5000; i=i+1 271 y(i) = 61.0000; t(i) = .7500; i=i+1 272 y(i) = 32.9000; t(i) = 1.5000; i=i+1 273 y(i) = 13.8700; t(i) = 3.0000; i=i+1 274 y(i) = 11.8100; t(i) = 3.0000; i=i+1 275 y(i) = 13.3100; t(i) = 3.0000; i=i+1 276 y(i) = 5.4400; t(i) = 6.0000; i=i+1 277 y(i) = 78.0000; t(i) = .5000; i=i+1 278 y(i) = 63.5000; t(i) = .7500; i=i+1 279 y(i) = 33.8000; t(i) = 1.5000; i=i+1 280 y(i) = 12.5600; t(i) = 3.0000; i=i+1 281 y(i) = 5.6300; t(i) = 6.0000; i=i+1 282 y(i) = 12.7500; t(i) = 3.0000; i=i+1 283 y(i) = 13.1200; t(i) = 3.0000; i=i+1 284 y(i) = 5.4400; t(i) = 6.0000; i=i+1 285 y(i) = 76.8000; t(i) = .5000; i=i+1 286 y(i) = 60.0000; t(i) = .7500; i=i+1 287 y(i) = 47.8000; t(i) = 1.0000; i=i+1 288 y(i) = 32.0000; t(i) = 1.5000; i=i+1 289 y(i) = 22.2000; t(i) = 2.0000; i=i+1 290 y(i) = 22.5700; t(i) = 2.0000; i=i+1 291 y(i) = 18.8200; t(i) = 2.5000; i=i+1 292 y(i) = 13.9500; t(i) = 3.0000; i=i+1 293 y(i) = 11.2500; t(i) = 4.0000; i=i+1 294 y(i) = 9.0000; t(i) = 5.0000; i=i+1 295 y(i) = 6.6700; t(i) = 6.0000; i=i+1 296 y(i) = 75.8000; t(i) = .5000; i=i+1 297 y(i) = 62.0000; t(i) = .7500; i=i+1 298 y(i) = 48.8000; t(i) = 1.0000; i=i+1 299 y(i) = 35.2000; t(i) = 1.5000; i=i+1 300 y(i) = 20.0000; t(i) = 2.0000; i=i+1 301 y(i) = 20.3200; t(i) = 2.0000; i=i+1 302 y(i) = 19.3100; t(i) = 2.5000; i=i+1 303 y(i) = 12.7500; t(i) = 3.0000; i=i+1 304 y(i) = 10.4200; t(i) = 4.0000; i=i+1 305 y(i) = 7.3100; t(i) = 5.0000; i=i+1 306 y(i) = 7.4200; t(i) = 6.0000; i=i+1 307 y(i) = 70.5000; t(i) = .5000; i=i+1 308 y(i) = 59.5000; t(i) = .7500; i=i+1 309 y(i) = 48.5000; t(i) = 1.0000; i=i+1 310 y(i) = 35.8000; t(i) = 1.5000; i=i+1 311 y(i) = 21.0000; t(i) = 2.0000; i=i+1 312 y(i) = 21.6700; t(i) = 2.0000; i=i+1 313 y(i) = 21.0000; t(i) = 2.5000; i=i+1 314 y(i) = 15.6400; t(i) = 3.0000; i=i+1 315 y(i) = 8.1700; t(i) = 4.0000; i=i+1 316 y(i) = 8.5500; t(i) = 5.0000; i=i+1 317 y(i) = 10.1200; t(i) = 6.0000; i=i+1 318 y(i) = 78.0000; t(i) = .5000; i=i+1 319 y(i) = 66.0000; t(i) = .6250; i=i+1 320 y(i) = 62.0000; t(i) = .7500; i=i+1 321 y(i) = 58.0000; t(i) = .8750; i=i+1 322 y(i) = 47.7000; t(i) = 1.0000; i=i+1 323 y(i) = 37.8000; t(i) = 1.2500; i=i+1 324 y(i) = 20.2000; t(i) = 2.2500; i=i+1 325 y(i) = 21.0700; t(i) = 2.2500; i=i+1 326 y(i) = 13.8700; t(i) = 2.7500; i=i+1 327 y(i) = 9.6700; t(i) = 3.2500; i=i+1 328 y(i) = 7.7600; t(i) = 3.7500; i=i+1 329 y(i) = 5.4400; t(i) = 4.2500; i=i+1 330 y(i) = 4.8700; t(i) = 4.7500; i=i+1 331 y(i) = 4.0100; t(i) = 5.2500; i=i+1 332 y(i) = 3.7500; t(i) = 5.7500; i=i+1 333 y(i) = 24.1900; t(i) = 3.0000; i=i+1 334 y(i) = 25.7600; t(i) = 3.0000; i=i+1 335 y(i) = 18.0700; t(i) = 3.0000; i=i+1 336 y(i) = 11.8100; t(i) = 3.0000; i=i+1 337 y(i) = 12.0700; t(i) = 3.0000; i=i+1 338 y(i) = 16.1200; t(i) = 3.0000; i=i+1 339 y(i) = 70.8000; t(i) = .5000; i=i+1 340 y(i) = 54.7000; t(i) = .7500; i=i+1 341 y(i) = 48.0000; t(i) = 1.0000; i=i+1 342 y(i) = 39.8000; t(i) = 1.5000; i=i+1 343 y(i) = 29.8000; t(i) = 2.0000; i=i+1 344 y(i) = 23.7000; t(i) = 2.5000; i=i+1 345 y(i) = 29.6200; t(i) = 2.0000; i=i+1 346 y(i) = 23.8100; t(i) = 2.5000; i=i+1 347 y(i) = 17.7000; t(i) = 3.0000; i=i+1 348 y(i) = 11.5500; t(i) = 4.0000; i=i+1 349 y(i) = 12.0700; t(i) = 5.0000; i=i+1 350 y(i) = 8.7400; t(i) = 6.0000; i=i+1 351 y(i) = 80.7000; t(i) = .5000; i=i+1 352 y(i) = 61.3000; t(i) = .7500; i=i+1 353 y(i) = 47.5000; t(i) = 1.0000; i=i+1 354 y(i) = 29.0000; t(i) = 1.5000; i=i+1 355 y(i) = 24.0000; t(i) = 2.0000; i=i+1 356 y(i) = 17.7000; t(i) = 2.5000; i=i+1 357 y(i) = 24.5600; t(i) = 2.0000; i=i+1 358 y(i) = 18.6700; t(i) = 2.5000; i=i+1 359 y(i) = 16.2400; t(i) = 3.0000; i=i+1 360 y(i) = 8.7400; t(i) = 4.0000; i=i+1 361 y(i) = 7.8700; t(i) = 5.0000; i=i+1 362 y(i) = 8.5100; t(i) = 6.0000; i=i+1 363 y(i) = 66.7000; t(i) = .5000; i=i+1 364 y(i) = 59.2000; t(i) = .7500; i=i+1 365 y(i) = 40.8000; t(i) = 1.0000; i=i+1 366 y(i) = 30.7000; t(i) = 1.5000; i=i+1 367 y(i) = 25.7000; t(i) = 2.0000; i=i+1 368 y(i) = 16.3000; t(i) = 2.5000; i=i+1 369 y(i) = 25.9900; t(i) = 2.0000; i=i+1 370 y(i) = 16.9500; t(i) = 2.5000; i=i+1 371 y(i) = 13.3500; t(i) = 3.0000; i=i+1 372 y(i) = 8.6200; t(i) = 4.0000; i=i+1 373 y(i) = 7.2000; t(i) = 5.0000; i=i+1 374 y(i) = 6.6400; t(i) = 6.0000; i=i+1 375 y(i) = 13.6900; t(i) = 3.0000; i=i+1 376 y(i) = 81.0000; t(i) = .5000; i=i+1 377 y(i) = 64.5000; t(i) = .7500; i=i+1 378 y(i) = 35.5000; t(i) = 1.5000; i=i+1 379 y(i) = 13.3100; t(i) = 3.0000; i=i+1 380 y(i) = 4.8700; t(i) = 6.0000; i=i+1 381 y(i) = 12.9400; t(i) = 3.0000; i=i+1 382 y(i) = 5.0600; t(i) = 6.0000; i=i+1 383 y(i) = 15.1900; t(i) = 3.0000; i=i+1 384 y(i) = 14.6200; t(i) = 3.0000; i=i+1 385 y(i) = 15.6400; t(i) = 3.0000; i=i+1 386 y(i) = 25.5000; t(i) = 1.7500; i=i+1 387 y(i) = 25.9500; t(i) = 1.7500; i=i+1 388 y(i) = 81.7000; t(i) = .5000; i=i+1 389 y(i) = 61.6000; t(i) = .7500; i=i+1 390 y(i) = 29.8000; t(i) = 1.7500; i=i+1 391 y(i) = 29.8100; t(i) = 1.7500; i=i+1 392 y(i) = 17.1700; t(i) = 2.7500; i=i+1 393 y(i) = 10.3900; t(i) = 3.7500; i=i+1 394 y(i) = 28.4000; t(i) = 1.7500; i=i+1 395 y(i) = 28.6900; t(i) = 1.7500; i=i+1 396 y(i) = 81.3000; t(i) = .5000; i=i+1 397 y(i) = 60.9000; t(i) = .7500; i=i+1 398 y(i) = 16.6500; t(i) = 2.7500; i=i+1 399 y(i) = 10.0500; t(i) = 3.7500; i=i+1 400 y(i) = 28.9000; t(i) = 1.7500; i=i+1 401 y(i) = 28.9500; t(i) = 1.7500; i=i+1 402 403 return 404 end 405 406!/*TEST 407! 408! build: 409! requires: !complex 410! 411! test: 412! args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-5 413! requires: !single 414! 415!TEST*/ 416