Actual source code: maros.c

petsc-dev 2014-07-21
Report Typos and Errors
  1: /* Program usage: mpirun -np 1 maros1 [-help] [all TAO options] */

  3: /* ----------------------------------------------------------------------
  4: TODO Explain maros example
  5: ---------------------------------------------------------------------- */

  7: #include <petsctao.h>

  9: static  char help[]="";

 11: /*T
 12:    Concepts: TAO^Solving an unconstrained minimization problem
 13:    Routines: TaoCreate(); TaoSetType();
 14:    Routines: TaoSetInitialVector();
 15:    Routines: TaoSetObjectiveAndGradientRoutine();
 16:    Routines: TaoSetEqualityConstraintsRoutine();
 17:    Routines: TaoSetInequalityConstraintsRoutine();
 18:    Routines: TaoSetEqualityJacobianRoutine();
 19:    Routines: TaoSetInequalityJacobianRoutine();
 20:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 21:    Routines: TaoGetKSP(); TaoSolve();
 22:    Routines: TaoGetConvergedReason(); TaoDestroy();
 23:    Processors: 1
 24: T*/

 26: /*
 27:    User-defined application context - contains data needed by the
 28:    application-provided call-back routines, FormFunction(),
 29:    FormGradient(), and FormHessian().
 30: */

 32: /*
 33:    x,d in R^n
 34:    f in R
 35:    bin in R^mi
 36:    beq in R^me
 37:    Aeq in R^(me x n)
 38:    Ain in R^(mi x n)
 39:    H in R^(n x n)
 40:    min f=(1/2)*x'*H*x + d'*x
 41:    s.t.  Aeq*x == beq
 42:          Ain*x >= bin
 43: */
 44: typedef struct {
 45:   char     name[32];
 46:   PetscInt n; /* Length x */
 47:   PetscInt me; /* number of equality constraints */
 48:   PetscInt mi; /* number of inequality constraints */
 49:   PetscInt m;  /* me+mi */
 50:   Mat      Aeq,Ain,H;
 51:   Vec      beq,bin,d;
 52: } AppCtx;

 54: /* -------- User-defined Routines --------- */

 56: PetscErrorCode InitializeProblem(AppCtx*);
 57: PetscErrorCode DestroyProblem(AppCtx *);
 58: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
 59: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
 60: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
 61: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
 62: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
 63: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);

 67: PetscErrorCode main(int argc,char **argv)
 68: {
 69:   PetscErrorCode     ierr;                /* used to check for functions returning nonzeros */
 70:   PetscMPIInt        size;
 71:   Vec                x;                   /* solution */
 72:   KSP                ksp;
 73:   PC                 pc;
 74:   Vec                ceq,cin;
 75:   PetscBool          flg;                 /* A return value when checking for use options */
 76:   Tao                tao;                 /* Tao solver context */
 77:   TaoConvergedReason reason;
 78:   AppCtx             user;                /* application context */

 80:   /* Initialize TAO,PETSc */
 81:   PetscInitialize(&argc,&argv,(char *)0,help);
 82:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 83:   /* Specify default parameters for the problem, check for command-line overrides */
 84:   PetscStrncpy(user.name,"HS21",8);
 85:   PetscOptionsGetString(NULL,"-cutername",user.name,24,&flg);

 87:   PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name);
 88:   InitializeProblem(&user);
 89:   VecDuplicate(user.d,&x);
 90:   VecDuplicate(user.beq,&ceq);
 91:   VecDuplicate(user.bin,&cin);
 92:   VecSet(x,1.0);

 94:   TaoCreate(PETSC_COMM_WORLD,&tao);
 95:   TaoSetType(tao,TAOIPM);
 96:   TaoSetInitialVector(tao,x);
 97:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);
 98:   TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user);
 99:   TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user);
100:   TaoSetInequalityBounds(tao,user.bin,NULL);
101:   TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user);
102:   TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user);
103:   TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
104:   TaoGetKSP(tao,&ksp);
105:   KSPGetPC(ksp,&pc);
106:   PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU);
107:   /* TODO -- why didn't that work? */
108:   if (size == 1) {
109:     PetscOptionsSetValue("-pc_factor_mat_solver_package","superlu");
110:   } else {
111:     PetscOptionsSetValue("-pc_factor_mat_solver_package","superlu_dist");
112:   }
113:   PCSetType(pc,PCLU);
114:   KSPSetType(ksp,KSPPREONLY);
115:   TaoSetTolerances(tao,1e-12,0,0,0,0);

117:   TaoSetFromOptions(tao);

119:   /* Solve */
120:   TaoSolve(tao);

122:   /* Analyze solution */
123:   TaoGetConvergedReason(tao,&reason);
124:   if (reason < 0) {
125:     PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");
126:   } else {
127:     PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);
128:   }

130:   DestroyProblem(&user);
131:   VecDestroy(&x);
132:   VecDestroy(&ceq);
133:   VecDestroy(&cin);
134:   TaoDestroy(&tao);

136:   PetscFinalize();
137:   return 0;
138: }

142: PetscErrorCode InitializeProblem(AppCtx *user)
143: {
145:   PetscViewer    loader;
146:   MPI_Comm       comm;
147:   PetscInt       nrows,ncols,i;
148:   PetscScalar    one=1.0;
149:   char           filebase[128];
150:   char           filename[128];

153:   comm = PETSC_COMM_WORLD;
154:   PetscStrncpy(filebase,user->name,128);
155:   PetscStrncat(filebase,"/",1);
156:   PetscStrncpy(filename,filebase,128);
157:   PetscStrncat(filename,"f",3);
158:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);

160:   VecCreate(comm,&user->d);
161:   VecLoad(user->d,loader);
162:   PetscViewerDestroy(&loader);
163:   VecGetSize(user->d,&nrows);
164:   VecSetFromOptions(user->d);
165:   user->n = nrows;

167:   PetscStrncpy(filename,filebase,128);
168:   PetscStrncat(filename,"H",3);
169:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);

171:   MatCreate(comm,&user->H);
172:   MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);
173:   MatLoad(user->H,loader);
174:   PetscViewerDestroy(&loader);
175:   MatGetSize(user->H,&nrows,&ncols);
176:   if (nrows != user->n) SETERRQ(comm,0,"H: nrows != n\n");
177:   if (ncols != user->n) SETERRQ(comm,0,"H: ncols != n\n");
178:   MatSetFromOptions(user->H);

180:   PetscStrncpy(filename,filebase,128);
181:   PetscStrncat(filename,"Aeq",3);
182:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
183:   if (ierr) {
184:     user->Aeq = NULL;
185:     user->me  = 0;
186:   } else {
187:     MatCreate(comm,&user->Aeq);
188:     MatLoad(user->Aeq,loader);
189:     PetscViewerDestroy(&loader);
190:     MatGetSize(user->Aeq,&nrows,&ncols);
191:     if (ncols != user->n) SETERRQ(comm,0,"Aeq ncols != H nrows\n");
192:     MatSetFromOptions(user->Aeq);
193:     user->me = nrows;
194:   }

196:   PetscStrncpy(filename,filebase,128);
197:   PetscStrncat(filename,"Beq",3);
198:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
199:   if (ierr) {
200:     user->beq = 0;
201:   } else {
202:     VecCreate(comm,&user->beq);
203:     VecLoad(user->beq,loader);
204:     PetscViewerDestroy(&loader);
205:     VecGetSize(user->beq,&nrows);
206:     if (nrows != user->me) SETERRQ(comm,0,"Aeq nrows != Beq n\n");
207:     VecSetFromOptions(user->beq);
208:   }

210:   user->mi = user->n;
211:   /* Ain = eye(n,n) */
212:   MatCreate(comm,&user->Ain);
213:   MatSetType(user->Ain,MATAIJ);
214:   MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);

216:   MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);
217:   MatSeqAIJSetPreallocation(user->Ain,1,NULL);

219:   for (i=0;i<user->mi;i++) {
220:     MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);
221:   }
222:   MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);
223:   MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);
224:   MatSetFromOptions(user->Ain);

226:   /* bin = [0,0 ... 0]' */
227:   VecCreate(comm,&user->bin);
228:   VecSetType(user->bin,VECMPI);
229:   VecSetSizes(user->bin,PETSC_DECIDE,user->mi);
230:   VecSet(user->bin,0.0);
231:   VecSetFromOptions(user->bin);
232:   user->m = user->me + user->mi;
233:   return(0);
234: }

238: PetscErrorCode DestroyProblem(AppCtx *user)
239: {

243:   MatDestroy(&user->H);
244:   MatDestroy(&user->Aeq);
245:   MatDestroy(&user->Ain);
246:   VecDestroy(&user->beq);
247:   VecDestroy(&user->bin);
248:   VecDestroy(&user->d);
249:   return(0);
250: }
253: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
254: {
255:   AppCtx         *user = (AppCtx*)ctx;
256:   PetscScalar    xtHx;

260:   MatMult(user->H,x,g);
261:   VecDot(x,g,&xtHx);
262:   VecDot(x,user->d,f);
263:   *f += 0.5*xtHx;
264:   VecAXPY(g,1.0,user->d);
265:   return(0);
266: }

270: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
271: {
273:   return(0);
274: }

278: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
279: {
280:   AppCtx         *user = (AppCtx*)ctx;

284:   MatMult(user->Ain,x,ci);
285:   return(0);
286: }

290: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx)
291: {
292:   AppCtx         *user = (AppCtx*)ctx;

296:   MatMult(user->Aeq,x,ce);
297:   VecAXPY(ce,-1.0,user->beq);
298:   return(0);
299: }

303: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre,  void *ctx)
304: {
306:   return(0);
307: }

311: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
312: {
314:   return(0);
315: }