#include "tao.h" /* Include "tao.h" so we can use the TAO solvers with PETSc support. */ static char help[] = "This example demonstrates use of the TAO package to\n\ solve an unconstrained system of equations. This example is based on a\n\ problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and\n\ boundary values along the edges of the domain, the objective is to find the\n\ surface with the minimal area that satisfies the boundary conditions.\n\ This application solves this problem using complimentarity -- We are actually\n\ solving the system (grad f)_i >= 0, if x_i == l_i \n\ (grad f)_i = 0, if l_i < x_i < u_i \n\ (grad f)_i <= 0, if x_i == u_i \n\ where f is the function to be minimized. \n\ \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\ -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n"; /*T Concepts: TAO - Solving a complementarity problem Routines: TaoInitialize(); TaoFinalize(); Routines: TaoCreate(); TaoDestroy(); Routines: TaoApplicationCreate(); TaoAppDestroy(); Routines: TaoSetConstraintRoutine(); TaoSetJacobianRoutine(); Routines: TaoGetVariableBounds(); TaoAppSetInitialSolutionVec(); Routines: TaoSetOptions(); Routines: TaoSolve(); Processors: 1 T*/ /* User-defined application context - contains data needed by the application-provided call-back routines, FormFunctionGradient(), FormHessian(). */ typedef struct { int mx, my; double *bottom, *top, *left, *right; } AppCtx; /* -------- User-defined Routines --------- */ static int MSA_BoundaryConditions(AppCtx *); static int MSA_InitialPoint(AppCtx *, Vec); int FormConstraints(TAO_APPLICATION, Vec, Vec, void *); int FormJacobian(TAO_APPLICATION, Vec, Mat *, Mat*, MatStructure*,void *); #undef __FUNCT__ #define __FUNCT__ "main" int main(int argc, char **argv) { int info; /* used to check for functions returning nonzeros */ Vec x; /* solution vector */ Vec c; /* Constraints function vector */ Vec xl,xu; /* Bounds on the variables */ PetscTruth flg; /* A return variable when checking for user options */ TAO_SOLVER tao; /* TAO_SOLVER solver context */ TAO_APPLICATION my_app; /* The PETSc application */ Mat J; /* Jacobian matrix */ int N; /* Number of elements in vector */ PetscScalar lb = -TAO_INFINITY; /* lower bound constant */ PetscScalar ub = TAO_INFINITY; /* upper bound constant */ AppCtx user; /* user-defined work context */ /* Initialize PETSc, TAO */ PetscInitialize(&argc, &argv, (char *)0, help ); TaoInitialize(&argc, &argv, (char *)0, help ); /* Specify default dimension of the problem */ user.mx = 4; user.my = 4; /* Check for any command line arguments that override defaults */ info = PetscOptionsGetInt(PETSC_NULL, "-mx", &user.mx, &flg); CHKERRQ(info); info = PetscOptionsGetInt(PETSC_NULL, "-my", &user.my, &flg); CHKERRQ(info); /* Calculate any derived values from parameters */ N = user.mx*user.my; PetscPrintf(PETSC_COMM_SELF,"\n---- Minimum Surface Area Problem -----\n"); PetscPrintf(PETSC_COMM_SELF,"mx:%d, my:%d\n", user.mx,user.my); /* Create appropriate vectors and matrices */ info = VecCreateSeq(MPI_COMM_SELF, N, &x); info = VecDuplicate(x, &c); CHKERRQ(info); info = MatCreateSeqAIJ(MPI_COMM_SELF, N, N, 7, PETSC_NULL, &J); CHKERRQ(info); /* The TAO code begins here */ /* Create TAO solver and set desired solution method */ info = TaoCreate(PETSC_COMM_SELF, "tao_ssils", &tao); CHKERRQ(info); info = TaoApplicationCreate(PETSC_COMM_SELF, &my_app); CHKERRQ(info); /* Set data structure */ info = TaoAppSetInitialSolutionVec(my_app, x); CHKERRQ(info); /* Set routines for constraints function and Jacobian evaluation */ info = TaoAppSetConstraintRoutine(my_app, FormConstraints, (void *)&user); CHKERRQ(info); info = TaoAppSetJacobianRoutine(my_app, FormJacobian, (void *)&user); CHKERRQ(info); info = TaoAppSetJacobianMat(my_app, J, J); CHKERRQ(info); /* Set the variable bounds */ info = MSA_BoundaryConditions(&user); CHKERRQ(info); /* Set initial solution guess */ info = MSA_InitialPoint(&user, x); CHKERRQ(info); /* Now that the PETSc application is all set, attach to tao context */ info = TaoSetupApplicationSolver(my_app,tao); CHKERRQ(info); /* Set Bounds on variables */ info = VecDuplicate(x, &xl); CHKERRQ(info); info = VecDuplicate(x, &xu); CHKERRQ(info); info = VecSet(xl, lb); CHKERRQ(info); info = VecSet(xu, ub); CHKERRQ(info); info = TaoAppSetVariableBounds(my_app,xl,xu); CHKERRQ(info); /* Check for any tao command line options */ info = TaoSetOptions(my_app,tao); CHKERRQ(info); /* Solve the application */ info = TaoSolveApplication(my_app,tao); CHKERRQ(info); /* Free Tao data structures */ info = TaoDestroy(tao); CHKERRQ(info); info = TaoAppDestroy(my_app); CHKERRQ(info); /* Free PETSc data structures */ info = VecDestroy(x); CHKERRQ(info); info = VecDestroy(xl); CHKERRQ(info); info = VecDestroy(xu); CHKERRQ(info); info = VecDestroy(c); CHKERRQ(info); info = MatDestroy(J); CHKERRQ(info); /* Free user-created data structures */ PetscFree(user.bottom); PetscFree(user.top); PetscFree(user.left); PetscFree(user.right); /* Finalize TAO and PETSc */ PetscFinalize(); TaoFinalize(); return 0; } /* -------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "FormConstraints" /* FormConstraints - Evaluates gradient of f. Input Parameters: . tao - the TAO_APPLICATION context . X - input vector . ptr - optional user-defined context, as set by TaoAppSetConstraintRoutine() Output Parameters: . G - vector containing the newly evaluated gradient */ int FormConstraints(TAO_APPLICATION tao, Vec X, Vec G, void *ptr){ AppCtx *user = (AppCtx *) ptr; int info,i,j,row; int mx=user->mx, my=user->my; double hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy; double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; double df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc; PetscScalar zero=0.0; PetscScalar *g, *x; /* Initialize vector to zero */ info = VecSet(G, zero); CHKERRQ(info); /* Get pointers to vector data */ info = VecGetArray(X, &x); CHKERRQ(info); info = VecGetArray(G, &g); CHKERRQ(info); /* Compute function over the locally owned part of the mesh */ for (j=0; j<my; j++){ for (i=0; i< mx; i++){ row= j*mx + i; xc = x[row]; xlt=xrb=xl=xr=xb=xt=xc; if (i==0){ /* left side */ xl= user->left[j+1]; xlt = user->left[j+2]; } else { xl = x[row-1]; } if (j==0){ /* bottom side */ xb=user->bottom[i+1]; xrb = user->bottom[i+2]; } else { xb = x[row-mx]; } if (i+1 == mx){ /* right side */ xr=user->right[j+1]; xrb = user->right[j]; } else { xr = x[row+1]; } if (j+1==0+my){ /* top side */ xt=user->top[i+1]; xlt = user->top[i]; }else { xt = x[row+mx]; } if (i>0 && j+1<my){ xlt = x[row-1+mx]; } if (j>0 && i+1<mx){ xrb = x[row+1-mx]; } d1 = (xc-xl); d2 = (xc-xr); d3 = (xc-xt); d4 = (xc-xb); d5 = (xr-xrb); d6 = (xrb-xb); d7 = (xlt-xl); d8 = (xt-xlt); df1dxc = d1*hydhx; df2dxc = ( d1*hydhx + d4*hxdhy ); df3dxc = d3*hxdhy; df4dxc = ( d2*hydhx + d3*hxdhy ); df5dxc = d2*hydhx; df6dxc = d4*hxdhy; d1 /= hx; d2 /= hx; d3 /= hy; d4 /= hy; d5 /= hy; d6 /= hx; d7 /= hy; d8 /= hx; f1 = sqrt( 1.0 + d1*d1 + d7*d7); f2 = sqrt( 1.0 + d1*d1 + d4*d4); f3 = sqrt( 1.0 + d3*d3 + d8*d8); f4 = sqrt( 1.0 + d3*d3 + d2*d2); f5 = sqrt( 1.0 + d2*d2 + d5*d5); f6 = sqrt( 1.0 + d4*d4 + d6*d6); df1dxc /= f1; df2dxc /= f2; df3dxc /= f3; df4dxc /= f4; df5dxc /= f5; df6dxc /= f6; g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc )/2.0; } } /* Restore vectors */ info = VecRestoreArray(X, &x); CHKERRQ(info); info = VecRestoreArray(G, &g); CHKERRQ(info); info = PetscLogFlops(67*mx*my); CHKERRQ(info); return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "FormJacobian" /* FormJacobian - Evaluates Jacobian matrix. Input Parameters: . tao - the TAO_APPLICATION context . X - input vector . ptr - optional user-defined context, as set by TaoSetJacobian() Output Parameters: . tH - Jacobian matrix */ int FormJacobian(TAO_APPLICATION tao, Vec X, Mat *tH, Mat* tHPre, MatStructure* flag, void *ptr) { AppCtx *user = (AppCtx *) ptr; Mat H = *tH; int i,j,k,row,info; int mx=user->mx, my=user->my; int col[7]; double hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy; double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; double hl,hr,ht,hb,hc,htl,hbr; PetscScalar *x, v[7]; PetscTruth assembled; /* Set various matrix options */ info = MatSetOption(H,MAT_IGNORE_OFF_PROC_ENTRIES); CHKERRQ(info); info = MatSetOption(H,MAT_COLUMNS_SORTED); CHKERRQ(info); info = MatSetOption(H,MAT_ROWS_SORTED); CHKERRQ(info); info = MatAssembled(H,&assembled); CHKERRQ(info); if (assembled){info = MatZeroEntries(H); CHKERRQ(info);} *flag=SAME_NONZERO_PATTERN; /* Get pointers to vector data */ info = VecGetArray(X, &x); CHKERRQ(info); /* Compute Jacobian over the locally owned part of the mesh */ for (i=0; i< mx; i++){ for (j=0; j<my; j++){ row= j*mx + i; xc = x[row]; xlt=xrb=xl=xr=xb=xt=xc; /* Left side */ if (i==0){ xl= user->left[j+1]; xlt = user->left[j+2]; } else { xl = x[row-1]; } if (j==0){ xb=user->bottom[i+1]; xrb = user->bottom[i+2]; } else { xb = x[row-mx]; } if (i+1 == mx){ xr=user->right[j+1]; xrb = user->right[j]; } else { xr = x[row+1]; } if (j+1==my){ xt=user->top[i+1]; xlt = user->top[i]; }else { xt = x[row+mx]; } if (i>0 && j+1<my){ xlt = x[row-1+mx]; } if (j>0 && i+1<mx){ xrb = x[row+1-mx]; } d1 = (xc-xl)/hx; d2 = (xc-xr)/hx; d3 = (xc-xt)/hy; d4 = (xc-xb)/hy; d5 = (xrb-xr)/hy; d6 = (xrb-xb)/hx; d7 = (xlt-xl)/hy; d8 = (xlt-xt)/hx; f1 = sqrt( 1.0 + d1*d1 + d7*d7); f2 = sqrt( 1.0 + d1*d1 + d4*d4); f3 = sqrt( 1.0 + d3*d3 + d8*d8); f4 = sqrt( 1.0 + d3*d3 + d2*d2); f5 = sqrt( 1.0 + d2*d2 + d5*d5); f6 = sqrt( 1.0 + d4*d4 + d6*d6); hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2); hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4); ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4); hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2); hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6); htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3); 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) + (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); hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0; k=0; if (j>0){ v[k]=hb; col[k]=row - mx; k++; } if (j>0 && i < mx -1){ v[k]=hbr; col[k]=row - mx+1; k++; } if (i>0){ v[k]= hl; col[k]=row - 1; k++; } v[k]= hc; col[k]=row; k++; if (i < mx-1 ){ v[k]= hr; col[k]=row+1; k++; } if (i>0 && j < my-1 ){ v[k]= htl; col[k] = row+mx-1; k++; } if (j < my-1 ){ v[k]= ht; col[k] = row+mx; k++; } /* Set matrix values using local numbering, which was defined earlier, in the main routine. */ info = MatSetValues(H,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info); } } /* Restore vectors */ info = VecRestoreArray(X,&x); CHKERRQ(info); /* Assemble the matrix */ info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info); info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info); info = PetscLogFlops(199*mx*my); CHKERRQ(info); return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "MSA_BoundaryConditions" /* MSA_BoundaryConditions - Calculates the boundary conditions for the region. Input Parameter: . user - user-defined application context Output Parameter: . user - user-defined application context */ static int MSA_BoundaryConditions(AppCtx * user) { int i,j,k,limit=0,info,maxits=5; int mx=user->mx,my=user->my; int bsize=0, lsize=0, tsize=0, rsize=0; 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; double *boundary; bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2; info = PetscMalloc(bsize*sizeof(double), &user->bottom);CHKERRQ(info); info = PetscMalloc(tsize*sizeof(double), &user->top);CHKERRQ(info); info = PetscMalloc(lsize*sizeof(double), &user->left);CHKERRQ(info); info = PetscMalloc(rsize*sizeof(double), &user->right);CHKERRQ(info); hx= (r-l)/(mx+1); hy=(t-b)/(my+1); for (j=0; j<4; j++){ if (j==0){ yt=b; xt=l; limit=bsize; boundary=user->bottom; } else if (j==1){ yt=t; xt=l; limit=tsize; boundary=user->top; } else if (j==2){ yt=b; xt=l; limit=lsize; boundary=user->left; } else { // if (j==3) yt=b; xt=r; limit=rsize; boundary=user->right; } 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; } boundary[i]=u1*u1-u2*u2; if (j==0 || j==1) { xt=xt+hx; } else { // if (j==2 || j==3) yt=yt+hy; } } } return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "MSA_InitialPoint" /* MSA_InitialPoint - Calculates the initial guess in one of three ways. Input Parameters: . user - user-defined application context . X - vector for initial guess Output Parameters: . X - newly computed initial guess */ static int MSA_InitialPoint(AppCtx * user, Vec X) { int start=-1,i,j,info; PetscScalar zero=0.0; PetscTruth flg; info = PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); CHKERRQ(info); if (flg && start==0){ /* The zero vector is reasonable */ info = VecSet(X, zero); CHKERRQ(info); /* PLogInfo(user,"Min. Surface Area Problem: Start with 0 vector \n"); */ } else { /* Take an average of the boundary conditions */ int row; int mx=user->mx,my=user->my; PetscScalar *x; /* Get pointers to vector data */ info = VecGetArray(X,&x); CHKERRQ(info); /* Perform local computations */ for (j=0; j<my; j++){ for (i=0; i< mx; i++){ row=(j)*mx + (i); 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; } } /* Restore vectors */ info = VecRestoreArray(X,&x); CHKERRQ(info); } return 0; }