Actual source code: maros.c
petsc-3.13.6 2020-09-29
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 Section 1.5 Writing Application Codes with PETSc context - contains data needed by the
28: Section 1.5 Writing Application Codes with PETSc-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: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
68: PetscMPIInt size;
69: Vec x; /* solution */
70: KSP ksp;
71: PC pc;
72: Vec ceq,cin;
73: PetscBool flg; /* A return value when checking for use options */
74: Tao tao; /* Tao solver context */
75: TaoConvergedReason reason;
76: AppCtx user; /* Section 1.5 Writing Application Codes with PETSc context */
78: /* Initialize TAO,PETSc */
79: PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
80: MPI_Comm_size(PETSC_COMM_WORLD,&size);
81: /* Specify default parameters for the problem, check for command-line overrides */
82: PetscStrncpy(user.name,"HS21",8);
83: PetscOptionsGetString(NULL,NULL,"-cutername",user.name,24,&flg);
85: PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name);
86: InitializeProblem(&user);
87: VecDuplicate(user.d,&x);
88: VecDuplicate(user.beq,&ceq);
89: VecDuplicate(user.bin,&cin);
90: VecSet(x,1.0);
92: TaoCreate(PETSC_COMM_WORLD,&tao);
93: TaoSetType(tao,TAOIPM);
94: TaoSetInitialVector(tao,x);
95: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);
96: TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user);
97: TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user);
98: TaoSetInequalityBounds(tao,user.bin,NULL);
99: TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user);
100: TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user);
101: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
102: TaoGetKSP(tao,&ksp);
103: KSPGetPC(ksp,&pc);
104: PCSetType(pc,PCLU);
105: /*
106: This algorithm produces matrices with zeros along the diagonal therefore we need to use
107: SuperLU which does partial pivoting
108: */
109: PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);
110: KSPSetType(ksp,KSPPREONLY);
111: TaoSetTolerances(tao,0,0,0);
113: TaoSetFromOptions(tao);
114: TaoSolve(tao);
115: TaoGetConvergedReason(tao,&reason);
116: if (reason < 0) {
117: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason]);
118: } else {
119: PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason]);
120: }
122: DestroyProblem(&user);
123: VecDestroy(&x);
124: VecDestroy(&ceq);
125: VecDestroy(&cin);
126: TaoDestroy(&tao);
128: PetscFinalize();
129: return ierr;
130: }
132: PetscErrorCode InitializeProblem(AppCtx *user)
133: {
135: PetscViewer loader;
136: MPI_Comm comm;
137: PetscInt nrows,ncols,i;
138: PetscScalar one=1.0;
139: char filebase[128];
140: char filename[128];
143: comm = PETSC_COMM_WORLD;
144: PetscStrncpy(filebase,user->name,sizeof(filebase));
145: PetscStrlcat(filebase,"/",sizeof(filebase));
146: PetscStrncpy(filename,filebase,sizeof(filename));
147: PetscStrlcat(filename,"f",sizeof(filename));
148: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
150: VecCreate(comm,&user->d);
151: VecLoad(user->d,loader);
152: PetscViewerDestroy(&loader);
153: VecGetSize(user->d,&nrows);
154: VecSetFromOptions(user->d);
155: user->n = nrows;
157: PetscStrncpy(filename,filebase,sizeof(filename));
158: PetscStrlcat(filename,"H",sizeof(filename));
159: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
161: MatCreate(comm,&user->H);
162: MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);
163: MatLoad(user->H,loader);
164: PetscViewerDestroy(&loader);
165: MatGetSize(user->H,&nrows,&ncols);
166: if (nrows != user->n) SETERRQ(comm,0,"H: nrows != n\n");
167: if (ncols != user->n) SETERRQ(comm,0,"H: ncols != n\n");
168: MatSetFromOptions(user->H);
170: PetscStrncpy(filename,filebase,sizeof(filename));
171: PetscStrlcat(filename,"Aeq",sizeof(filename));
172: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
173: if (ierr) {
174: user->Aeq = NULL;
175: user->me = 0;
176: } else {
177: MatCreate(comm,&user->Aeq);
178: MatLoad(user->Aeq,loader);
179: PetscViewerDestroy(&loader);
180: MatGetSize(user->Aeq,&nrows,&ncols);
181: if (ncols != user->n) SETERRQ(comm,0,"Aeq ncols != H nrows\n");
182: MatSetFromOptions(user->Aeq);
183: user->me = nrows;
184: }
186: PetscStrncpy(filename,filebase,sizeof(filename));
187: PetscStrlcat(filename,"Beq",sizeof(filename));
188: PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
189: if (ierr) {
190: user->beq = 0;
191: } else {
192: VecCreate(comm,&user->beq);
193: VecLoad(user->beq,loader);
194: PetscViewerDestroy(&loader);
195: VecGetSize(user->beq,&nrows);
196: if (nrows != user->me) SETERRQ(comm,0,"Aeq nrows != Beq n\n");
197: VecSetFromOptions(user->beq);
198: }
200: user->mi = user->n;
201: /* Ain = eye(n,n) */
202: MatCreate(comm,&user->Ain);
203: MatSetType(user->Ain,MATAIJ);
204: MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);
206: MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);
207: MatSeqAIJSetPreallocation(user->Ain,1,NULL);
209: for (i=0;i<user->mi;i++) {
210: MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);
211: }
212: MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);
213: MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);
214: MatSetFromOptions(user->Ain);
216: /* bin = [0,0 ... 0]' */
217: VecCreate(comm,&user->bin);
218: VecSetType(user->bin,VECMPI);
219: VecSetSizes(user->bin,PETSC_DECIDE,user->mi);
220: VecSet(user->bin,0.0);
221: VecSetFromOptions(user->bin);
222: user->m = user->me + user->mi;
223: return(0);
224: }
226: PetscErrorCode DestroyProblem(AppCtx *user)
227: {
231: MatDestroy(&user->H);
232: MatDestroy(&user->Aeq);
233: MatDestroy(&user->Ain);
234: VecDestroy(&user->beq);
235: VecDestroy(&user->bin);
236: VecDestroy(&user->d);
237: return(0);
238: }
239: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
240: {
241: AppCtx *user = (AppCtx*)ctx;
242: PetscScalar xtHx;
246: MatMult(user->H,x,g);
247: VecDot(x,g,&xtHx);
248: VecDot(x,user->d,f);
249: *f += 0.5*xtHx;
250: VecAXPY(g,1.0,user->d);
251: return(0);
252: }
254: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
255: {
257: return(0);
258: }
260: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
261: {
262: AppCtx *user = (AppCtx*)ctx;
266: MatMult(user->Ain,x,ci);
267: return(0);
268: }
270: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx)
271: {
272: AppCtx *user = (AppCtx*)ctx;
276: MatMult(user->Aeq,x,ce);
277: VecAXPY(ce,-1.0,user->beq);
278: return(0);
279: }
281: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
282: {
284: return(0);
285: }
287: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
288: {
290: return(0);
291: }
294: /*TEST
296: build:
297: requires: !complex
299: test:
300: requires: superlu
301: localrunfiles: HS21
303: TEST*/