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: /*
12: User-defined application context - contains data needed by the
13: application-provided call-back routines, FormFunction(),
14: FormGradient(), and FormHessian().
15: */
17: /*
18: x,d in R^n
19: f in R
20: bin in R^mi
21: beq in R^me
22: Aeq in R^(me x n)
23: Ain in R^(mi x n)
24: H in R^(n x n)
25: min f=(1/2)*x'*H*x + d'*x
26: s.t. Aeq*x == beq
27: Ain*x >= bin
28: */
29: typedef struct {
30: char name[32];
31: PetscInt n; /* Length x */
32: PetscInt me; /* number of equality constraints */
33: PetscInt mi; /* number of inequality constraints */
34: PetscInt m; /* me+mi */
35: Mat Aeq,Ain,H;
36: Vec beq,bin,d;
37: } AppCtx;
39: /* -------- User-defined Routines --------- */
41: PetscErrorCode InitializeProblem(AppCtx*);
42: PetscErrorCode DestroyProblem(AppCtx *);
43: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
44: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
45: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
46: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
47: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
48: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
50: PetscErrorCode main(int argc,char **argv)
51: {
52: PetscMPIInt size;
53: Vec x; /* solution */
54: KSP ksp;
55: PC pc;
56: Vec ceq,cin;
57: PetscBool flg; /* A return value when checking for use options */
58: Tao tao; /* Tao solver context */
59: TaoConvergedReason reason;
60: AppCtx user; /* application context */
62: /* Initialize TAO,PETSc */
63: PetscInitialize(&argc,&argv,(char *)0,help);
64: MPI_Comm_size(PETSC_COMM_WORLD,&size);
65: /* Specify default parameters for the problem, check for command-line overrides */
66: PetscStrncpy(user.name,"HS21",sizeof(user.name));
67: PetscOptionsGetString(NULL,NULL,"-cutername",user.name,sizeof(user.name),&flg);
69: PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name);
70: InitializeProblem(&user);
71: VecDuplicate(user.d,&x);
72: VecDuplicate(user.beq,&ceq);
73: VecDuplicate(user.bin,&cin);
74: VecSet(x,1.0);
76: TaoCreate(PETSC_COMM_WORLD,&tao);
77: TaoSetType(tao,TAOIPM);
78: TaoSetSolution(tao,x);
79: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user);
80: TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user);
81: TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user);
82: TaoSetInequalityBounds(tao,user.bin,NULL);
83: TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user);
84: TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user);
85: TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user);
86: TaoGetKSP(tao,&ksp);
87: KSPGetPC(ksp,&pc);
88: PCSetType(pc,PCLU);
89: /*
90: This algorithm produces matrices with zeros along the diagonal therefore we need to use
91: SuperLU which does partial pivoting
92: */
93: PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);
94: KSPSetType(ksp,KSPPREONLY);
95: TaoSetTolerances(tao,0,0,0);
97: TaoSetFromOptions(tao);
98: TaoSolve(tao);
99: TaoGetConvergedReason(tao,&reason);
100: if (reason < 0) {
101: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason]);
102: } else {
103: PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason]);
104: }
106: DestroyProblem(&user);
107: VecDestroy(&x);
108: VecDestroy(&ceq);
109: VecDestroy(&cin);
110: TaoDestroy(&tao);
112: PetscFinalize();
113: return 0;
114: }
116: PetscErrorCode InitializeProblem(AppCtx *user)
117: {
118: PetscViewer loader;
119: MPI_Comm comm;
120: PetscInt nrows,ncols,i;
121: PetscScalar one = 1.0;
122: char filebase[128];
123: char filename[128];
125: comm = PETSC_COMM_WORLD;
126: PetscStrncpy(filebase,user->name,sizeof(filebase));
127: PetscStrlcat(filebase,"/",sizeof(filebase));
128: PetscStrncpy(filename,filebase,sizeof(filename));
129: PetscStrlcat(filename,"f",sizeof(filename));
130: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
132: VecCreate(comm,&user->d);
133: VecLoad(user->d,loader);
134: PetscViewerDestroy(&loader);
135: VecGetSize(user->d,&nrows);
136: VecSetFromOptions(user->d);
137: user->n = nrows;
139: PetscStrncpy(filename,filebase,sizeof(filename));
140: PetscStrlcat(filename,"H",sizeof(filename));
141: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
143: MatCreate(comm,&user->H);
144: MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);
145: MatLoad(user->H,loader);
146: PetscViewerDestroy(&loader);
147: MatGetSize(user->H,&nrows,&ncols);
150: MatSetFromOptions(user->H);
152: PetscStrncpy(filename,filebase,sizeof(filename));
153: PetscStrlcat(filename,"Aeq",sizeof(filename));
154: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
155: MatCreate(comm,&user->Aeq);
156: MatLoad(user->Aeq,loader);
157: PetscViewerDestroy(&loader);
158: MatGetSize(user->Aeq,&nrows,&ncols);
160: MatSetFromOptions(user->Aeq);
161: user->me = nrows;
163: PetscStrncpy(filename,filebase,sizeof(filename));
164: PetscStrlcat(filename,"Beq",sizeof(filename));
165: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
166: VecCreate(comm,&user->beq);
167: VecLoad(user->beq,loader);
168: PetscViewerDestroy(&loader);
169: VecGetSize(user->beq,&nrows);
171: VecSetFromOptions(user->beq);
173: user->mi = user->n;
174: /* Ain = eye(n,n) */
175: MatCreate(comm,&user->Ain);
176: MatSetType(user->Ain,MATAIJ);
177: MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);
179: MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);
180: MatSeqAIJSetPreallocation(user->Ain,1,NULL);
182: for (i=0;i<user->mi;i++) MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);
183: MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);
184: MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);
185: MatSetFromOptions(user->Ain);
187: /* bin = [0,0 ... 0]' */
188: VecCreate(comm,&user->bin);
189: VecSetType(user->bin,VECMPI);
190: VecSetSizes(user->bin,PETSC_DECIDE,user->mi);
191: VecSet(user->bin,0.0);
192: VecSetFromOptions(user->bin);
193: user->m = user->me + user->mi;
194: return 0;
195: }
197: PetscErrorCode DestroyProblem(AppCtx *user)
198: {
199: MatDestroy(&user->H);
200: MatDestroy(&user->Aeq);
201: MatDestroy(&user->Ain);
202: VecDestroy(&user->beq);
203: VecDestroy(&user->bin);
204: VecDestroy(&user->d);
205: return 0;
206: }
208: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
209: {
210: AppCtx *user = (AppCtx*)ctx;
211: PetscScalar xtHx;
213: MatMult(user->H,x,g);
214: VecDot(x,g,&xtHx);
215: VecDot(x,user->d,f);
216: *f += 0.5*xtHx;
217: VecAXPY(g,1.0,user->d);
218: return 0;
219: }
221: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
222: {
223: return 0;
224: }
226: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
227: {
228: AppCtx *user = (AppCtx*)ctx;
230: MatMult(user->Ain,x,ci);
231: return 0;
232: }
234: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx)
235: {
236: AppCtx *user = (AppCtx*)ctx;
238: MatMult(user->Aeq,x,ce);
239: VecAXPY(ce,-1.0,user->beq);
240: return 0;
241: }
243: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
244: {
245: return 0;
246: }
248: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
249: {
250: return 0;
251: }
253: /*TEST
255: build:
256: requires: !complex
258: test:
259: requires: superlu
260: localrunfiles: HS21
262: TEST*/