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