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