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