xref: /petsc/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c (revision 58d68138c660dfb4e9f5b03334792cd4f2ffd7cc)
1 static char help[] = "TAO/Pounders Matlab Testing on the More'-Wild Benchmark Problems\n\
2 The interface calls:\n\
3     TestingInitialize.m to initialize the problem set\n\
4     ProblemInitialize.m to initialize each instance\n\
5     ProblemFinalize.m to store the performance data for the instance solved\n\
6     TestingFinalize.m to store the entire set of performance data\n\
7 \n\
8 TestingPlot.m is called outside of TAO/Pounders to produce a performance profile\n\
9 of the results compared to the Matlab fminsearch algorithm.\n";
10 
11 #include <petsctao.h>
12 #include <petscmatlab.h>
13 
14 typedef struct {
15   PetscMatlabEngine mengine;
16 
17   double delta; /* Initial trust region radius */
18 
19   int n;     /* Number of inputs */
20   int m;     /* Number of outputs */
21   int nfmax; /* Maximum function evaluations */
22   int npmax; /* Maximum interpolation points */
23 } AppCtx;
24 
25 static PetscErrorCode EvaluateResidual(Tao tao, Vec X, Vec F, void *ptr) {
26   AppCtx *user = (AppCtx *)ptr;
27 
28   PetscFunctionBegin;
29   PetscCall(PetscObjectSetName((PetscObject)X, "X"));
30   PetscCall(PetscMatlabEnginePut(user->mengine, (PetscObject)X));
31   PetscCall(PetscMatlabEngineEvaluate(user->mengine, "F = func(X);"));
32   PetscCall(PetscObjectSetName((PetscObject)F, "F"));
33   PetscCall(PetscMatlabEngineGet(user->mengine, (PetscObject)F));
34   PetscFunctionReturn(0);
35 }
36 
37 static PetscErrorCode EvaluateJacobian(Tao tao, Vec X, Mat J, Mat JPre, void *ptr) {
38   AppCtx *user = (AppCtx *)ptr;
39 
40   PetscFunctionBegin;
41   PetscCall(PetscObjectSetName((PetscObject)X, "X"));
42   PetscCall(PetscMatlabEnginePut(user->mengine, (PetscObject)X));
43   PetscCall(PetscMatlabEngineEvaluate(user->mengine, "J = jac(X);"));
44   PetscCall(PetscObjectSetName((PetscObject)J, "J"));
45   PetscCall(PetscMatlabEngineGet(user->mengine, (PetscObject)J));
46   PetscFunctionReturn(0);
47 }
48 
49 static PetscErrorCode TaoPounders(AppCtx *user) {
50   Tao  tao;
51   Vec  X, F;
52   Mat  J;
53   char buf[1024];
54 
55   PetscFunctionBegin;
56 
57   /* Set the values for the algorithm options we want to use */
58   sprintf(buf, "%d", user->npmax);
59   PetscCall(PetscOptionsSetValue(NULL, "-tao_pounders_npmax", buf));
60   sprintf(buf, "%5.4e", user->delta);
61   PetscCall(PetscOptionsSetValue(NULL, "-tao_pounders_delta", buf));
62 
63   /* Create the TAO objects and set the type */
64   PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
65 
66   /* Create starting point and initialize */
67   PetscCall(VecCreateSeq(PETSC_COMM_SELF, user->n, &X));
68   PetscCall(PetscObjectSetName((PetscObject)X, "X0"));
69   PetscCall(PetscMatlabEngineGet(user->mengine, (PetscObject)X));
70   PetscCall(TaoSetSolution(tao, X));
71 
72   /* Create residuals vector and set residual function */
73   PetscCall(VecCreateSeq(PETSC_COMM_SELF, user->m, &F));
74   PetscCall(PetscObjectSetName((PetscObject)F, "F"));
75   PetscCall(TaoSetResidualRoutine(tao, F, EvaluateResidual, (void *)user));
76 
77   /* Create Jacobian matrix and set residual Jacobian routine */
78   PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, user->m, user->n, user->n, NULL, &J));
79   PetscCall(PetscObjectSetName((PetscObject)J, "J"));
80   PetscCall(TaoSetJacobianResidualRoutine(tao, J, J, EvaluateJacobian, (void *)user));
81 
82   /* Solve the problem */
83   PetscCall(TaoSetType(tao, TAOPOUNDERS));
84   PetscCall(TaoSetMaximumFunctionEvaluations(tao, user->nfmax));
85   PetscCall(TaoSetFromOptions(tao));
86   PetscCall(TaoSolve(tao));
87 
88   /* Finish the problem */
89   PetscCall(MatDestroy(&J));
90   PetscCall(VecDestroy(&X));
91   PetscCall(VecDestroy(&F));
92   PetscCall(TaoDestroy(&tao));
93   PetscFunctionReturn(0);
94 }
95 
96 int main(int argc, char **argv) {
97   AppCtx      user;
98   PetscScalar tmp;
99   PetscInt    prob_id = 0;
100   PetscBool   flg, testall = PETSC_FALSE;
101   int         i, i0, imax;
102 
103   PetscFunctionBeginUser;
104   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
105   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_all", &testall, NULL));
106   PetscCall(PetscOptionsGetInt(NULL, NULL, "-prob_id", &prob_id, &flg));
107   if (!testall) {
108     if (!flg) {
109       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Problem number must be specified with -prob_id");
110     } else if ((prob_id < 1) || (prob_id > 53)) {
111       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Problem number must be between 1 and 53!");
112     } else {
113       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Running problem %d\n", prob_id));
114     }
115   } else {
116     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Running all problems\n"));
117   }
118 
119   PetscCall(PetscMatlabEngineCreate(PETSC_COMM_SELF, NULL, &user.mengine));
120   PetscCall(PetscMatlabEngineEvaluate(user.mengine, "TestingInitialize"));
121 
122   if (testall) {
123     i0   = 1;
124     imax = 53;
125   } else {
126     i0   = (int)prob_id;
127     imax = (int)prob_id;
128   }
129 
130   for (i = i0; i <= imax; ++i) {
131     PetscCall(PetscPrintf(PETSC_COMM_SELF, "%d\n", i));
132     PetscCall(PetscMatlabEngineEvaluate(user.mengine, "np = %d; ProblemInitialize", i));
133     PetscCall(PetscMatlabEngineGetArray(user.mengine, 1, 1, &tmp, "n"));
134     user.n = (int)tmp;
135     PetscCall(PetscMatlabEngineGetArray(user.mengine, 1, 1, &tmp, "m"));
136     user.m = (int)tmp;
137     PetscCall(PetscMatlabEngineGetArray(user.mengine, 1, 1, &tmp, "nfmax"));
138     user.nfmax = (int)tmp;
139     PetscCall(PetscMatlabEngineGetArray(user.mengine, 1, 1, &tmp, "npmax"));
140     user.npmax = (int)tmp;
141     PetscCall(PetscMatlabEngineGetArray(user.mengine, 1, 1, &tmp, "delta"));
142     user.delta = (double)tmp;
143 
144     /* Ignore return code for now -- do not stop testing on inf or nan errors */
145     PetscCall(TaoPounders(&user));
146 
147     PetscCall(PetscMatlabEngineEvaluate(user.mengine, "ProblemFinalize"));
148   }
149 
150   PetscCall(PetscMatlabEngineEvaluate(user.mengine, "TestingFinalize"));
151   PetscCall(PetscMatlabEngineDestroy(&user.mengine));
152   PetscCall(PetscFinalize());
153   return 0;
154 }
155 
156 /*TEST
157 
158    build:
159       requires: matlab_engine
160 
161    test:
162       localrunfiles: more_wild_probs TestingInitialize.m TestingFinalize.m ProblemInitialize.m ProblemFinalize.m
163       args: -tao_smonitor -prob_id 5
164 
165 TEST*/
166