1 /* 2 Include "petsctao.h" so that we can use TAO solvers. Note that this 3 file automatically includes libraries such as: 4 petsc.h - base PETSc routines petscvec.h - vectors 5 petscsys.h - sysem routines petscmat.h - matrices 6 petscis.h - index sets petscksp.h - Krylov subspace methods 7 petscviewer.h - viewers petscpc.h - preconditioners 8 9 This version tests correlated terms using both vector and listed forms 10 */ 11 12 #include <petsctao.h> 13 14 /* 15 Description: These data are the result of a NIST study involving 16 ultrasonic calibration. The response variable is 17 ultrasonic response, and the predictor variable is 18 metal distance. 19 20 Reference: Chwirut, D., NIST (197?). 21 Ultrasonic Reference Block Study. 22 */ 23 24 static char help[]="Finds the nonlinear least-squares solution to the model \n\ 25 y = exp[-b1*x]/(b2+b3*x) + e \n"; 26 27 /*T 28 Concepts: TAO^Solving a system of nonlinear equations, nonlinear least squares 29 Routines: TaoCreate(); 30 Routines: TaoSetType(); 31 Routines: TaoSetResidualRoutine(); 32 Routines: TaoSetJacobianRoutine(); 33 Routines: TaoSetInitialVector(); 34 Routines: TaoSetFromOptions(); 35 Routines: TaoSetConvergenceHistory(); TaoGetConvergenceHistory(); 36 Routines: TaoSolve(); 37 Routines: TaoView(); TaoDestroy(); 38 Processors: 1 39 T*/ 40 41 #define NOBSERVATIONS 214 42 #define NPARAMETERS 3 43 44 /* User-defined application context */ 45 typedef struct { 46 /* Working space */ 47 PetscReal t[NOBSERVATIONS]; /* array of independent variables of observation */ 48 PetscReal y[NOBSERVATIONS]; /* array of dependent variables */ 49 PetscReal j[NOBSERVATIONS][NPARAMETERS]; /* dense jacobian matrix array*/ 50 PetscInt idm[NOBSERVATIONS]; /* Matrix indices for jacobian */ 51 PetscInt idn[NPARAMETERS]; 52 } AppCtx; 53 54 /* User provided Routines */ 55 PetscErrorCode InitializeData(AppCtx *user); 56 PetscErrorCode FormStartingPoint(Vec); 57 PetscErrorCode EvaluateFunction(Tao, Vec, Vec, void *); 58 PetscErrorCode EvaluateJacobian(Tao, Vec, Mat, Mat, void *); 59 60 /*--------------------------------------------------------------------*/ 61 int main(int argc,char **argv) 62 { 63 PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 64 PetscInt wtype=0; 65 Vec x, f; /* solution, function */ 66 Vec w; /* weights */ 67 Mat J; /* Jacobian matrix */ 68 Tao tao; /* Tao solver context */ 69 PetscInt i; /* iteration information */ 70 PetscReal hist[100],resid[100]; 71 PetscInt lits[100]; 72 PetscInt w_row[NOBSERVATIONS]; /* explicit weights */ 73 PetscInt w_col[NOBSERVATIONS]; 74 PetscReal w_vals[NOBSERVATIONS]; 75 PetscBool flg; 76 AppCtx user; /* user-defined work context */ 77 78 ierr = PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr; 79 ierr = PetscOptionsGetInt(NULL,NULL,"-wtype",&wtype,&flg);CHKERRQ(ierr); 80 ierr = PetscPrintf(PETSC_COMM_WORLD,"wtype=%d\n",wtype);CHKERRQ(ierr); 81 /* Allocate vectors */ 82 ierr = VecCreateSeq(MPI_COMM_SELF,NPARAMETERS,&x);CHKERRQ(ierr); 83 ierr = VecCreateSeq(MPI_COMM_SELF,NOBSERVATIONS,&f);CHKERRQ(ierr); 84 85 ierr = VecDuplicate(f,&w);CHKERRQ(ierr); 86 87 /* no correlation, but set in different ways */ 88 ierr = VecSet(w,1.0);CHKERRQ(ierr); 89 for (i=0;i<NOBSERVATIONS;i++) { 90 w_row[i]=i; w_col[i]=i; w_vals[i]=1.0; 91 } 92 93 /* Create the Jacobian matrix. */ 94 ierr = MatCreateSeqDense(MPI_COMM_SELF,NOBSERVATIONS,NPARAMETERS,NULL,&J);CHKERRQ(ierr); 95 96 for (i=0;i<NOBSERVATIONS;i++) user.idm[i] = i; 97 98 for (i=0;i<NPARAMETERS;i++) user.idn[i] = i; 99 100 /* Create TAO solver and set desired solution method */ 101 ierr = TaoCreate(PETSC_COMM_SELF,&tao);CHKERRQ(ierr); 102 ierr = TaoSetType(tao,TAOPOUNDERS);CHKERRQ(ierr); 103 104 /* Set the function and Jacobian routines. */ 105 ierr = InitializeData(&user);CHKERRQ(ierr); 106 ierr = FormStartingPoint(x);CHKERRQ(ierr); 107 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr); 108 ierr = TaoSetResidualRoutine(tao,f,EvaluateFunction,(void*)&user);CHKERRQ(ierr); 109 if (wtype == 1) { 110 ierr = TaoSetResidualWeights(tao,w,0,NULL,NULL,NULL);CHKERRQ(ierr); 111 } else if (wtype == 2) { 112 ierr = TaoSetResidualWeights(tao,NULL,NOBSERVATIONS,w_row,w_col,w_vals);CHKERRQ(ierr); 113 } 114 ierr = TaoSetJacobianResidualRoutine(tao, J, J, EvaluateJacobian, (void*)&user);CHKERRQ(ierr); 115 ierr = TaoSetTolerances(tao,1e-5,0.0,PETSC_DEFAULT);CHKERRQ(ierr); 116 117 /* Check for any TAO command line arguments */ 118 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 119 120 ierr = TaoSetConvergenceHistory(tao,hist,resid,0,lits,100,PETSC_TRUE);CHKERRQ(ierr); 121 /* Perform the Solve */ 122 ierr = TaoSolve(tao);CHKERRQ(ierr); 123 124 /* Free TAO data structures */ 125 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 126 127 /* Free PETSc data structures */ 128 ierr = VecDestroy(&x);CHKERRQ(ierr); 129 ierr = VecDestroy(&w);CHKERRQ(ierr); 130 ierr = VecDestroy(&f);CHKERRQ(ierr); 131 ierr = MatDestroy(&J);CHKERRQ(ierr); 132 133 ierr = PetscFinalize(); 134 return ierr; 135 } 136 137 /*--------------------------------------------------------------------*/ 138 PetscErrorCode EvaluateFunction(Tao tao, Vec X, Vec F, void *ptr) 139 { 140 AppCtx *user = (AppCtx *)ptr; 141 PetscInt i; 142 PetscReal *y=user->y,*f,*t=user->t; 143 const PetscReal *x; 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 148 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 149 150 for (i=0;i<NOBSERVATIONS;i++) { 151 f[i] = y[i] - PetscExpScalar(-x[0]*t[i])/(x[1] + x[2]*t[i]); 152 } 153 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 154 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 155 PetscLogFlops(6*NOBSERVATIONS); 156 PetscFunctionReturn(0); 157 } 158 159 /*------------------------------------------------------------*/ 160 /* J[i][j] = df[i]/dt[j] */ 161 PetscErrorCode EvaluateJacobian(Tao tao, Vec X, Mat J, Mat Jpre, void *ptr) 162 { 163 AppCtx *user = (AppCtx *)ptr; 164 PetscInt i; 165 PetscReal *t=user->t; 166 const PetscReal *x; 167 PetscReal base; 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 172 for (i=0;i<NOBSERVATIONS;i++) { 173 base = PetscExpScalar(-x[0]*t[i])/(x[1] + x[2]*t[i]); 174 175 user->j[i][0] = t[i]*base; 176 user->j[i][1] = base/(x[1] + x[2]*t[i]); 177 user->j[i][2] = base*t[i]/(x[1] + x[2]*t[i]); 178 } 179 180 /* Assemble the matrix */ 181 ierr = MatSetValues(J,NOBSERVATIONS,user->idm, NPARAMETERS, user->idn,(PetscReal *)user->j,INSERT_VALUES);CHKERRQ(ierr); 182 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 183 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 184 185 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 186 PetscLogFlops(NOBSERVATIONS * 13); 187 PetscFunctionReturn(0); 188 } 189 190 /* ------------------------------------------------------------ */ 191 PetscErrorCode FormStartingPoint(Vec X) 192 { 193 PetscReal *x; 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 198 x[0] = 1.19; 199 x[1] = -1.86; 200 x[2] = 1.08; 201 VecRestoreArray(X,&x);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204 205 /* ---------------------------------------------------------------------- */ 206 PetscErrorCode InitializeData(AppCtx *user) 207 { 208 PetscReal *t=user->t,*y=user->y; 209 PetscInt i=0; 210 211 PetscFunctionBegin; 212 y[i] = 92.9000; t[i++] = 0.5000; 213 y[i] = 78.7000; t[i++] = 0.6250; 214 y[i] = 64.2000; t[i++] = 0.7500; 215 y[i] = 64.9000; t[i++] = 0.8750; 216 y[i] = 57.1000; t[i++] = 1.0000; 217 y[i] = 43.3000; t[i++] = 1.2500; 218 y[i] = 31.1000; t[i++] = 1.7500; 219 y[i] = 23.6000; t[i++] = 2.2500; 220 y[i] = 31.0500; t[i++] = 1.7500; 221 y[i] = 23.7750; t[i++] = 2.2500; 222 y[i] = 17.7375; t[i++] = 2.7500; 223 y[i] = 13.8000; t[i++] = 3.2500; 224 y[i] = 11.5875; t[i++] = 3.7500; 225 y[i] = 9.4125; t[i++] = 4.2500; 226 y[i] = 7.7250; t[i++] = 4.7500; 227 y[i] = 7.3500; t[i++] = 5.2500; 228 y[i] = 8.0250; t[i++] = 5.7500; 229 y[i] = 90.6000; t[i++] = 0.5000; 230 y[i] = 76.9000; t[i++] = 0.6250; 231 y[i] = 71.6000; t[i++] = 0.7500; 232 y[i] = 63.6000; t[i++] = 0.8750; 233 y[i] = 54.0000; t[i++] = 1.0000; 234 y[i] = 39.2000; t[i++] = 1.2500; 235 y[i] = 29.3000; t[i++] = 1.7500; 236 y[i] = 21.4000; t[i++] = 2.2500; 237 y[i] = 29.1750; t[i++] = 1.7500; 238 y[i] = 22.1250; t[i++] = 2.2500; 239 y[i] = 17.5125; t[i++] = 2.7500; 240 y[i] = 14.2500; t[i++] = 3.2500; 241 y[i] = 9.4500; t[i++] = 3.7500; 242 y[i] = 9.1500; t[i++] = 4.2500; 243 y[i] = 7.9125; t[i++] = 4.7500; 244 y[i] = 8.4750; t[i++] = 5.2500; 245 y[i] = 6.1125; t[i++] = 5.7500; 246 y[i] = 80.0000; t[i++] = 0.5000; 247 y[i] = 79.0000; t[i++] = 0.6250; 248 y[i] = 63.8000; t[i++] = 0.7500; 249 y[i] = 57.2000; t[i++] = 0.8750; 250 y[i] = 53.2000; t[i++] = 1.0000; 251 y[i] = 42.5000; t[i++] = 1.2500; 252 y[i] = 26.8000; t[i++] = 1.7500; 253 y[i] = 20.4000; t[i++] = 2.2500; 254 y[i] = 26.8500; t[i++] = 1.7500; 255 y[i] = 21.0000; t[i++] = 2.2500; 256 y[i] = 16.4625; t[i++] = 2.7500; 257 y[i] = 12.5250; t[i++] = 3.2500; 258 y[i] = 10.5375; t[i++] = 3.7500; 259 y[i] = 8.5875; t[i++] = 4.2500; 260 y[i] = 7.1250; t[i++] = 4.7500; 261 y[i] = 6.1125; t[i++] = 5.2500; 262 y[i] = 5.9625; t[i++] = 5.7500; 263 y[i] = 74.1000; t[i++] = 0.5000; 264 y[i] = 67.3000; t[i++] = 0.6250; 265 y[i] = 60.8000; t[i++] = 0.7500; 266 y[i] = 55.5000; t[i++] = 0.8750; 267 y[i] = 50.3000; t[i++] = 1.0000; 268 y[i] = 41.0000; t[i++] = 1.2500; 269 y[i] = 29.4000; t[i++] = 1.7500; 270 y[i] = 20.4000; t[i++] = 2.2500; 271 y[i] = 29.3625; t[i++] = 1.7500; 272 y[i] = 21.1500; t[i++] = 2.2500; 273 y[i] = 16.7625; t[i++] = 2.7500; 274 y[i] = 13.2000; t[i++] = 3.2500; 275 y[i] = 10.8750; t[i++] = 3.7500; 276 y[i] = 8.1750; t[i++] = 4.2500; 277 y[i] = 7.3500; t[i++] = 4.7500; 278 y[i] = 5.9625; t[i++] = 5.2500; 279 y[i] = 5.6250; t[i++] = 5.7500; 280 y[i] = 81.5000; t[i++] = .5000; 281 y[i] = 62.4000; t[i++] = .7500; 282 y[i] = 32.5000; t[i++] = 1.5000; 283 y[i] = 12.4100; t[i++] = 3.0000; 284 y[i] = 13.1200; t[i++] = 3.0000; 285 y[i] = 15.5600; t[i++] = 3.0000; 286 y[i] = 5.6300; t[i++] = 6.0000; 287 y[i] = 78.0000; t[i++] = .5000; 288 y[i] = 59.9000; t[i++] = .7500; 289 y[i] = 33.2000; t[i++] = 1.5000; 290 y[i] = 13.8400; t[i++] = 3.0000; 291 y[i] = 12.7500; t[i++] = 3.0000; 292 y[i] = 14.6200; t[i++] = 3.0000; 293 y[i] = 3.9400; t[i++] = 6.0000; 294 y[i] = 76.8000; t[i++] = .5000; 295 y[i] = 61.0000; t[i++] = .7500; 296 y[i] = 32.9000; t[i++] = 1.5000; 297 y[i] = 13.8700; t[i++] = 3.0000; 298 y[i] = 11.8100; t[i++] = 3.0000; 299 y[i] = 13.3100; t[i++] = 3.0000; 300 y[i] = 5.4400; t[i++] = 6.0000; 301 y[i] = 78.0000; t[i++] = .5000; 302 y[i] = 63.5000; t[i++] = .7500; 303 y[i] = 33.8000; t[i++] = 1.5000; 304 y[i] = 12.5600; t[i++] = 3.0000; 305 y[i] = 5.6300; t[i++] = 6.0000; 306 y[i] = 12.7500; t[i++] = 3.0000; 307 y[i] = 13.1200; t[i++] = 3.0000; 308 y[i] = 5.4400; t[i++] = 6.0000; 309 y[i] = 76.8000; t[i++] = .5000; 310 y[i] = 60.0000; t[i++] = .7500; 311 y[i] = 47.8000; t[i++] = 1.0000; 312 y[i] = 32.0000; t[i++] = 1.5000; 313 y[i] = 22.2000; t[i++] = 2.0000; 314 y[i] = 22.5700; t[i++] = 2.0000; 315 y[i] = 18.8200; t[i++] = 2.5000; 316 y[i] = 13.9500; t[i++] = 3.0000; 317 y[i] = 11.2500; t[i++] = 4.0000; 318 y[i] = 9.0000; t[i++] = 5.0000; 319 y[i] = 6.6700; t[i++] = 6.0000; 320 y[i] = 75.8000; t[i++] = .5000; 321 y[i] = 62.0000; t[i++] = .7500; 322 y[i] = 48.8000; t[i++] = 1.0000; 323 y[i] = 35.2000; t[i++] = 1.5000; 324 y[i] = 20.0000; t[i++] = 2.0000; 325 y[i] = 20.3200; t[i++] = 2.0000; 326 y[i] = 19.3100; t[i++] = 2.5000; 327 y[i] = 12.7500; t[i++] = 3.0000; 328 y[i] = 10.4200; t[i++] = 4.0000; 329 y[i] = 7.3100; t[i++] = 5.0000; 330 y[i] = 7.4200; t[i++] = 6.0000; 331 y[i] = 70.5000; t[i++] = .5000; 332 y[i] = 59.5000; t[i++] = .7500; 333 y[i] = 48.5000; t[i++] = 1.0000; 334 y[i] = 35.8000; t[i++] = 1.5000; 335 y[i] = 21.0000; t[i++] = 2.0000; 336 y[i] = 21.6700; t[i++] = 2.0000; 337 y[i] = 21.0000; t[i++] = 2.5000; 338 y[i] = 15.6400; t[i++] = 3.0000; 339 y[i] = 8.1700; t[i++] = 4.0000; 340 y[i] = 8.5500; t[i++] = 5.0000; 341 y[i] = 10.1200; t[i++] = 6.0000; 342 y[i] = 78.0000; t[i++] = .5000; 343 y[i] = 66.0000; t[i++] = .6250; 344 y[i] = 62.0000; t[i++] = .7500; 345 y[i] = 58.0000; t[i++] = .8750; 346 y[i] = 47.7000; t[i++] = 1.0000; 347 y[i] = 37.8000; t[i++] = 1.2500; 348 y[i] = 20.2000; t[i++] = 2.2500; 349 y[i] = 21.0700; t[i++] = 2.2500; 350 y[i] = 13.8700; t[i++] = 2.7500; 351 y[i] = 9.6700; t[i++] = 3.2500; 352 y[i] = 7.7600; t[i++] = 3.7500; 353 y[i] = 5.4400; t[i++] = 4.2500; 354 y[i] = 4.8700; t[i++] = 4.7500; 355 y[i] = 4.0100; t[i++] = 5.2500; 356 y[i] = 3.7500; t[i++] = 5.7500; 357 y[i] = 24.1900; t[i++] = 3.0000; 358 y[i] = 25.7600; t[i++] = 3.0000; 359 y[i] = 18.0700; t[i++] = 3.0000; 360 y[i] = 11.8100; t[i++] = 3.0000; 361 y[i] = 12.0700; t[i++] = 3.0000; 362 y[i] = 16.1200; t[i++] = 3.0000; 363 y[i] = 70.8000; t[i++] = .5000; 364 y[i] = 54.7000; t[i++] = .7500; 365 y[i] = 48.0000; t[i++] = 1.0000; 366 y[i] = 39.8000; t[i++] = 1.5000; 367 y[i] = 29.8000; t[i++] = 2.0000; 368 y[i] = 23.7000; t[i++] = 2.5000; 369 y[i] = 29.6200; t[i++] = 2.0000; 370 y[i] = 23.8100; t[i++] = 2.5000; 371 y[i] = 17.7000; t[i++] = 3.0000; 372 y[i] = 11.5500; t[i++] = 4.0000; 373 y[i] = 12.0700; t[i++] = 5.0000; 374 y[i] = 8.7400; t[i++] = 6.0000; 375 y[i] = 80.7000; t[i++] = .5000; 376 y[i] = 61.3000; t[i++] = .7500; 377 y[i] = 47.5000; t[i++] = 1.0000; 378 y[i] = 29.0000; t[i++] = 1.5000; 379 y[i] = 24.0000; t[i++] = 2.0000; 380 y[i] = 17.7000; t[i++] = 2.5000; 381 y[i] = 24.5600; t[i++] = 2.0000; 382 y[i] = 18.6700; t[i++] = 2.5000; 383 y[i] = 16.2400; t[i++] = 3.0000; 384 y[i] = 8.7400; t[i++] = 4.0000; 385 y[i] = 7.8700; t[i++] = 5.0000; 386 y[i] = 8.5100; t[i++] = 6.0000; 387 y[i] = 66.7000; t[i++] = .5000; 388 y[i] = 59.2000; t[i++] = .7500; 389 y[i] = 40.8000; t[i++] = 1.0000; 390 y[i] = 30.7000; t[i++] = 1.5000; 391 y[i] = 25.7000; t[i++] = 2.0000; 392 y[i] = 16.3000; t[i++] = 2.5000; 393 y[i] = 25.9900; t[i++] = 2.0000; 394 y[i] = 16.9500; t[i++] = 2.5000; 395 y[i] = 13.3500; t[i++] = 3.0000; 396 y[i] = 8.6200; t[i++] = 4.0000; 397 y[i] = 7.2000; t[i++] = 5.0000; 398 y[i] = 6.6400; t[i++] = 6.0000; 399 y[i] = 13.6900; t[i++] = 3.0000; 400 y[i] = 81.0000; t[i++] = .5000; 401 y[i] = 64.5000; t[i++] = .7500; 402 y[i] = 35.5000; t[i++] = 1.5000; 403 y[i] = 13.3100; t[i++] = 3.0000; 404 y[i] = 4.8700; t[i++] = 6.0000; 405 y[i] = 12.9400; t[i++] = 3.0000; 406 y[i] = 5.0600; t[i++] = 6.0000; 407 y[i] = 15.1900; t[i++] = 3.0000; 408 y[i] = 14.6200; t[i++] = 3.0000; 409 y[i] = 15.6400; t[i++] = 3.0000; 410 y[i] = 25.5000; t[i++] = 1.7500; 411 y[i] = 25.9500; t[i++] = 1.7500; 412 y[i] = 81.7000; t[i++] = .5000; 413 y[i] = 61.6000; t[i++] = .7500; 414 y[i] = 29.8000; t[i++] = 1.7500; 415 y[i] = 29.8100; t[i++] = 1.7500; 416 y[i] = 17.1700; t[i++] = 2.7500; 417 y[i] = 10.3900; t[i++] = 3.7500; 418 y[i] = 28.4000; t[i++] = 1.7500; 419 y[i] = 28.6900; t[i++] = 1.7500; 420 y[i] = 81.3000; t[i++] = .5000; 421 y[i] = 60.9000; t[i++] = .7500; 422 y[i] = 16.6500; t[i++] = 2.7500; 423 y[i] = 10.0500; t[i++] = 3.7500; 424 y[i] = 28.9000; t[i++] = 1.7500; 425 y[i] = 28.9500; t[i++] = 1.7500; 426 PetscFunctionReturn(0); 427 } 428 429 /*TEST 430 431 build: 432 requires: !complex 433 434 test: 435 args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_pounders_delta 0.05 -tao_gatol 1.e-5 436 requires: !single 437 TODO: produces different output for many different systems 438 439 test: 440 suffix: 2 441 args: -tao_smonitor -tao_max_it 100 -wtype 1 -tao_type pounders -tao_pounders_delta 0.05 -tao_gatol 1.e-5 442 requires: !single 443 TODO: produces different output for many different systems 444 445 test: 446 suffix: 3 447 args: -tao_smonitor -tao_max_it 100 -wtype 2 -tao_type pounders -tao_pounders_delta 0.05 -tao_gatol 1.e-5 448 requires: !single 449 TODO: produces different output for many different systems 450 451 test: 452 suffix: 4 453 args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_pounders_delta 0.05 -pounders_subsolver_tao_type blmvm -tao_gatol 1.e-5 454 requires: !single 455 TODO: produces different output for many different systems 456 457 TEST*/ 458