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) {
200: PetscCall(VecDuplicate(tao->solution, &tao->gradient));
201: PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
202: PetscCall(VecDuplicate(tao->solution, &ipmP->rd));
203: PetscCall(VecDuplicate(tao->solution, &ipmP->rhs_x));
204: PetscCall(VecDuplicate(tao->solution, &ipmP->work));
205: PetscCall(VecDuplicate(tao->solution, &ipmP->save_x));
206: }
207: if (tao->constraints_equality) {
208: PetscCall(VecGetSize(tao->constraints_equality, &ipmP->me));
209: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->lambdae));
210: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->dlambdae));
211: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->rhs_lambdae));
212: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->save_lambdae));
213: PetscCall(VecDuplicate(tao->constraints_equality, &ipmP->rpe));
214: PetscCall(VecDuplicate(tao->constraints_equality, &tao->DE));
215: }
216: if (tao->constraints_inequality) PetscCall(VecDuplicate(tao->constraints_inequality, &tao->DI));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: static PetscErrorCode IPMInitializeBounds(Tao tao)
221: {
222: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
223: Vec xtmp;
224: PetscInt xstart, xend;
225: PetscInt ucstart, ucend; /* user ci */
226: PetscInt ucestart, uceend; /* user ce */
227: PetscInt sstart = 0, send = 0;
228: PetscInt bigsize;
229: PetscInt i, counter, nloc;
230: PetscInt *cind, *xind, *ucind, *uceind, *stepind;
231: VecType vtype;
232: const PetscInt *xli, *xui;
233: PetscInt xl_offset, xu_offset;
234: IS bigxl, bigxu, isuc, isc, isx, sis, is1;
235: MPI_Comm comm;
237: PetscFunctionBegin;
238: cind = xind = ucind = uceind = stepind = NULL;
239: ipmP->mi = 0;
240: ipmP->nxlb = 0;
241: ipmP->nxub = 0;
242: ipmP->nb = 0;
243: ipmP->nslack = 0;
245: PetscCall(VecDuplicate(tao->solution, &xtmp));
246: PetscCall(TaoComputeVariableBounds(tao));
247: if (tao->XL) {
248: PetscCall(VecSet(xtmp, PETSC_NINFINITY));
249: PetscCall(VecWhichGreaterThan(tao->XL, xtmp, &ipmP->isxl));
250: PetscCall(ISGetSize(ipmP->isxl, &ipmP->nxlb));
251: } else {
252: ipmP->nxlb = 0;
253: }
254: if (tao->XU) {
255: PetscCall(VecSet(xtmp, PETSC_INFINITY));
256: PetscCall(VecWhichLessThan(tao->XU, xtmp, &ipmP->isxu));
257: PetscCall(ISGetSize(ipmP->isxu, &ipmP->nxub));
258: } else {
259: ipmP->nxub = 0;
260: }
261: PetscCall(VecDestroy(&xtmp));
262: if (tao->constraints_inequality) {
263: PetscCall(VecGetSize(tao->constraints_inequality, &ipmP->mi));
264: } else {
265: ipmP->mi = 0;
266: }
267: ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi;
269: PetscCall(PetscObjectGetComm((PetscObject)tao->solution, &comm));
271: bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
272: PetscCall(PetscMalloc1(bigsize, &stepind));
273: PetscCall(PetscMalloc1(ipmP->n, &xind));
274: PetscCall(PetscMalloc1(ipmP->me, &uceind));
275: PetscCall(VecGetOwnershipRange(tao->solution, &xstart, &xend));
277: if (ipmP->nb > 0) {
278: PetscCall(VecCreate(comm, &ipmP->s));
279: PetscCall(VecSetSizes(ipmP->s, PETSC_DECIDE, ipmP->nb));
280: PetscCall(VecSetFromOptions(ipmP->s));
281: PetscCall(VecDuplicate(ipmP->s, &ipmP->ds));
282: PetscCall(VecDuplicate(ipmP->s, &ipmP->rhs_s));
283: PetscCall(VecDuplicate(ipmP->s, &ipmP->complementarity));
284: PetscCall(VecDuplicate(ipmP->s, &ipmP->ci));
286: PetscCall(VecDuplicate(ipmP->s, &ipmP->lambdai));
287: PetscCall(VecDuplicate(ipmP->s, &ipmP->dlambdai));
288: PetscCall(VecDuplicate(ipmP->s, &ipmP->rhs_lambdai));
289: PetscCall(VecDuplicate(ipmP->s, &ipmP->save_lambdai));
291: PetscCall(VecDuplicate(ipmP->s, &ipmP->save_s));
292: PetscCall(VecDuplicate(ipmP->s, &ipmP->rpi));
293: PetscCall(VecDuplicate(ipmP->s, &ipmP->Zero_nb));
294: PetscCall(VecSet(ipmP->Zero_nb, 0.0));
295: PetscCall(VecDuplicate(ipmP->s, &ipmP->One_nb));
296: PetscCall(VecSet(ipmP->One_nb, 1.0));
297: PetscCall(VecDuplicate(ipmP->s, &ipmP->Inf_nb));
298: PetscCall(VecSet(ipmP->Inf_nb, PETSC_INFINITY));
300: PetscCall(PetscMalloc1(ipmP->nb, &cind));
301: PetscCall(PetscMalloc1(ipmP->mi, &ucind));
302: PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send));
304: if (ipmP->mi > 0) {
305: PetscCall(VecGetOwnershipRange(tao->constraints_inequality, &ucstart, &ucend));
306: counter = 0;
307: for (i = ucstart; i < ucend; i++) cind[counter++] = i;
308: PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isuc));
309: PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isc));
310: PetscCall(VecScatterCreate(tao->constraints_inequality, isuc, ipmP->ci, isc, &ipmP->ci_scat));
312: PetscCall(ISDestroy(&isuc));
313: PetscCall(ISDestroy(&isc));
314: }
315: /* need to know how may xbound indices are on each process */
316: /* TODO better way */
317: if (ipmP->nxlb) {
318: PetscCall(ISAllGather(ipmP->isxl, &bigxl));
319: PetscCall(ISGetIndices(bigxl, &xli));
320: /* find offsets for this processor */
321: xl_offset = ipmP->mi;
322: for (i = 0; i < ipmP->nxlb; i++) {
323: if (xli[i] < xstart) {
324: xl_offset++;
325: } else break;
326: }
327: PetscCall(ISRestoreIndices(bigxl, &xli));
329: PetscCall(ISGetIndices(ipmP->isxl, &xli));
330: PetscCall(ISGetLocalSize(ipmP->isxl, &nloc));
331: for (i = 0; i < nloc; i++) {
332: xind[i] = xli[i];
333: cind[i] = xl_offset + i;
334: }
336: PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx));
337: PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc));
338: PetscCall(VecScatterCreate(tao->XL, isx, ipmP->ci, isc, &ipmP->xl_scat));
339: PetscCall(ISDestroy(&isx));
340: PetscCall(ISDestroy(&isc));
341: PetscCall(ISDestroy(&bigxl));
342: }
344: if (ipmP->nxub) {
345: PetscCall(ISAllGather(ipmP->isxu, &bigxu));
346: PetscCall(ISGetIndices(bigxu, &xui));
347: /* find offsets for this processor */
348: xu_offset = ipmP->mi + ipmP->nxlb;
349: for (i = 0; i < ipmP->nxub; i++) {
350: if (xui[i] < xstart) {
351: xu_offset++;
352: } else break;
353: }
354: PetscCall(ISRestoreIndices(bigxu, &xui));
356: PetscCall(ISGetIndices(ipmP->isxu, &xui));
357: PetscCall(ISGetLocalSize(ipmP->isxu, &nloc));
358: for (i = 0; i < nloc; i++) {
359: xind[i] = xui[i];
360: cind[i] = xu_offset + i;
361: }
363: PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx));
364: PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc));
365: PetscCall(VecScatterCreate(tao->XU, isx, ipmP->ci, isc, &ipmP->xu_scat));
366: PetscCall(ISDestroy(&isx));
367: PetscCall(ISDestroy(&isc));
368: PetscCall(ISDestroy(&bigxu));
369: }
370: }
371: PetscCall(VecCreate(comm, &ipmP->bigrhs));
372: PetscCall(VecGetType(tao->solution, &vtype));
373: PetscCall(VecSetType(ipmP->bigrhs, vtype));
374: PetscCall(VecSetSizes(ipmP->bigrhs, PETSC_DECIDE, bigsize));
375: PetscCall(VecSetFromOptions(ipmP->bigrhs));
376: PetscCall(VecDuplicate(ipmP->bigrhs, &ipmP->bigstep));
378: /* create scatters for step->x and x->rhs */
379: for (i = xstart; i < xend; i++) {
380: stepind[i - xstart] = i;
381: xind[i - xstart] = i;
382: }
383: PetscCall(ISCreateGeneral(comm, xend - xstart, stepind, PETSC_COPY_VALUES, &sis));
384: PetscCall(ISCreateGeneral(comm, xend - xstart, xind, PETSC_COPY_VALUES, &is1));
385: PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->solution, is1, &ipmP->step1));
386: PetscCall(VecScatterCreate(tao->solution, is1, ipmP->bigrhs, sis, &ipmP->rhs1));
387: PetscCall(ISDestroy(&sis));
388: PetscCall(ISDestroy(&is1));
390: if (ipmP->nb > 0) {
391: for (i = sstart; i < send; i++) {
392: stepind[i - sstart] = i + ipmP->n;
393: cind[i - sstart] = i;
394: }
395: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
396: PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1));
397: PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step2));
398: PetscCall(ISDestroy(&sis));
400: for (i = sstart; i < send; i++) {
401: stepind[i - sstart] = i + ipmP->n + ipmP->me;
402: cind[i - sstart] = i;
403: }
404: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
405: PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs3));
406: PetscCall(ISDestroy(&sis));
407: PetscCall(ISDestroy(&is1));
408: }
410: if (ipmP->me > 0) {
411: PetscCall(VecGetOwnershipRange(tao->constraints_equality, &ucestart, &uceend));
412: for (i = ucestart; i < uceend; i++) {
413: stepind[i - ucestart] = i + ipmP->n + ipmP->nb;
414: uceind[i - ucestart] = i;
415: }
417: PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis));
418: PetscCall(ISCreateGeneral(comm, uceend - ucestart, uceind, PETSC_COPY_VALUES, &is1));
419: PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->constraints_equality, is1, &ipmP->step3));
420: PetscCall(ISDestroy(&sis));
422: for (i = ucestart; i < uceend; i++) stepind[i - ucestart] = i + ipmP->n;
424: PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis));
425: PetscCall(VecScatterCreate(tao->constraints_equality, is1, ipmP->bigrhs, sis, &ipmP->rhs2));
426: PetscCall(ISDestroy(&sis));
427: PetscCall(ISDestroy(&is1));
428: }
430: if (ipmP->nb > 0) {
431: for (i = sstart; i < send; i++) {
432: stepind[i - sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
433: cind[i - sstart] = i;
434: }
435: PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1));
436: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
437: PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step4));
438: PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs4));
439: PetscCall(ISDestroy(&sis));
440: PetscCall(ISDestroy(&is1));
441: }
443: PetscCall(PetscFree(stepind));
444: PetscCall(PetscFree(cind));
445: PetscCall(PetscFree(ucind));
446: PetscCall(PetscFree(uceind));
447: PetscCall(PetscFree(xind));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: static PetscErrorCode TaoDestroy_IPM(Tao tao)
452: {
453: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
455: PetscFunctionBegin;
456: PetscCall(VecDestroy(&ipmP->rd));
457: PetscCall(VecDestroy(&ipmP->rpe));
458: PetscCall(VecDestroy(&ipmP->rpi));
459: PetscCall(VecDestroy(&ipmP->work));
460: PetscCall(VecDestroy(&ipmP->lambdae));
461: PetscCall(VecDestroy(&ipmP->lambdai));
462: PetscCall(VecDestroy(&ipmP->s));
463: PetscCall(VecDestroy(&ipmP->ds));
464: PetscCall(VecDestroy(&ipmP->ci));
466: PetscCall(VecDestroy(&ipmP->rhs_x));
467: PetscCall(VecDestroy(&ipmP->rhs_lambdae));
468: PetscCall(VecDestroy(&ipmP->rhs_lambdai));
469: PetscCall(VecDestroy(&ipmP->rhs_s));
471: PetscCall(VecDestroy(&ipmP->save_x));
472: PetscCall(VecDestroy(&ipmP->save_lambdae));
473: PetscCall(VecDestroy(&ipmP->save_lambdai));
474: PetscCall(VecDestroy(&ipmP->save_s));
476: PetscCall(VecScatterDestroy(&ipmP->step1));
477: PetscCall(VecScatterDestroy(&ipmP->step2));
478: PetscCall(VecScatterDestroy(&ipmP->step3));
479: PetscCall(VecScatterDestroy(&ipmP->step4));
481: PetscCall(VecScatterDestroy(&ipmP->rhs1));
482: PetscCall(VecScatterDestroy(&ipmP->rhs2));
483: PetscCall(VecScatterDestroy(&ipmP->rhs3));
484: PetscCall(VecScatterDestroy(&ipmP->rhs4));
486: PetscCall(VecScatterDestroy(&ipmP->ci_scat));
487: PetscCall(VecScatterDestroy(&ipmP->xl_scat));
488: PetscCall(VecScatterDestroy(&ipmP->xu_scat));
490: PetscCall(VecDestroy(&ipmP->dlambdai));
491: PetscCall(VecDestroy(&ipmP->dlambdae));
492: PetscCall(VecDestroy(&ipmP->Zero_nb));
493: PetscCall(VecDestroy(&ipmP->One_nb));
494: PetscCall(VecDestroy(&ipmP->Inf_nb));
495: PetscCall(VecDestroy(&ipmP->complementarity));
497: PetscCall(VecDestroy(&ipmP->bigrhs));
498: PetscCall(VecDestroy(&ipmP->bigstep));
499: PetscCall(MatDestroy(&ipmP->Ai));
500: PetscCall(MatDestroy(&ipmP->K));
501: PetscCall(ISDestroy(&ipmP->isxu));
502: PetscCall(ISDestroy(&ipmP->isxl));
503: PetscCall(KSPDestroy(&tao->ksp));
504: PetscCall(PetscFree(tao->data));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode TaoSetFromOptions_IPM(Tao tao, PetscOptionItems *PetscOptionsObject)
509: {
510: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
512: PetscFunctionBegin;
513: PetscOptionsHeadBegin(PetscOptionsObject, "IPM method for constrained optimization");
514: PetscCall(PetscOptionsBool("-tao_ipm_monitorkkt", "monitor kkt status", NULL, ipmP->monitorkkt, &ipmP->monitorkkt, NULL));
515: PetscCall(PetscOptionsReal("-tao_ipm_pushs", "parameter to push initial slack variables away from bounds", NULL, ipmP->pushs, &ipmP->pushs, NULL));
516: PetscCall(PetscOptionsReal("-tao_ipm_pushnu", "parameter to push initial (inequality) dual variables away from bounds", NULL, ipmP->pushnu, &ipmP->pushnu, NULL));
517: PetscOptionsHeadEnd();
518: PetscCall(KSPSetFromOptions(tao->ksp));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
523: {
524: return PETSC_SUCCESS;
525: }
527: /* IPMObjectiveAndGradient()
528: f = d'x + 0.5 * x' * H * x
529: rd = H*x + d + Ae'*lame - Ai'*lami
530: rpe = Ae*x - be
531: rpi = Ai*x - yi - bi
532: mu = yi' * lami/mi;
533: com = yi.*lami
535: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
536: */
537: /*
538: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
539: {
540: Tao tao = (Tao)tptr;
541: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
542: PetscFunctionBegin;
543: PetscCall(IPMComputeKKT(tao));
544: *f = ipmP->phi;
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
547: */
549: /*
550: f = d'x + 0.5 * x' * H * x
551: rd = H*x + d + Ae'*lame - Ai'*lami
552: Ai = jac_ineq
553: I (w/lb)
554: -I (w/ub)
556: rpe = ce
557: rpi = ci - s;
558: com = s.*lami
559: mu = yi' * lami/mi;
561: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
562: */
563: static PetscErrorCode IPMComputeKKT(Tao tao)
564: {
565: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
566: PetscScalar norm;
568: PetscFunctionBegin;
569: PetscCall(VecCopy(tao->gradient, ipmP->rd));
571: if (ipmP->me > 0) {
572: /* rd = gradient + Ae'*lambdae */
573: PetscCall(MatMultTranspose(tao->jacobian_equality, ipmP->lambdae, ipmP->work));
574: PetscCall(VecAXPY(ipmP->rd, 1.0, ipmP->work));
576: /* rpe = ce(x) */
577: PetscCall(VecCopy(tao->constraints_equality, ipmP->rpe));
578: }
579: if (ipmP->nb > 0) {
580: /* rd = rd - Ai'*lambdai */
581: PetscCall(MatMultTranspose(ipmP->Ai, ipmP->lambdai, ipmP->work));
582: PetscCall(VecAXPY(ipmP->rd, -1.0, ipmP->work));
584: /* rpi = cin - s */
585: PetscCall(VecCopy(ipmP->ci, ipmP->rpi));
586: PetscCall(VecAXPY(ipmP->rpi, -1.0, ipmP->s));
588: /* com = s .* lami */
589: PetscCall(VecPointwiseMult(ipmP->complementarity, ipmP->s, ipmP->lambdai));
590: }
591: /* phi = ||rd; rpe; rpi; com|| */
592: PetscCall(VecDot(ipmP->rd, ipmP->rd, &norm));
593: ipmP->phi = norm;
594: if (ipmP->me > 0) {
595: PetscCall(VecDot(ipmP->rpe, ipmP->rpe, &norm));
596: ipmP->phi += norm;
597: }
598: if (ipmP->nb > 0) {
599: PetscCall(VecDot(ipmP->rpi, ipmP->rpi, &norm));
600: ipmP->phi += norm;
601: PetscCall(VecDot(ipmP->complementarity, ipmP->complementarity, &norm));
602: ipmP->phi += norm;
603: /* mu = s'*lami/nb */
604: PetscCall(VecDot(ipmP->s, ipmP->lambdai, &ipmP->mu));
605: ipmP->mu /= ipmP->nb;
606: } else {
607: ipmP->mu = 1.0;
608: }
610: ipmP->phi = PetscSqrtScalar(ipmP->phi);
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: /* evaluate user info at current point */
615: static PetscErrorCode IPMEvaluate(Tao tao)
616: {
617: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
619: PetscFunctionBegin;
620: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &ipmP->kkt_f, tao->gradient));
621: PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
622: if (ipmP->me > 0) {
623: PetscCall(TaoComputeEqualityConstraints(tao, tao->solution, tao->constraints_equality));
624: PetscCall(TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre));
625: }
626: if (ipmP->mi > 0) {
627: PetscCall(TaoComputeInequalityConstraints(tao, tao->solution, tao->constraints_inequality));
628: PetscCall(TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre));
629: }
630: if (ipmP->nb > 0) {
631: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
632: PetscCall(IPMUpdateAi(tao));
633: }
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: /* Push initial point away from bounds */
638: static PetscErrorCode IPMPushInitialPoint(Tao tao)
639: {
640: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
642: PetscFunctionBegin;
643: PetscCall(TaoComputeVariableBounds(tao));
644: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
645: if (ipmP->nb > 0) {
646: PetscCall(VecSet(ipmP->s, ipmP->pushs));
647: PetscCall(VecSet(ipmP->lambdai, ipmP->pushnu));
648: if (ipmP->mi > 0) PetscCall(VecSet(tao->DI, ipmP->pushnu));
649: }
650: if (ipmP->me > 0) {
651: PetscCall(VecSet(tao->DE, 1.0));
652: PetscCall(VecSet(ipmP->lambdae, 1.0));
653: }
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: static PetscErrorCode IPMUpdateAi(Tao tao)
658: {
659: /* Ai = Ji
660: I (w/lb)
661: -I (w/ub) */
663: /* Ci = user->ci
664: Xi - lb (w/lb)
665: -Xi + ub (w/ub) */
667: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
668: MPI_Comm comm;
669: PetscInt i;
670: PetscScalar newval;
671: PetscInt newrow, newcol, ncols;
672: const PetscScalar *vals;
673: const PetscInt *cols;
674: PetscInt astart, aend, jstart, jend;
675: PetscInt *nonzeros;
676: PetscInt r2, r3, r4;
677: PetscMPIInt size;
678: Vec solu;
679: PetscInt nloc;
681: PetscFunctionBegin;
682: r2 = ipmP->mi;
683: r3 = r2 + ipmP->nxlb;
684: r4 = r3 + ipmP->nxub;
686: if (!ipmP->nb) PetscFunctionReturn(PETSC_SUCCESS);
688: /* Create Ai matrix if it doesn't exist yet */
689: if (!ipmP->Ai) {
690: comm = ((PetscObject)(tao->solution))->comm;
691: PetscCallMPI(MPI_Comm_size(comm, &size));
692: if (size == 1) {
693: PetscCall(PetscMalloc1(ipmP->nb, &nonzeros));
694: for (i = 0; i < ipmP->mi; i++) {
695: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
696: nonzeros[i] = ncols;
697: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
698: }
699: for (i = r2; i < r4; i++) nonzeros[i] = 1;
700: }
701: PetscCall(MatCreate(comm, &ipmP->Ai));
702: PetscCall(MatSetType(ipmP->Ai, MATAIJ));
704: PetscCall(TaoGetSolution(tao, &solu));
705: PetscCall(VecGetLocalSize(solu, &nloc));
706: PetscCall(MatSetSizes(ipmP->Ai, PETSC_DECIDE, nloc, ipmP->nb, PETSC_DECIDE));
707: PetscCall(MatSetFromOptions(ipmP->Ai));
708: PetscCall(MatMPIAIJSetPreallocation(ipmP->Ai, ipmP->nb, NULL, ipmP->nb, NULL));
709: PetscCall(MatSeqAIJSetPreallocation(ipmP->Ai, PETSC_DEFAULT, nonzeros));
710: if (size == 1) PetscCall(PetscFree(nonzeros));
711: }
713: /* Copy values from user jacobian to Ai */
714: PetscCall(MatGetOwnershipRange(ipmP->Ai, &astart, &aend));
716: /* Ai w/lb */
717: if (ipmP->mi) {
718: PetscCall(MatZeroEntries(ipmP->Ai));
719: PetscCall(MatGetOwnershipRange(tao->jacobian_inequality, &jstart, &jend));
720: for (i = jstart; i < jend; i++) {
721: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
722: newrow = i;
723: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, ncols, cols, vals, INSERT_VALUES));
724: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
725: }
726: }
728: /* I w/ xlb */
729: if (ipmP->nxlb) {
730: for (i = 0; i < ipmP->nxlb; i++) {
731: if (i >= astart && i < aend) {
732: newrow = i + r2;
733: newcol = i;
734: newval = 1.0;
735: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
736: }
737: }
738: }
739: if (ipmP->nxub) {
740: /* I w/ xub */
741: for (i = 0; i < ipmP->nxub; i++) {
742: if (i >= astart && i < aend) {
743: newrow = i + r3;
744: newcol = i;
745: newval = -1.0;
746: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
747: }
748: }
749: }
751: PetscCall(MatAssemblyBegin(ipmP->Ai, MAT_FINAL_ASSEMBLY));
752: PetscCall(MatAssemblyEnd(ipmP->Ai, MAT_FINAL_ASSEMBLY));
753: CHKMEMQ;
755: PetscCall(VecSet(ipmP->ci, 0.0));
757: /* user ci */
758: if (ipmP->mi > 0) {
759: PetscCall(VecScatterBegin(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
760: PetscCall(VecScatterEnd(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
761: }
762: if (!ipmP->work) PetscCall(VecDuplicate(tao->solution, &ipmP->work));
763: PetscCall(VecCopy(tao->solution, ipmP->work));
764: if (tao->XL) {
765: PetscCall(VecAXPY(ipmP->work, -1.0, tao->XL));
767: /* lower bounds on variables */
768: if (ipmP->nxlb > 0) {
769: PetscCall(VecScatterBegin(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
770: PetscCall(VecScatterEnd(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
771: }
772: }
773: if (tao->XU) {
774: /* upper bounds on variables */
775: PetscCall(VecCopy(tao->solution, ipmP->work));
776: PetscCall(VecScale(ipmP->work, -1.0));
777: PetscCall(VecAXPY(ipmP->work, 1.0, tao->XU));
778: if (ipmP->nxub > 0) {
779: PetscCall(VecScatterBegin(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
780: PetscCall(VecScatterEnd(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
781: }
782: }
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: /* create K = [ Hlag , 0 , Ae', -Ai'];
787: [Ae , 0, 0 , 0];
788: [Ai ,-I, 0 , 0];
789: [ 0 , S , 0, Y ]; */
790: static PetscErrorCode IPMUpdateK(Tao tao)
791: {
792: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
793: MPI_Comm comm;
794: PetscMPIInt size;
795: PetscInt i, j, row;
796: PetscInt ncols, newcol, newcols[2], newrow;
797: const PetscInt *cols;
798: const PetscReal *vals;
799: const PetscReal *l, *y;
800: PetscReal *newvals;
801: PetscReal newval;
802: PetscInt subsize;
803: const PetscInt *indices;
804: PetscInt *nonzeros, *d_nonzeros, *o_nonzeros;
805: PetscInt bigsize;
806: PetscInt r1, r2, r3;
807: PetscInt c1, c2, c3;
808: PetscInt klocalsize;
809: PetscInt hstart, hend, kstart, kend;
810: PetscInt aistart, aiend, aestart, aeend;
811: PetscInt sstart, send;
813: PetscFunctionBegin;
814: comm = ((PetscObject)(tao->solution))->comm;
815: PetscCallMPI(MPI_Comm_size(comm, &size));
816: PetscCall(IPMUpdateAi(tao));
818: /* allocate workspace */
819: subsize = PetscMax(ipmP->n, ipmP->nb);
820: subsize = PetscMax(ipmP->me, subsize);
821: subsize = PetscMax(2, subsize);
822: PetscCall(PetscMalloc1(subsize, (PetscInt **)&indices));
823: PetscCall(PetscMalloc1(subsize, &newvals));
825: r1 = c1 = ipmP->n;
826: r2 = r1 + ipmP->me;
827: c2 = c1 + ipmP->nb;
828: r3 = c3 = r2 + ipmP->nb;
830: bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
831: PetscCall(VecGetOwnershipRange(ipmP->bigrhs, &kstart, &kend));
832: PetscCall(MatGetOwnershipRange(tao->hessian, &hstart, &hend));
833: klocalsize = kend - kstart;
834: if (!ipmP->K) {
835: if (size == 1) {
836: PetscCall(PetscMalloc1(kend - kstart, &nonzeros));
837: for (i = 0; i < bigsize; i++) {
838: if (i < r1) {
839: PetscCall(MatGetRow(tao->hessian, i, &ncols, NULL, NULL));
840: nonzeros[i] = ncols;
841: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, NULL, NULL));
842: nonzeros[i] += ipmP->me + ipmP->nb;
843: } else if (i < r2) {
844: nonzeros[i - kstart] = ipmP->n;
845: } else if (i < r3) {
846: nonzeros[i - kstart] = ipmP->n + 1;
847: } else if (i < bigsize) {
848: nonzeros[i - kstart] = 2;
849: }
850: }
851: PetscCall(MatCreate(comm, &ipmP->K));
852: PetscCall(MatSetType(ipmP->K, MATSEQAIJ));
853: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
854: PetscCall(MatSeqAIJSetPreallocation(ipmP->K, 0, nonzeros));
855: PetscCall(MatSetFromOptions(ipmP->K));
856: PetscCall(PetscFree(nonzeros));
857: } else {
858: PetscCall(PetscMalloc1(kend - kstart, &d_nonzeros));
859: PetscCall(PetscMalloc1(kend - kstart, &o_nonzeros));
860: for (i = kstart; i < kend; i++) {
861: if (i < r1) {
862: /* TODO fix preallocation for mpi mats */
863: d_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, kend - kstart);
864: o_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, bigsize - (kend - kstart));
865: } else if (i < r2) {
866: d_nonzeros[i - kstart] = PetscMin(ipmP->n, kend - kstart);
867: o_nonzeros[i - kstart] = PetscMin(ipmP->n, bigsize - (kend - kstart));
868: } else if (i < r3) {
869: d_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, kend - kstart);
870: o_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, bigsize - (kend - kstart));
871: } else {
872: d_nonzeros[i - kstart] = PetscMin(2, kend - kstart);
873: o_nonzeros[i - kstart] = PetscMin(2, bigsize - (kend - kstart));
874: }
875: }
876: PetscCall(MatCreate(comm, &ipmP->K));
877: PetscCall(MatSetType(ipmP->K, MATMPIAIJ));
878: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
879: PetscCall(MatMPIAIJSetPreallocation(ipmP->K, 0, d_nonzeros, 0, o_nonzeros));
880: PetscCall(PetscFree(d_nonzeros));
881: PetscCall(PetscFree(o_nonzeros));
882: PetscCall(MatSetFromOptions(ipmP->K));
883: }
884: }
886: PetscCall(MatZeroEntries(ipmP->K));
887: /* Copy H */
888: for (i = hstart; i < hend; i++) {
889: PetscCall(MatGetRow(tao->hessian, i, &ncols, &cols, &vals));
890: if (ncols > 0) PetscCall(MatSetValues(ipmP->K, 1, &i, ncols, cols, vals, INSERT_VALUES));
891: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, &cols, &vals));
892: }
894: /* Copy Ae and Ae' */
895: if (ipmP->me > 0) {
896: PetscCall(MatGetOwnershipRange(tao->jacobian_equality, &aestart, &aeend));
897: for (i = aestart; i < aeend; i++) {
898: PetscCall(MatGetRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
899: if (ncols > 0) {
900: /*Ae*/
901: row = i + r1;
902: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
903: /*Ae'*/
904: for (j = 0; j < ncols; j++) {
905: newcol = i + c2;
906: newrow = cols[j];
907: newval = vals[j];
908: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
909: }
910: }
911: PetscCall(MatRestoreRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
912: }
913: }
915: if (ipmP->nb > 0) {
916: PetscCall(MatGetOwnershipRange(ipmP->Ai, &aistart, &aiend));
917: /* Copy Ai,and Ai' */
918: for (i = aistart; i < aiend; i++) {
919: row = i + r2;
920: PetscCall(MatGetRow(ipmP->Ai, i, &ncols, &cols, &vals));
921: if (ncols > 0) {
922: /*Ai*/
923: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
924: /*-Ai'*/
925: for (j = 0; j < ncols; j++) {
926: newcol = i + c3;
927: newrow = cols[j];
928: newval = -vals[j];
929: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
930: }
931: }
932: PetscCall(MatRestoreRow(ipmP->Ai, i, &ncols, &cols, &vals));
933: }
935: /* -I */
936: for (i = kstart; i < kend; i++) {
937: if (i >= r2 && i < r3) {
938: newrow = i;
939: newcol = i - r2 + c1;
940: newval = -1.0;
941: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
942: }
943: }
945: /* Copy L,Y */
946: PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send));
947: PetscCall(VecGetArrayRead(ipmP->lambdai, &l));
948: PetscCall(VecGetArrayRead(ipmP->s, &y));
950: for (i = sstart; i < send; i++) {
951: newcols[0] = c1 + i;
952: newcols[1] = c3 + i;
953: newvals[0] = l[i - sstart];
954: newvals[1] = y[i - sstart];
955: newrow = r3 + i;
956: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 2, newcols, newvals, INSERT_VALUES));
957: }
959: PetscCall(VecRestoreArrayRead(ipmP->lambdai, &l));
960: PetscCall(VecRestoreArrayRead(ipmP->s, &y));
961: }
963: PetscCall(PetscFree(indices));
964: PetscCall(PetscFree(newvals));
965: PetscCall(MatAssemblyBegin(ipmP->K, MAT_FINAL_ASSEMBLY));
966: PetscCall(MatAssemblyEnd(ipmP->K, MAT_FINAL_ASSEMBLY));
967: PetscFunctionReturn(PETSC_SUCCESS);
968: }
970: static PetscErrorCode IPMGatherRHS(Tao tao, Vec RHS, Vec X1, Vec X2, Vec X3, Vec X4)
971: {
972: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
974: PetscFunctionBegin;
975: /* rhs = [x1 (n)
976: x2 (me)
977: x3 (nb)
978: x4 (nb)] */
979: if (X1) {
980: PetscCall(VecScatterBegin(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
981: PetscCall(VecScatterEnd(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
982: }
983: if (ipmP->me > 0 && X2) {
984: PetscCall(VecScatterBegin(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
985: PetscCall(VecScatterEnd(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
986: }
987: if (ipmP->nb > 0) {
988: if (X3) {
989: PetscCall(VecScatterBegin(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
990: PetscCall(VecScatterEnd(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
991: }
992: if (X4) {
993: PetscCall(VecScatterBegin(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
994: PetscCall(VecScatterEnd(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
995: }
996: }
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: static PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1001: {
1002: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1004: PetscFunctionBegin;
1005: CHKMEMQ;
1006: /* [x1 (n)
1007: x2 (nb) may be 0
1008: x3 (me) may be 0
1009: x4 (nb) may be 0 */
1010: if (X1) {
1011: PetscCall(VecScatterBegin(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1012: PetscCall(VecScatterEnd(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1013: }
1014: if (X2 && ipmP->nb > 0) {
1015: PetscCall(VecScatterBegin(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1016: PetscCall(VecScatterEnd(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1017: }
1018: if (X3 && ipmP->me > 0) {
1019: PetscCall(VecScatterBegin(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1020: PetscCall(VecScatterEnd(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1021: }
1022: if (X4 && ipmP->nb > 0) {
1023: PetscCall(VecScatterBegin(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1024: PetscCall(VecScatterEnd(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1025: }
1026: CHKMEMQ;
1027: PetscFunctionReturn(PETSC_SUCCESS);
1028: }
1030: /*MC
1031: TAOIPM - Interior point algorithm for generally constrained optimization.
1033: Options Database Keys:
1034: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1035: - -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1037: Notes:
1038: 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.
1039: Level: beginner
1041: .seealso: `Tao`, `TAOPDIPM`, `TaoType`
1042: M*/
1044: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1045: {
1046: TAO_IPM *ipmP;
1048: PetscFunctionBegin;
1049: tao->ops->setup = TaoSetup_IPM;
1050: tao->ops->solve = TaoSolve_IPM;
1051: tao->ops->view = TaoView_IPM;
1052: tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1053: tao->ops->destroy = TaoDestroy_IPM;
1054: /* tao->ops->computedual = TaoComputeDual_IPM; */
1056: PetscCall(PetscNew(&ipmP));
1057: tao->data = (void *)ipmP;
1059: /* Override default settings (unless already changed) */
1060: if (!tao->max_it_changed) tao->max_it = 200;
1061: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1063: ipmP->dec = 10000; /* line search criteria */
1064: ipmP->taumin = 0.995;
1065: ipmP->monitorkkt = PETSC_FALSE;
1066: ipmP->pushs = 100;
1067: ipmP->pushnu = 100;
1068: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
1069: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
1070: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }