Actual source code: maros.c
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: TaoSetSolution();
15: Routines: TaoSetObjectiveAndGradient();
16: Routines: TaoSetEqualityConstraintsRoutine();
17: Routines: TaoSetInequalityConstraintsRoutine();
18: Routines: TaoSetEqualityJacobianRoutine();
19: Routines: TaoSetInequalityJacobianRoutine();
20: Routines: TaoSetHessian(); 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*);
65: PetscErrorCode main(int argc,char **argv)
66: {
67: PetscMPIInt size;
68: Vec x; /* solution */
69: KSP ksp;
70: PC pc;
71: Vec ceq,cin;
72: PetscBool flg; /* A return value when checking for use options */
73: Tao tao; /* Tao solver context */
74: TaoConvergedReason reason;
75: AppCtx user; /* application context */
77: /* Initialize TAO,PETSc */
78: PetscInitialize(&argc,&argv,(char *)0,help);
79: MPI_Comm_size(PETSC_COMM_WORLD,&size);
80: /* Specify default parameters for the problem, check for command-line overrides */
81: PetscStrncpy(user.name,"HS21",sizeof(user.name));
82: PetscOptionsGetString(NULL,NULL,"-cutername",user.name,sizeof(user.name),&flg);
84: PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name);
85: InitializeProblem(&user);
86: VecDuplicate(user.d,&x);
87: VecDuplicate(user.beq,&ceq);
88: VecDuplicate(user.bin,&cin);
89: VecSet(x,1.0);
91: TaoCreate(PETSC_COMM_WORLD,&tao);
92: TaoSetType(tao,TAOIPM);
93: TaoSetSolution(tao,x);
94: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user);
95: TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user);
96: TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user);
97: TaoSetInequalityBounds(tao,user.bin,NULL);
98: TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user);
99: TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user);
100: TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user);
101: TaoGetKSP(tao,&ksp);
102: KSPGetPC(ksp,&pc);
103: PCSetType(pc,PCLU);
104: /*
105: This algorithm produces matrices with zeros along the diagonal therefore we need to use
106: SuperLU which does partial pivoting
107: */
108: PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);
109: KSPSetType(ksp,KSPPREONLY);
110: TaoSetTolerances(tao,0,0,0);
112: TaoSetFromOptions(tao);
113: TaoSolve(tao);
114: TaoGetConvergedReason(tao,&reason);
115: if (reason < 0) {
116: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason]);
117: } else {
118: PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason]);
119: }
121: DestroyProblem(&user);
122: VecDestroy(&x);
123: VecDestroy(&ceq);
124: VecDestroy(&cin);
125: TaoDestroy(&tao);
127: PetscFinalize();
128: return 0;
129: }
131: PetscErrorCode InitializeProblem(AppCtx *user)
132: {
133: PetscViewer loader;
134: MPI_Comm comm;
135: PetscInt nrows,ncols,i;
136: PetscScalar one = 1.0;
137: char filebase[128];
138: char filename[128];
140: comm = PETSC_COMM_WORLD;
141: PetscStrncpy(filebase,user->name,sizeof(filebase));
142: PetscStrlcat(filebase,"/",sizeof(filebase));
143: PetscStrncpy(filename,filebase,sizeof(filename));
144: PetscStrlcat(filename,"f",sizeof(filename));
145: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
147: VecCreate(comm,&user->d);
148: VecLoad(user->d,loader);
149: PetscViewerDestroy(&loader);
150: VecGetSize(user->d,&nrows);
151: VecSetFromOptions(user->d);
152: user->n = nrows;
154: PetscStrncpy(filename,filebase,sizeof(filename));
155: PetscStrlcat(filename,"H",sizeof(filename));
156: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
158: MatCreate(comm,&user->H);
159: MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);
160: MatLoad(user->H,loader);
161: PetscViewerDestroy(&loader);
162: MatGetSize(user->H,&nrows,&ncols);
165: MatSetFromOptions(user->H);
167: PetscStrncpy(filename,filebase,sizeof(filename));
168: PetscStrlcat(filename,"Aeq",sizeof(filename));
169: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
170: MatCreate(comm,&user->Aeq);
171: MatLoad(user->Aeq,loader);
172: PetscViewerDestroy(&loader);
173: MatGetSize(user->Aeq,&nrows,&ncols);
175: MatSetFromOptions(user->Aeq);
176: user->me = nrows;
178: PetscStrncpy(filename,filebase,sizeof(filename));
179: PetscStrlcat(filename,"Beq",sizeof(filename));
180: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
181: VecCreate(comm,&user->beq);
182: VecLoad(user->beq,loader);
183: PetscViewerDestroy(&loader);
184: VecGetSize(user->beq,&nrows);
186: VecSetFromOptions(user->beq);
188: user->mi = user->n;
189: /* Ain = eye(n,n) */
190: MatCreate(comm,&user->Ain);
191: MatSetType(user->Ain,MATAIJ);
192: MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);
194: MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);
195: MatSeqAIJSetPreallocation(user->Ain,1,NULL);
197: for (i=0;i<user->mi;i++) MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);
198: MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);
199: MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);
200: MatSetFromOptions(user->Ain);
202: /* bin = [0,0 ... 0]' */
203: VecCreate(comm,&user->bin);
204: VecSetType(user->bin,VECMPI);
205: VecSetSizes(user->bin,PETSC_DECIDE,user->mi);
206: VecSet(user->bin,0.0);
207: VecSetFromOptions(user->bin);
208: user->m = user->me + user->mi;
209: return 0;
210: }
212: PetscErrorCode DestroyProblem(AppCtx *user)
213: {
214: MatDestroy(&user->H);
215: MatDestroy(&user->Aeq);
216: MatDestroy(&user->Ain);
217: VecDestroy(&user->beq);
218: VecDestroy(&user->bin);
219: VecDestroy(&user->d);
220: return 0;
221: }
223: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
224: {
225: AppCtx *user = (AppCtx*)ctx;
226: PetscScalar xtHx;
228: MatMult(user->H,x,g);
229: VecDot(x,g,&xtHx);
230: VecDot(x,user->d,f);
231: *f += 0.5*xtHx;
232: VecAXPY(g,1.0,user->d);
233: return 0;
234: }
236: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
237: {
238: return 0;
239: }
241: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
242: {
243: AppCtx *user = (AppCtx*)ctx;
245: MatMult(user->Ain,x,ci);
246: return 0;
247: }
249: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx)
250: {
251: AppCtx *user = (AppCtx*)ctx;
253: MatMult(user->Aeq,x,ce);
254: VecAXPY(ce,-1.0,user->beq);
255: return 0;
256: }
258: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
259: {
260: return 0;
261: }
263: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
264: {
265: return 0;
266: }
268: /*TEST
270: build:
271: requires: !complex
273: test:
274: requires: superlu
275: localrunfiles: HS21
277: TEST*/