/*$Id: minsurf3.c,v 1.1 2002/05/13 19:32:57 benson Exp $*/ /* Program usage: mpirun -np <proc> minsurf3 [-help] [all TAO options] */ /* minsurf1: boundary values like in COPS example: - parabola for top and bottom boundaries - zero o.w. */ /* Include "tao.h" so we can use TAO solvers. petscda.h for distributed array ad_deriv.h for AD gradient */ #include "petscda.h" #include "tao.h" #include "taodaapplication.h" static char help[] = "Given a rectangular 2-D domain and boundary values along \n\ the edges of the domain, the objective is to find the surface with \n\ the minimal area that satisfies the boundary conditions.\n\ The command line options are:\n\ -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\ -byelement, if computation is made by functions on rectangular elements\n\ -adic, if AD is used (AD is not used by default)\n\ -bottom <b>, where <b> = bottom bound on the rectangular domain\n\ -top <t>, where <t> = bottom bound on the rectangular domain\n\ -left <l>, where <l> = left bound on the rectangular domain\n\ -right <r>, where <r> = right bound on the rectangular domain\n\ -cops <r>, if COPS boundaries should be use (default MINPACK bounds)\n\ -obst <obst>, where <obst> = 1 is obstacle is present, zero otherwise \n\n"; /*T Concepts: TAO - Solving a bounded minimization problem Routines: TaoInitialize(); TaoFinalize(); Routines: TaoCreate(); TaoDestroy(); Routines: DAAppSetVariableBoundsRoutine(); Routines: DAAppSetElementObjectiveAndGradientRoutine(); Routines: DAAppSetElementHessianRoutine(); Routines: DAAppSetObjectiveAndGradientRoutine(); Routines: DAAppSetADElementFunctionGradient(); Routines: DAAppSetHessianRoutine(); Routines: TaoSetOptions(); Routines: TaoAppGetSolutionStatus(); TaoDAAppSolve(); Routines: DAAppSetBeforeMonitor(); TaoView(); Routines: DAAppGetSolution(); Routines: DAAppGetInterpolationMatrix(); Processors: n T*/ /* User-defined application context - contains data needed by the application-provided call-back routines. */ /* This structure is used only when an ADIC generated gradient is used. An InactiveDouble type is a double */ typedef struct { InactiveDouble hx, hy; /* increment size in both directions */ InactiveDouble area; /* area of the triangles */ } ADFGCtx; int ad_MSurfLocalFunction(int[2], DERIV_TYPE[4], DERIV_TYPE*, void*); typedef struct { PetscReal b, t, l, r; /* domain boundaries */ PetscReal bheight; /* height of obstacle */ PetscReal fx, fy; /* relative size of obstacle */ double hx, hy; /* increment size in both directions */ double area; /* area of the triangles */ int mx, my; /* discretization including boundaries */ int bmx, bmy; /* discretization of obstacle */ ADFGCtx fgctx; /* Used only when an ADIC generated gradient is used */ } AppCtx; /* User-defined routines */ static int AppCtxInitialize(void *); static int MSurfLocalFunctionGradient(int[2], double[4], double *, double[4], void *); static int WholeMSurfFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*); static int MSurfLocalHessian(int[2], double[4], double[4][4], void *); static int WholeMSurfHessian(TAO_APPLICATION,DA,Vec,Mat,void*); static int COPS_Bounds(TAO_APPLICATION, DA, Vec, Vec, void*); static int MINPACK_Bounds(TAO_APPLICATION, DA, Vec, Vec, void *); static int MyGridMonitorBefore(TAO_APPLICATION, DA, int, void *); #undef __FUNCT__ #define __FUNCT__ "main" int main( int argc, char **argv ) { int info; /* used to check for functions returning nonzeros */ int iter,nlevels; int Nx,Ny; double ff,gnorm; DA DAarray[20]; Vec X; PetscTruth flg, PreLoad = PETSC_TRUE; /* flags */ TaoMethod method = "tao_bnls"; /* minimization method */ TaoTerminateReason reason; TAO_SOLVER tao; /* TAO_SOLVER solver context */ TAO_APPLICATION MSurfApp; /* The PETSc application */ AppCtx user; /* user-defined work context */ /* Initialize TAO */ PetscInitialize(&argc, &argv, (char *)0, help); TaoInitialize(&argc, &argv, (char *)0, help); PreLoadBegin(PreLoad,"Solve"); info = AppCtxInitialize((void*)&user); CHKERRQ(info); nlevels=4; info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info); if (PreLoadIt == 0) { nlevels = 1; user.mx = 11; user.my = 11; } PetscPrintf(MPI_COMM_WORLD,"\n---- Minimal Surface Problem (simple boundary) -----\n\n"); /* Let PETSc determine the vector distribution */ Nx = PETSC_DECIDE; Ny = PETSC_DECIDE; /* Create distributed array (DA) to manage parallel grid and vectors */ info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx, user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info); for (iter=1;iter<nlevels;iter++){ info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info); } /* Create TAO solver and set desired solution method */ info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info); info = TaoApplicationCreate(PETSC_COMM_WORLD, &MSurfApp); CHKERRQ(info); info = TaoAppSetDAApp(MSurfApp,DAarray,nlevels); CHKERRQ(info); /* Sets routines for function, gradient and bounds evaluation */ info = PetscOptionsHasName(TAO_NULL, "-cops", &flg); CHKERRQ(info); if (flg){ info = DAAppSetVariableBoundsRoutine(MSurfApp,COPS_Bounds,(void *)&user); CHKERRQ(info); } else { info = DAAppSetVariableBoundsRoutine(MSurfApp,MINPACK_Bounds,(void *)&user); CHKERRQ(info); } info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info); if (flg) { /* Sets routines for function and gradient evaluation, element by element */ info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info); if (flg) { info = DAAppSetADElementFunctionGradient(MSurfApp,ad_MSurfLocalFunction,150,(void *)&user.fgctx); CHKERRQ(info); } else { info = DAAppSetElementObjectiveAndGradientRoutine(MSurfApp,MSurfLocalFunctionGradient,36,(void *)&user); CHKERRQ(info); } /* Sets routines for Hessian evaluation, element by element */ info = DAAppSetElementHessianRoutine(MSurfApp,MSurfLocalHessian,87,(void*)&user); CHKERRQ(info); } else { /* Sets routines for function and gradient evaluation, all in one routine */ info = DAAppSetObjectiveAndGradientRoutine(MSurfApp,WholeMSurfFunctionGradient,(void *)&user); CHKERRQ(info); /* Sets routines for Hessian evaluation, all in one routine */ info = DAAppSetHessianRoutine(MSurfApp,WholeMSurfHessian,(void*)&user); CHKERRQ(info); } info = DAAppSetBeforeMonitor(MSurfApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info); info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info); if (flg){ info = DAAppPrintStageTimes(MSurfApp); CHKERRQ(info); info = DAAppPrintInterpolationError(MSurfApp); CHKERRQ(info); } info = TaoAppSetRelativeTolerance(MSurfApp,1.0e-8); CHKERRQ(info); info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info); info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info); /* Check for any tao command line options */ info = TaoSetOptions(MSurfApp,tao); CHKERRQ(info); /* SOLVE THE APPLICATION */ info = TaoDAAppSolve(MSurfApp,tao); CHKERRQ(info); /* Get information on termination */ info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info); if (reason <= 0 ){ PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n"); PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm); } info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info); if (flg){ info = DAAppGetSolution(MSurfApp,nlevels-1,&X); CHKERRQ(info); info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info); } /* To View TAO solver information */ // info = TaoView(tao); CHKERRQ(info); /* Free TAO data structures */ info = TaoDestroy(tao); CHKERRQ(info); info = TaoAppDestroy(MSurfApp); CHKERRQ(info); /* Free PETSc data structures */ for (iter=0;iter<nlevels;iter++){ info = DADestroy(DAarray[iter]); CHKERRQ(info); } PreLoadEnd(); /* Finalize TAO */ TaoFinalize(); PetscFinalize(); return 0; } /* main */ /*----- The following two routines MyGridMonitorBefore MyGridMonitorAfter help diplay info of iterations at every grid level */ #undef __FUNCT__ #define __FUNCT__ "MyGridMonitorBefore" static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, int level, void *ctx) { AppCtx *user = (AppCtx*)ctx; int info,mx,my; info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL, PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info); user->mx = mx; user->my = my; user->hx = (user->r - user->l) / (user->mx - 1); user->hy = (user->t - user->b) / (user->my - 1); user->area = 0.5 * user->hx * user->hy; user->fgctx.hx = user->hx; user->fgctx.hy = user->hy; user->fgctx.area = user->area; user->bmx=(int)((mx+1)*user->fx); user->bmy=(int)((my+1)*user->fy); PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my); return 0; } #undef __FUNCT__ #define __FUNCT__ "AppCtxInitialize" /* AppCtxInitialize - Sets initial values for the application context parameters Input: ptr - void user-defined application context Output: ptr - user-defined application context with the default or user-provided parameters */ static int AppCtxInitialize(void *ptr) { AppCtx *user = (AppCtx*)ptr; PetscTruth flg; int info; /* Specify default parameters */ user->mx = user->my = 11; user->b = -0.5; user->t = 0.5; user->l = -0.5; user->r = 0.5; user->fx=0.5; user->fy=0.5; user->bheight=0.0; /* Check for command line arguments that override defaults */ info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info); info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info); info = PetscOptionsGetReal(TAO_NULL, "-bottom", &user->b, &flg); CHKERRQ(info); info = PetscOptionsGetReal(TAO_NULL, "-top", &user->t, &flg); CHKERRQ(info); info = PetscOptionsGetReal(TAO_NULL, "-left", &user->l, &flg); CHKERRQ(info); info = PetscOptionsGetReal(TAO_NULL, "-right", &user->r, &flg); CHKERRQ(info); info = PetscOptionsGetReal(PETSC_NULL,"-bmx",&user->fx,&flg); CHKERRQ(info); info = PetscOptionsGetReal(PETSC_NULL,"-bmy",&user->fy,&flg); CHKERRQ(info); info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user->bheight,&flg); CHKERRQ(info); user->hx = (user->r - user->l) / (user->mx - 1); user->hy = (user->t - user->b) / (user->my - 1); user->area = 0.5 * user->hx * user->hy; info = PetscLogFlops(8); CHKERRQ(info); return 0; } /* AppCtxInitialize */ /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/ #undef __FUNCT__ #define __FUNCT__ "COPS_Bounds" /* FormBounds - Forms bounds on the variables Input: user - user-defined application context Output: XL - vector of lower bounds XU - vector of upper bounds */ static int COPS_Bounds(TAO_APPLICATION tao, DA da, Vec XL, Vec XU, void *ptr) { AppCtx *user = (AppCtx*)ptr; PetscTruth flg; int i, j, mx, my, info, xs, xm, ys, ym; double lb = -TAO_INFINITY; double ub = TAO_INFINITY; double **xl, **xu; double xi, xi1, xi2; double cx, cy, radx, rady, height = 1.0; info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL, PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info); user->hx = (user->r - user->l) / (mx - 1); user->hy = (user->t - user->b) / (my - 1); user->area = 0.5 * user->hx * user->hy; info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info); info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info); info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info); /* Compute default bounds */ for (j = ys; j < ys+ym; j++){ for (i = xs; i < xs+xm; i++){ if (j == 0 || j == my - 1) { xi = user->l + i * user->hx; xl[j][i] = xu[j][i] = -4 * (xi - user->l) * (xi - user->r); } else if (i == 0 || i == mx - 1) { xl[j][i] = xu[j][i] = 0.0; } else { xl[j][i] = lb; xu[j][i] = ub; } } } /* Adjust lower bound if obstacle is present */ info = PetscOptionsHasName(PETSC_NULL, "-obst", &flg); CHKERRQ(info); if (flg) { radx = (user->r - user->l) * 0.25; cx = user->l + 2.0 * radx; rady = (user->t - user->b) * 0.25; cy = user->b + 2.0 * rady; for (j = ys; j < ys+ym; j++){ for (i = xs; i < xs+xm; i++){ xi1 = user->l + i * user->hx; xi2 = user->b + j * user->hy; if ( fabs(xi1 - cx) <= radx && fabs(xi2 - cy) <= rady ) { xl[j][i] = height; } } } info = PetscLogFlops(8 + xm * ym * 6); CHKERRQ(info); } info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info); info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info); info = PetscLogFlops(12 * ym + 6); CHKERRQ(info); return 0; } /* DAGetBounds2d */ #undef __FUNCT__ #define __FUNCT__ "MINPACK_Bounds" static int MINPACK_Bounds(TAO_APPLICATION tao, DA da, Vec XL,Vec XU, void *user1){ AppCtx *user=(AppCtx*)user1; int i,j,k,limit=0,info,maxits=5; int xs,ys,xm,ym; int mx, my, bmy, bmx; int row=0, bsize=0, lsize=0, tsize=0, rsize=0; double bheight=user->bheight; double one=1.0, two=2.0, three=3.0, tol=1e-10; double fnorm,det,hx,hy,xt=0,yt=0; double u1,u2,nf1,nf2,njac11,njac12,njac21,njac22; double b=-0.5, t=0.5, l=-0.5, r=0.5; PetscScalar scl,*xl,*xu, lb=TAO_NINFINITY, ub=TAO_INFINITY; PetscTruth flg; double **xll; info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL, PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info); bheight=user->bheight,lb=-1000; ub=1000; bmx=user->bmx; bmy=user->bmy; info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info); info = VecSet(XL, lb); CHKERRQ(info); info = VecSet(XU, ub); CHKERRQ(info); /* Get pointers to vector data */ info = DAVecGetArray(da,XL,(void**)&xll); CHKERRQ(info); if (ys==0) bsize=xm; if (xs==0) lsize=ym; if (xs+xm==mx) rsize=ym; if (ys+ym==my) tsize=xm; hx= (r-l)/(mx-1); hy=(t-b)/(my-1); /* Compute the optional lower box */ for (i=xs; i< xs+xm; i++){ for (j=ys; j<ys+ym; j++){ if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 && j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){ xll[j][i] = bheight; } } } info = DAVecRestoreArray(da,XL,(void**)&xll); CHKERRQ(info); /* Boundary Values */ info = VecGetArray(XL,&xl); CHKERRQ(info); info = VecGetArray(XU,&xu); CHKERRQ(info); for (j=0; j<4; j++){ if (j==0){ yt=b; xt=l+hx*xs; limit=bsize; scl=1.0; info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg); } else if (j==1){ yt=t; xt=l+hx*xs; limit=tsize; scl=1.0; info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg); } else if (j==2){ yt=b+hy*ys; xt=l; limit=lsize; scl=1.0; info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg); } else if (j==3){ yt=b+hy*ys; xt=r; limit=rsize; scl=1.0; info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg); } for (i=0; i<limit; i++){ u1=xt; u2=-yt; for (k=0; k<maxits; k++){ nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt; nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt; fnorm=sqrt(nf1*nf1+nf2*nf2); if (fnorm <= tol) break; njac11=one+u2*u2-u1*u1; njac12=two*u1*u2; njac21=-two*u1*u2; njac22=-one - u1*u1 + u2*u2; det = njac11*njac22-njac21*njac12; u1 = u1-(njac22*nf1-njac12*nf2)/det; u2 = u2-(njac11*nf2-njac21*nf1)/det; } if (j==0){ row = i; } else if (j==1){ row = (ym-1)*xm+i; } else if (j==2){ row = (xm*i); } else if (j==3){ row = xm*(i+1)-1; } xl[row]=(u1*u1-u2*u2)*scl; xu[row]=(u1*u1-u2*u2)*scl; if (j==0 || j==1) { xt=xt+hx; } else if (j==2 || j==3){ yt=yt+hy; } } } info = VecRestoreArray(XL,&xl); CHKERRQ(info); info = VecRestoreArray(XU,&xu); CHKERRQ(info); return 0; } /*------- USER-DEFINED: routine to evaluate the function and gradient at a local (rectangular element) level -------*/ #undef __FUNCT__ #define __FUNCT__ "MSurfLocalFunctionGradient" /* MSurfLocalFunctionGradient - Evaluates function and gradient over the local rectangular element Input: coor - vector with the indices of the position of current element in the first, second and third directions x - current point (values over the current rectangular element) ptr - user-defined application context Output: f - value of the objective funtion at the local rectangular element g - gradient of the local function */ static int MSurfLocalFunctionGradient(int coor[2], double x[4], double *f, double g[4], void *ptr) { AppCtx *user = (AppCtx*)ptr; double hx, hy, area; double dvdx, dvdy, flow, fup; double d1,d2; hx = user->hx; hy = user->hy; area = user->area; /* lower triangle contribution */ dvdx = (x[0] - x[1]); dvdy = (x[0] - x[2]); g[1] = (dvdx * (hy/hx))/(-2); g[2] = (dvdy * (hx/hy))/(-2); dvdx /= hx; dvdy /= hy; flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); *f=flow; g[1] /= flow; g[2] /= flow; g[0] = -(g[1]+g[2]); /* upper triangle contribution */ dvdx = (x[3] - x[2]); dvdy = (x[3] - x[1]); d1 = (dvdy*(hx/hy))/(-2); d2 = (dvdx*(hy/hx))/(-2); dvdx /= hx; dvdy /= hy; fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); *f += fup; g[1] += d1/fup; g[2] += d2/fup; g[3] = -(d1+d2)/fup; *f *= area; return 0; } /* MSurfLocalFunctionGradient */ #undef __FUNCT__ #define __FUNCT__ "MSurfLocalHessian" /* MSurfLocalHessian - Computes the Hessian of the local (partial) function defined over the current rectangle Input: coor - vector with the indices of the position of current element in the first, second and third directions x - current local solution (over the rectangle only) ptr - user-defined application context Output: H - Hessian matrix of the local function (wrt the four points of the rectangle only) */ static int MSurfLocalHessian(int coor[2], double x[4], double H[4][4],void *ptr) { AppCtx *user = (AppCtx*)ptr; double hx, hy, area, byhxhx, byhyhy; double dvdx, dvdy, flow, fup; double areadivf, areadivf3; hx = user->hx; hy = user->hy; area = user->area; byhxhx = 1.0 / (hx * hx); byhyhy = 1.0 / (hy * hy); /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */ dvdx = (x[0] - x[1]) / hx; /* lower triangle contribution */ dvdy = (x[0] - x[2]) / hy; flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); dvdx = dvdx / hx; dvdy = dvdy / hy; areadivf = area / flow; areadivf3 = areadivf / (flow * flow); H[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy); H[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx); H[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy); H[0][3] = 0.0; H[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx; H[1][2] = areadivf3 * (-dvdx) * dvdy; H[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy; /* upper triangle contribution */ dvdx = (x[3] - x[2]) / hx; dvdy = (x[3] - x[1]) / hy; fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); dvdx = dvdx / hx; dvdy = dvdy / hy; areadivf = area / fup; areadivf3 = areadivf / (fup * fup); H[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy; H[1][2] += areadivf3 * (-dvdy) * dvdx; H[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx); H[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy; H[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx; H[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy); H[1][0] = H[0][1]; H[2][0] = H[0][2]; H[3][0] = H[0][3]; H[2][1] = H[1][2]; H[3][1] = H[1][3]; H[3][2] = H[2][3]; return 0; } /* MSurfLocalHessian */ /*------- USER-DEFINED: routine to evaluate the function and gradient at the whole grid -------*/ #undef __FUNCT__ #define __FUNCT__ "WholeMSurfFunctionGradient" /* WholeMSurfFunctionGradient - Evaluates function and gradient over the whole grid Input: daapplication - TAO application object da - distributed array X - the current point, at which the function and gradient are evaluated ptr - user-defined application context Output: f - value of the objective funtion at X G - gradient at X */ static int WholeMSurfFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) { AppCtx *user = (AppCtx*)ptr; Vec localX, localG; int info, i, j; int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye; double **x, **g; double floc = 0.0; PetscScalar zero = 0.0; double hx, hy, area; double dvdx, dvdy, flow, fup; double areadivf; hx = user->hx; hy = user->hy; area = user->area; info = DAGetLocalVector(da, &localX); CHKERRQ(info); info = DAGetLocalVector(da, &localG); CHKERRQ(info); info = VecSet(G, zero); CHKERRQ(info); info = VecSet(localG, zero); CHKERRQ(info); info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info); info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info); info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info); info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info); info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info); info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info); xe = gxs + gxm - 1; ye = gys + gym - 1; for (j = ys; j < ye; j++) { for (i = xs; i < xe; i++) { /* lower triangle contribution */ dvdx = (x[j][i] - x[j][i+1]) / hx; dvdy = (x[j][i] - x[j+1][i]) / hy; flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); areadivf = area / flow; g[j][i] += (dvdx / hx + dvdy / hy) * areadivf; g[j][i+1] += (-dvdx / hx) * areadivf; g[j+1][i] += (-dvdy / hy) * areadivf; /* upper triangle contribution */ dvdx = (x[j+1][i+1] - x[j+1][i]) / hx; dvdy = (x[j+1][i+1] - x[j][i+1]) / hy; fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); areadivf = area / fup; g[j][i+1] += (-dvdy / hy) * areadivf; g[j+1][i] += (-dvdx / hx) * areadivf; g[j+1][i+1] += (dvdx / hx + dvdy / hy) * areadivf; floc += area * (flow + fup); } } info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info); info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info); info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info); info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info); info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info); info = DARestoreLocalVector(da, &localX); CHKERRQ(info); info = DARestoreLocalVector(da, &localG); CHKERRQ(info); info = PetscLogFlops((xe-xs) * (ye-ys) * 42); CHKERRQ(info); return 0; } /* WholeMSurfFunctionGradient */ /*------- USER-DEFINED: routine to evaluate the Hessian at the whole grid -------*/ #undef __FUNCT__ #define __FUNCT__ "WholeMSurfHessian" /* WholeMSurfHessian - Evaluates Hessian over the whole grid Input: daapplication - TAO application object da - distributed array X - the current point, at which the function and gradient are evaluated ptr - user-defined application context Output: H - Hessian at X */ static int WholeMSurfHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) { AppCtx *user = (AppCtx*)ptr; Vec localX; int info, i, j, ind[4]; int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye; double smallH[4][4]; double **x; double hx, hy, area, byhxhx, byhyhy; double dvdx, dvdy, flow, fup; double areadivf, areadivf3; PetscTruth assembled; hx = user->hx; hy = user->hy; area = user->area; byhxhx = 1.0 / (hx * hx); byhyhy = 1.0 / (hy * hy); info = DAGetLocalVector(da, &localX); CHKERRQ(info); info = MatAssembled(H,&assembled); CHKERRQ(info); if (assembled){info = MatZeroEntries(H); CHKERRQ(info);} info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info); info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info); info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info); info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info); info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info); xe = gxs + gxm - 1; ye = gys + gym - 1; for (j = ys; j < ye; j++) { for (i = xs; i < xe; i++) { /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */ dvdx = (x[j][i] - x[j][i+1]) / hx; /* lower triangle contribution */ dvdy = (x[j][i] - x[j+1][i]) / hy; flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); dvdx = dvdx / hx; dvdy = dvdy / hy; areadivf = area / flow; areadivf3 = areadivf / (flow * flow); smallH[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy); smallH[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx); smallH[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy); smallH[0][3] = 0.0; smallH[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx; smallH[1][2] = areadivf3 * (-dvdx) * dvdy; smallH[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy; /* upper triangle contribution */ dvdx = (x[j+1][i+1] - x[j+1][i]) / hx; dvdy = (x[j+1][i+1] - x[j][i+1]) / hy; fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy ); dvdx = dvdx / hx; dvdy = dvdy / hy; areadivf = area / fup; areadivf3 = areadivf / (fup * fup); smallH[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy; smallH[1][2] += areadivf3 * (-dvdy) * dvdx; smallH[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx); smallH[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy; smallH[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx; smallH[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy); smallH[1][0] = smallH[0][1]; smallH[2][0] = smallH[0][2]; smallH[3][0] = smallH[0][3]; smallH[2][1] = smallH[1][2]; smallH[3][1] = smallH[1][3]; smallH[3][2] = smallH[2][3]; ind[0] = (j-gys) * gxm + (i-gxs); ind[1] = ind[0] + 1; ind[2] = ind[0] + gxm; ind[3] = ind[2] + 1; info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info); } } info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info); info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info); info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info); info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info); info = DARestoreLocalVector(da, &localX); CHKERRQ(info); info = PetscLogFlops((xe-xs) * (ye-ys) * 83 + 4); CHKERRQ(info); return 0; } /* WholeMSurfHessian */