Actual source code: ipm.c
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/constrained/impls/ipm/ipm.h>
4: /*
5: x,d in R^n
6: f in R
7: nb = mi + nlb+nub
8: s in R^nb is slack vector CI(x) / x-XL / -x+XU
9: bin in R^mi (tao->constraints_inequality)
10: beq in R^me (tao->constraints_equality)
11: lambdai in R^nb (ipmP->lambdai)
12: lambdae in R^me (ipmP->lambdae)
13: Jeq in R^(me x n) (tao->jacobian_equality)
14: Jin in R^(mi x n) (tao->jacobian_inequality)
15: Ai in R^(nb x n) (ipmP->Ai)
16: H in R^(n x n) (tao->hessian)
17: min f=(1/2)*x'*H*x + d'*x
18: s.t. CE(x) == 0
19: CI(x) >= 0
20: x >= tao->XL
21: -x >= -tao->XU
22: */
24: static PetscErrorCode IPMComputeKKT(Tao tao);
25: static PetscErrorCode IPMPushInitialPoint(Tao tao);
26: static PetscErrorCode IPMEvaluate(Tao tao);
27: static PetscErrorCode IPMUpdateK(Tao tao);
28: static PetscErrorCode IPMUpdateAi(Tao tao);
29: static PetscErrorCode IPMGatherRHS(Tao tao, Vec, Vec, Vec, Vec, Vec);
30: static PetscErrorCode IPMScatterStep(Tao tao, Vec, Vec, Vec, Vec, Vec);
31: static PetscErrorCode IPMInitializeBounds(Tao tao);
33: static PetscErrorCode TaoSolve_IPM(Tao tao)
34: {
35: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
36: PetscInt its, i;
37: PetscScalar stepsize = 1.0;
38: PetscScalar step_s, step_l, alpha, tau, sigma, phi_target;
40: PetscFunctionBegin;
41: /* Push initial point away from bounds */
42: PetscCall(IPMInitializeBounds(tao));
43: PetscCall(IPMPushInitialPoint(tao));
44: PetscCall(VecCopy(tao->solution, ipmP->rhs_x));
45: PetscCall(IPMEvaluate(tao));
46: PetscCall(IPMComputeKKT(tao));
48: tao->reason = TAO_CONTINUE_ITERATING;
49: PetscCall(TaoLogConvergenceHistory(tao, ipmP->kkt_f, ipmP->phi, 0.0, tao->ksp_its));
50: PetscCall(TaoMonitor(tao, tao->niter, ipmP->kkt_f, ipmP->phi, 0.0, 1.0));
51: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
53: while (tao->reason == TAO_CONTINUE_ITERATING) {
54: /* Call general purpose update function */
55: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
57: tao->ksp_its = 0;
58: PetscCall(IPMUpdateK(tao));
59: /*
60: rhs.x = -rd
61: rhs.lame = -rpe
62: rhs.lami = -rpi
63: rhs.com = -com
64: */
66: PetscCall(VecCopy(ipmP->rd, ipmP->rhs_x));
67: if (ipmP->me > 0) PetscCall(VecCopy(ipmP->rpe, ipmP->rhs_lambdae));
68: if (ipmP->nb > 0) {
69: PetscCall(VecCopy(ipmP->rpi, ipmP->rhs_lambdai));
70: PetscCall(VecCopy(ipmP->complementarity, ipmP->rhs_s));
71: }
72: PetscCall(IPMGatherRHS(tao, ipmP->bigrhs, ipmP->rhs_x, ipmP->rhs_lambdae, ipmP->rhs_lambdai, ipmP->rhs_s));
73: PetscCall(VecScale(ipmP->bigrhs, -1.0));
75: /* solve K * step = rhs */
76: PetscCall(KSPSetOperators(tao->ksp, ipmP->K, ipmP->K));
77: PetscCall(KSPSolve(tao->ksp, ipmP->bigrhs, ipmP->bigstep));
79: PetscCall(IPMScatterStep(tao, ipmP->bigstep, tao->stepdirection, ipmP->ds, ipmP->dlambdae, ipmP->dlambdai));
80: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
81: tao->ksp_its += its;
82: tao->ksp_tot_its += its;
83: /* Find distance along step direction to closest bound */
84: if (ipmP->nb > 0) {
85: PetscCall(VecStepBoundInfo(ipmP->s, ipmP->ds, ipmP->Zero_nb, ipmP->Inf_nb, &step_s, NULL, NULL));
86: PetscCall(VecStepBoundInfo(ipmP->lambdai, ipmP->dlambdai, ipmP->Zero_nb, ipmP->Inf_nb, &step_l, NULL, NULL));
87: alpha = PetscMin(step_s, step_l);
88: alpha = PetscMin(alpha, 1.0);
89: ipmP->alpha1 = alpha;
90: } else {
91: ipmP->alpha1 = alpha = 1.0;
92: }
94: /* x_aff = x + alpha*d */
95: PetscCall(VecCopy(tao->solution, ipmP->save_x));
96: if (ipmP->me > 0) PetscCall(VecCopy(ipmP->lambdae, ipmP->save_lambdae));
97: if (ipmP->nb > 0) {
98: PetscCall(VecCopy(ipmP->lambdai, ipmP->save_lambdai));
99: PetscCall(VecCopy(ipmP->s, ipmP->save_s));
100: }
102: PetscCall(VecAXPY(tao->solution, alpha, tao->stepdirection));
103: if (ipmP->me > 0) PetscCall(VecAXPY(ipmP->lambdae, alpha, ipmP->dlambdae));
104: if (ipmP->nb > 0) {
105: PetscCall(VecAXPY(ipmP->lambdai, alpha, ipmP->dlambdai));
106: PetscCall(VecAXPY(ipmP->s, alpha, ipmP->ds));
107: }
109: /* Recompute kkt to find centering parameter sigma = (new_mu/old_mu)^3 */
110: if (ipmP->mu == 0.0) {
111: sigma = 0.0;
112: } else {
113: sigma = 1.0 / ipmP->mu;
114: }
115: PetscCall(IPMComputeKKT(tao));
116: sigma *= ipmP->mu;
117: sigma *= sigma * sigma;
119: /* revert kkt info */
120: PetscCall(VecCopy(ipmP->save_x, tao->solution));
121: if (ipmP->me > 0) PetscCall(VecCopy(ipmP->save_lambdae, ipmP->lambdae));
122: if (ipmP->nb > 0) {
123: PetscCall(VecCopy(ipmP->save_lambdai, ipmP->lambdai));
124: PetscCall(VecCopy(ipmP->save_s, ipmP->s));
125: }
126: PetscCall(IPMComputeKKT(tao));
128: /* update rhs with new complementarity vector */
129: if (ipmP->nb > 0) {
130: PetscCall(VecCopy(ipmP->complementarity, ipmP->rhs_s));
131: PetscCall(VecScale(ipmP->rhs_s, -1.0));
132: PetscCall(VecShift(ipmP->rhs_s, sigma * ipmP->mu));
133: }
134: PetscCall(IPMGatherRHS(tao, ipmP->bigrhs, NULL, NULL, NULL, ipmP->rhs_s));
136: /* solve K * step = rhs */
137: PetscCall(KSPSetOperators(tao->ksp, ipmP->K, ipmP->K));
138: PetscCall(KSPSolve(tao->ksp, ipmP->bigrhs, ipmP->bigstep));
140: PetscCall(IPMScatterStep(tao, ipmP->bigstep, tao->stepdirection, ipmP->ds, ipmP->dlambdae, ipmP->dlambdai));
141: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
142: tao->ksp_its += its;
143: tao->ksp_tot_its += its;
144: if (ipmP->nb > 0) {
145: /* Get max step size and apply frac-to-boundary */
146: tau = PetscMax(ipmP->taumin, 1.0 - ipmP->mu);
147: tau = PetscMin(tau, 1.0);
148: if (tau != 1.0) {
149: PetscCall(VecScale(ipmP->s, tau));
150: PetscCall(VecScale(ipmP->lambdai, tau));
151: }
152: PetscCall(VecStepBoundInfo(ipmP->s, ipmP->ds, ipmP->Zero_nb, ipmP->Inf_nb, &step_s, NULL, NULL));
153: PetscCall(VecStepBoundInfo(ipmP->lambdai, ipmP->dlambdai, ipmP->Zero_nb, ipmP->Inf_nb, &step_l, NULL, NULL));
154: if (tau != 1.0) {
155: PetscCall(VecCopy(ipmP->save_s, ipmP->s));
156: PetscCall(VecCopy(ipmP->save_lambdai, ipmP->lambdai));
157: }
158: alpha = PetscMin(step_s, step_l);
159: alpha = PetscMin(alpha, 1.0);
160: } else {
161: alpha = 1.0;
162: }
163: ipmP->alpha2 = alpha;
164: /* TODO make phi_target meaningful */
165: phi_target = ipmP->dec * ipmP->phi;
166: for (i = 0; i < 11; i++) {
167: PetscCall(VecAXPY(tao->solution, alpha, tao->stepdirection));
168: if (ipmP->nb > 0) {
169: PetscCall(VecAXPY(ipmP->s, alpha, ipmP->ds));
170: PetscCall(VecAXPY(ipmP->lambdai, alpha, ipmP->dlambdai));
171: }
172: if (ipmP->me > 0) PetscCall(VecAXPY(ipmP->lambdae, alpha, ipmP->dlambdae));
174: /* update dual variables */
175: if (ipmP->me > 0) PetscCall(VecCopy(ipmP->lambdae, tao->DE));
177: PetscCall(IPMEvaluate(tao));
178: PetscCall(IPMComputeKKT(tao));
179: if (ipmP->phi <= phi_target) break;
180: alpha /= 2.0;
181: }
183: PetscCall(TaoLogConvergenceHistory(tao, ipmP->kkt_f, ipmP->phi, 0.0, tao->ksp_its));
184: PetscCall(TaoMonitor(tao, tao->niter, ipmP->kkt_f, ipmP->phi, 0.0, stepsize));
185: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
186: tao->niter++;
187: }
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: static PetscErrorCode TaoSetup_IPM(Tao tao)
192: {
193: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
195: PetscFunctionBegin;
196: ipmP->nb = ipmP->mi = ipmP->me = 0;
197: ipmP->K = NULL;
198: PetscCall(VecGetSize(tao->solution, &ipmP->n));
199: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
200: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
201: if (!ipmP->rd) PetscCall(VecDuplicate(tao->solution, &ipmP->rd));
202: if (!ipmP->rhs_x) PetscCall(VecDuplicate(tao->solution, &ipmP->rhs_x));
203: if (!ipmP->work) PetscCall(VecDuplicate(tao->solution, &ipmP->work));
204: if (!ipmP->save_x) PetscCall(VecDuplicate(tao->solution, &ipmP->save_x));
205: if (tao->constraints_equality) {
206: PetscCall(VecGetSize(tao->constraints_equality, &ipmP->me));
207: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->lambdae));
208: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->dlambdae));
209: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->rhs_lambdae));
210: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->save_lambdae));
211: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->rpe));
212: PetscCall(VecDuplicate(tao->constraints_equality, &tao->DE));
213: }
214: if (tao->constraints_inequality) PetscCall(VecDuplicate(tao->constraints_inequality, &tao->DI));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: static PetscErrorCode IPMInitializeBounds(Tao tao)
219: {
220: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
221: Vec xtmp;
222: PetscInt xstart, xend;
223: PetscInt ucstart, ucend; /* user ci */
224: PetscInt ucestart, uceend; /* user ce */
225: PetscInt sstart = 0, send = 0;
226: PetscInt bigsize;
227: PetscInt i, counter, nloc;
228: PetscInt *cind, *xind, *ucind, *uceind, *stepind;
229: VecType vtype;
230: const PetscInt *xli, *xui;
231: PetscInt xl_offset, xu_offset;
232: IS bigxl, bigxu, isuc, isc, isx, sis, is1;
233: MPI_Comm comm;
235: PetscFunctionBegin;
236: cind = xind = ucind = uceind = stepind = NULL;
237: ipmP->mi = 0;
238: ipmP->nxlb = 0;
239: ipmP->nxub = 0;
240: ipmP->nb = 0;
241: ipmP->nslack = 0;
243: PetscCall(VecDuplicate(tao->solution, &xtmp));
244: PetscCall(TaoComputeVariableBounds(tao));
245: if (tao->XL) {
246: PetscCall(VecSet(xtmp, PETSC_NINFINITY));
247: PetscCall(VecWhichGreaterThan(tao->XL, xtmp, &ipmP->isxl));
248: PetscCall(ISGetSize(ipmP->isxl, &ipmP->nxlb));
249: } else {
250: ipmP->nxlb = 0;
251: }
252: if (tao->XU) {
253: PetscCall(VecSet(xtmp, PETSC_INFINITY));
254: PetscCall(VecWhichLessThan(tao->XU, xtmp, &ipmP->isxu));
255: PetscCall(ISGetSize(ipmP->isxu, &ipmP->nxub));
256: } else {
257: ipmP->nxub = 0;
258: }
259: PetscCall(VecDestroy(&xtmp));
260: if (tao->constraints_inequality) {
261: PetscCall(VecGetSize(tao->constraints_inequality, &ipmP->mi));
262: } else {
263: ipmP->mi = 0;
264: }
265: ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi;
267: PetscCall(PetscObjectGetComm((PetscObject)tao->solution, &comm));
269: bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
270: PetscCall(PetscMalloc1(bigsize, &stepind));
271: PetscCall(PetscMalloc1(ipmP->n, &xind));
272: PetscCall(PetscMalloc1(ipmP->me, &uceind));
273: PetscCall(VecGetOwnershipRange(tao->solution, &xstart, &xend));
275: if (ipmP->nb > 0) {
276: PetscCall(VecCreate(comm, &ipmP->s));
277: PetscCall(VecSetSizes(ipmP->s, PETSC_DECIDE, ipmP->nb));
278: PetscCall(VecSetFromOptions(ipmP->s));
279: PetscCall(VecDuplicate(ipmP->s, &ipmP->ds));
280: PetscCall(VecDuplicate(ipmP->s, &ipmP->rhs_s));
281: PetscCall(VecDuplicate(ipmP->s, &ipmP->complementarity));
282: PetscCall(VecDuplicate(ipmP->s, &ipmP->ci));
284: PetscCall(VecDuplicate(ipmP->s, &ipmP->lambdai));
285: PetscCall(VecDuplicate(ipmP->s, &ipmP->dlambdai));
286: PetscCall(VecDuplicate(ipmP->s, &ipmP->rhs_lambdai));
287: PetscCall(VecDuplicate(ipmP->s, &ipmP->save_lambdai));
289: PetscCall(VecDuplicate(ipmP->s, &ipmP->save_s));
290: PetscCall(VecDuplicate(ipmP->s, &ipmP->rpi));
291: PetscCall(VecDuplicate(ipmP->s, &ipmP->Zero_nb));
292: PetscCall(VecDuplicate(ipmP->s, &ipmP->One_nb));
293: PetscCall(VecSet(ipmP->One_nb, 1.0));
294: PetscCall(VecDuplicate(ipmP->s, &ipmP->Inf_nb));
295: PetscCall(VecSet(ipmP->Inf_nb, PETSC_INFINITY));
297: PetscCall(PetscMalloc1(ipmP->nb, &cind));
298: PetscCall(PetscMalloc1(ipmP->mi, &ucind));
299: PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send));
301: if (ipmP->mi > 0) {
302: PetscCall(VecGetOwnershipRange(tao->constraints_inequality, &ucstart, &ucend));
303: counter = 0;
304: for (i = ucstart; i < ucend; i++) cind[counter++] = i;
305: PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isuc));
306: PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isc));
307: PetscCall(VecScatterCreate(tao->constraints_inequality, isuc, ipmP->ci, isc, &ipmP->ci_scat));
309: PetscCall(ISDestroy(&isuc));
310: PetscCall(ISDestroy(&isc));
311: }
312: /* need to know how may xbound indices are on each process */
313: /* TODO better way */
314: if (ipmP->nxlb) {
315: PetscCall(ISAllGather(ipmP->isxl, &bigxl));
316: PetscCall(ISGetIndices(bigxl, &xli));
317: /* find offsets for this processor */
318: xl_offset = ipmP->mi;
319: for (i = 0; i < ipmP->nxlb; i++) {
320: if (xli[i] < xstart) {
321: xl_offset++;
322: } else break;
323: }
324: PetscCall(ISRestoreIndices(bigxl, &xli));
326: PetscCall(ISGetIndices(ipmP->isxl, &xli));
327: PetscCall(ISGetLocalSize(ipmP->isxl, &nloc));
328: for (i = 0; i < nloc; i++) {
329: xind[i] = xli[i];
330: cind[i] = xl_offset + i;
331: }
333: PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx));
334: PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc));
335: PetscCall(VecScatterCreate(tao->XL, isx, ipmP->ci, isc, &ipmP->xl_scat));
336: PetscCall(ISDestroy(&isx));
337: PetscCall(ISDestroy(&isc));
338: PetscCall(ISDestroy(&bigxl));
339: }
341: if (ipmP->nxub) {
342: PetscCall(ISAllGather(ipmP->isxu, &bigxu));
343: PetscCall(ISGetIndices(bigxu, &xui));
344: /* find offsets for this processor */
345: xu_offset = ipmP->mi + ipmP->nxlb;
346: for (i = 0; i < ipmP->nxub; i++) {
347: if (xui[i] < xstart) {
348: xu_offset++;
349: } else break;
350: }
351: PetscCall(ISRestoreIndices(bigxu, &xui));
353: PetscCall(ISGetIndices(ipmP->isxu, &xui));
354: PetscCall(ISGetLocalSize(ipmP->isxu, &nloc));
355: for (i = 0; i < nloc; i++) {
356: xind[i] = xui[i];
357: cind[i] = xu_offset + i;
358: }
360: PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx));
361: PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc));
362: PetscCall(VecScatterCreate(tao->XU, isx, ipmP->ci, isc, &ipmP->xu_scat));
363: PetscCall(ISDestroy(&isx));
364: PetscCall(ISDestroy(&isc));
365: PetscCall(ISDestroy(&bigxu));
366: }
367: }
368: PetscCall(VecCreate(comm, &ipmP->bigrhs));
369: PetscCall(VecGetType(tao->solution, &vtype));
370: PetscCall(VecSetType(ipmP->bigrhs, vtype));
371: PetscCall(VecSetSizes(ipmP->bigrhs, PETSC_DECIDE, bigsize));
372: PetscCall(VecSetFromOptions(ipmP->bigrhs));
373: PetscCall(VecDuplicate(ipmP->bigrhs, &ipmP->bigstep));
375: /* create scatters for step->x and x->rhs */
376: for (i = xstart; i < xend; i++) {
377: stepind[i - xstart] = i;
378: xind[i - xstart] = i;
379: }
380: PetscCall(ISCreateGeneral(comm, xend - xstart, stepind, PETSC_COPY_VALUES, &sis));
381: PetscCall(ISCreateGeneral(comm, xend - xstart, xind, PETSC_COPY_VALUES, &is1));
382: PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->solution, is1, &ipmP->step1));
383: PetscCall(VecScatterCreate(tao->solution, is1, ipmP->bigrhs, sis, &ipmP->rhs1));
384: PetscCall(ISDestroy(&sis));
385: PetscCall(ISDestroy(&is1));
387: if (ipmP->nb > 0) {
388: for (i = sstart; i < send; i++) {
389: stepind[i - sstart] = i + ipmP->n;
390: cind[i - sstart] = i;
391: }
392: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
393: PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1));
394: PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step2));
395: PetscCall(ISDestroy(&sis));
397: for (i = sstart; i < send; i++) {
398: stepind[i - sstart] = i + ipmP->n + ipmP->me;
399: cind[i - sstart] = i;
400: }
401: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
402: PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs3));
403: PetscCall(ISDestroy(&sis));
404: PetscCall(ISDestroy(&is1));
405: }
407: if (ipmP->me > 0) {
408: PetscCall(VecGetOwnershipRange(tao->constraints_equality, &ucestart, &uceend));
409: for (i = ucestart; i < uceend; i++) {
410: stepind[i - ucestart] = i + ipmP->n + ipmP->nb;
411: uceind[i - ucestart] = i;
412: }
414: PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis));
415: PetscCall(ISCreateGeneral(comm, uceend - ucestart, uceind, PETSC_COPY_VALUES, &is1));
416: PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->constraints_equality, is1, &ipmP->step3));
417: PetscCall(ISDestroy(&sis));
419: for (i = ucestart; i < uceend; i++) stepind[i - ucestart] = i + ipmP->n;
421: PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis));
422: PetscCall(VecScatterCreate(tao->constraints_equality, is1, ipmP->bigrhs, sis, &ipmP->rhs2));
423: PetscCall(ISDestroy(&sis));
424: PetscCall(ISDestroy(&is1));
425: }
427: if (ipmP->nb > 0) {
428: for (i = sstart; i < send; i++) {
429: stepind[i - sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
430: cind[i - sstart] = i;
431: }
432: PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1));
433: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
434: PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step4));
435: PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs4));
436: PetscCall(ISDestroy(&sis));
437: PetscCall(ISDestroy(&is1));
438: }
440: PetscCall(PetscFree(stepind));
441: PetscCall(PetscFree(cind));
442: PetscCall(PetscFree(ucind));
443: PetscCall(PetscFree(uceind));
444: PetscCall(PetscFree(xind));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode TaoDestroy_IPM(Tao tao)
449: {
450: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
452: PetscFunctionBegin;
453: PetscCall(VecDestroy(&ipmP->rd));
454: PetscCall(VecDestroy(&ipmP->rpe));
455: PetscCall(VecDestroy(&ipmP->rpi));
456: PetscCall(VecDestroy(&ipmP->work));
457: PetscCall(VecDestroy(&ipmP->lambdae));
458: PetscCall(VecDestroy(&ipmP->lambdai));
459: PetscCall(VecDestroy(&ipmP->s));
460: PetscCall(VecDestroy(&ipmP->ds));
461: PetscCall(VecDestroy(&ipmP->ci));
463: PetscCall(VecDestroy(&ipmP->rhs_x));
464: PetscCall(VecDestroy(&ipmP->rhs_lambdae));
465: PetscCall(VecDestroy(&ipmP->rhs_lambdai));
466: PetscCall(VecDestroy(&ipmP->rhs_s));
468: PetscCall(VecDestroy(&ipmP->save_x));
469: PetscCall(VecDestroy(&ipmP->save_lambdae));
470: PetscCall(VecDestroy(&ipmP->save_lambdai));
471: PetscCall(VecDestroy(&ipmP->save_s));
473: PetscCall(VecScatterDestroy(&ipmP->step1));
474: PetscCall(VecScatterDestroy(&ipmP->step2));
475: PetscCall(VecScatterDestroy(&ipmP->step3));
476: PetscCall(VecScatterDestroy(&ipmP->step4));
478: PetscCall(VecScatterDestroy(&ipmP->rhs1));
479: PetscCall(VecScatterDestroy(&ipmP->rhs2));
480: PetscCall(VecScatterDestroy(&ipmP->rhs3));
481: PetscCall(VecScatterDestroy(&ipmP->rhs4));
483: PetscCall(VecScatterDestroy(&ipmP->ci_scat));
484: PetscCall(VecScatterDestroy(&ipmP->xl_scat));
485: PetscCall(VecScatterDestroy(&ipmP->xu_scat));
487: PetscCall(VecDestroy(&ipmP->dlambdai));
488: PetscCall(VecDestroy(&ipmP->dlambdae));
489: PetscCall(VecDestroy(&ipmP->Zero_nb));
490: PetscCall(VecDestroy(&ipmP->One_nb));
491: PetscCall(VecDestroy(&ipmP->Inf_nb));
492: PetscCall(VecDestroy(&ipmP->complementarity));
494: PetscCall(VecDestroy(&ipmP->bigrhs));
495: PetscCall(VecDestroy(&ipmP->bigstep));
496: PetscCall(MatDestroy(&ipmP->Ai));
497: PetscCall(MatDestroy(&ipmP->K));
498: PetscCall(ISDestroy(&ipmP->isxu));
499: PetscCall(ISDestroy(&ipmP->isxl));
500: PetscCall(KSPDestroy(&tao->ksp));
501: PetscCall(PetscFree(tao->data));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: static PetscErrorCode TaoSetFromOptions_IPM(Tao tao, PetscOptionItems PetscOptionsObject)
506: {
507: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
509: PetscFunctionBegin;
510: PetscOptionsHeadBegin(PetscOptionsObject, "IPM method for constrained optimization");
511: PetscCall(PetscOptionsBool("-tao_ipm_monitorkkt", "monitor kkt status", NULL, ipmP->monitorkkt, &ipmP->monitorkkt, NULL));
512: PetscCall(PetscOptionsReal("-tao_ipm_pushs", "parameter to push initial slack variables away from bounds", NULL, ipmP->pushs, &ipmP->pushs, NULL));
513: PetscCall(PetscOptionsReal("-tao_ipm_pushnu", "parameter to push initial (inequality) dual variables away from bounds", NULL, ipmP->pushnu, &ipmP->pushnu, NULL));
514: PetscOptionsHeadEnd();
515: PetscCall(KSPSetFromOptions(tao->ksp));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
520: {
521: return PETSC_SUCCESS;
522: }
524: /* IPMObjectiveAndGradient()
525: f = d'x + 0.5 * x' * H * x
526: rd = H*x + d + Ae'*lame - Ai'*lami
527: rpe = Ae*x - be
528: rpi = Ai*x - yi - bi
529: mu = yi' * lami/mi;
530: com = yi.*lami
532: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
533: */
534: /*
535: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
536: {
537: Tao tao = (Tao)tptr;
538: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
540: PetscFunctionBegin;
541: PetscCall(IPMComputeKKT(tao));
542: *f = ipmP->phi;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
545: */
547: /*
548: f = d'x + 0.5 * x' * H * x
549: rd = H*x + d + Ae'*lame - Ai'*lami
550: Ai = jac_ineq
551: I (w/lb)
552: -I (w/ub)
554: rpe = ce
555: rpi = ci - s;
556: com = s.*lami
557: mu = yi' * lami/mi;
559: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
560: */
561: static PetscErrorCode IPMComputeKKT(Tao tao)
562: {
563: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
564: PetscScalar norm;
566: PetscFunctionBegin;
567: PetscCall(VecCopy(tao->gradient, ipmP->rd));
569: if (ipmP->me > 0) {
570: /* rd = gradient + Ae'*lambdae */
571: PetscCall(MatMultTranspose(tao->jacobian_equality, ipmP->lambdae, ipmP->work));
572: PetscCall(VecAXPY(ipmP->rd, 1.0, ipmP->work));
574: /* rpe = ce(x) */
575: PetscCall(VecCopy(tao->constraints_equality, ipmP->rpe));
576: }
577: if (ipmP->nb > 0) {
578: /* rd = rd - Ai'*lambdai */
579: PetscCall(MatMultTranspose(ipmP->Ai, ipmP->lambdai, ipmP->work));
580: PetscCall(VecAXPY(ipmP->rd, -1.0, ipmP->work));
582: /* rpi = cin - s */
583: PetscCall(VecCopy(ipmP->ci, ipmP->rpi));
584: PetscCall(VecAXPY(ipmP->rpi, -1.0, ipmP->s));
586: /* com = s .* lami */
587: PetscCall(VecPointwiseMult(ipmP->complementarity, ipmP->s, ipmP->lambdai));
588: }
589: /* phi = ||rd; rpe; rpi; com|| */
590: PetscCall(VecDot(ipmP->rd, ipmP->rd, &norm));
591: ipmP->phi = norm;
592: if (ipmP->me > 0) {
593: PetscCall(VecDot(ipmP->rpe, ipmP->rpe, &norm));
594: ipmP->phi += norm;
595: }
596: if (ipmP->nb > 0) {
597: PetscCall(VecDot(ipmP->rpi, ipmP->rpi, &norm));
598: ipmP->phi += norm;
599: PetscCall(VecDot(ipmP->complementarity, ipmP->complementarity, &norm));
600: ipmP->phi += norm;
601: /* mu = s'*lami/nb */
602: PetscCall(VecDot(ipmP->s, ipmP->lambdai, &ipmP->mu));
603: ipmP->mu /= ipmP->nb;
604: } else {
605: ipmP->mu = 1.0;
606: }
608: ipmP->phi = PetscSqrtScalar(ipmP->phi);
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: /* evaluate user info at current point */
613: static PetscErrorCode IPMEvaluate(Tao tao)
614: {
615: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
617: PetscFunctionBegin;
618: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &ipmP->kkt_f, tao->gradient));
619: PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
620: if (ipmP->me > 0) {
621: PetscCall(TaoComputeEqualityConstraints(tao, tao->solution, tao->constraints_equality));
622: PetscCall(TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre));
623: }
624: if (ipmP->mi > 0) {
625: PetscCall(TaoComputeInequalityConstraints(tao, tao->solution, tao->constraints_inequality));
626: PetscCall(TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre));
627: }
628: if (ipmP->nb > 0) {
629: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
630: PetscCall(IPMUpdateAi(tao));
631: }
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: /* Push initial point away from bounds */
636: static PetscErrorCode IPMPushInitialPoint(Tao tao)
637: {
638: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
640: PetscFunctionBegin;
641: PetscCall(TaoComputeVariableBounds(tao));
642: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
643: if (ipmP->nb > 0) {
644: PetscCall(VecSet(ipmP->s, ipmP->pushs));
645: PetscCall(VecSet(ipmP->lambdai, ipmP->pushnu));
646: if (ipmP->mi > 0) PetscCall(VecSet(tao->DI, ipmP->pushnu));
647: }
648: if (ipmP->me > 0) {
649: PetscCall(VecSet(tao->DE, 1.0));
650: PetscCall(VecSet(ipmP->lambdae, 1.0));
651: }
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: static PetscErrorCode IPMUpdateAi(Tao tao)
656: {
657: /* Ai = Ji
658: I (w/lb)
659: -I (w/ub) */
661: /* Ci = user->ci
662: Xi - lb (w/lb)
663: -Xi + ub (w/ub) */
665: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
666: MPI_Comm comm;
667: PetscInt i;
668: PetscScalar newval;
669: PetscInt newrow, newcol, ncols;
670: const PetscScalar *vals;
671: const PetscInt *cols;
672: PetscInt astart, aend, jstart, jend;
673: PetscInt *nonzeros;
674: PetscInt r2, r3, r4;
675: PetscMPIInt size;
676: Vec solu;
677: PetscInt nloc;
679: PetscFunctionBegin;
680: r2 = ipmP->mi;
681: r3 = r2 + ipmP->nxlb;
682: r4 = r3 + ipmP->nxub;
684: if (!ipmP->nb) PetscFunctionReturn(PETSC_SUCCESS);
686: /* Create Ai matrix if it doesn't exist yet */
687: if (!ipmP->Ai) {
688: comm = ((PetscObject)tao->solution)->comm;
689: PetscCallMPI(MPI_Comm_size(comm, &size));
690: if (size == 1) {
691: PetscCall(PetscMalloc1(ipmP->nb, &nonzeros));
692: for (i = 0; i < ipmP->mi; i++) {
693: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
694: nonzeros[i] = ncols;
695: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
696: }
697: for (i = r2; i < r4; i++) nonzeros[i] = 1;
698: }
699: PetscCall(MatCreate(comm, &ipmP->Ai));
700: PetscCall(MatSetType(ipmP->Ai, MATAIJ));
702: PetscCall(TaoGetSolution(tao, &solu));
703: PetscCall(VecGetLocalSize(solu, &nloc));
704: PetscCall(MatSetSizes(ipmP->Ai, PETSC_DECIDE, nloc, ipmP->nb, PETSC_DECIDE));
705: PetscCall(MatSetFromOptions(ipmP->Ai));
706: PetscCall(MatMPIAIJSetPreallocation(ipmP->Ai, ipmP->nb, NULL, ipmP->nb, NULL));
707: PetscCall(MatSeqAIJSetPreallocation(ipmP->Ai, PETSC_DEFAULT, nonzeros));
708: if (size == 1) PetscCall(PetscFree(nonzeros));
709: }
711: /* Copy values from user jacobian to Ai */
712: PetscCall(MatGetOwnershipRange(ipmP->Ai, &astart, &aend));
714: /* Ai w/lb */
715: if (ipmP->mi) {
716: PetscCall(MatZeroEntries(ipmP->Ai));
717: PetscCall(MatGetOwnershipRange(tao->jacobian_inequality, &jstart, &jend));
718: for (i = jstart; i < jend; i++) {
719: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
720: newrow = i;
721: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, ncols, cols, vals, INSERT_VALUES));
722: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
723: }
724: }
726: /* I w/ xlb */
727: if (ipmP->nxlb) {
728: for (i = 0; i < ipmP->nxlb; i++) {
729: if (i >= astart && i < aend) {
730: newrow = i + r2;
731: newcol = i;
732: newval = 1.0;
733: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
734: }
735: }
736: }
737: if (ipmP->nxub) {
738: /* I w/ xub */
739: for (i = 0; i < ipmP->nxub; i++) {
740: if (i >= astart && i < aend) {
741: newrow = i + r3;
742: newcol = i;
743: newval = -1.0;
744: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
745: }
746: }
747: }
749: PetscCall(MatAssemblyBegin(ipmP->Ai, MAT_FINAL_ASSEMBLY));
750: PetscCall(MatAssemblyEnd(ipmP->Ai, MAT_FINAL_ASSEMBLY));
751: CHKMEMQ;
753: PetscCall(VecSet(ipmP->ci, 0.0));
755: /* user ci */
756: if (ipmP->mi > 0) {
757: PetscCall(VecScatterBegin(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
758: PetscCall(VecScatterEnd(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
759: }
760: if (!ipmP->work) PetscCall(VecDuplicate(tao->solution, &ipmP->work));
761: PetscCall(VecCopy(tao->solution, ipmP->work));
762: if (tao->XL) {
763: PetscCall(VecAXPY(ipmP->work, -1.0, tao->XL));
765: /* lower bounds on variables */
766: if (ipmP->nxlb > 0) {
767: PetscCall(VecScatterBegin(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
768: PetscCall(VecScatterEnd(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
769: }
770: }
771: if (tao->XU) {
772: /* upper bounds on variables */
773: PetscCall(VecCopy(tao->solution, ipmP->work));
774: PetscCall(VecScale(ipmP->work, -1.0));
775: PetscCall(VecAXPY(ipmP->work, 1.0, tao->XU));
776: if (ipmP->nxub > 0) {
777: PetscCall(VecScatterBegin(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
778: PetscCall(VecScatterEnd(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
779: }
780: }
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: /* create K = [ Hlag , 0 , Ae', -Ai'];
785: [Ae , 0, 0 , 0];
786: [Ai ,-I, 0 , 0];
787: [ 0 , S , 0, Y ]; */
788: static PetscErrorCode IPMUpdateK(Tao tao)
789: {
790: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
791: MPI_Comm comm;
792: PetscMPIInt size;
793: PetscInt i, j, row;
794: PetscInt ncols, newcol, newcols[2], newrow;
795: const PetscInt *cols;
796: const PetscReal *vals;
797: const PetscReal *l, *y;
798: PetscReal *newvals;
799: PetscReal newval;
800: PetscInt subsize;
801: const PetscInt *indices;
802: PetscInt *nonzeros, *d_nonzeros, *o_nonzeros;
803: PetscInt bigsize;
804: PetscInt r1, r2, r3;
805: PetscInt c1, c2, c3;
806: PetscInt klocalsize;
807: PetscInt hstart, hend, kstart, kend;
808: PetscInt aistart, aiend, aestart, aeend;
809: PetscInt sstart, send;
811: PetscFunctionBegin;
812: comm = ((PetscObject)tao->solution)->comm;
813: PetscCallMPI(MPI_Comm_size(comm, &size));
814: PetscCall(IPMUpdateAi(tao));
816: /* allocate workspace */
817: subsize = PetscMax(ipmP->n, ipmP->nb);
818: subsize = PetscMax(ipmP->me, subsize);
819: subsize = PetscMax(2, subsize);
820: PetscCall(PetscMalloc1(subsize, &indices));
821: PetscCall(PetscMalloc1(subsize, &newvals));
823: r1 = c1 = ipmP->n;
824: r2 = r1 + ipmP->me;
825: c2 = c1 + ipmP->nb;
826: r3 = c3 = r2 + ipmP->nb;
828: bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
829: PetscCall(VecGetOwnershipRange(ipmP->bigrhs, &kstart, &kend));
830: PetscCall(MatGetOwnershipRange(tao->hessian, &hstart, &hend));
831: klocalsize = kend - kstart;
832: if (!ipmP->K) {
833: if (size == 1) {
834: PetscCall(PetscMalloc1(kend - kstart, &nonzeros));
835: for (i = 0; i < bigsize; i++) {
836: if (i < r1) {
837: PetscCall(MatGetRow(tao->hessian, i, &ncols, NULL, NULL));
838: nonzeros[i] = ncols;
839: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, NULL, NULL));
840: nonzeros[i] += ipmP->me + ipmP->nb;
841: } else if (i < r2) {
842: nonzeros[i - kstart] = ipmP->n;
843: } else if (i < r3) {
844: nonzeros[i - kstart] = ipmP->n + 1;
845: } else if (i < bigsize) {
846: nonzeros[i - kstart] = 2;
847: }
848: }
849: PetscCall(MatCreate(comm, &ipmP->K));
850: PetscCall(MatSetType(ipmP->K, MATSEQAIJ));
851: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
852: PetscCall(MatSeqAIJSetPreallocation(ipmP->K, 0, nonzeros));
853: PetscCall(MatSetFromOptions(ipmP->K));
854: PetscCall(PetscFree(nonzeros));
855: } else {
856: PetscCall(PetscMalloc1(kend - kstart, &d_nonzeros));
857: PetscCall(PetscMalloc1(kend - kstart, &o_nonzeros));
858: for (i = kstart; i < kend; i++) {
859: if (i < r1) {
860: /* TODO fix preallocation for mpi mats */
861: d_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, kend - kstart);
862: o_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, bigsize - (kend - kstart));
863: } else if (i < r2) {
864: d_nonzeros[i - kstart] = PetscMin(ipmP->n, kend - kstart);
865: o_nonzeros[i - kstart] = PetscMin(ipmP->n, bigsize - (kend - kstart));
866: } else if (i < r3) {
867: d_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, kend - kstart);
868: o_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, bigsize - (kend - kstart));
869: } else {
870: d_nonzeros[i - kstart] = PetscMin(2, kend - kstart);
871: o_nonzeros[i - kstart] = PetscMin(2, bigsize - (kend - kstart));
872: }
873: }
874: PetscCall(MatCreate(comm, &ipmP->K));
875: PetscCall(MatSetType(ipmP->K, MATMPIAIJ));
876: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
877: PetscCall(MatMPIAIJSetPreallocation(ipmP->K, 0, d_nonzeros, 0, o_nonzeros));
878: PetscCall(PetscFree(d_nonzeros));
879: PetscCall(PetscFree(o_nonzeros));
880: PetscCall(MatSetFromOptions(ipmP->K));
881: }
882: }
884: PetscCall(MatZeroEntries(ipmP->K));
885: /* Copy H */
886: for (i = hstart; i < hend; i++) {
887: PetscCall(MatGetRow(tao->hessian, i, &ncols, &cols, &vals));
888: if (ncols > 0) PetscCall(MatSetValues(ipmP->K, 1, &i, ncols, cols, vals, INSERT_VALUES));
889: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, &cols, &vals));
890: }
892: /* Copy Ae and Ae' */
893: if (ipmP->me > 0) {
894: PetscCall(MatGetOwnershipRange(tao->jacobian_equality, &aestart, &aeend));
895: for (i = aestart; i < aeend; i++) {
896: PetscCall(MatGetRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
897: if (ncols > 0) {
898: /*Ae*/
899: row = i + r1;
900: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
901: /*Ae'*/
902: for (j = 0; j < ncols; j++) {
903: newcol = i + c2;
904: newrow = cols[j];
905: newval = vals[j];
906: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
907: }
908: }
909: PetscCall(MatRestoreRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
910: }
911: }
913: if (ipmP->nb > 0) {
914: PetscCall(MatGetOwnershipRange(ipmP->Ai, &aistart, &aiend));
915: /* Copy Ai,and Ai' */
916: for (i = aistart; i < aiend; i++) {
917: row = i + r2;
918: PetscCall(MatGetRow(ipmP->Ai, i, &ncols, &cols, &vals));
919: if (ncols > 0) {
920: /*Ai*/
921: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
922: /*-Ai'*/
923: for (j = 0; j < ncols; j++) {
924: newcol = i + c3;
925: newrow = cols[j];
926: newval = -vals[j];
927: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
928: }
929: }
930: PetscCall(MatRestoreRow(ipmP->Ai, i, &ncols, &cols, &vals));
931: }
933: /* -I */
934: for (i = kstart; i < kend; i++) {
935: if (i >= r2 && i < r3) {
936: newrow = i;
937: newcol = i - r2 + c1;
938: newval = -1.0;
939: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
940: }
941: }
943: /* Copy L,Y */
944: PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send));
945: PetscCall(VecGetArrayRead(ipmP->lambdai, &l));
946: PetscCall(VecGetArrayRead(ipmP->s, &y));
948: for (i = sstart; i < send; i++) {
949: newcols[0] = c1 + i;
950: newcols[1] = c3 + i;
951: newvals[0] = l[i - sstart];
952: newvals[1] = y[i - sstart];
953: newrow = r3 + i;
954: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 2, newcols, newvals, INSERT_VALUES));
955: }
957: PetscCall(VecRestoreArrayRead(ipmP->lambdai, &l));
958: PetscCall(VecRestoreArrayRead(ipmP->s, &y));
959: }
961: PetscCall(PetscFree(indices));
962: PetscCall(PetscFree(newvals));
963: PetscCall(MatAssemblyBegin(ipmP->K, MAT_FINAL_ASSEMBLY));
964: PetscCall(MatAssemblyEnd(ipmP->K, MAT_FINAL_ASSEMBLY));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: static PetscErrorCode IPMGatherRHS(Tao tao, Vec RHS, Vec X1, Vec X2, Vec X3, Vec X4)
969: {
970: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
972: PetscFunctionBegin;
973: /* rhs = [x1 (n)
974: x2 (me)
975: x3 (nb)
976: x4 (nb)] */
977: if (X1) {
978: PetscCall(VecScatterBegin(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
979: PetscCall(VecScatterEnd(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
980: }
981: if (ipmP->me > 0 && X2) {
982: PetscCall(VecScatterBegin(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
983: PetscCall(VecScatterEnd(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
984: }
985: if (ipmP->nb > 0) {
986: if (X3) {
987: PetscCall(VecScatterBegin(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
988: PetscCall(VecScatterEnd(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
989: }
990: if (X4) {
991: PetscCall(VecScatterBegin(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
992: PetscCall(VecScatterEnd(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
993: }
994: }
995: PetscFunctionReturn(PETSC_SUCCESS);
996: }
998: static PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
999: {
1000: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1002: PetscFunctionBegin;
1003: CHKMEMQ;
1004: /* [x1 (n)
1005: x2 (nb) may be 0
1006: x3 (me) may be 0
1007: x4 (nb) may be 0 */
1008: if (X1) {
1009: PetscCall(VecScatterBegin(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1010: PetscCall(VecScatterEnd(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1011: }
1012: if (X2 && ipmP->nb > 0) {
1013: PetscCall(VecScatterBegin(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1014: PetscCall(VecScatterEnd(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1015: }
1016: if (X3 && ipmP->me > 0) {
1017: PetscCall(VecScatterBegin(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1018: PetscCall(VecScatterEnd(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1019: }
1020: if (X4 && ipmP->nb > 0) {
1021: PetscCall(VecScatterBegin(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1022: PetscCall(VecScatterEnd(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1023: }
1024: CHKMEMQ;
1025: PetscFunctionReturn(PETSC_SUCCESS);
1026: }
1028: /*MC
1029: TAOIPM - Interior point algorithm for generally constrained optimization.
1031: Options Database Keys:
1032: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1033: - -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1035: Notes:
1036: This algorithm is more of a place-holder for future constrained optimization algorithms and should not yet be used for large problems or production code.
1037: Level: beginner
1039: .seealso: `Tao`, `TAOPDIPM`, `TaoType`
1040: M*/
1042: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1043: {
1044: TAO_IPM *ipmP;
1046: PetscFunctionBegin;
1047: tao->ops->setup = TaoSetup_IPM;
1048: tao->ops->solve = TaoSolve_IPM;
1049: tao->ops->view = TaoView_IPM;
1050: tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1051: tao->ops->destroy = TaoDestroy_IPM;
1052: tao->uses_gradient = PETSC_TRUE;
1053: tao->uses_hessian_matrices = PETSC_TRUE;
1054: /* tao->ops->computedual = TaoComputeDual_IPM; */
1056: PetscCall(PetscNew(&ipmP));
1057: tao->data = (void *)ipmP;
1059: /* Override default settings (unless already changed) */
1060: PetscCall(TaoParametersInitialize(tao));
1061: PetscObjectParameterSetDefault(tao, max_it, 200);
1062: PetscObjectParameterSetDefault(tao, max_funcs, 500);
1064: ipmP->dec = 10000; /* line search criteria */
1065: ipmP->taumin = 0.995;
1066: ipmP->monitorkkt = PETSC_FALSE;
1067: ipmP->pushs = 100;
1068: ipmP->pushnu = 100;
1069: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
1070: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
1071: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
1072: PetscFunctionReturn(PETSC_SUCCESS);
1073: }