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(VecSet(ipmP->Zero_nb, 0.0));
293: PetscCall(VecDuplicate(ipmP->s, &ipmP->One_nb));
294: PetscCall(VecSet(ipmP->One_nb, 1.0));
295: PetscCall(VecDuplicate(ipmP->s, &ipmP->Inf_nb));
296: PetscCall(VecSet(ipmP->Inf_nb, PETSC_INFINITY));
298: PetscCall(PetscMalloc1(ipmP->nb, &cind));
299: PetscCall(PetscMalloc1(ipmP->mi, &ucind));
300: PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send));
302: if (ipmP->mi > 0) {
303: PetscCall(VecGetOwnershipRange(tao->constraints_inequality, &ucstart, &ucend));
304: counter = 0;
305: for (i = ucstart; i < ucend; i++) cind[counter++] = i;
306: PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isuc));
307: PetscCall(ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isc));
308: PetscCall(VecScatterCreate(tao->constraints_inequality, isuc, ipmP->ci, isc, &ipmP->ci_scat));
310: PetscCall(ISDestroy(&isuc));
311: PetscCall(ISDestroy(&isc));
312: }
313: /* need to know how may xbound indices are on each process */
314: /* TODO better way */
315: if (ipmP->nxlb) {
316: PetscCall(ISAllGather(ipmP->isxl, &bigxl));
317: PetscCall(ISGetIndices(bigxl, &xli));
318: /* find offsets for this processor */
319: xl_offset = ipmP->mi;
320: for (i = 0; i < ipmP->nxlb; i++) {
321: if (xli[i] < xstart) {
322: xl_offset++;
323: } else break;
324: }
325: PetscCall(ISRestoreIndices(bigxl, &xli));
327: PetscCall(ISGetIndices(ipmP->isxl, &xli));
328: PetscCall(ISGetLocalSize(ipmP->isxl, &nloc));
329: for (i = 0; i < nloc; i++) {
330: xind[i] = xli[i];
331: cind[i] = xl_offset + i;
332: }
334: PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx));
335: PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc));
336: PetscCall(VecScatterCreate(tao->XL, isx, ipmP->ci, isc, &ipmP->xl_scat));
337: PetscCall(ISDestroy(&isx));
338: PetscCall(ISDestroy(&isc));
339: PetscCall(ISDestroy(&bigxl));
340: }
342: if (ipmP->nxub) {
343: PetscCall(ISAllGather(ipmP->isxu, &bigxu));
344: PetscCall(ISGetIndices(bigxu, &xui));
345: /* find offsets for this processor */
346: xu_offset = ipmP->mi + ipmP->nxlb;
347: for (i = 0; i < ipmP->nxub; i++) {
348: if (xui[i] < xstart) {
349: xu_offset++;
350: } else break;
351: }
352: PetscCall(ISRestoreIndices(bigxu, &xui));
354: PetscCall(ISGetIndices(ipmP->isxu, &xui));
355: PetscCall(ISGetLocalSize(ipmP->isxu, &nloc));
356: for (i = 0; i < nloc; i++) {
357: xind[i] = xui[i];
358: cind[i] = xu_offset + i;
359: }
361: PetscCall(ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx));
362: PetscCall(ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc));
363: PetscCall(VecScatterCreate(tao->XU, isx, ipmP->ci, isc, &ipmP->xu_scat));
364: PetscCall(ISDestroy(&isx));
365: PetscCall(ISDestroy(&isc));
366: PetscCall(ISDestroy(&bigxu));
367: }
368: }
369: PetscCall(VecCreate(comm, &ipmP->bigrhs));
370: PetscCall(VecGetType(tao->solution, &vtype));
371: PetscCall(VecSetType(ipmP->bigrhs, vtype));
372: PetscCall(VecSetSizes(ipmP->bigrhs, PETSC_DECIDE, bigsize));
373: PetscCall(VecSetFromOptions(ipmP->bigrhs));
374: PetscCall(VecDuplicate(ipmP->bigrhs, &ipmP->bigstep));
376: /* create scatters for step->x and x->rhs */
377: for (i = xstart; i < xend; i++) {
378: stepind[i - xstart] = i;
379: xind[i - xstart] = i;
380: }
381: PetscCall(ISCreateGeneral(comm, xend - xstart, stepind, PETSC_COPY_VALUES, &sis));
382: PetscCall(ISCreateGeneral(comm, xend - xstart, xind, PETSC_COPY_VALUES, &is1));
383: PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->solution, is1, &ipmP->step1));
384: PetscCall(VecScatterCreate(tao->solution, is1, ipmP->bigrhs, sis, &ipmP->rhs1));
385: PetscCall(ISDestroy(&sis));
386: PetscCall(ISDestroy(&is1));
388: if (ipmP->nb > 0) {
389: for (i = sstart; i < send; i++) {
390: stepind[i - sstart] = i + ipmP->n;
391: cind[i - sstart] = i;
392: }
393: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
394: PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1));
395: PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step2));
396: PetscCall(ISDestroy(&sis));
398: for (i = sstart; i < send; i++) {
399: stepind[i - sstart] = i + ipmP->n + ipmP->me;
400: cind[i - sstart] = i;
401: }
402: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
403: PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs3));
404: PetscCall(ISDestroy(&sis));
405: PetscCall(ISDestroy(&is1));
406: }
408: if (ipmP->me > 0) {
409: PetscCall(VecGetOwnershipRange(tao->constraints_equality, &ucestart, &uceend));
410: for (i = ucestart; i < uceend; i++) {
411: stepind[i - ucestart] = i + ipmP->n + ipmP->nb;
412: uceind[i - ucestart] = i;
413: }
415: PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis));
416: PetscCall(ISCreateGeneral(comm, uceend - ucestart, uceind, PETSC_COPY_VALUES, &is1));
417: PetscCall(VecScatterCreate(ipmP->bigstep, sis, tao->constraints_equality, is1, &ipmP->step3));
418: PetscCall(ISDestroy(&sis));
420: for (i = ucestart; i < uceend; i++) stepind[i - ucestart] = i + ipmP->n;
422: PetscCall(ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis));
423: PetscCall(VecScatterCreate(tao->constraints_equality, is1, ipmP->bigrhs, sis, &ipmP->rhs2));
424: PetscCall(ISDestroy(&sis));
425: PetscCall(ISDestroy(&is1));
426: }
428: if (ipmP->nb > 0) {
429: for (i = sstart; i < send; i++) {
430: stepind[i - sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
431: cind[i - sstart] = i;
432: }
433: PetscCall(ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1));
434: PetscCall(ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis));
435: PetscCall(VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step4));
436: PetscCall(VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs4));
437: PetscCall(ISDestroy(&sis));
438: PetscCall(ISDestroy(&is1));
439: }
441: PetscCall(PetscFree(stepind));
442: PetscCall(PetscFree(cind));
443: PetscCall(PetscFree(ucind));
444: PetscCall(PetscFree(uceind));
445: PetscCall(PetscFree(xind));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode TaoDestroy_IPM(Tao tao)
450: {
451: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
453: PetscFunctionBegin;
454: PetscCall(VecDestroy(&ipmP->rd));
455: PetscCall(VecDestroy(&ipmP->rpe));
456: PetscCall(VecDestroy(&ipmP->rpi));
457: PetscCall(VecDestroy(&ipmP->work));
458: PetscCall(VecDestroy(&ipmP->lambdae));
459: PetscCall(VecDestroy(&ipmP->lambdai));
460: PetscCall(VecDestroy(&ipmP->s));
461: PetscCall(VecDestroy(&ipmP->ds));
462: PetscCall(VecDestroy(&ipmP->ci));
464: PetscCall(VecDestroy(&ipmP->rhs_x));
465: PetscCall(VecDestroy(&ipmP->rhs_lambdae));
466: PetscCall(VecDestroy(&ipmP->rhs_lambdai));
467: PetscCall(VecDestroy(&ipmP->rhs_s));
469: PetscCall(VecDestroy(&ipmP->save_x));
470: PetscCall(VecDestroy(&ipmP->save_lambdae));
471: PetscCall(VecDestroy(&ipmP->save_lambdai));
472: PetscCall(VecDestroy(&ipmP->save_s));
474: PetscCall(VecScatterDestroy(&ipmP->step1));
475: PetscCall(VecScatterDestroy(&ipmP->step2));
476: PetscCall(VecScatterDestroy(&ipmP->step3));
477: PetscCall(VecScatterDestroy(&ipmP->step4));
479: PetscCall(VecScatterDestroy(&ipmP->rhs1));
480: PetscCall(VecScatterDestroy(&ipmP->rhs2));
481: PetscCall(VecScatterDestroy(&ipmP->rhs3));
482: PetscCall(VecScatterDestroy(&ipmP->rhs4));
484: PetscCall(VecScatterDestroy(&ipmP->ci_scat));
485: PetscCall(VecScatterDestroy(&ipmP->xl_scat));
486: PetscCall(VecScatterDestroy(&ipmP->xu_scat));
488: PetscCall(VecDestroy(&ipmP->dlambdai));
489: PetscCall(VecDestroy(&ipmP->dlambdae));
490: PetscCall(VecDestroy(&ipmP->Zero_nb));
491: PetscCall(VecDestroy(&ipmP->One_nb));
492: PetscCall(VecDestroy(&ipmP->Inf_nb));
493: PetscCall(VecDestroy(&ipmP->complementarity));
495: PetscCall(VecDestroy(&ipmP->bigrhs));
496: PetscCall(VecDestroy(&ipmP->bigstep));
497: PetscCall(MatDestroy(&ipmP->Ai));
498: PetscCall(MatDestroy(&ipmP->K));
499: PetscCall(ISDestroy(&ipmP->isxu));
500: PetscCall(ISDestroy(&ipmP->isxl));
501: PetscCall(KSPDestroy(&tao->ksp));
502: PetscCall(PetscFree(tao->data));
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode TaoSetFromOptions_IPM(Tao tao, PetscOptionItems PetscOptionsObject)
507: {
508: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
510: PetscFunctionBegin;
511: PetscOptionsHeadBegin(PetscOptionsObject, "IPM method for constrained optimization");
512: PetscCall(PetscOptionsBool("-tao_ipm_monitorkkt", "monitor kkt status", NULL, ipmP->monitorkkt, &ipmP->monitorkkt, NULL));
513: PetscCall(PetscOptionsReal("-tao_ipm_pushs", "parameter to push initial slack variables away from bounds", NULL, ipmP->pushs, &ipmP->pushs, NULL));
514: PetscCall(PetscOptionsReal("-tao_ipm_pushnu", "parameter to push initial (inequality) dual variables away from bounds", NULL, ipmP->pushnu, &ipmP->pushnu, NULL));
515: PetscOptionsHeadEnd();
516: PetscCall(KSPSetFromOptions(tao->ksp));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
521: {
522: return PETSC_SUCCESS;
523: }
525: /* IPMObjectiveAndGradient()
526: f = d'x + 0.5 * x' * H * x
527: rd = H*x + d + Ae'*lame - Ai'*lami
528: rpe = Ae*x - be
529: rpi = Ai*x - yi - bi
530: mu = yi' * lami/mi;
531: com = yi.*lami
533: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
534: */
535: /*
536: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
537: {
538: Tao tao = (Tao)tptr;
539: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
541: PetscFunctionBegin;
542: PetscCall(IPMComputeKKT(tao));
543: *f = ipmP->phi;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
546: */
548: /*
549: f = d'x + 0.5 * x' * H * x
550: rd = H*x + d + Ae'*lame - Ai'*lami
551: Ai = jac_ineq
552: I (w/lb)
553: -I (w/ub)
555: rpe = ce
556: rpi = ci - s;
557: com = s.*lami
558: mu = yi' * lami/mi;
560: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
561: */
562: static PetscErrorCode IPMComputeKKT(Tao tao)
563: {
564: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
565: PetscScalar norm;
567: PetscFunctionBegin;
568: PetscCall(VecCopy(tao->gradient, ipmP->rd));
570: if (ipmP->me > 0) {
571: /* rd = gradient + Ae'*lambdae */
572: PetscCall(MatMultTranspose(tao->jacobian_equality, ipmP->lambdae, ipmP->work));
573: PetscCall(VecAXPY(ipmP->rd, 1.0, ipmP->work));
575: /* rpe = ce(x) */
576: PetscCall(VecCopy(tao->constraints_equality, ipmP->rpe));
577: }
578: if (ipmP->nb > 0) {
579: /* rd = rd - Ai'*lambdai */
580: PetscCall(MatMultTranspose(ipmP->Ai, ipmP->lambdai, ipmP->work));
581: PetscCall(VecAXPY(ipmP->rd, -1.0, ipmP->work));
583: /* rpi = cin - s */
584: PetscCall(VecCopy(ipmP->ci, ipmP->rpi));
585: PetscCall(VecAXPY(ipmP->rpi, -1.0, ipmP->s));
587: /* com = s .* lami */
588: PetscCall(VecPointwiseMult(ipmP->complementarity, ipmP->s, ipmP->lambdai));
589: }
590: /* phi = ||rd; rpe; rpi; com|| */
591: PetscCall(VecDot(ipmP->rd, ipmP->rd, &norm));
592: ipmP->phi = norm;
593: if (ipmP->me > 0) {
594: PetscCall(VecDot(ipmP->rpe, ipmP->rpe, &norm));
595: ipmP->phi += norm;
596: }
597: if (ipmP->nb > 0) {
598: PetscCall(VecDot(ipmP->rpi, ipmP->rpi, &norm));
599: ipmP->phi += norm;
600: PetscCall(VecDot(ipmP->complementarity, ipmP->complementarity, &norm));
601: ipmP->phi += norm;
602: /* mu = s'*lami/nb */
603: PetscCall(VecDot(ipmP->s, ipmP->lambdai, &ipmP->mu));
604: ipmP->mu /= ipmP->nb;
605: } else {
606: ipmP->mu = 1.0;
607: }
609: ipmP->phi = PetscSqrtScalar(ipmP->phi);
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: /* evaluate user info at current point */
614: static PetscErrorCode IPMEvaluate(Tao tao)
615: {
616: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
618: PetscFunctionBegin;
619: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &ipmP->kkt_f, tao->gradient));
620: PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
621: if (ipmP->me > 0) {
622: PetscCall(TaoComputeEqualityConstraints(tao, tao->solution, tao->constraints_equality));
623: PetscCall(TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre));
624: }
625: if (ipmP->mi > 0) {
626: PetscCall(TaoComputeInequalityConstraints(tao, tao->solution, tao->constraints_inequality));
627: PetscCall(TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre));
628: }
629: if (ipmP->nb > 0) {
630: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
631: PetscCall(IPMUpdateAi(tao));
632: }
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /* Push initial point away from bounds */
637: static PetscErrorCode IPMPushInitialPoint(Tao tao)
638: {
639: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
641: PetscFunctionBegin;
642: PetscCall(TaoComputeVariableBounds(tao));
643: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
644: if (ipmP->nb > 0) {
645: PetscCall(VecSet(ipmP->s, ipmP->pushs));
646: PetscCall(VecSet(ipmP->lambdai, ipmP->pushnu));
647: if (ipmP->mi > 0) PetscCall(VecSet(tao->DI, ipmP->pushnu));
648: }
649: if (ipmP->me > 0) {
650: PetscCall(VecSet(tao->DE, 1.0));
651: PetscCall(VecSet(ipmP->lambdae, 1.0));
652: }
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: static PetscErrorCode IPMUpdateAi(Tao tao)
657: {
658: /* Ai = Ji
659: I (w/lb)
660: -I (w/ub) */
662: /* Ci = user->ci
663: Xi - lb (w/lb)
664: -Xi + ub (w/ub) */
666: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
667: MPI_Comm comm;
668: PetscInt i;
669: PetscScalar newval;
670: PetscInt newrow, newcol, ncols;
671: const PetscScalar *vals;
672: const PetscInt *cols;
673: PetscInt astart, aend, jstart, jend;
674: PetscInt *nonzeros;
675: PetscInt r2, r3, r4;
676: PetscMPIInt size;
677: Vec solu;
678: PetscInt nloc;
680: PetscFunctionBegin;
681: r2 = ipmP->mi;
682: r3 = r2 + ipmP->nxlb;
683: r4 = r3 + ipmP->nxub;
685: if (!ipmP->nb) PetscFunctionReturn(PETSC_SUCCESS);
687: /* Create Ai matrix if it doesn't exist yet */
688: if (!ipmP->Ai) {
689: comm = ((PetscObject)tao->solution)->comm;
690: PetscCallMPI(MPI_Comm_size(comm, &size));
691: if (size == 1) {
692: PetscCall(PetscMalloc1(ipmP->nb, &nonzeros));
693: for (i = 0; i < ipmP->mi; i++) {
694: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
695: nonzeros[i] = ncols;
696: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, NULL, NULL));
697: }
698: for (i = r2; i < r4; i++) nonzeros[i] = 1;
699: }
700: PetscCall(MatCreate(comm, &ipmP->Ai));
701: PetscCall(MatSetType(ipmP->Ai, MATAIJ));
703: PetscCall(TaoGetSolution(tao, &solu));
704: PetscCall(VecGetLocalSize(solu, &nloc));
705: PetscCall(MatSetSizes(ipmP->Ai, PETSC_DECIDE, nloc, ipmP->nb, PETSC_DECIDE));
706: PetscCall(MatSetFromOptions(ipmP->Ai));
707: PetscCall(MatMPIAIJSetPreallocation(ipmP->Ai, ipmP->nb, NULL, ipmP->nb, NULL));
708: PetscCall(MatSeqAIJSetPreallocation(ipmP->Ai, PETSC_DEFAULT, nonzeros));
709: if (size == 1) PetscCall(PetscFree(nonzeros));
710: }
712: /* Copy values from user jacobian to Ai */
713: PetscCall(MatGetOwnershipRange(ipmP->Ai, &astart, &aend));
715: /* Ai w/lb */
716: if (ipmP->mi) {
717: PetscCall(MatZeroEntries(ipmP->Ai));
718: PetscCall(MatGetOwnershipRange(tao->jacobian_inequality, &jstart, &jend));
719: for (i = jstart; i < jend; i++) {
720: PetscCall(MatGetRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
721: newrow = i;
722: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, ncols, cols, vals, INSERT_VALUES));
723: PetscCall(MatRestoreRow(tao->jacobian_inequality, i, &ncols, &cols, &vals));
724: }
725: }
727: /* I w/ xlb */
728: if (ipmP->nxlb) {
729: for (i = 0; i < ipmP->nxlb; i++) {
730: if (i >= astart && i < aend) {
731: newrow = i + r2;
732: newcol = i;
733: newval = 1.0;
734: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
735: }
736: }
737: }
738: if (ipmP->nxub) {
739: /* I w/ xub */
740: for (i = 0; i < ipmP->nxub; i++) {
741: if (i >= astart && i < aend) {
742: newrow = i + r3;
743: newcol = i;
744: newval = -1.0;
745: PetscCall(MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
746: }
747: }
748: }
750: PetscCall(MatAssemblyBegin(ipmP->Ai, MAT_FINAL_ASSEMBLY));
751: PetscCall(MatAssemblyEnd(ipmP->Ai, MAT_FINAL_ASSEMBLY));
752: CHKMEMQ;
754: PetscCall(VecSet(ipmP->ci, 0.0));
756: /* user ci */
757: if (ipmP->mi > 0) {
758: PetscCall(VecScatterBegin(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
759: PetscCall(VecScatterEnd(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
760: }
761: if (!ipmP->work) PetscCall(VecDuplicate(tao->solution, &ipmP->work));
762: PetscCall(VecCopy(tao->solution, ipmP->work));
763: if (tao->XL) {
764: PetscCall(VecAXPY(ipmP->work, -1.0, tao->XL));
766: /* lower bounds on variables */
767: if (ipmP->nxlb > 0) {
768: PetscCall(VecScatterBegin(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
769: PetscCall(VecScatterEnd(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
770: }
771: }
772: if (tao->XU) {
773: /* upper bounds on variables */
774: PetscCall(VecCopy(tao->solution, ipmP->work));
775: PetscCall(VecScale(ipmP->work, -1.0));
776: PetscCall(VecAXPY(ipmP->work, 1.0, tao->XU));
777: if (ipmP->nxub > 0) {
778: PetscCall(VecScatterBegin(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
779: PetscCall(VecScatterEnd(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD));
780: }
781: }
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
785: /* create K = [ Hlag , 0 , Ae', -Ai'];
786: [Ae , 0, 0 , 0];
787: [Ai ,-I, 0 , 0];
788: [ 0 , S , 0, Y ]; */
789: static PetscErrorCode IPMUpdateK(Tao tao)
790: {
791: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
792: MPI_Comm comm;
793: PetscMPIInt size;
794: PetscInt i, j, row;
795: PetscInt ncols, newcol, newcols[2], newrow;
796: const PetscInt *cols;
797: const PetscReal *vals;
798: const PetscReal *l, *y;
799: PetscReal *newvals;
800: PetscReal newval;
801: PetscInt subsize;
802: const PetscInt *indices;
803: PetscInt *nonzeros, *d_nonzeros, *o_nonzeros;
804: PetscInt bigsize;
805: PetscInt r1, r2, r3;
806: PetscInt c1, c2, c3;
807: PetscInt klocalsize;
808: PetscInt hstart, hend, kstart, kend;
809: PetscInt aistart, aiend, aestart, aeend;
810: PetscInt sstart, send;
812: PetscFunctionBegin;
813: comm = ((PetscObject)tao->solution)->comm;
814: PetscCallMPI(MPI_Comm_size(comm, &size));
815: PetscCall(IPMUpdateAi(tao));
817: /* allocate workspace */
818: subsize = PetscMax(ipmP->n, ipmP->nb);
819: subsize = PetscMax(ipmP->me, subsize);
820: subsize = PetscMax(2, subsize);
821: PetscCall(PetscMalloc1(subsize, &indices));
822: PetscCall(PetscMalloc1(subsize, &newvals));
824: r1 = c1 = ipmP->n;
825: r2 = r1 + ipmP->me;
826: c2 = c1 + ipmP->nb;
827: r3 = c3 = r2 + ipmP->nb;
829: bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
830: PetscCall(VecGetOwnershipRange(ipmP->bigrhs, &kstart, &kend));
831: PetscCall(MatGetOwnershipRange(tao->hessian, &hstart, &hend));
832: klocalsize = kend - kstart;
833: if (!ipmP->K) {
834: if (size == 1) {
835: PetscCall(PetscMalloc1(kend - kstart, &nonzeros));
836: for (i = 0; i < bigsize; i++) {
837: if (i < r1) {
838: PetscCall(MatGetRow(tao->hessian, i, &ncols, NULL, NULL));
839: nonzeros[i] = ncols;
840: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, NULL, NULL));
841: nonzeros[i] += ipmP->me + ipmP->nb;
842: } else if (i < r2) {
843: nonzeros[i - kstart] = ipmP->n;
844: } else if (i < r3) {
845: nonzeros[i - kstart] = ipmP->n + 1;
846: } else if (i < bigsize) {
847: nonzeros[i - kstart] = 2;
848: }
849: }
850: PetscCall(MatCreate(comm, &ipmP->K));
851: PetscCall(MatSetType(ipmP->K, MATSEQAIJ));
852: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
853: PetscCall(MatSeqAIJSetPreallocation(ipmP->K, 0, nonzeros));
854: PetscCall(MatSetFromOptions(ipmP->K));
855: PetscCall(PetscFree(nonzeros));
856: } else {
857: PetscCall(PetscMalloc1(kend - kstart, &d_nonzeros));
858: PetscCall(PetscMalloc1(kend - kstart, &o_nonzeros));
859: for (i = kstart; i < kend; i++) {
860: if (i < r1) {
861: /* TODO fix preallocation for mpi mats */
862: d_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, kend - kstart);
863: o_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, bigsize - (kend - kstart));
864: } else if (i < r2) {
865: d_nonzeros[i - kstart] = PetscMin(ipmP->n, kend - kstart);
866: o_nonzeros[i - kstart] = PetscMin(ipmP->n, bigsize - (kend - kstart));
867: } else if (i < r3) {
868: d_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, kend - kstart);
869: o_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, bigsize - (kend - kstart));
870: } else {
871: d_nonzeros[i - kstart] = PetscMin(2, kend - kstart);
872: o_nonzeros[i - kstart] = PetscMin(2, bigsize - (kend - kstart));
873: }
874: }
875: PetscCall(MatCreate(comm, &ipmP->K));
876: PetscCall(MatSetType(ipmP->K, MATMPIAIJ));
877: PetscCall(MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE));
878: PetscCall(MatMPIAIJSetPreallocation(ipmP->K, 0, d_nonzeros, 0, o_nonzeros));
879: PetscCall(PetscFree(d_nonzeros));
880: PetscCall(PetscFree(o_nonzeros));
881: PetscCall(MatSetFromOptions(ipmP->K));
882: }
883: }
885: PetscCall(MatZeroEntries(ipmP->K));
886: /* Copy H */
887: for (i = hstart; i < hend; i++) {
888: PetscCall(MatGetRow(tao->hessian, i, &ncols, &cols, &vals));
889: if (ncols > 0) PetscCall(MatSetValues(ipmP->K, 1, &i, ncols, cols, vals, INSERT_VALUES));
890: PetscCall(MatRestoreRow(tao->hessian, i, &ncols, &cols, &vals));
891: }
893: /* Copy Ae and Ae' */
894: if (ipmP->me > 0) {
895: PetscCall(MatGetOwnershipRange(tao->jacobian_equality, &aestart, &aeend));
896: for (i = aestart; i < aeend; i++) {
897: PetscCall(MatGetRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
898: if (ncols > 0) {
899: /*Ae*/
900: row = i + r1;
901: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
902: /*Ae'*/
903: for (j = 0; j < ncols; j++) {
904: newcol = i + c2;
905: newrow = cols[j];
906: newval = vals[j];
907: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
908: }
909: }
910: PetscCall(MatRestoreRow(tao->jacobian_equality, i, &ncols, &cols, &vals));
911: }
912: }
914: if (ipmP->nb > 0) {
915: PetscCall(MatGetOwnershipRange(ipmP->Ai, &aistart, &aiend));
916: /* Copy Ai,and Ai' */
917: for (i = aistart; i < aiend; i++) {
918: row = i + r2;
919: PetscCall(MatGetRow(ipmP->Ai, i, &ncols, &cols, &vals));
920: if (ncols > 0) {
921: /*Ai*/
922: PetscCall(MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES));
923: /*-Ai'*/
924: for (j = 0; j < ncols; j++) {
925: newcol = i + c3;
926: newrow = cols[j];
927: newval = -vals[j];
928: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
929: }
930: }
931: PetscCall(MatRestoreRow(ipmP->Ai, i, &ncols, &cols, &vals));
932: }
934: /* -I */
935: for (i = kstart; i < kend; i++) {
936: if (i >= r2 && i < r3) {
937: newrow = i;
938: newcol = i - r2 + c1;
939: newval = -1.0;
940: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES));
941: }
942: }
944: /* Copy L,Y */
945: PetscCall(VecGetOwnershipRange(ipmP->s, &sstart, &send));
946: PetscCall(VecGetArrayRead(ipmP->lambdai, &l));
947: PetscCall(VecGetArrayRead(ipmP->s, &y));
949: for (i = sstart; i < send; i++) {
950: newcols[0] = c1 + i;
951: newcols[1] = c3 + i;
952: newvals[0] = l[i - sstart];
953: newvals[1] = y[i - sstart];
954: newrow = r3 + i;
955: PetscCall(MatSetValues(ipmP->K, 1, &newrow, 2, newcols, newvals, INSERT_VALUES));
956: }
958: PetscCall(VecRestoreArrayRead(ipmP->lambdai, &l));
959: PetscCall(VecRestoreArrayRead(ipmP->s, &y));
960: }
962: PetscCall(PetscFree(indices));
963: PetscCall(PetscFree(newvals));
964: PetscCall(MatAssemblyBegin(ipmP->K, MAT_FINAL_ASSEMBLY));
965: PetscCall(MatAssemblyEnd(ipmP->K, MAT_FINAL_ASSEMBLY));
966: PetscFunctionReturn(PETSC_SUCCESS);
967: }
969: static PetscErrorCode IPMGatherRHS(Tao tao, Vec RHS, Vec X1, Vec X2, Vec X3, Vec X4)
970: {
971: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
973: PetscFunctionBegin;
974: /* rhs = [x1 (n)
975: x2 (me)
976: x3 (nb)
977: x4 (nb)] */
978: if (X1) {
979: PetscCall(VecScatterBegin(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
980: PetscCall(VecScatterEnd(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD));
981: }
982: if (ipmP->me > 0 && X2) {
983: PetscCall(VecScatterBegin(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
984: PetscCall(VecScatterEnd(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD));
985: }
986: if (ipmP->nb > 0) {
987: if (X3) {
988: PetscCall(VecScatterBegin(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
989: PetscCall(VecScatterEnd(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD));
990: }
991: if (X4) {
992: PetscCall(VecScatterBegin(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
993: PetscCall(VecScatterEnd(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD));
994: }
995: }
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: static PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1000: {
1001: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1003: PetscFunctionBegin;
1004: CHKMEMQ;
1005: /* [x1 (n)
1006: x2 (nb) may be 0
1007: x3 (me) may be 0
1008: x4 (nb) may be 0 */
1009: if (X1) {
1010: PetscCall(VecScatterBegin(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1011: PetscCall(VecScatterEnd(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD));
1012: }
1013: if (X2 && ipmP->nb > 0) {
1014: PetscCall(VecScatterBegin(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1015: PetscCall(VecScatterEnd(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD));
1016: }
1017: if (X3 && ipmP->me > 0) {
1018: PetscCall(VecScatterBegin(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1019: PetscCall(VecScatterEnd(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD));
1020: }
1021: if (X4 && ipmP->nb > 0) {
1022: PetscCall(VecScatterBegin(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1023: PetscCall(VecScatterEnd(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD));
1024: }
1025: CHKMEMQ;
1026: PetscFunctionReturn(PETSC_SUCCESS);
1027: }
1029: /*MC
1030: TAOIPM - Interior point algorithm for generally constrained optimization.
1032: Options Database Keys:
1033: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1034: - -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1036: Notes:
1037: 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.
1038: Level: beginner
1040: .seealso: `Tao`, `TAOPDIPM`, `TaoType`
1041: M*/
1043: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1044: {
1045: TAO_IPM *ipmP;
1047: PetscFunctionBegin;
1048: tao->ops->setup = TaoSetup_IPM;
1049: tao->ops->solve = TaoSolve_IPM;
1050: tao->ops->view = TaoView_IPM;
1051: tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1052: tao->ops->destroy = TaoDestroy_IPM;
1053: tao->uses_gradient = PETSC_TRUE;
1054: tao->uses_hessian_matrices = PETSC_TRUE;
1055: /* tao->ops->computedual = TaoComputeDual_IPM; */
1057: PetscCall(PetscNew(&ipmP));
1058: tao->data = (void *)ipmP;
1060: /* Override default settings (unless already changed) */
1061: PetscCall(TaoParametersInitialize(tao));
1062: PetscObjectParameterSetDefault(tao, max_it, 200);
1063: PetscObjectParameterSetDefault(tao, max_funcs, 500);
1065: ipmP->dec = 10000; /* line search criteria */
1066: ipmP->taumin = 0.995;
1067: ipmP->monitorkkt = PETSC_FALSE;
1068: ipmP->pushs = 100;
1069: ipmP->pushnu = 100;
1070: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
1071: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
1072: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
1073: PetscFunctionReturn(PETSC_SUCCESS);
1074: }