Actual source code: maros.c
petsc-3.7.5 2017-01-01
1: /* Program usage: mpiexec -n 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,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: PCSetType(pc,PCLU);
107: /*
108: This algorithm produces matrices with zeros along the diagonal therefore we need to use
109: SuperLU which does partial pivoting
110: */
111: PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU);
112: KSPSetType(ksp,KSPPREONLY);
113: TaoSetTolerances(tao,0,0,0);
115: TaoSetFromOptions(tao);
116: TaoSolve(tao);
117: TaoGetConvergedReason(tao,&reason);
118: if (reason < 0) {
119: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason]);
120: } else {
121: PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason]);
122: }
124: DestroyProblem(&user);
125: VecDestroy(&x);
126: VecDestroy(&ceq);
127: VecDestroy(&cin);
128: TaoDestroy(&tao);
130: PetscFinalize();
131: return 0;
132: }
136: PetscErrorCode InitializeProblem(AppCtx *user)
137: {
139: PetscViewer loader;
140: MPI_Comm comm;
141: PetscInt nrows,ncols,i;
142: PetscScalar one=1.0;
143: char filebase[128];
144: char filename[128];
147: comm = PETSC_COMM_WORLD;
148: PetscStrncpy(filebase,user->name,128);
149: PetscStrncat(filebase,"/",1);
150: PetscStrncpy(filename,filebase,128);
151: PetscStrncat(filename,"f",3);
152: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
154: VecCreate(comm,&user->d);
155: VecLoad(user->d,loader);
156: PetscViewerDestroy(&loader);
157: VecGetSize(user->d,&nrows);
158: VecSetFromOptions(user->d);
159: user->n = nrows;
161: PetscStrncpy(filename,filebase,128);
162: PetscStrncat(filename,"H",3);
163: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
165: MatCreate(comm,&user->H);
166: MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);
167: MatLoad(user->H,loader);
168: PetscViewerDestroy(&loader);
169: MatGetSize(user->H,&nrows,&ncols);
170: if (nrows != user->n) SETERRQ(comm,0,"H: nrows != n\n");
171: if (ncols != user->n) SETERRQ(comm,0,"H: ncols != n\n");
172: MatSetFromOptions(user->H);
174: PetscStrncpy(filename,filebase,128);
175: PetscStrncat(filename,"Aeq",3);
176: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
177: if (ierr) {
178: user->Aeq = NULL;
179: user->me = 0;
180: } else {
181: MatCreate(comm,&user->Aeq);
182: MatLoad(user->Aeq,loader);
183: PetscViewerDestroy(&loader);
184: MatGetSize(user->Aeq,&nrows,&ncols);
185: if (ncols != user->n) SETERRQ(comm,0,"Aeq ncols != H nrows\n");
186: MatSetFromOptions(user->Aeq);
187: user->me = nrows;
188: }
190: PetscStrncpy(filename,filebase,128);
191: PetscStrncat(filename,"Beq",3);
192: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
193: if (ierr) {
194: user->beq = 0;
195: } else {
196: VecCreate(comm,&user->beq);
197: VecLoad(user->beq,loader);
198: PetscViewerDestroy(&loader);
199: VecGetSize(user->beq,&nrows);
200: if (nrows != user->me) SETERRQ(comm,0,"Aeq nrows != Beq n\n");
201: VecSetFromOptions(user->beq);
202: }
204: user->mi = user->n;
205: /* Ain = eye(n,n) */
206: MatCreate(comm,&user->Ain);
207: MatSetType(user->Ain,MATAIJ);
208: MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);
210: MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);
211: MatSeqAIJSetPreallocation(user->Ain,1,NULL);
213: for (i=0;i<user->mi;i++) {
214: MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);
215: }
216: MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);
217: MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);
218: MatSetFromOptions(user->Ain);
220: /* bin = [0,0 ... 0]' */
221: VecCreate(comm,&user->bin);
222: VecSetType(user->bin,VECMPI);
223: VecSetSizes(user->bin,PETSC_DECIDE,user->mi);
224: VecSet(user->bin,0.0);
225: VecSetFromOptions(user->bin);
226: user->m = user->me + user->mi;
227: return(0);
228: }
232: PetscErrorCode DestroyProblem(AppCtx *user)
233: {
237: MatDestroy(&user->H);
238: MatDestroy(&user->Aeq);
239: MatDestroy(&user->Ain);
240: VecDestroy(&user->beq);
241: VecDestroy(&user->bin);
242: VecDestroy(&user->d);
243: return(0);
244: }
247: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
248: {
249: AppCtx *user = (AppCtx*)ctx;
250: PetscScalar xtHx;
254: MatMult(user->H,x,g);
255: VecDot(x,g,&xtHx);
256: VecDot(x,user->d,f);
257: *f += 0.5*xtHx;
258: VecAXPY(g,1.0,user->d);
259: return(0);
260: }
264: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
265: {
267: return(0);
268: }
272: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
273: {
274: AppCtx *user = (AppCtx*)ctx;
278: MatMult(user->Ain,x,ci);
279: return(0);
280: }
284: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx)
285: {
286: AppCtx *user = (AppCtx*)ctx;
290: MatMult(user->Aeq,x,ce);
291: VecAXPY(ce,-1.0,user->beq);
292: return(0);
293: }
297: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
298: {
300: return(0);
301: }
305: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
306: {
308: return(0);
309: }