Actual source code: maros.c
petsc-3.5.4 2015-05-23
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: }