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