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;
543: PetscFunctionBegin;
544: PetscCall(IPMComputeKKT(tao));
545: *f = ipmP->phi;
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
548: */
550: /*
551: f = d'x + 0.5 * x' * H * x
552: rd = H*x + d + Ae'*lame - Ai'*lami
553: Ai = jac_ineq
554: I (w/lb)
555: -I (w/ub)
557: rpe = ce
558: rpi = ci - s;
559: com = s.*lami
560: mu = yi' * lami/mi;
562: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
563: */
564: static PetscErrorCode IPMComputeKKT(Tao tao)
565: {
566: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
567: PetscScalar norm;
569: PetscFunctionBegin;
570: PetscCall(VecCopy(tao->gradient, ipmP->rd));
572: if (ipmP->me > 0) {
573: /* rd = gradient + Ae'*lambdae */
574: PetscCall(MatMultTranspose(tao->jacobian_equality, ipmP->lambdae, ipmP->work));
575: PetscCall(VecAXPY(ipmP->rd, 1.0, ipmP->work));
577: /* rpe = ce(x) */
578: PetscCall(VecCopy(tao->constraints_equality, ipmP->rpe));
579: }
580: if (ipmP->nb > 0) {
581: /* rd = rd - Ai'*lambdai */
582: PetscCall(MatMultTranspose(ipmP->Ai, ipmP->lambdai, ipmP->work));
583: PetscCall(VecAXPY(ipmP->rd, -1.0, ipmP->work));
585: /* rpi = cin - s */
586: PetscCall(VecCopy(ipmP->ci, ipmP->rpi));
587: PetscCall(VecAXPY(ipmP->rpi, -1.0, ipmP->s));
589: /* com = s .* lami */
590: PetscCall(VecPointwiseMult(ipmP->complementarity, ipmP->s, ipmP->lambdai));
591: }
592: /* phi = ||rd; rpe; rpi; com|| */
593: PetscCall(VecDot(ipmP->rd, ipmP->rd, &norm));
594: ipmP->phi = norm;
595: if (ipmP->me > 0) {
596: PetscCall(VecDot(ipmP->rpe, ipmP->rpe, &norm));
597: ipmP->phi += norm;
598: }
599: if (ipmP->nb > 0) {
600: PetscCall(VecDot(ipmP->rpi, ipmP->rpi, &norm));
601: ipmP->phi += norm;
602: PetscCall(VecDot(ipmP->complementarity, ipmP->complementarity, &norm));
603: ipmP->phi += norm;
604: /* mu = s'*lami/nb */
605: PetscCall(VecDot(ipmP->s, ipmP->lambdai, &ipmP->mu));
606: ipmP->mu /= ipmP->nb;
607: } else {
608: ipmP->mu = 1.0;
609: }
611: ipmP->phi = PetscSqrtScalar(ipmP->phi);
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /* evaluate user info at current point */
616: static PetscErrorCode IPMEvaluate(Tao tao)
617: {
618: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
620: PetscFunctionBegin;
621: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &ipmP->kkt_f, tao->gradient));
622: PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
623: if (ipmP->me > 0) {
624: PetscCall(TaoComputeEqualityConstraints(tao, tao->solution, tao->constraints_equality));
625: PetscCall(TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre));
626: }
627: if (ipmP->mi > 0) {
628: PetscCall(TaoComputeInequalityConstraints(tao, tao->solution, tao->constraints_inequality));
629: PetscCall(TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre));
630: }
631: if (ipmP->nb > 0) {
632: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
633: PetscCall(IPMUpdateAi(tao));
634: }
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /* Push initial point away from bounds */
639: static PetscErrorCode IPMPushInitialPoint(Tao tao)
640: {
641: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
643: PetscFunctionBegin;
644: PetscCall(TaoComputeVariableBounds(tao));
645: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
646: if (ipmP->nb > 0) {
647: PetscCall(VecSet(ipmP->s, ipmP->pushs));
648: PetscCall(VecSet(ipmP->lambdai, ipmP->pushnu));
649: if (ipmP->mi > 0) PetscCall(VecSet(tao->DI, ipmP->pushnu));
650: }
651: if (ipmP->me > 0) {
652: PetscCall(VecSet(tao->DE, 1.0));
653: PetscCall(VecSet(ipmP->lambdae, 1.0));
654: }
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: static PetscErrorCode IPMUpdateAi(Tao tao)
659: {
660: /* Ai = Ji
661: I (w/lb)
662: -I (w/ub) */
664: /* Ci = user->ci
665: Xi - lb (w/lb)
666: -Xi + ub (w/ub) */
668: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
669: MPI_Comm comm;
670: PetscInt i;
671: PetscScalar newval;
672: PetscInt newrow, newcol, ncols;
673: const PetscScalar *vals;
674: const PetscInt *cols;
675: PetscInt astart, aend, jstart, jend;
676: PetscInt *nonzeros;
677: PetscInt r2, r3, r4;
678: PetscMPIInt size;
679: Vec solu;
680: PetscInt nloc;
682: PetscFunctionBegin;
683: r2 = ipmP->mi;
684: r3 = r2 + ipmP->nxlb;
685: r4 = r3 + ipmP->nxub;
687: if (!ipmP->nb) PetscFunctionReturn(PETSC_SUCCESS);
689: /* Create Ai matrix if it doesn't exist yet */
690: if (!ipmP->Ai) {
691: comm = ((PetscObject)tao->solution)->comm;
692: PetscCallMPI(MPI_Comm_size(comm, &size));
693: if (size == 1) {
694: PetscCall(PetscMalloc1(ipmP->nb, &nonzeros));
695: for (i = 0; i < ipmP->mi; i++) {
696: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
697: nonzeros[i] = ncols;
698: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
699: }
700: for (i = r2; i < r4; i++) nonzeros[i] = 1;
701: }
702: PetscCall(MatCreate(comm, &ipmP->Ai));
703: PetscCall(MatSetType(ipmP->Ai, MATAIJ));
705: PetscCall(TaoGetSolution(tao, &solu));
706: PetscCall(VecGetLocalSize(solu, &nloc));
707: PetscCall(MatSetSizes(ipmP->Ai, PETSC_DECIDE, nloc, ipmP->nb, PETSC_DECIDE));
708: PetscCall(MatSetFromOptions(ipmP->Ai));
709: PetscCall(MatMPIAIJSetPreallocation(ipmP->Ai, ipmP->nb, NULL, ipmP->nb, NULL));
710: PetscCall(MatSeqAIJSetPreallocation(ipmP->Ai, PETSC_DEFAULT, nonzeros));
711: if (size == 1) PetscCall(PetscFree(nonzeros));
712: }
714: /* Copy values from user jacobian to Ai */
715: PetscCall(MatGetOwnershipRange(ipmP->Ai, &astart, &aend));
717: /* Ai w/lb */
718: if (ipmP->mi) {
719: PetscCall(MatZeroEntries(ipmP->Ai));
720: PetscCall(MatGetOwnershipRange(tao->jacobian_inequality, &jstart, &jend));
721: for (i = jstart; i < jend; i++) {
722: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
723: newrow = i;
724: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, ncols, cols, vals, INSERT_VALUES));
725: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
726: }
727: }
729: /* I w/ xlb */
730: if (ipmP->nxlb) {
731: for (i = 0; i < ipmP->nxlb; i++) {
732: if (i >= astart && i < aend) {
733: newrow = i + r2;
734: newcol = i;
735: newval = 1.0;
736: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
737: }
738: }
739: }
740: if (ipmP->nxub) {
741: /* I w/ xub */
742: for (i = 0; i < ipmP->nxub; i++) {
743: if (i >= astart && i < aend) {
744: newrow = i + r3;
745: newcol = i;
746: newval = -1.0;
747: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
748: }
749: }
750: }
752: PetscCall(MatAssemblyBegin(ipmP->Ai, MAT_FINAL_ASSEMBLY));
753: PetscCall(MatAssemblyEnd(ipmP->Ai, MAT_FINAL_ASSEMBLY));
754: CHKMEMQ;
756: PetscCall(VecSet(ipmP->ci, 0.0));
758: /* user ci */
759: if (ipmP->mi > 0) {
760: PetscCall(VecScatterBegin(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
761: PetscCall(VecScatterEnd(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
762: }
763: if (!ipmP->work) PetscCall(VecDuplicate(tao->solution, &ipmP->work));
764: PetscCall(VecCopy(tao->solution, ipmP->work));
765: if (tao->XL) {
766: PetscCall(VecAXPY(ipmP->work, -1.0, tao->XL));
768: /* lower bounds on variables */
769: if (ipmP->nxlb > 0) {
770: PetscCall(VecScatterBegin(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
771: PetscCall(VecScatterEnd(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
772: }
773: }
774: if (tao->XU) {
775: /* upper bounds on variables */
776: PetscCall(VecCopy(tao->solution, ipmP->work));
777: PetscCall(VecScale(ipmP->work, -1.0));
778: PetscCall(VecAXPY(ipmP->work, 1.0, tao->XU));
779: if (ipmP->nxub > 0) {
780: PetscCall(VecScatterBegin(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
781: PetscCall(VecScatterEnd(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
782: }
783: }
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: /* create K = [ Hlag , 0 , Ae', -Ai'];
788: [Ae , 0, 0 , 0];
789: [Ai ,-I, 0 , 0];
790: [ 0 , S , 0, Y ]; */
791: static PetscErrorCode IPMUpdateK(Tao tao)
792: {
793: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
794: MPI_Comm comm;
795: PetscMPIInt size;
796: PetscInt i, j, row;
797: PetscInt ncols, newcol, newcols[2], newrow;
798: const PetscInt *cols;
799: const PetscReal *vals;
800: const PetscReal *l, *y;
801: PetscReal *newvals;
802: PetscReal newval;
803: PetscInt subsize;
804: const PetscInt *indices;
805: PetscInt *nonzeros, *d_nonzeros, *o_nonzeros;
806: PetscInt bigsize;
807: PetscInt r1, r2, r3;
808: PetscInt c1, c2, c3;
809: PetscInt klocalsize;
810: PetscInt hstart, hend, kstart, kend;
811: PetscInt aistart, aiend, aestart, aeend;
812: PetscInt sstart, send;
814: PetscFunctionBegin;
815: comm = ((PetscObject)tao->solution)->comm;
816: PetscCallMPI(MPI_Comm_size(comm, &size));
817: PetscCall(IPMUpdateAi(tao));
819: /* allocate workspace */
820: subsize = PetscMax(ipmP->n, ipmP->nb);
821: subsize = PetscMax(ipmP->me, subsize);
822: subsize = PetscMax(2, subsize);
823: PetscCall(PetscMalloc1(subsize, (PetscInt **)&indices));
824: PetscCall(PetscMalloc1(subsize, &newvals));
826: r1 = c1 = ipmP->n;
827: r2 = r1 + ipmP->me;
828: c2 = c1 + ipmP->nb;
829: r3 = c3 = r2 + ipmP->nb;
831: bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
832: PetscCall(VecGetOwnershipRange(ipmP->bigrhs, &kstart, &kend));
833: PetscCall(MatGetOwnershipRange(tao->hessian, &hstart, &hend));
834: klocalsize = kend - kstart;
835: if (!ipmP->K) {
836: if (size == 1) {
837: PetscCall(PetscMalloc1(kend - kstart, &nonzeros));
838: for (i = 0; i < bigsize; i++) {
839: if (i < r1) {
840: PetscCall(MatGetRow(tao->hessian, i, &ncols, NULL, NULL));
841: nonzeros[i] = ncols;
842: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, NULL, NULL));
843: nonzeros[i] += ipmP->me + ipmP->nb;
844: } else if (i < r2) {
845: nonzeros[i - kstart] = ipmP->n;
846: } else if (i < r3) {
847: nonzeros[i - kstart] = ipmP->n + 1;
848: } else if (i < bigsize) {
849: nonzeros[i - kstart] = 2;
850: }
851: }
852: PetscCall(MatCreate(comm, &ipmP->K));
853: PetscCall(MatSetType(ipmP->K, MATSEQAIJ));
854: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
855: PetscCall(MatSeqAIJSetPreallocation(ipmP->K, 0, nonzeros));
856: PetscCall(MatSetFromOptions(ipmP->K));
857: PetscCall(PetscFree(nonzeros));
858: } else {
859: PetscCall(PetscMalloc1(kend - kstart, &d_nonzeros));
860: PetscCall(PetscMalloc1(kend - kstart, &o_nonzeros));
861: for (i = kstart; i < kend; i++) {
862: if (i < r1) {
863: /* TODO fix preallocation for mpi mats */
864: d_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, kend - kstart);
865: o_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, bigsize - (kend - kstart));
866: } else if (i < r2) {
867: d_nonzeros[i - kstart] = PetscMin(ipmP->n, kend - kstart);
868: o_nonzeros[i - kstart] = PetscMin(ipmP->n, bigsize - (kend - kstart));
869: } else if (i < r3) {
870: d_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, kend - kstart);
871: o_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, bigsize - (kend - kstart));
872: } else {
873: d_nonzeros[i - kstart] = PetscMin(2, kend - kstart);
874: o_nonzeros[i - kstart] = PetscMin(2, bigsize - (kend - kstart));
875: }
876: }
877: PetscCall(MatCreate(comm, &ipmP->K));
878: PetscCall(MatSetType(ipmP->K, MATMPIAIJ));
879: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
880: PetscCall(MatMPIAIJSetPreallocation(ipmP->K, 0, d_nonzeros, 0, o_nonzeros));
881: PetscCall(PetscFree(d_nonzeros));
882: PetscCall(PetscFree(o_nonzeros));
883: PetscCall(MatSetFromOptions(ipmP->K));
884: }
885: }
887: PetscCall(MatZeroEntries(ipmP->K));
888: /* Copy H */
889: for (i = hstart; i < hend; i++) {
890: PetscCall(MatGetRow(tao->hessian, i, &ncols, &cols, &vals));
891: if (ncols > 0) PetscCall(MatSetValues(ipmP->K, 1, &i, ncols, cols, vals, INSERT_VALUES));
892: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, &cols, &vals));
893: }
895: /* Copy Ae and Ae' */
896: if (ipmP->me > 0) {
897: PetscCall(MatGetOwnershipRange(tao->jacobian_equality, &aestart, &aeend));
898: for (i = aestart; i < aeend; i++) {
899: PetscCall(MatGetRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
900: if (ncols > 0) {
901: /*Ae*/
902: row = i + r1;
903: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
904: /*Ae'*/
905: for (j = 0; j < ncols; j++) {
906: newcol = i + c2;
907: newrow = cols[j];
908: newval = vals[j];
909: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
910: }
911: }
912: PetscCall(MatRestoreRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
913: }
914: }
916: if (ipmP->nb > 0) {
917: PetscCall(MatGetOwnershipRange(ipmP->Ai, &aistart, &aiend));
918: /* Copy Ai,and Ai' */
919: for (i = aistart; i < aiend; i++) {
920: row = i + r2;
921: PetscCall(MatGetRow(ipmP->Ai, i, &ncols, &cols, &vals));
922: if (ncols > 0) {
923: /*Ai*/
924: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
925: /*-Ai'*/
926: for (j = 0; j < ncols; j++) {
927: newcol = i + c3;
928: newrow = cols[j];
929: newval = -vals[j];
930: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
931: }
932: }
933: PetscCall(MatRestoreRow(ipmP->Ai, i, &ncols, &cols, &vals));
934: }
936: /* -I */
937: for (i = kstart; i < kend; i++) {
938: if (i >= r2 && i < r3) {
939: newrow = i;
940: newcol = i - r2 + c1;
941: newval = -1.0;
942: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
943: }
944: }
946: /* Copy L,Y */
947: PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send));
948: PetscCall(VecGetArrayRead(ipmP->lambdai, &l));
949: PetscCall(VecGetArrayRead(ipmP->s, &y));
951: for (i = sstart; i < send; i++) {
952: newcols[0] = c1 + i;
953: newcols[1] = c3 + i;
954: newvals[0] = l[i - sstart];
955: newvals[1] = y[i - sstart];
956: newrow = r3 + i;
957: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 2, newcols, newvals, INSERT_VALUES));
958: }
960: PetscCall(VecRestoreArrayRead(ipmP->lambdai, &l));
961: PetscCall(VecRestoreArrayRead(ipmP->s, &y));
962: }
964: PetscCall(PetscFree(indices));
965: PetscCall(PetscFree(newvals));
966: PetscCall(MatAssemblyBegin(ipmP->K, MAT_FINAL_ASSEMBLY));
967: PetscCall(MatAssemblyEnd(ipmP->K, MAT_FINAL_ASSEMBLY));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: static PetscErrorCode IPMGatherRHS(Tao tao, Vec RHS, Vec X1, Vec X2, Vec X3, Vec X4)
972: {
973: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
975: PetscFunctionBegin;
976: /* rhs = [x1 (n)
977: x2 (me)
978: x3 (nb)
979: x4 (nb)] */
980: if (X1) {
981: PetscCall(VecScatterBegin(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
982: PetscCall(VecScatterEnd(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
983: }
984: if (ipmP->me > 0 && X2) {
985: PetscCall(VecScatterBegin(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
986: PetscCall(VecScatterEnd(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
987: }
988: if (ipmP->nb > 0) {
989: if (X3) {
990: PetscCall(VecScatterBegin(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
991: PetscCall(VecScatterEnd(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
992: }
993: if (X4) {
994: PetscCall(VecScatterBegin(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
995: PetscCall(VecScatterEnd(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
996: }
997: }
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: static PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1002: {
1003: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1005: PetscFunctionBegin;
1006: CHKMEMQ;
1007: /* [x1 (n)
1008: x2 (nb) may be 0
1009: x3 (me) may be 0
1010: x4 (nb) may be 0 */
1011: if (X1) {
1012: PetscCall(VecScatterBegin(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1013: PetscCall(VecScatterEnd(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1014: }
1015: if (X2 && ipmP->nb > 0) {
1016: PetscCall(VecScatterBegin(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1017: PetscCall(VecScatterEnd(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1018: }
1019: if (X3 && ipmP->me > 0) {
1020: PetscCall(VecScatterBegin(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1021: PetscCall(VecScatterEnd(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1022: }
1023: if (X4 && ipmP->nb > 0) {
1024: PetscCall(VecScatterBegin(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1025: PetscCall(VecScatterEnd(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1026: }
1027: CHKMEMQ;
1028: PetscFunctionReturn(PETSC_SUCCESS);
1029: }
1031: /*MC
1032: TAOIPM - Interior point algorithm for generally constrained optimization.
1034: Options Database Keys:
1035: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1036: - -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1038: Notes:
1039: 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.
1040: Level: beginner
1042: .seealso: `Tao`, `TAOPDIPM`, `TaoType`
1043: M*/
1045: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1046: {
1047: TAO_IPM *ipmP;
1049: PetscFunctionBegin;
1050: tao->ops->setup = TaoSetup_IPM;
1051: tao->ops->solve = TaoSolve_IPM;
1052: tao->ops->view = TaoView_IPM;
1053: tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1054: tao->ops->destroy = TaoDestroy_IPM;
1055: /* tao->ops->computedual = TaoComputeDual_IPM; */
1057: PetscCall(PetscNew(&ipmP));
1058: tao->data = (void *)ipmP;
1060: /* Override default settings (unless already changed) */
1061: if (!tao->max_it_changed) tao->max_it = 200;
1062: if (!tao->max_funcs_changed) 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: }