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