1 #include <petsctaolinesearch.h>
2 #include <../src/tao/constrained/impls/ipm/pdipm.h>
3 #include <petscsnes.h>
4
5 /*
6 TaoPDIPMEvaluateFunctionsAndJacobians - Evaluate the objective function f, gradient fx, constraints, and all the Jacobians at current vector
7
8 Collective on tao
9
10 Input Parameter:
11 + tao - solver context
12 - x - vector at which all objects to be evaluated
13
14 Level: beginner
15
16 .seealso: TaoPDIPMUpdateConstraints(), TaoPDIPMSetUpBounds()
17 */
TaoPDIPMEvaluateFunctionsAndJacobians(Tao tao,Vec x)18 PetscErrorCode TaoPDIPMEvaluateFunctionsAndJacobians(Tao tao,Vec x)
19 {
20 PetscErrorCode ierr;
21 TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
22
23 PetscFunctionBegin;
24 /* Compute user objective function and gradient */
25 ierr = TaoComputeObjectiveAndGradient(tao,x,&pdipm->obj,tao->gradient);CHKERRQ(ierr);
26
27 /* Equality constraints and Jacobian */
28 if (pdipm->Ng) {
29 ierr = TaoComputeEqualityConstraints(tao,x,tao->constraints_equality);CHKERRQ(ierr);
30 ierr = TaoComputeJacobianEquality(tao,x,tao->jacobian_equality,tao->jacobian_equality_pre);CHKERRQ(ierr);
31 }
32
33 /* Inequality constraints and Jacobian */
34 if (pdipm->Nh) {
35 ierr = TaoComputeInequalityConstraints(tao,x,tao->constraints_inequality);CHKERRQ(ierr);
36 ierr = TaoComputeJacobianInequality(tao,x,tao->jacobian_inequality,tao->jacobian_inequality_pre);CHKERRQ(ierr);
37 }
38 PetscFunctionReturn(0);
39 }
40
41 /*
42 TaoPDIPMUpdateConstraints - Update the vectors ce and ci at x
43
44 Collective
45
46 Input Parameter:
47 + tao - Tao context
48 - x - vector at which constraints to be evaluted
49
50 Level: beginner
51
52 .seealso: TaoPDIPMEvaluateFunctionsAndJacobians()
53 */
TaoPDIPMUpdateConstraints(Tao tao,Vec x)54 PetscErrorCode TaoPDIPMUpdateConstraints(Tao tao,Vec x)
55 {
56 PetscErrorCode ierr;
57 TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
58 PetscInt i,offset,offset1,k,xstart;
59 PetscScalar *carr;
60 const PetscInt *ubptr,*lbptr,*bxptr,*fxptr;
61 const PetscScalar *xarr,*xuarr,*xlarr,*garr,*harr;
62
63 PetscFunctionBegin;
64 ierr = VecGetOwnershipRange(x,&xstart,NULL);CHKERRQ(ierr);
65
66 ierr = VecGetArrayRead(x,&xarr);CHKERRQ(ierr);
67 ierr = VecGetArrayRead(tao->XU,&xuarr);CHKERRQ(ierr);
68 ierr = VecGetArrayRead(tao->XL,&xlarr);CHKERRQ(ierr);
69
70 /* (1) Update ce vector */
71 ierr = VecGetArray(pdipm->ce,&carr);CHKERRQ(ierr);
72
73 if (pdipm->Ng) {
74 /* (1.a) Inserting updated g(x) */
75 ierr = VecGetArrayRead(tao->constraints_equality,&garr);CHKERRQ(ierr);
76 ierr = PetscMemcpy(carr,garr,pdipm->ng*sizeof(PetscScalar));CHKERRQ(ierr);
77 ierr = VecRestoreArrayRead(tao->constraints_equality,&garr);CHKERRQ(ierr);
78 }
79
80 /* (1.b) Update xfixed */
81 if (pdipm->Nxfixed) {
82 offset = pdipm->ng;
83 ierr = ISGetIndices(pdipm->isxfixed,&fxptr);CHKERRQ(ierr); /* global indices in x */
84 for (k=0;k < pdipm->nxfixed;k++){
85 i = fxptr[k]-xstart;
86 carr[offset + k] = xarr[i] - xuarr[i];
87 }
88 }
89 ierr = VecRestoreArray(pdipm->ce,&carr);CHKERRQ(ierr);
90
91 /* (2) Update ci vector */
92 ierr = VecGetArray(pdipm->ci,&carr);CHKERRQ(ierr);
93
94 if (pdipm->Nh) {
95 /* (2.a) Inserting updated h(x) */
96 ierr = VecGetArrayRead(tao->constraints_inequality,&harr);CHKERRQ(ierr);
97 ierr = PetscMemcpy(carr,harr,pdipm->nh*sizeof(PetscScalar));CHKERRQ(ierr);
98 ierr = VecRestoreArrayRead(tao->constraints_inequality,&harr);CHKERRQ(ierr);
99 }
100
101 /* (2.b) Update xub */
102 offset = pdipm->nh;
103 if (pdipm->Nxub) {
104 ierr = ISGetIndices(pdipm->isxub,&ubptr);CHKERRQ(ierr);
105 for (k=0; k<pdipm->nxub; k++){
106 i = ubptr[k]-xstart;
107 carr[offset + k] = xuarr[i] - xarr[i];
108 }
109 }
110
111 if (pdipm->Nxlb) {
112 /* (2.c) Update xlb */
113 offset += pdipm->nxub;
114 ierr = ISGetIndices(pdipm->isxlb,&lbptr);CHKERRQ(ierr); /* global indices in x */
115 for (k=0; k<pdipm->nxlb; k++){
116 i = lbptr[k]-xstart;
117 carr[offset + k] = xarr[i] - xlarr[i];
118 }
119 }
120
121 if (pdipm->Nxbox) {
122 /* (2.d) Update xbox */
123 offset += pdipm->nxlb;
124 offset1 = offset + pdipm->nxbox;
125 ierr = ISGetIndices(pdipm->isxbox,&bxptr);CHKERRQ(ierr); /* global indices in x */
126 for (k=0; k<pdipm->nxbox; k++){
127 i = bxptr[k]-xstart; /* local indices in x */
128 carr[offset+k] = xuarr[i] - xarr[i];
129 carr[offset1+k] = xarr[i] - xlarr[i];
130 }
131 }
132 ierr = VecRestoreArray(pdipm->ci,&carr);CHKERRQ(ierr);
133
134 /* Restoring Vectors */
135 ierr = VecRestoreArrayRead(x,&xarr);CHKERRQ(ierr);
136 ierr = VecRestoreArrayRead(tao->XU,&xuarr);CHKERRQ(ierr);
137 ierr = VecRestoreArrayRead(tao->XL,&xlarr);CHKERRQ(ierr);
138 PetscFunctionReturn(0);
139 }
140
141 /*
142 TaoPDIPMSetUpBounds - Create upper and lower bound vectors of x
143
144 Collective
145
146 Input Parameter:
147 . tao - holds pdipm and XL & XU
148
149 Level: beginner
150
151 .seealso: TaoPDIPMUpdateConstraints
152 */
TaoPDIPMSetUpBounds(Tao tao)153 PetscErrorCode TaoPDIPMSetUpBounds(Tao tao)
154 {
155 PetscErrorCode ierr;
156 TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
157 const PetscScalar *xl,*xu;
158 PetscInt n,*ixlb,*ixub,*ixfixed,*ixfree,*ixbox,i,low,high,idx;
159 MPI_Comm comm;
160 PetscInt sendbuf[5],recvbuf[5];
161
162 PetscFunctionBegin;
163 /* Creates upper and lower bounds vectors on x, if not created already */
164 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr);
165
166 ierr = VecGetLocalSize(tao->XL,&n);CHKERRQ(ierr);
167 ierr = PetscMalloc5(n,&ixlb,n,&ixub,n,&ixfree,n,&ixfixed,n,&ixbox);CHKERRQ(ierr);
168
169 ierr = VecGetOwnershipRange(tao->XL,&low,&high);CHKERRQ(ierr);
170 ierr = VecGetArrayRead(tao->XL,&xl);CHKERRQ(ierr);
171 ierr = VecGetArrayRead(tao->XU,&xu);CHKERRQ(ierr);
172 for (i=0; i<n; i++) {
173 idx = low + i;
174 if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
175 if (PetscRealPart(xl[i]) == PetscRealPart(xu[i])) {
176 ixfixed[pdipm->nxfixed++] = idx;
177 } else ixbox[pdipm->nxbox++] = idx;
178 } else {
179 if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) >= PETSC_INFINITY)) {
180 ixlb[pdipm->nxlb++] = idx;
181 } else if ((PetscRealPart(xl[i]) <= PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
182 ixub[pdipm->nxlb++] = idx;
183 } else ixfree[pdipm->nxfree++] = idx;
184 }
185 }
186 ierr = VecRestoreArrayRead(tao->XL,&xl);CHKERRQ(ierr);
187 ierr = VecRestoreArrayRead(tao->XU,&xu);CHKERRQ(ierr);
188
189 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
190 sendbuf[0] = pdipm->nxlb;
191 sendbuf[1] = pdipm->nxub;
192 sendbuf[2] = pdipm->nxfixed;
193 sendbuf[3] = pdipm->nxbox;
194 sendbuf[4] = pdipm->nxfree;
195
196 ierr = MPI_Allreduce(sendbuf,recvbuf,5,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr);
197 pdipm->Nxlb = recvbuf[0];
198 pdipm->Nxub = recvbuf[1];
199 pdipm->Nxfixed = recvbuf[2];
200 pdipm->Nxbox = recvbuf[3];
201 pdipm->Nxfree = recvbuf[4];
202
203 if (pdipm->Nxlb) {
204 ierr = ISCreateGeneral(comm,pdipm->nxlb,ixlb,PETSC_COPY_VALUES,&pdipm->isxlb);CHKERRQ(ierr);
205 }
206 if (pdipm->Nxub) {
207 ierr = ISCreateGeneral(comm,pdipm->nxub,ixub,PETSC_COPY_VALUES,&pdipm->isxub);CHKERRQ(ierr);
208 }
209 if (pdipm->Nxfixed) {
210 ierr = ISCreateGeneral(comm,pdipm->nxfixed,ixfixed,PETSC_COPY_VALUES,&pdipm->isxfixed);CHKERRQ(ierr);
211 }
212 if (pdipm->Nxbox) {
213 ierr = ISCreateGeneral(comm,pdipm->nxbox,ixbox,PETSC_COPY_VALUES,&pdipm->isxbox);CHKERRQ(ierr);
214 }
215 if (pdipm->Nxfree) {
216 ierr = ISCreateGeneral(comm,pdipm->nxfree,ixfree,PETSC_COPY_VALUES,&pdipm->isxfree);CHKERRQ(ierr);
217 }
218 ierr = PetscFree5(ixlb,ixub,ixfixed,ixbox,ixfree);CHKERRQ(ierr);
219 PetscFunctionReturn(0);
220 }
221
222 /*
223 TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
224 X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
225 four subvectors need to be initialized and its values copied over to X. Instead
226 of copying, we use VecPlace/ResetArray functions to share the memory locations for
227 X and the subvectors
228
229 Collective
230
231 Input Parameter:
232 . tao - Tao context
233
234 Level: beginner
235 */
TaoPDIPMInitializeSolution(Tao tao)236 PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
237 {
238 PetscErrorCode ierr;
239 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
240 PetscScalar *Xarr,*z,*lambdai;
241 PetscInt i;
242 const PetscScalar *xarr,*h;
243
244 PetscFunctionBegin;
245 ierr = VecGetArray(pdipm->X,&Xarr);CHKERRQ(ierr);
246
247 /* Set Initialize X.x = tao->solution */
248 ierr = VecGetArrayRead(tao->solution,&xarr);CHKERRQ(ierr);
249 ierr = PetscMemcpy(Xarr,xarr,pdipm->nx*sizeof(PetscScalar));CHKERRQ(ierr);
250 ierr = VecRestoreArrayRead(tao->solution,&xarr);CHKERRQ(ierr);
251
252 /* Initialize X.lambdae = 0.0 */
253 if (pdipm->lambdae) {
254 ierr = VecSet(pdipm->lambdae,0.0);CHKERRQ(ierr);
255 }
256 /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
257 if (pdipm->lambdai) {
258 ierr = VecSet(pdipm->lambdai,pdipm->push_init_lambdai);CHKERRQ(ierr);
259 }
260 if (pdipm->z) {
261 ierr = VecSet(pdipm->z,pdipm->push_init_slack);CHKERRQ(ierr);
262 }
263
264 /* Additional modification for X.lambdai and X.z */
265 if (pdipm->lambdai) {
266 ierr = VecGetArray(pdipm->lambdai,&lambdai);CHKERRQ(ierr);
267 }
268 if (pdipm->z) {
269 ierr = VecGetArray(pdipm->z,&z);CHKERRQ(ierr);
270 }
271 if (pdipm->Nh) {
272 ierr = VecGetArrayRead(tao->constraints_inequality,&h);CHKERRQ(ierr);
273 for (i=0; i < pdipm->nh; i++) {
274 if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
275 if (pdipm->mu/z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu/z[i];
276 }
277 ierr = VecRestoreArrayRead(tao->constraints_inequality,&h);CHKERRQ(ierr);
278 }
279 if (pdipm->lambdai) {
280 ierr = VecRestoreArray(pdipm->lambdai,&lambdai);CHKERRQ(ierr);
281 }
282 if (pdipm->z) {
283 ierr = VecRestoreArray(pdipm->z,&z);CHKERRQ(ierr);
284 }
285
286 ierr = VecRestoreArray(pdipm->X,&Xarr);CHKERRQ(ierr);
287 PetscFunctionReturn(0);
288 }
289
290 /*
291 TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X
292
293 Input Parameter:
294 snes - SNES context
295 X - KKT Vector
296 *ctx - pdipm context
297
298 Output Parameter:
299 J - Hessian matrix
300 Jpre - Preconditioner
301 */
TaoSNESJacobian_PDIPM(SNES snes,Vec X,Mat J,Mat Jpre,void * ctx)302 PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes,Vec X, Mat J, Mat Jpre, void *ctx)
303 {
304 PetscErrorCode ierr;
305 Tao tao=(Tao)ctx;
306 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
307 PetscInt i,row,cols[2],Jrstart,rjstart,nc,j;
308 const PetscInt *aj,*ranges,*Jranges,*rranges,*cranges;
309 const PetscScalar *Xarr,*aa;
310 PetscScalar vals[2];
311 PetscInt proc,nx_all,*nce_all=pdipm->nce_all;
312 MPI_Comm comm;
313 PetscMPIInt rank,size;
314 Mat jac_equality_trans=pdipm->jac_equality_trans,jac_inequality_trans=pdipm->jac_inequality_trans;
315
316 PetscFunctionBegin;
317 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr);
318 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
319 ierr = MPI_Comm_rank(comm,&size);CHKERRQ(ierr);
320
321 ierr = MatGetOwnershipRanges(Jpre,&Jranges);CHKERRQ(ierr);
322 ierr = MatGetOwnershipRange(Jpre,&Jrstart,NULL);CHKERRQ(ierr);
323 ierr = MatGetOwnershipRangesColumn(tao->hessian,&rranges);CHKERRQ(ierr);
324 ierr = MatGetOwnershipRangesColumn(tao->hessian,&cranges);CHKERRQ(ierr);
325
326 ierr = VecGetArrayRead(X,&Xarr);CHKERRQ(ierr);
327
328 /* (2) insert Z and Ci to Jpre -- overwrite existing values */
329 for (i=0; i < pdipm->nci; i++) {
330 row = Jrstart + pdipm->off_z + i;
331 cols[0] = Jrstart + pdipm->off_lambdai + i;
332 cols[1] = row;
333 vals[0] = Xarr[pdipm->off_z + i];
334 vals[1] = Xarr[pdipm->off_lambdai + i];
335 ierr = MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
336 }
337
338 /* (3) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
339 if (pdipm->Ng) {
340 ierr = MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);CHKERRQ(ierr);
341 for (i=0; i<pdipm->ng; i++){
342 row = Jrstart + pdipm->off_lambdae + i;
343
344 ierr = MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
345 proc = 0;
346 for (j=0; j < nc; j++) {
347 while (aj[j] >= cranges[proc+1]) proc++;
348 cols[0] = aj[j] - cranges[proc] + Jranges[proc];
349 ierr = MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);CHKERRQ(ierr);
350 }
351 ierr = MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
352 }
353 }
354
355 if (pdipm->Nh) {
356 /* (4) insert 3nd row block of Jpre: [ grad h, 0, 0, 0] */
357 ierr = MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);CHKERRQ(ierr);
358 for (i=0; i < pdipm->nh; i++){
359 row = Jrstart + pdipm->off_lambdai + i;
360
361 ierr = MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
362 proc = 0;
363 for (j=0; j < nc; j++) {
364 while (aj[j] >= cranges[proc+1]) proc++;
365 cols[0] = aj[j] - cranges[proc] + Jranges[proc];
366 ierr = MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);CHKERRQ(ierr);
367 }
368 ierr = MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
369 }
370 }
371
372 /* (5) insert Wxx, grad g' and -grad h' to Jpre */
373 if (pdipm->Ng) {
374 ierr = MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);CHKERRQ(ierr);
375 }
376 if (pdipm->Nh) {
377 ierr = MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);CHKERRQ(ierr);
378 }
379
380 ierr = VecPlaceArray(pdipm->x,Xarr);CHKERRQ(ierr);
381 ierr = TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
382 ierr = VecResetArray(pdipm->x);CHKERRQ(ierr);
383
384 ierr = MatGetOwnershipRange(tao->hessian,&rjstart,NULL);CHKERRQ(ierr);
385 for (i=0; i<pdipm->nx; i++){
386 row = Jrstart + i;
387
388 /* insert Wxx */
389 ierr = MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
390 proc = 0;
391 for (j=0; j < nc; j++) {
392 while (aj[j] >= cranges[proc+1]) proc++;
393 cols[0] = aj[j] - cranges[proc] + Jranges[proc];
394 ierr = MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);CHKERRQ(ierr);
395 }
396 ierr = MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
397
398 if (pdipm->ng) {
399 /* insert grad g' */
400 ierr = MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
401 ierr = MatGetOwnershipRanges(tao->jacobian_equality,&ranges);CHKERRQ(ierr);
402 proc = 0;
403 for (j=0; j < nc; j++) {
404 /* find row ownership of */
405 while (aj[j] >= ranges[proc+1]) proc++;
406 nx_all = rranges[proc+1] - rranges[proc];
407 cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
408 ierr = MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);CHKERRQ(ierr);
409 }
410 ierr = MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
411 }
412
413 if (pdipm->nh) {
414 /* insert -grad h' */
415 ierr = MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
416 ierr = MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);CHKERRQ(ierr);
417 proc = 0;
418 for (j=0; j < nc; j++) {
419 /* find row ownership of */
420 while (aj[j] >= ranges[proc+1]) proc++;
421 nx_all = rranges[proc+1] - rranges[proc];
422 cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
423 ierr = MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);CHKERRQ(ierr);
424 }
425 ierr = MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);CHKERRQ(ierr);
426 }
427 }
428 ierr = VecRestoreArrayRead(X,&Xarr);CHKERRQ(ierr);
429
430 /* (6) assemble Jpre and J */
431 ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
432 ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
433
434 if (J != Jpre) {
435 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
436 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
437 }
438 PetscFunctionReturn(0);
439 }
440
441 /*
442 TaoSnesFunction_PDIPM - Evaluate KKT function at X
443
444 Input Parameter:
445 snes - SNES context
446 X - KKT Vector
447 *ctx - pdipm
448
449 Output Parameter:
450 F - Updated Lagrangian vector
451 */
TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void * ctx)452 PetscErrorCode TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void *ctx)
453 {
454 PetscErrorCode ierr;
455 Tao tao=(Tao)ctx;
456 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
457 PetscScalar *Farr;
458 Vec x,L1;
459 PetscInt i;
460 PetscReal res[2],cnorm[2];
461 const PetscScalar *Xarr,*carr,*zarr,*larr;
462
463 PetscFunctionBegin;
464 ierr = VecSet(F,0.0);CHKERRQ(ierr);
465
466 ierr = VecGetArrayRead(X,&Xarr);CHKERRQ(ierr);
467 ierr = VecGetArray(F,&Farr);CHKERRQ(ierr);
468
469 /* (0) Evaluate f, fx, Gx, Hx at X.x Note: pdipm->x is not changed below */
470 x = pdipm->x;
471 ierr = VecPlaceArray(x,Xarr);CHKERRQ(ierr);
472 ierr = TaoPDIPMEvaluateFunctionsAndJacobians(tao,x);CHKERRQ(ierr);
473
474 /* Update ce, ci, and Jci at X.x */
475 ierr = TaoPDIPMUpdateConstraints(tao,x);CHKERRQ(ierr);
476 ierr = VecResetArray(x);CHKERRQ(ierr);
477
478 /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
479 L1 = pdipm->x;
480 ierr = VecPlaceArray(L1,Farr);CHKERRQ(ierr);
481 if (pdipm->Nci) {
482 if (pdipm->Nh) {
483 /* L1 += gradH'*DI. Note: tao->DI is not changed below */
484 ierr = VecPlaceArray(tao->DI,Xarr+pdipm->off_lambdai);CHKERRQ(ierr);
485 ierr = MatMultTransposeAdd(tao->jacobian_inequality,tao->DI,L1,L1);CHKERRQ(ierr);
486 ierr = VecResetArray(tao->DI);CHKERRQ(ierr);
487 }
488
489 /* L1 += Jci_xb'*lambdai_xb */
490 ierr = VecPlaceArray(pdipm->lambdai_xb,Xarr+pdipm->off_lambdai+pdipm->nh);CHKERRQ(ierr);
491 ierr = MatMultTransposeAdd(pdipm->Jci_xb,pdipm->lambdai_xb,L1,L1);CHKERRQ(ierr);
492 ierr = VecResetArray(pdipm->lambdai_xb);CHKERRQ(ierr);
493
494 /* (1.4) L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
495 ierr = VecScale(L1,-1.0);CHKERRQ(ierr);
496 }
497
498 /* L1 += fx */
499 ierr = VecAXPY(L1,1.0,tao->gradient);CHKERRQ(ierr);
500
501 if (pdipm->Nce) {
502 if (pdipm->Ng) {
503 /* L1 += gradG'*DE. Note: tao->DE is not changed below */
504 ierr = VecPlaceArray(tao->DE,Xarr+pdipm->off_lambdae);CHKERRQ(ierr);
505 ierr = MatMultTransposeAdd(tao->jacobian_equality,tao->DE,L1,L1);CHKERRQ(ierr);
506 ierr = VecResetArray(tao->DE);CHKERRQ(ierr);
507 }
508 if (pdipm->Nxfixed) {
509 /* L1 += Jce_xfixed'*lambdae_xfixed */
510 ierr = VecPlaceArray(pdipm->lambdae_xfixed,Xarr+pdipm->off_lambdae+pdipm->ng);CHKERRQ(ierr);
511 ierr = MatMultTransposeAdd(pdipm->Jce_xfixed,pdipm->lambdae_xfixed,L1,L1);CHKERRQ(ierr);
512 ierr = VecResetArray(pdipm->lambdae_xfixed);CHKERRQ(ierr);
513 }
514 }
515 ierr = VecNorm(L1,NORM_2,&res[0]);CHKERRQ(ierr);
516 ierr = VecResetArray(L1);CHKERRQ(ierr);
517
518 /* (2) L2 = ce(x) */
519 if (pdipm->Nce) {
520 ierr = VecGetArrayRead(pdipm->ce,&carr);CHKERRQ(ierr);
521 for (i=0; i<pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
522 ierr = VecRestoreArrayRead(pdipm->ce,&carr);CHKERRQ(ierr);
523 }
524 ierr = VecNorm(pdipm->ce,NORM_2,&cnorm[0]);CHKERRQ(ierr);
525
526 if (pdipm->Nci) {
527 /* (3) L3 = ci(x) - z;
528 (4) L4 = Z * Lambdai * e - mu * e
529 */
530 ierr = VecGetArrayRead(pdipm->ci,&carr);CHKERRQ(ierr);
531 larr = Xarr+pdipm->off_lambdai;
532 zarr = Xarr+pdipm->off_z;
533 for (i=0; i<pdipm->nci; i++) {
534 Farr[pdipm->off_lambdai + i] = carr[i] - zarr[i];
535 Farr[pdipm->off_z + i] = zarr[i]*larr[i] - pdipm->mu;
536 }
537 ierr = VecRestoreArrayRead(pdipm->ci,&carr);CHKERRQ(ierr);
538 }
539
540 ierr = VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);CHKERRQ(ierr);
541 ierr = VecNorm(pdipm->ci,NORM_2,&cnorm[1]);CHKERRQ(ierr);
542 ierr = VecResetArray(pdipm->ci);CHKERRQ(ierr);
543
544 /* note: pdipm->z is not changed below */
545 if (pdipm->z) {
546 ierr = VecPlaceArray(pdipm->z,Farr+pdipm->off_z);CHKERRQ(ierr);
547 ierr = VecNorm(pdipm->z,NORM_2,&res[1]);CHKERRQ(ierr);
548 ierr = VecResetArray(pdipm->z);CHKERRQ(ierr);
549 } else res[1] = 0.0;
550
551 tao->residual = PetscSqrtReal(res[0]*res[0] + res[1]*res[1]);
552 tao->cnorm = PetscSqrtReal(cnorm[0]*cnorm[0] + cnorm[1]*cnorm[1]);
553
554 ierr = VecRestoreArrayRead(X,&Xarr);CHKERRQ(ierr);
555 ierr = VecRestoreArray(F,&Farr);CHKERRQ(ierr);
556 PetscFunctionReturn(0);
557 }
558
559 /*
560 PDIPMLineSearch - Custom line search used with PDIPM.
561
562 Collective on TAO
563
564 Notes:
565 PDIPMLineSearch employs a simple backtracking line-search to keep
566 the slack variables (z) and inequality constraints lagrange multipliers
567 (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
568 alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
569 slack variables are updated as X = X + alpha_d*dx. The constraint multipliers
570 are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
571 is also updated as mu = mu + z'lambdai/Nci
572 */
PDIPMLineSearch(SNESLineSearch linesearch,void * ctx)573 PetscErrorCode PDIPMLineSearch(SNESLineSearch linesearch,void *ctx)
574 {
575 PetscErrorCode ierr;
576 Tao tao=(Tao)ctx;
577 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
578 SNES snes;
579 Vec X,F,Y,W,G;
580 PetscInt i,iter;
581 PetscReal alpha_p=1.0,alpha_d=1.0,alpha[4];
582 PetscScalar *Xarr,*z,*lambdai,dot;
583 const PetscScalar *dXarr,*dz,*dlambdai;
584 PetscScalar *taosolarr;
585
586 PetscFunctionBegin;
587 ierr = SNESLineSearchGetSNES(linesearch,&snes);CHKERRQ(ierr);
588 ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr);
589
590 ierr = SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);CHKERRQ(ierr);
591 ierr = SNESLineSearchGetVecs(linesearch,&X,&F,&Y,&W,&G);CHKERRQ(ierr);
592
593 ierr = VecGetArray(X,&Xarr);CHKERRQ(ierr);
594 ierr = VecGetArrayRead(Y,&dXarr);CHKERRQ(ierr);
595 z = Xarr + pdipm->off_z;
596 dz = dXarr + pdipm->off_z;
597 for (i=0; i < pdipm->nci; i++) {
598 if (z[i] - dz[i] < 0.0) {
599 alpha_p = PetscMin(alpha_p,0.9999*z[i]/dz[i]);
600 }
601 }
602
603 lambdai = Xarr + pdipm->off_lambdai;
604 dlambdai = dXarr + pdipm->off_lambdai;
605
606 for (i=0; i<pdipm->nci; i++) {
607 if (lambdai[i] - dlambdai[i] < 0.0) {
608 alpha_d = PetscMin(0.9999*lambdai[i]/dlambdai[i],alpha_d);
609 }
610 }
611
612 alpha[0] = alpha_p;
613 alpha[1] = alpha_d;
614 ierr = VecRestoreArrayRead(Y,&dXarr);CHKERRQ(ierr);
615 ierr = VecRestoreArray(X,&Xarr);CHKERRQ(ierr);
616
617 /* alpha = min(alpha) over all processes */
618 ierr = MPI_Allreduce(alpha,alpha+2,2,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)tao));CHKERRQ(ierr);
619
620 alpha_p = alpha[2];
621 alpha_d = alpha[3];
622
623 ierr = VecGetArray(X,&Xarr);CHKERRQ(ierr);
624 ierr = VecGetArrayRead(Y,&dXarr);CHKERRQ(ierr);
625 for (i=0; i<pdipm->nx; i++) {
626 Xarr[i] = Xarr[i] - alpha_p * dXarr[i];
627 }
628
629 for (i=0; i<pdipm->nce; i++) {
630 Xarr[i+pdipm->off_lambdae] = Xarr[i+pdipm->off_lambdae] - alpha_d * dXarr[i+pdipm->off_lambdae];
631 }
632
633 for (i=0; i<pdipm->nci; i++) {
634 Xarr[i+pdipm->off_lambdai] = Xarr[i+pdipm->off_lambdai] - alpha_d * dXarr[i+pdipm->off_lambdai];
635 }
636
637 for (i=0; i<pdipm->nci; i++) {
638 Xarr[i+pdipm->off_z] = Xarr[i+pdipm->off_z] - alpha_p * dXarr[i+pdipm->off_z];
639 }
640
641 ierr = VecGetArray(tao->solution,&taosolarr);CHKERRQ(ierr);
642 ierr = PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));CHKERRQ(ierr);
643 ierr = VecRestoreArray(tao->solution,&taosolarr);CHKERRQ(ierr);
644
645
646 ierr = VecRestoreArray(X,&Xarr);CHKERRQ(ierr);
647 ierr = VecRestoreArrayRead(Y,&dXarr);CHKERRQ(ierr);
648
649 /* Evaluate F at X */
650 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
651 ierr = SNESLineSearchComputeNorms(linesearch);CHKERRQ(ierr); /* must call this func, do not know why */
652
653 /* update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
654 if (pdipm->z) {
655 ierr = VecDot(pdipm->z,pdipm->lambdai,&dot);CHKERRQ(ierr);
656 } else dot = 0.0;
657
658 /* if (PetscAbsReal(pdipm->gradL) < 0.9*pdipm->mu) */
659 pdipm->mu = pdipm->mu_update_factor * dot/pdipm->Nci;
660
661 /* Update F; get tao->residual and tao->cnorm */
662 ierr = TaoSNESFunction_PDIPM(snes,X,F,(void*)tao);CHKERRQ(ierr);
663
664 tao->niter++;
665 ierr = TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);CHKERRQ(ierr);
666 ierr = TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);CHKERRQ(ierr);
667
668 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
669 if (tao->reason) {
670 ierr = SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);CHKERRQ(ierr);
671 }
672 PetscFunctionReturn(0);
673 }
674
675 /*
676 TaoSolve_PDIPM
677
678 Input Parameter:
679 tao - TAO context
680
681 Output Parameter:
682 tao - TAO context
683 */
TaoSolve_PDIPM(Tao tao)684 PetscErrorCode TaoSolve_PDIPM(Tao tao)
685 {
686 PetscErrorCode ierr;
687 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
688 SNESLineSearch linesearch; /* SNESLineSearch context */
689 Vec dummy;
690
691 PetscFunctionBegin;
692 if (!tao->constraints_equality && !tao->constraints_inequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_NULL,"Equality and inequality constraints are not set. Either set them or switch to a different algorithm");
693
694 /* Initialize all variables */
695 ierr = TaoPDIPMInitializeSolution(tao);CHKERRQ(ierr);
696
697 /* Set linesearch */
698 ierr = SNESGetLineSearch(pdipm->snes,&linesearch);CHKERRQ(ierr);
699 ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);CHKERRQ(ierr);
700 ierr = SNESLineSearchShellSetUserFunc(linesearch,PDIPMLineSearch,tao);CHKERRQ(ierr);
701 ierr = SNESLineSearchSetFromOptions(linesearch);CHKERRQ(ierr);
702
703 tao->reason = TAO_CONTINUE_ITERATING;
704
705 /* -tao_monitor for iteration 0 and check convergence */
706 ierr = VecDuplicate(pdipm->X,&dummy);CHKERRQ(ierr);
707 ierr = TaoSNESFunction_PDIPM(pdipm->snes,pdipm->X,dummy,(void*)tao);CHKERRQ(ierr);
708
709 ierr = TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);CHKERRQ(ierr);
710 ierr = TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);CHKERRQ(ierr);
711 ierr = VecDestroy(&dummy);CHKERRQ(ierr);
712 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
713 if (tao->reason) {
714 ierr = SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);CHKERRQ(ierr);
715 }
716
717 while (tao->reason == TAO_CONTINUE_ITERATING) {
718 SNESConvergedReason reason;
719 ierr = SNESSolve(pdipm->snes,NULL,pdipm->X);CHKERRQ(ierr);
720
721 /* Check SNES convergence */
722 ierr = SNESGetConvergedReason(pdipm->snes,&reason);CHKERRQ(ierr);
723 if (reason < 0) {
724 ierr = PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes),"SNES solve did not converged due to reason %D\n",reason);CHKERRQ(ierr);
725 }
726
727 /* Check TAO convergence */
728 if (PetscIsInfOrNanReal(pdipm->obj)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"User-provided compute function generated Inf or NaN");
729 }
730 PetscFunctionReturn(0);
731 }
732
733 /*
734 TaoSetup_PDIPM - Sets up tao and pdipm
735
736 Input Parameter:
737 tao - TAO object
738
739 Output: pdipm - initialized object
740 */
TaoSetup_PDIPM(Tao tao)741 PetscErrorCode TaoSetup_PDIPM(Tao tao)
742 {
743 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
744 PetscErrorCode ierr;
745 MPI_Comm comm;
746 PetscMPIInt rank,size;
747 PetscInt row,col,Jcrstart,Jcrend,k,tmp,nc,proc,*nh_all,*ng_all;
748 PetscInt offset,*xa,*xb,i,j,rstart,rend;
749 PetscScalar one=1.0,neg_one=-1.0,*Xarr;
750 const PetscInt *cols,*rranges,*cranges,*aj,*ranges;
751 const PetscScalar *aa;
752 Mat J,jac_equality_trans,jac_inequality_trans;
753 Mat Jce_xfixed_trans,Jci_xb_trans;
754 PetscInt *dnz,*onz,rjstart,nx_all,*nce_all,*Jranges,cols1[2];
755
756 PetscFunctionBegin;
757 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr);
758 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
759 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
760
761 /* (1) Setup Bounds and create Tao vectors */
762 ierr = TaoPDIPMSetUpBounds(tao);CHKERRQ(ierr);
763
764 if (!tao->gradient) {
765 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);
766 ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr);
767 }
768
769 /* (2) Get sizes */
770 /* Size of vector x - This is set by TaoSetInitialVector */
771 ierr = VecGetSize(tao->solution,&pdipm->Nx);CHKERRQ(ierr);
772 ierr = VecGetLocalSize(tao->solution,&pdipm->nx);CHKERRQ(ierr);
773
774 /* Size of equality constraints and vectors */
775 if (tao->constraints_equality) {
776 ierr = VecGetSize(tao->constraints_equality,&pdipm->Ng);CHKERRQ(ierr);
777 ierr = VecGetLocalSize(tao->constraints_equality,&pdipm->ng);CHKERRQ(ierr);
778 } else {
779 pdipm->ng = pdipm->Ng = 0;
780 }
781
782 pdipm->nce = pdipm->ng + pdipm->nxfixed;
783 pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;
784
785 /* Size of inequality constraints and vectors */
786 if (tao->constraints_inequality) {
787 ierr = VecGetSize(tao->constraints_inequality,&pdipm->Nh);CHKERRQ(ierr);
788 ierr = VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);CHKERRQ(ierr);
789 } else {
790 pdipm->nh = pdipm->Nh = 0;
791 }
792
793 pdipm->nci = pdipm->nh + pdipm->nxlb + pdipm->nxub + 2*pdipm->nxbox;
794 pdipm->Nci = pdipm->Nh + pdipm->Nxlb + pdipm->Nxub + 2*pdipm->Nxbox;
795
796 /* Full size of the KKT system to be solved */
797 pdipm->n = pdipm->nx + pdipm->nce + 2*pdipm->nci;
798 pdipm->N = pdipm->Nx + pdipm->Nce + 2*pdipm->Nci;
799
800 /* list below to TaoView_PDIPM()? */
801 /* ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] nce %d = ng %d + nxfixed %d\n",rank,pdipm->nce,pdipm->ng,pdipm->nxfixed); */
802 /* ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] nci %d = nh %d + nxlb %d + nxub %d + 2*nxbox %d\n",rank,pdipm->nci,pdipm->nh,pdipm->nxlb,pdipm->nxub,pdipm->nxbox); */
803 /* ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] n %d = nx %d + nce %d + 2*nci %d\n",rank,pdipm->n,pdipm->nx,pdipm->nce,pdipm->nci); */
804
805 /* (3) Offsets for subvectors */
806 pdipm->off_lambdae = pdipm->nx;
807 pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
808 pdipm->off_z = pdipm->off_lambdai + pdipm->nci;
809
810 /* (4) Create vectors and subvectors */
811 /* Ce and Ci vectors */
812 ierr = VecCreate(comm,&pdipm->ce);CHKERRQ(ierr);
813 ierr = VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);CHKERRQ(ierr);
814 ierr = VecSetFromOptions(pdipm->ce);CHKERRQ(ierr);
815
816 ierr = VecCreate(comm,&pdipm->ci);CHKERRQ(ierr);
817 ierr = VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);CHKERRQ(ierr);
818 ierr = VecSetFromOptions(pdipm->ci);CHKERRQ(ierr);
819
820 /* X=[x; lambdae; lambdai; z] for the big KKT system */
821 ierr = VecCreate(comm,&pdipm->X);CHKERRQ(ierr);
822 ierr = VecSetSizes(pdipm->X,pdipm->n,pdipm->N);CHKERRQ(ierr);
823 ierr = VecSetFromOptions(pdipm->X);CHKERRQ(ierr);
824
825 /* Subvectors; they share local arrays with X */
826 ierr = VecGetArray(pdipm->X,&Xarr);CHKERRQ(ierr);
827 /* x shares local array with X.x */
828 if (pdipm->Nx) {
829 ierr = VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);CHKERRQ(ierr);
830 }
831
832 /* lambdae shares local array with X.lambdae */
833 if (pdipm->Nce) {
834 ierr = VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);CHKERRQ(ierr);
835 }
836
837 /* tao->DE shares local array with X.lambdae_g */
838 if (pdipm->Ng) {
839 ierr = VecCreateMPIWithArray(comm,1,pdipm->ng,pdipm->Ng,Xarr+pdipm->off_lambdae,&tao->DE);CHKERRQ(ierr);
840
841 ierr = VecCreate(comm,&pdipm->lambdae_xfixed);CHKERRQ(ierr);
842 ierr = VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);CHKERRQ(ierr);
843 ierr = VecSetFromOptions(pdipm->lambdae_xfixed);CHKERRQ(ierr);
844 }
845
846 if (pdipm->Nci) {
847 /* lambdai shares local array with X.lambdai */
848 ierr = VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_lambdai,&pdipm->lambdai);CHKERRQ(ierr);
849
850 /* z for slack variables; it shares local array with X.z */
851 ierr = VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_z,&pdipm->z);CHKERRQ(ierr);
852 }
853
854 /* tao->DI which shares local array with X.lambdai_h */
855 if (pdipm->Nh) {
856 ierr = VecCreateMPIWithArray(comm,1,pdipm->nh,pdipm->Nh,Xarr+pdipm->off_lambdai,&tao->DI);CHKERRQ(ierr);
857 }
858
859 ierr = VecCreate(comm,&pdipm->lambdai_xb);CHKERRQ(ierr);
860 ierr = VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);CHKERRQ(ierr);
861 ierr = VecSetFromOptions(pdipm->lambdai_xb);CHKERRQ(ierr);
862
863 ierr = VecRestoreArray(pdipm->X,&Xarr);CHKERRQ(ierr);
864
865 /* (5) Create Jacobians Jce_xfixed and Jci */
866 /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
867 if (pdipm->Nxfixed) {
868 /* Create Jce_xfixed */
869 ierr = MatCreate(comm,&pdipm->Jce_xfixed);CHKERRQ(ierr);
870 ierr = MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);CHKERRQ(ierr);
871 ierr = MatSetFromOptions(pdipm->Jce_xfixed);CHKERRQ(ierr);
872 ierr = MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);CHKERRQ(ierr);
873 ierr = MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);CHKERRQ(ierr);
874
875 ierr = MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);CHKERRQ(ierr);
876 ierr = ISGetIndices(pdipm->isxfixed,&cols);CHKERRQ(ierr);
877 k = 0;
878 for (row = Jcrstart; row < Jcrend; row++) {
879 ierr = MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);CHKERRQ(ierr);
880 k++;
881 }
882 ierr = ISRestoreIndices(pdipm->isxfixed, &cols);CHKERRQ(ierr);
883 ierr = MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
884 ierr = MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
885 }
886
887 /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
888 ierr = MatCreate(comm,&pdipm->Jci_xb);CHKERRQ(ierr);
889 ierr = MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);CHKERRQ(ierr);
890 ierr = MatSetFromOptions(pdipm->Jci_xb);CHKERRQ(ierr);
891 ierr = MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);CHKERRQ(ierr);
892 ierr = MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);CHKERRQ(ierr);
893
894 ierr = MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);CHKERRQ(ierr);
895 offset = Jcrstart;
896 if (pdipm->Nxub) {
897 /* Add xub to Jci_xb */
898 ierr = ISGetIndices(pdipm->isxub,&cols);CHKERRQ(ierr);
899 k = 0;
900 for (row = offset; row < offset + pdipm->nxub; row++) {
901 ierr = MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
902 k++;
903 }
904 ierr = ISRestoreIndices(pdipm->isxub, &cols);CHKERRQ(ierr);
905 }
906
907 if (pdipm->Nxlb) {
908 /* Add xlb to Jci_xb */
909 ierr = ISGetIndices(pdipm->isxlb,&cols);CHKERRQ(ierr);
910 k = 0;
911 offset += pdipm->nxub;
912 for (row = offset; row < offset + pdipm->nxlb; row++) {
913 ierr = MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);CHKERRQ(ierr);
914 k++;
915 }
916 ierr = ISRestoreIndices(pdipm->isxlb, &cols);CHKERRQ(ierr);
917 }
918
919 /* Add xbox to Jci_xb */
920 if (pdipm->Nxbox) {
921 ierr = ISGetIndices(pdipm->isxbox,&cols);CHKERRQ(ierr);
922 k = 0;
923 offset += pdipm->nxlb;
924 for (row = offset; row < offset + pdipm->nxbox; row++) {
925 ierr = MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
926 tmp = row + pdipm->nxbox;
927 ierr = MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);CHKERRQ(ierr);
928 k++;
929 }
930 ierr = ISRestoreIndices(pdipm->isxbox, &cols);CHKERRQ(ierr);
931 }
932
933 ierr = MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
934 ierr = MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
935 /* ierr = MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
936
937 /* (6) Set up ISs for PC Fieldsplit */
938 if (pdipm->solve_reduced_kkt) {
939 ierr = PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);CHKERRQ(ierr);
940 for (i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
941 for (i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;
942
943 ierr = ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);CHKERRQ(ierr);
944 ierr = ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);CHKERRQ(ierr);
945 }
946
947 /* (7) Gather offsets from all processes */
948 ierr = PetscMalloc1(size,&pdipm->nce_all);CHKERRQ(ierr);
949
950 /* Get rstart of KKT matrix */
951 ierr = MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr);
952 rstart -= pdipm->n;
953
954 ierr = MPI_Allgather(&pdipm->nce,1,MPIU_INT,pdipm->nce_all,1,MPIU_INT,comm);CHKERRQ(ierr);
955
956 ierr = PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);CHKERRQ(ierr);
957 ierr = MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);CHKERRQ(ierr);
958 ierr = MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);CHKERRQ(ierr);
959 ierr = MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);CHKERRQ(ierr);
960
961 ierr = MatGetOwnershipRanges(tao->hessian,&rranges);CHKERRQ(ierr);
962 ierr = MatGetOwnershipRangesColumn(tao->hessian,&cranges);CHKERRQ(ierr);
963
964 if (pdipm->Ng) {
965 ierr = TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);CHKERRQ(ierr);
966 ierr = MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);CHKERRQ(ierr);
967 }
968 if (pdipm->Nh) {
969 ierr = TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);CHKERRQ(ierr);
970 ierr = MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);CHKERRQ(ierr);
971 }
972
973 /* Count dnz,onz for preallocation of KKT matrix */
974 jac_equality_trans = pdipm->jac_equality_trans;
975 jac_inequality_trans = pdipm->jac_inequality_trans;
976 nce_all = pdipm->nce_all;
977
978 if (pdipm->Nxfixed) {
979 ierr = MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);CHKERRQ(ierr);
980 }
981 ierr = MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);CHKERRQ(ierr);
982
983 ierr = MatPreallocateInitialize(comm,pdipm->n,pdipm->n,dnz,onz);CHKERRQ(ierr);
984
985 /* 1st row block of KKT matrix: [Wxx; gradCe'; -gradCi'; 0] */
986 ierr = TaoPDIPMEvaluateFunctionsAndJacobians(tao,pdipm->x);CHKERRQ(ierr);
987 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
988
989 /* Insert tao->hessian */
990 ierr = MatGetOwnershipRange(tao->hessian,&rjstart,NULL);CHKERRQ(ierr);
991 for (i=0; i<pdipm->nx; i++){
992 row = rstart + i;
993
994 ierr = MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
995 proc = 0;
996 for (j=0; j < nc; j++) {
997 while (aj[j] >= cranges[proc+1]) proc++;
998 col = aj[j] - cranges[proc] + Jranges[proc];
999 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1000 }
1001 ierr = MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1002
1003 if (pdipm->ng) {
1004 /* Insert grad g' */
1005 ierr = MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1006 ierr = MatGetOwnershipRanges(tao->jacobian_equality,&ranges);CHKERRQ(ierr);
1007 proc = 0;
1008 for (j=0; j < nc; j++) {
1009 /* find row ownership of */
1010 while (aj[j] >= ranges[proc+1]) proc++;
1011 nx_all = rranges[proc+1] - rranges[proc];
1012 col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
1013 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1014 }
1015 ierr = MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1016 }
1017
1018 /* Insert Jce_xfixed^T' */
1019 if (pdipm->nxfixed) {
1020 ierr = MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1021 ierr = MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);CHKERRQ(ierr);
1022 proc = 0;
1023 for (j=0; j < nc; j++) {
1024 /* find row ownership of */
1025 while (aj[j] >= ranges[proc+1]) proc++;
1026 nx_all = rranges[proc+1] - rranges[proc];
1027 col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1028 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1029 }
1030 ierr = MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1031 }
1032
1033 if (pdipm->nh) {
1034 /* Insert -grad h' */
1035 ierr = MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1036 ierr = MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);CHKERRQ(ierr);
1037 proc = 0;
1038 for (j=0; j < nc; j++) {
1039 /* find row ownership of */
1040 while (aj[j] >= ranges[proc+1]) proc++;
1041 nx_all = rranges[proc+1] - rranges[proc];
1042 col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1043 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1044 }
1045 ierr = MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1046 }
1047
1048 /* Insert Jci_xb^T' */
1049 ierr = MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1050 ierr = MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);CHKERRQ(ierr);
1051 proc = 0;
1052 for (j=0; j < nc; j++) {
1053 /* find row ownership of */
1054 while (aj[j] >= ranges[proc+1]) proc++;
1055 nx_all = rranges[proc+1] - rranges[proc];
1056 col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1057 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1058 }
1059 ierr = MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1060 }
1061
1062 /* 2nd Row block of KKT matrix: [grad Ce, 0, 0, 0] */
1063 if (pdipm->Ng) {
1064 ierr = MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);CHKERRQ(ierr);
1065 for (i=0; i < pdipm->ng; i++){
1066 row = rstart + pdipm->off_lambdae + i;
1067
1068 ierr = MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1069 proc = 0;
1070 for (j=0; j < nc; j++) {
1071 while (aj[j] >= cranges[proc+1]) proc++;
1072 col = aj[j] - cranges[proc] + Jranges[proc];
1073 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr); /* grad g */
1074 }
1075 ierr = MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1076 }
1077 }
1078 /* Jce_xfixed */
1079 if (pdipm->Nxfixed) {
1080 ierr = MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);CHKERRQ(ierr);
1081 for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1082 row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1083
1084 ierr = MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
1085 if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1086
1087 proc = 0;
1088 j = 0;
1089 while (cols[j] >= cranges[proc+1]) proc++;
1090 col = cols[j] - cranges[proc] + Jranges[proc];
1091 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1092 ierr = MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
1093 }
1094 }
1095
1096 /* 3rd Row block of KKT matrix: [ gradCi, 0, 0, -I] */
1097 if (pdipm->Nh) {
1098 ierr = MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);CHKERRQ(ierr);
1099 for (i=0; i < pdipm->nh; i++){
1100 row = rstart + pdipm->off_lambdai + i;
1101
1102 ierr = MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1103 proc = 0;
1104 for (j=0; j < nc; j++) {
1105 while (aj[j] >= cranges[proc+1]) proc++;
1106 col = aj[j] - cranges[proc] + Jranges[proc];
1107 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr); /* grad h */
1108 }
1109 ierr = MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);
1110 }
1111 /* -I */
1112 for (i=0; i < pdipm->nh; i++){
1113 row = rstart + pdipm->off_lambdai + i;
1114 col = rstart + pdipm->off_z + i;
1115 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1116 }
1117 }
1118
1119 /* Jci_xb */
1120 ierr = MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);CHKERRQ(ierr);
1121 for (i=0; i < (pdipm->nci - pdipm->nh); i++){
1122 row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1123
1124 ierr = MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
1125 if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1126 proc = 0;
1127 for (j=0; j < nc; j++) {
1128 while (cols[j] >= cranges[proc+1]) proc++;
1129 col = cols[j] - cranges[proc] + Jranges[proc];
1130 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1131 }
1132 ierr = MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);
1133 /* -I */
1134 col = rstart + pdipm->off_z + pdipm->nh + i;
1135 ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);
1136 }
1137
1138 /* 4-th Row block of KKT matrix: Z and Ci */
1139 for (i=0; i < pdipm->nci; i++) {
1140 row = rstart + pdipm->off_z + i;
1141 cols1[0] = rstart + pdipm->off_lambdai + i;
1142 cols1[1] = row;
1143 ierr = MatPreallocateSet(row,2,cols1,dnz,onz);CHKERRQ(ierr);
1144 }
1145
1146 /* diagonal entry */
1147 for (i=0; i<pdipm->n; i++) dnz[i]++; /* diagonal entry */
1148
1149 /* Create KKT matrix */
1150 ierr = MatCreate(comm,&J);CHKERRQ(ierr);
1151 ierr = MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
1152 ierr = MatSetFromOptions(J);CHKERRQ(ierr);
1153 ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr);
1154 ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr);
1155 /* ierr = MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); */
1156 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
1157 pdipm->K = J;
1158
1159 /* (8) Set up nonlinear solver SNES */
1160 ierr = SNESSetFunction(pdipm->snes,NULL,TaoSNESFunction_PDIPM,(void*)tao);CHKERRQ(ierr);
1161 ierr = SNESSetJacobian(pdipm->snes,J,J,TaoSNESJacobian_PDIPM,(void*)tao);CHKERRQ(ierr);
1162
1163 if (pdipm->solve_reduced_kkt) {
1164 PC pc;
1165 ierr = KSPGetPC(tao->ksp,&pc);CHKERRQ(ierr);
1166 ierr = PCSetType(pc,PCFIELDSPLIT);CHKERRQ(ierr);
1167 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr);
1168 ierr = PCFieldSplitSetIS(pc,"2",pdipm->is2);CHKERRQ(ierr);
1169 ierr = PCFieldSplitSetIS(pc,"1",pdipm->is1);CHKERRQ(ierr);
1170 }
1171 ierr = SNESSetFromOptions(pdipm->snes);CHKERRQ(ierr);
1172
1173 /* (9) Insert constant entries to K */
1174 /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1175 ierr = MatGetOwnershipRange(J,&rstart,&rend);CHKERRQ(ierr);
1176 for (i=rstart; i<rend; i++){
1177 ierr = MatSetValue(J,i,i,0.0,INSERT_VALUES);CHKERRQ(ierr);
1178 }
1179
1180 /* Row block of K: [ grad Ce, 0, 0, 0] */
1181 if (pdipm->Nxfixed) {
1182 ierr = MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);CHKERRQ(ierr);
1183 for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1184 row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1185
1186 ierr = MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);CHKERRQ(ierr);
1187 proc = 0;
1188 for (j=0; j < nc; j++) {
1189 while (cols[j] >= cranges[proc+1]) proc++;
1190 col = cols[j] - cranges[proc] + Jranges[proc];
1191 ierr = MatSetValue(J,row,col,aa[j],INSERT_VALUES);CHKERRQ(ierr); /* grad Ce */
1192 ierr = MatSetValue(J,col,row,aa[j],INSERT_VALUES);CHKERRQ(ierr); /* grad Ce' */
1193 }
1194 ierr = MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);CHKERRQ(ierr);
1195 }
1196 }
1197
1198 /* Row block of K: [ grad Ci, 0, 0, -I] */
1199 ierr = MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);CHKERRQ(ierr);
1200 for (i=0; i < pdipm->nci - pdipm->nh; i++){
1201 row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1202
1203 ierr = MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);CHKERRQ(ierr);
1204 proc = 0;
1205 for (j=0; j < nc; j++) {
1206 while (cols[j] >= cranges[proc+1]) proc++;
1207 col = cols[j] - cranges[proc] + Jranges[proc];
1208 ierr = MatSetValue(J,col,row,-aa[j],INSERT_VALUES);CHKERRQ(ierr);
1209 ierr = MatSetValue(J,row,col,aa[j],INSERT_VALUES);CHKERRQ(ierr);
1210 }
1211 ierr = MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);CHKERRQ(ierr);
1212
1213 col = rstart + pdipm->off_z + pdipm->nh + i;
1214 ierr = MatSetValue(J,row,col,-1,INSERT_VALUES);CHKERRQ(ierr);
1215 }
1216
1217 for (i=0; i < pdipm->nh; i++){
1218 row = rstart + pdipm->off_lambdai + i;
1219 col = rstart + pdipm->off_z + i;
1220 ierr = MatSetValue(J,row,col,-1,INSERT_VALUES);CHKERRQ(ierr);
1221 }
1222
1223 if (pdipm->Nxfixed) {
1224 ierr = MatDestroy(&Jce_xfixed_trans);CHKERRQ(ierr);
1225 }
1226 ierr = MatDestroy(&Jci_xb_trans);CHKERRQ(ierr);
1227 ierr = PetscFree3(ng_all,nh_all,Jranges);CHKERRQ(ierr);
1228 PetscFunctionReturn(0);
1229 }
1230
1231 /*
1232 TaoDestroy_PDIPM - Destroys the pdipm object
1233
1234 Input:
1235 full pdipm
1236
1237 Output:
1238 Destroyed pdipm
1239 */
TaoDestroy_PDIPM(Tao tao)1240 PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1241 {
1242 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1243 PetscErrorCode ierr;
1244
1245 PetscFunctionBegin;
1246 /* Freeing Vectors assocaiated with KKT (X) */
1247 ierr = VecDestroy(&pdipm->x);CHKERRQ(ierr); /* Solution x */
1248 ierr = VecDestroy(&pdipm->lambdae);CHKERRQ(ierr); /* Equality constraints lagrangian multiplier*/
1249 ierr = VecDestroy(&pdipm->lambdai);CHKERRQ(ierr); /* Inequality constraints lagrangian multiplier*/
1250 ierr = VecDestroy(&pdipm->z);CHKERRQ(ierr); /* Slack variables */
1251 ierr = VecDestroy(&pdipm->X);CHKERRQ(ierr); /* Big KKT system vector [x; lambdae; lambdai; z] */
1252
1253 /* work vectors */
1254 ierr = VecDestroy(&pdipm->lambdae_xfixed);CHKERRQ(ierr);
1255 ierr = VecDestroy(&pdipm->lambdai_xb);CHKERRQ(ierr);
1256
1257 /* Legrangian equality and inequality Vec */
1258 ierr = VecDestroy(&pdipm->ce);CHKERRQ(ierr); /* Vec of equality constraints */
1259 ierr = VecDestroy(&pdipm->ci);CHKERRQ(ierr); /* Vec of inequality constraints */
1260
1261 /* Matrices */
1262 ierr = MatDestroy(&pdipm->Jce_xfixed);CHKERRQ(ierr);
1263 ierr = MatDestroy(&pdipm->Jci_xb);CHKERRQ(ierr); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb); J(nxbx)] */
1264 ierr = MatDestroy(&pdipm->K);CHKERRQ(ierr);
1265
1266 /* Index Sets */
1267 if (pdipm->Nxub) {
1268 ierr = ISDestroy(&pdipm->isxub);CHKERRQ(ierr); /* Finite upper bound only -inf < x < ub */
1269 }
1270
1271 if (pdipm->Nxlb) {
1272 ierr = ISDestroy(&pdipm->isxlb);CHKERRQ(ierr); /* Finite lower bound only lb <= x < inf */
1273 }
1274
1275 if (pdipm->Nxfixed) {
1276 ierr = ISDestroy(&pdipm->isxfixed);CHKERRQ(ierr); /* Fixed variables lb = x = ub */
1277 }
1278
1279 if (pdipm->Nxbox) {
1280 ierr = ISDestroy(&pdipm->isxbox);CHKERRQ(ierr); /* Boxed variables lb <= x <= ub */
1281 }
1282
1283 if (pdipm->Nxfree) {
1284 ierr = ISDestroy(&pdipm->isxfree);CHKERRQ(ierr); /* Free variables -inf <= x <= inf */
1285 }
1286
1287 if (pdipm->solve_reduced_kkt) {
1288 ierr = ISDestroy(&pdipm->is1);CHKERRQ(ierr);
1289 ierr = ISDestroy(&pdipm->is2);CHKERRQ(ierr);
1290 }
1291
1292 /* SNES */
1293 ierr = SNESDestroy(&pdipm->snes);CHKERRQ(ierr); /* Nonlinear solver */
1294 ierr = PetscFree(pdipm->nce_all);CHKERRQ(ierr);
1295 ierr = MatDestroy(&pdipm->jac_equality_trans);CHKERRQ(ierr);
1296 ierr = MatDestroy(&pdipm->jac_inequality_trans);CHKERRQ(ierr);
1297
1298 /* Destroy pdipm */
1299 ierr = PetscFree(tao->data);CHKERRQ(ierr); /* Holding locations of pdipm */
1300
1301 /* Destroy Dual */
1302 ierr = VecDestroy(&tao->DE);CHKERRQ(ierr); /* equality dual */
1303 ierr = VecDestroy(&tao->DI);CHKERRQ(ierr); /* dinequality dual */
1304 PetscFunctionReturn(0);
1305 }
1306
TaoSetFromOptions_PDIPM(PetscOptionItems * PetscOptionsObject,Tao tao)1307 PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1308 {
1309 TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1310 PetscErrorCode ierr;
1311
1312 PetscFunctionBegin;
1313 ierr = PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");CHKERRQ(ierr);
1314 ierr = PetscOptionsReal("-tao_pdipm_push_init_slack","parameter to push initial slack variables away from bounds",NULL,pdipm->push_init_slack,&pdipm->push_init_slack,NULL);CHKERRQ(ierr);
1315 ierr = PetscOptionsReal("-tao_pdipm_push_init_lambdai","parameter to push initial (inequality) dual variables away from bounds",NULL,pdipm->push_init_lambdai,&pdipm->push_init_lambdai,NULL);CHKERRQ(ierr);
1316 ierr = PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);CHKERRQ(ierr);
1317 ierr = PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);CHKERRQ(ierr);
1318 ierr = PetscOptionsTail();CHKERRQ(ierr);
1319 PetscFunctionReturn(0);
1320 }
1321
1322 /*MC
1323 TAOPDIPM - Barrier-based primal-dual interior point algorithm for generally constrained optimization.
1324
1325 Option Database Keys:
1326 + -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1327 . -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1328 - -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1329
1330 Level: beginner
1331 M*/
TaoCreate_PDIPM(Tao tao)1332 PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1333 {
1334 TAO_PDIPM *pdipm;
1335 PetscErrorCode ierr;
1336
1337 PetscFunctionBegin;
1338 tao->ops->setup = TaoSetup_PDIPM;
1339 tao->ops->solve = TaoSolve_PDIPM;
1340 tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1341 tao->ops->destroy = TaoDestroy_PDIPM;
1342
1343 ierr = PetscNewLog(tao,&pdipm);CHKERRQ(ierr);
1344 tao->data = (void*)pdipm;
1345
1346 pdipm->nx = pdipm->Nx = 0;
1347 pdipm->nxfixed = pdipm->Nxfixed = 0;
1348 pdipm->nxlb = pdipm->Nxlb = 0;
1349 pdipm->nxub = pdipm->Nxub = 0;
1350 pdipm->nxbox = pdipm->Nxbox = 0;
1351 pdipm->nxfree = pdipm->Nxfree = 0;
1352
1353 pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1354 pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1355 pdipm->n = pdipm->N = 0;
1356 pdipm->mu = 1.0;
1357 pdipm->mu_update_factor = 0.1;
1358
1359 pdipm->push_init_slack = 1.0;
1360 pdipm->push_init_lambdai = 1.0;
1361 pdipm->solve_reduced_kkt = PETSC_FALSE;
1362
1363 /* Override default settings (unless already changed) */
1364 if (!tao->max_it_changed) tao->max_it = 200;
1365 if (!tao->max_funcs_changed) tao->max_funcs = 500;
1366
1367 ierr = SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);CHKERRQ(ierr);
1368 ierr = SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);CHKERRQ(ierr);
1369 ierr = SNESGetKSP(pdipm->snes,&tao->ksp);CHKERRQ(ierr);
1370 ierr = PetscObjectReference((PetscObject)tao->ksp);CHKERRQ(ierr);
1371 PetscFunctionReturn(0);
1372 }
1373