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