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