1 /* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */
2
3 /* Include "petsctao.h" so we can use TAO solvers. */
4 #include <petsctao.h>
5
6 static char help[] =
7 "This example demonstrates use of the TAO package to \n\
8 solve an unconstrained minimization problem. This example is based on a \n\
9 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
10 boundary values along the edges of the domain, the objective is to find the\n\
11 surface with the minimal area that satisfies the boundary conditions.\n\
12 The command line options are:\n\
13 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
14 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
15 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n";
16
17 /*T
18 Concepts: TAO^Solving an unconstrained minimization problem
19 Routines: TaoCreate(); TaoSetType();
20 Routines: TaoSetInitialVector();
21 Routines: TaoSetObjectiveAndGradientRoutine();
22 Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
23 Routines: TaoGetKSP(); TaoSolve();
24 Routines: TaoDestroy();
25 Processors: 1
26 T*/
27
28
29
30 /*
31 User-defined application context - contains data needed by the
32 application-provided call-back routines, FormFunctionGradient()
33 and FormHessian().
34 */
35 typedef struct {
36 PetscInt mx, my; /* discretization in x, y directions */
37 PetscReal *bottom, *top, *left, *right; /* boundary values */
38 Mat H;
39 } AppCtx;
40
41 /* -------- User-defined Routines --------- */
42
43 static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
44 static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
45 static PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
46 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
47 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
48
main(int argc,char ** argv)49 int main(int argc, char **argv)
50 {
51 PetscErrorCode ierr; /* used to check for functions returning nonzeros */
52 PetscInt N; /* Size of vector */
53 PetscMPIInt size; /* Number of processors */
54 Vec x; /* solution, gradient vectors */
55 KSP ksp; /* PETSc Krylov subspace method */
56 PetscBool flg; /* A return value when checking for user options */
57 Tao tao; /* Tao solver context */
58 AppCtx user; /* user-defined work context */
59
60 /* Initialize TAO,PETSc */
61 ierr = PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;
62
63 ierr = MPI_Comm_size(MPI_COMM_WORLD,&size);CHKERRQ(ierr);
64 if (size >1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors");
65
66 /* Specify default dimension of the problem */
67 user.mx = 4; user.my = 4;
68
69 /* Check for any command line arguments that override defaults */
70 ierr = PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);CHKERRQ(ierr);
71 ierr = PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);CHKERRQ(ierr);
72
73 ierr = PetscPrintf(PETSC_COMM_SELF,"\n---- Minimum Surface Area Problem -----\n");CHKERRQ(ierr);
74 ierr = PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",user.mx,user.my);CHKERRQ(ierr);
75
76 /* Calculate any derived values from parameters */
77 N = user.mx*user.my;
78
79 /* Create TAO solver and set desired solution method */
80 ierr = TaoCreate(PETSC_COMM_SELF,&tao);CHKERRQ(ierr);
81 ierr = TaoSetType(tao,TAOLMVM);CHKERRQ(ierr);
82
83 /* Initialize minsurf application data structure for use in the function evaluations */
84 ierr = MSA_BoundaryConditions(&user);CHKERRQ(ierr); /* Application specific routine */
85
86 /*
87 Create a vector to hold the variables. Compute an initial solution.
88 Set this vector, which will be used by TAO.
89 */
90 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr);
91 ierr = MSA_InitialPoint(&user,x);CHKERRQ(ierr); /* Application specific routine */
92 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr); /* A TAO routine */
93
94 /* Provide TAO routines for function, gradient, and Hessian evaluation */
95 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);CHKERRQ(ierr);
96
97 /* Create a matrix data structure to store the Hessian. This structure will be used by TAO */
98 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,7,NULL,&(user.H));CHKERRQ(ierr);
99 ierr = MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
100 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);CHKERRQ(ierr);
101
102 /* Check for any TAO command line options */
103 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
104
105 /* Limit the number of iterations in the KSP linear solver */
106 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
107 if (ksp) {
108 ierr = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,user.mx*user.my);CHKERRQ(ierr);
109 }
110
111 /* SOLVE THE APPLICATION */
112 ierr = TaoSolve(tao);CHKERRQ(ierr);
113
114 ierr = TaoDestroy(&tao);CHKERRQ(ierr);
115 ierr = VecDestroy(&x);CHKERRQ(ierr);
116 ierr = MatDestroy(&user.H);CHKERRQ(ierr);
117 ierr = PetscFree(user.bottom);CHKERRQ(ierr);
118 ierr = PetscFree(user.top);CHKERRQ(ierr);
119 ierr = PetscFree(user.left);CHKERRQ(ierr);
120 ierr = PetscFree(user.right);CHKERRQ(ierr);
121
122 ierr = PetscFinalize();
123 return ierr;
124 }
125
126 /* -------------------------------------------------------------------- */
127
128 /* FormFunctionGradient - Evaluates function and corresponding gradient.
129
130 Input Parameters:
131 . tao - the Tao context
132 . X - input vector
133 . userCtx - optional user-defined context, as set by TaoSetFunctionGradient()
134
135 Output Parameters:
136 . fcn - the newly evaluated function
137 . G - vector containing the newly evaluated gradient
138 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * fcn,Vec G,void * userCtx)139 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *fcn,Vec G,void *userCtx)
140 {
141 AppCtx *user = (AppCtx *) userCtx;
142 PetscErrorCode ierr;
143 PetscInt i,j,row;
144 PetscInt mx=user->mx, my=user->my;
145 PetscReal rhx=mx+1, rhy=my+1;
146 PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy, ft=0;
147 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
148 PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
149 PetscReal zero=0.0;
150 PetscScalar *g;
151 const PetscScalar *x;
152
153 PetscFunctionBeginUser;
154 ierr = VecSet(G, zero);CHKERRQ(ierr);
155
156 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
157 ierr = VecGetArray(G,&g);CHKERRQ(ierr);
158
159 /* Compute function over the locally owned part of the mesh */
160 for (j=0; j<my; j++){
161 for (i=0; i< mx; i++){
162 row=(j)*mx + (i);
163 xc = x[row];
164 xlt=xrb=xl=xr=xb=xt=xc;
165 if (i==0){ /* left side */
166 xl = user->left[j+1];
167 xlt = user->left[j+2];
168 } else {
169 xl = x[row-1];
170 }
171
172 if (j==0){ /* bottom side */
173 xb = user->bottom[i+1];
174 xrb = user->bottom[i+2];
175 } else {
176 xb = x[row-mx];
177 }
178
179 if (i+1 == mx){ /* right side */
180 xr = user->right[j+1];
181 xrb = user->right[j];
182 } else {
183 xr = x[row+1];
184 }
185
186 if (j+1==0+my){ /* top side */
187 xt = user->top[i+1];
188 xlt = user->top[i];
189 }else {
190 xt = x[row+mx];
191 }
192
193 if (i>0 && j+1<my){
194 xlt = x[row-1+mx];
195 }
196 if (j>0 && i+1<mx){
197 xrb = x[row+1-mx];
198 }
199
200 d1 = (xc-xl);
201 d2 = (xc-xr);
202 d3 = (xc-xt);
203 d4 = (xc-xb);
204 d5 = (xr-xrb);
205 d6 = (xrb-xb);
206 d7 = (xlt-xl);
207 d8 = (xt-xlt);
208
209 df1dxc = d1*hydhx;
210 df2dxc = (d1*hydhx + d4*hxdhy);
211 df3dxc = d3*hxdhy;
212 df4dxc = (d2*hydhx + d3*hxdhy);
213 df5dxc = d2*hydhx;
214 df6dxc = d4*hxdhy;
215
216 d1 *= rhx;
217 d2 *= rhx;
218 d3 *= rhy;
219 d4 *= rhy;
220 d5 *= rhy;
221 d6 *= rhx;
222 d7 *= rhy;
223 d8 *= rhx;
224
225 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
226 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
227 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
228 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
229 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
230 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
231
232 ft = ft + (f2 + f4);
233
234 df1dxc /= f1;
235 df2dxc /= f2;
236 df3dxc /= f3;
237 df4dxc /= f4;
238 df5dxc /= f5;
239 df6dxc /= f6;
240
241 g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc)/2.0;
242 }
243 }
244
245 for (j=0; j<my; j++){ /* left side */
246 d3 = (user->left[j+1] - user->left[j+2])*rhy;
247 d2 = (user->left[j+1] - x[j*mx])*rhx;
248 ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
249 }
250
251 for (i=0; i<mx; i++){ /* bottom */
252 d2 = (user->bottom[i+1]-user->bottom[i+2])*rhx;
253 d3 = (user->bottom[i+1]-x[i])*rhy;
254 ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
255 }
256
257 for (j=0; j< my; j++){ /* right side */
258 d1 = (x[(j+1)*mx-1]-user->right[j+1])*rhx;
259 d4 = (user->right[j]-user->right[j+1])*rhy;
260 ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
261 }
262
263 for (i=0; i<mx; i++){ /* top side */
264 d1 = (x[(my-1)*mx + i] - user->top[i+1])*rhy;
265 d4 = (user->top[i+1] - user->top[i])*rhx;
266 ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
267 }
268
269 /* Bottom left corner */
270 d1 = (user->left[0]-user->left[1])*rhy;
271 d2 = (user->bottom[0]-user->bottom[1])*rhx;
272 ft += PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
273
274 /* Top right corner */
275 d1 = (user->right[my+1] - user->right[my])*rhy;
276 d2 = (user->top[mx+1] - user->top[mx])*rhx;
277 ft += PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
278
279 (*fcn)=ft*area;
280
281 /* Restore vectors */
282 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
283 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr);
284 ierr = PetscLogFlops(67.0*mx*my);CHKERRQ(ierr);
285 PetscFunctionReturn(0);
286 }
287
288 /* ------------------------------------------------------------------- */
289 /*
290 FormHessian - Evaluates the Hessian matrix.
291
292 Input Parameters:
293 . tao - the Tao context
294 . x - input vector
295 . ptr - optional user-defined context, as set by TaoSetHessian()
296
297 Output Parameters:
298 . H - Hessian matrix
299 . Hpre - optionally different preconditioning matrix
300 . flg - flag indicating matrix structure
301
302 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)303 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
304 {
305 PetscErrorCode ierr;
306 AppCtx *user = (AppCtx *) ptr;
307
308 PetscFunctionBeginUser;
309 /* Evaluate the Hessian entries*/
310 ierr = QuadraticH(user,X,H);CHKERRQ(ierr);
311 PetscFunctionReturn(0);
312 }
313
314 /* ------------------------------------------------------------------- */
315 /*
316 QuadraticH - Evaluates the Hessian matrix.
317
318 Input Parameters:
319 . user - user-defined context, as set by TaoSetHessian()
320 . X - input vector
321
322 Output Parameter:
323 . H - Hessian matrix
324 */
QuadraticH(AppCtx * user,Vec X,Mat Hessian)325 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
326 {
327 PetscErrorCode ierr;
328 PetscInt i,j,k,row;
329 PetscInt mx=user->mx, my=user->my;
330 PetscInt col[7];
331 PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
332 PetscReal rhx=mx+1, rhy=my+1;
333 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
334 PetscReal hl,hr,ht,hb,hc,htl,hbr;
335 const PetscScalar *x;
336 PetscReal v[7];
337
338 /* Get pointers to vector data */
339 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
340
341 /* Initialize matrix entries to zero */
342 ierr = MatZeroEntries(Hessian);CHKERRQ(ierr);
343
344 /* Set various matrix options */
345 ierr = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
346
347 /* Compute Hessian over the locally owned part of the mesh */
348 for (i=0; i< mx; i++){
349 for (j=0; j<my; j++){
350
351 row=(j)*mx + (i);
352
353 xc = x[row];
354 xlt=xrb=xl=xr=xb=xt=xc;
355
356 /* Left side */
357 if (i==0){
358 xl= user->left[j+1];
359 xlt = user->left[j+2];
360 } else {
361 xl = x[row-1];
362 }
363
364 if (j==0){
365 xb=user->bottom[i+1];
366 xrb = user->bottom[i+2];
367 } else {
368 xb = x[row-mx];
369 }
370
371 if (i+1 == mx){
372 xr=user->right[j+1];
373 xrb = user->right[j];
374 } else {
375 xr = x[row+1];
376 }
377
378 if (j+1==my){
379 xt=user->top[i+1];
380 xlt = user->top[i];
381 }else {
382 xt = x[row+mx];
383 }
384
385 if (i>0 && j+1<my){
386 xlt = x[row-1+mx];
387 }
388 if (j>0 && i+1<mx){
389 xrb = x[row+1-mx];
390 }
391
392
393 d1 = (xc-xl)*rhx;
394 d2 = (xc-xr)*rhx;
395 d3 = (xc-xt)*rhy;
396 d4 = (xc-xb)*rhy;
397 d5 = (xrb-xr)*rhy;
398 d6 = (xrb-xb)*rhx;
399 d7 = (xlt-xl)*rhy;
400 d8 = (xlt-xt)*rhx;
401
402 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
403 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
404 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
405 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
406 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
407 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
408
409
410 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
411 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
412 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
413 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
414
415 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
416 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
417
418 hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
419 (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
420
421 hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
422
423 k=0;
424 if (j>0){
425 v[k]=hb; col[k]=row - mx; k++;
426 }
427
428 if (j>0 && i < mx -1){
429 v[k]=hbr; col[k]=row - mx+1; k++;
430 }
431
432 if (i>0){
433 v[k]= hl; col[k]=row - 1; k++;
434 }
435
436 v[k]= hc; col[k]=row; k++;
437
438 if (i < mx-1){
439 v[k]= hr; col[k]=row+1; k++;
440 }
441
442 if (i>0 && j < my-1){
443 v[k]= htl; col[k] = row+mx-1; k++;
444 }
445
446 if (j < my-1){
447 v[k]= ht; col[k] = row+mx; k++;
448 }
449
450 /*
451 Set matrix values using local numbering, which was defined
452 earlier, in the main routine.
453 */
454 ierr = MatSetValues(Hessian,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(ierr);
455 }
456 }
457
458 /* Restore vectors */
459 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
460
461 /* Assemble the matrix */
462 ierr = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
463 ierr = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
464
465 ierr = PetscLogFlops(199.0*mx*my);CHKERRQ(ierr);
466 PetscFunctionReturn(0);
467 }
468
469 /* ------------------------------------------------------------------- */
470 /*
471 MSA_BoundaryConditions - Calculates the boundary conditions for
472 the region.
473
474 Input Parameter:
475 . user - user-defined application context
476
477 Output Parameter:
478 . user - user-defined application context
479 */
MSA_BoundaryConditions(AppCtx * user)480 static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
481 {
482 PetscErrorCode ierr;
483 PetscInt i,j,k,limit=0;
484 PetscInt maxits=5;
485 PetscInt mx=user->mx,my=user->my;
486 PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
487 PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
488 PetscReal fnorm,det,hx,hy,xt=0,yt=0;
489 PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
490 PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
491 PetscReal *boundary;
492
493 PetscFunctionBeginUser;
494 bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2;
495
496 ierr = PetscMalloc1(bsize,&user->bottom);CHKERRQ(ierr);
497 ierr = PetscMalloc1(tsize,&user->top);CHKERRQ(ierr);
498 ierr = PetscMalloc1(lsize,&user->left);CHKERRQ(ierr);
499 ierr = PetscMalloc1(rsize,&user->right);CHKERRQ(ierr);
500
501 hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
502
503 for (j=0; j<4; j++){
504 if (j==0){
505 yt=b;
506 xt=l;
507 limit=bsize;
508 boundary=user->bottom;
509 } else if (j==1){
510 yt=t;
511 xt=l;
512 limit=tsize;
513 boundary=user->top;
514 } else if (j==2){
515 yt=b;
516 xt=l;
517 limit=lsize;
518 boundary=user->left;
519 } else { /* if (j==3) */
520 yt=b;
521 xt=r;
522 limit=rsize;
523 boundary=user->right;
524 }
525
526 for (i=0; i<limit; i++){
527 u1=xt;
528 u2=-yt;
529 for (k=0; k<maxits; k++){
530 nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
531 nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
532 fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
533 if (fnorm <= tol) break;
534 njac11=one+u2*u2-u1*u1;
535 njac12=two*u1*u2;
536 njac21=-two*u1*u2;
537 njac22=-one - u1*u1 + u2*u2;
538 det = njac11*njac22-njac21*njac12;
539 u1 = u1-(njac22*nf1-njac12*nf2)/det;
540 u2 = u2-(njac11*nf2-njac21*nf1)/det;
541 }
542
543 boundary[i]=u1*u1-u2*u2;
544 if (j==0 || j==1) {
545 xt=xt+hx;
546 } else { /* if (j==2 || j==3) */
547 yt=yt+hy;
548 }
549 }
550 }
551 PetscFunctionReturn(0);
552 }
553
554 /* ------------------------------------------------------------------- */
555 /*
556 MSA_InitialPoint - Calculates the initial guess in one of three ways.
557
558 Input Parameters:
559 . user - user-defined application context
560 . X - vector for initial guess
561
562 Output Parameters:
563 . X - newly computed initial guess
564 */
MSA_InitialPoint(AppCtx * user,Vec X)565 static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
566 {
567 PetscInt start=-1,i,j;
568 PetscErrorCode ierr;
569 PetscReal zero=0.0;
570 PetscBool flg;
571
572 PetscFunctionBeginUser;
573 ierr = VecSet(X, zero);CHKERRQ(ierr);
574 ierr = PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);CHKERRQ(ierr);
575
576 if (flg && start==0){ /* The zero vector is reasonable */
577 ierr = VecSet(X, zero);CHKERRQ(ierr);
578 } else { /* Take an average of the boundary conditions */
579 PetscInt row;
580 PetscInt mx=user->mx,my=user->my;
581 PetscReal *x;
582
583 /* Get pointers to vector data */
584 ierr = VecGetArray(X,&x);CHKERRQ(ierr);
585 /* Perform local computations */
586 for (j=0; j<my; j++){
587 for (i=0; i< mx; i++){
588 row=(j)*mx + (i);
589 x[row] = (((j+1)*user->bottom[i+1]+(my-j+1)*user->top[i+1])/(my+2)+ ((i+1)*user->left[j+1]+(mx-i+1)*user->right[j+1])/(mx+2))/2.0;
590 }
591 }
592 /* Restore vectors */
593 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
594 }
595 PetscFunctionReturn(0);
596 }
597
598
599 /*TEST
600
601 build:
602 requires: !complex
603
604 test:
605 args: -tao_smonitor -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4
606 requires: !single
607
608 test:
609 suffix: 2
610 args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
611 requires: !single
612
613 test:
614 suffix: 3
615 args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
616 requires: !single
617
618 test:
619 suffix: 4
620 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4
621 requires: !single
622
623 test:
624 suffix: 5
625 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4
626 requires: !single
627
628 test:
629 suffix: 6
630 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4
631 requires: !single
632
633 test:
634 suffix: 7
635 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
636 requires: !single
637
638 test:
639 suffix: 8
640 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
641 requires: !single
642
643 test:
644 suffix: 9
645 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
646 requires: !single
647
648 TEST*/
649