Actual source code: pounders.c
1: #include "pounders.h"
5: static PetscErrorCode pounders_h(TaoSolver subtao, Vec v, Mat *H, Mat *Hpre, MatStructure *flag, void *ctx)
6: {
8: *flag = SAME_NONZERO_PATTERN;
9: return(0);
10: }
13: static PetscErrorCode pounders_fg(TaoSolver subtao, Vec x, PetscReal *f, Vec g, void *ctx)
14: {
15: TAO_POUNDERS *mfqP = (TAO_POUNDERS*)ctx;
16: PetscReal d1,d2;
19: /* g = A*x (add b later)*/
20: MatMult(mfqP->Hs,x,g);
23: /* f = 1/2 * x'*(Ax) + b'*x */
24: VecDot(x,g,&d1);
25: VecDot(mfqP->b,x,&d2);
26: *f = 0.5 *d1 + d2;
28: /* now g = g + b */
29: VecAXPY(g, 1.0, mfqP->b);
31: return(0);
32: }
35: PetscErrorCode gqtwrap(TaoSolver tao,PetscReal *gnorm, PetscReal *qmin)
36: {
38: PetscReal atol=1.0e-10;
39: PetscInt info,its;
40: TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data;
41:
44: if (! mfqP->usegqt) {
45: Vec x,xl,xu,pdel,ndel;
46: PetscInt i,j;
47: VecCreateSeqWithArray(PETSC_COMM_SELF,mfqP->n,mfqP->Xsubproblem,&x);
48: VecCreateSeqWithArray(PETSC_COMM_SELF,mfqP->n,mfqP->Gres,&mfqP->b);
49: VecCreateSeq(PETSC_COMM_SELF,mfqP->n,&xl);
50: VecDuplicate(xl,&xu);
51: VecDuplicate(xl,&pdel);
52: VecDuplicate(xl,&ndel);
53: VecSet(x,0.0);
54: VecSet(ndel,-mfqP->delta);
55: VecSet(pdel,mfqP->delta);
56: for (i=0;i<mfqP->n;i++) {
57: for (j=i;j<mfqP->n;j++) {
58: mfqP->Hres[j+mfqP->n*i] = mfqP->Hres[mfqP->n*j+i];
59: }
60: }
61: MatCreateSeqDense(PETSC_COMM_SELF,mfqP->n,mfqP->n,mfqP->Hres,&mfqP->Hs);
62:
63: TaoResetStatistics(mfqP->subtao);
64: TaoSetInitialVector(mfqP->subtao,x);
65: TaoSetHessianRoutine(mfqP->subtao,mfqP->Hs,mfqP->Hs,pounders_h,(void*)mfqP);
66: TaoSetTolerances(mfqP->subtao,PETSC_NULL,PETSC_NULL,*gnorm,*gnorm,PETSC_NULL);
67: /* enforce bound constraints -- experimental */
68: if (tao->XU && tao->XL) {
69: VecCopy(tao->XU,xu);
70: VecAXPY(xu,-1.0,tao->solution);
71: VecScale(xu,1.0/mfqP->delta);
72: VecCopy(tao->XL,xl);
73: VecAXPY(xl,-1.0,tao->solution);
74: VecScale(xl,1.0/mfqP->delta);
75:
76: VecPointwiseMin(xu,xu,pdel);
77: VecPointwiseMax(xl,xl,ndel);
78: }
79: TaoSetVariableBounds(mfqP->subtao,xl,xu);
80: TaoSolve(mfqP->subtao);
81: TaoGetSolutionStatus(mfqP->subtao,PETSC_NULL,qmin,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
82: VecDestroy(&x);
83: VecDestroy(&xl);
84: VecDestroy(&xu);
85: VecDestroy(&pdel);
86: VecDestroy(&ndel);
87: VecDestroy(&mfqP->b);
88: MatDestroy(&mfqP->Hs);
89:
91: } else {
92: gqt(mfqP->n,mfqP->Hres,mfqP->n,mfqP->Gres,1.0,mfqP->gqt_rtol,atol,
93: mfqP->gqt_maxits,gnorm,qmin,mfqP->Xsubproblem,&info,&its,
94: mfqP->work,mfqP->work2, mfqP->work3);
95: /*dgqt_(&mfqP->n,mfqP->Hres,&mfqP->n,mfqP->Gres,&one,&mfqP->gqt_rtol,&atol,
96: &mfqP->gqt_maxits,gnorm,qmin,mfqP->Xsubproblem,&info,&its,
97: mfqP->work,mfqP->work2, mfqP->work3);*/
98: }
99: *qmin *= -1;
100: return(0);
101: }
106: PetscErrorCode phi2eval(PetscReal *x, PetscInt n, PetscReal *phi) {
107: /* Phi = .5*[x(1)^2 sqrt(2)*x(1)*x(2) ... sqrt(2)*x(1)*x(n) ... x(2)^2 sqrt(2)*x(2)*x(3) .. x(n)^2] */
108: PetscInt i,j,k;
109: PetscReal sqrt2 = PetscSqrtReal(2.0);
111: j=0;
113: for (i=0;i<n;i++) {
114: phi[j] = 0.5 * x[i]*x[i];
115: j++;
116: for (k=i+1;k<n;k++) {
117: phi[j] = x[i]*x[k]/sqrt2;
118: j++;
119: }
120:
121: }
123: return(0);
124: }
128: PetscErrorCode getquadpounders(TAO_POUNDERS *mfqP) {
129: /* Computes the parameters of the quadratic Q(x) = c + g'*x + 0.5*x*G*x'
130: that satisfies the interpolation conditions Q(X[:,j]) = f(j)
131: for j=1,...,m and with a Hessian matrix of least Frobenius norm */
133: /* NB --we are ignoring c */
134: PetscInt i,j,k,num,np = mfqP->nmodelpoints;
135: PetscReal one = 1.0,zero=0.0,negone=-1.0;
136: PetscBLASInt blasnpmax = mfqP->npmax;
137: PetscBLASInt blasnplus1 = mfqP->n+1;
138: PetscBLASInt blasnp = np;
139: PetscBLASInt blasint = mfqP->n*(mfqP->n+1) / 2;
140: PetscBLASInt blasint2 = np - mfqP->n-1;
141: PetscBLASInt blasinfo,ione=1;
142: PetscReal sqrt2 = PetscSqrtReal(2.0);
143:
146: for (i=0;i<mfqP->n*mfqP->m;i++) {
147: mfqP->Gdel[i] = 0;
148: }
149: for (i=0;i<mfqP->n*mfqP->n*mfqP->m;i++) {
150: mfqP->Hdel[i] = 0;
151: }
153: /* Let Ltmp = (L'*L) */
154: BLASgemm_("T","N",&blasint2,&blasint2,&blasint,&one,&mfqP->L[(mfqP->n+1)*blasint],&blasint,&mfqP->L[(mfqP->n+1)*blasint],&blasint,&zero,mfqP->L_tmp,&blasint);
155:
156: /* factor Ltmp */
157: LAPACKpotrf_("L",&blasint2,mfqP->L_tmp,&blasint,&blasinfo);
158: if (blasinfo != 0) {
159: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine potrf returned with value %D\n",blasinfo);
160: }
162: /* factor M */
163: LAPACKgetrf_(&blasnplus1,&blasnpmax,mfqP->M,&blasnplus1,mfqP->npmaxiwork,&blasinfo);
164: if (blasinfo != 0) {
165: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine getrf returned with value %D\n",blasinfo);
166: }
167:
168: for (k=0;k<mfqP->m;k++) {
169: /* Solve L'*L*Omega = Z' * RESk*/
170: BLASgemv_("T",&blasnp,&blasint2,&one,mfqP->Z,&blasnpmax,&mfqP->RES[mfqP->npmax*k],&ione,&zero,mfqP->omega,&ione);
172: LAPACKpotrs_("L",&blasint2,&ione,mfqP->L_tmp,&blasint,mfqP->omega,&blasint2,&blasinfo);
173: if (blasinfo != 0) {
174: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine potrs returned with value %D\n",blasinfo);
175: }
176:
177:
178: /* Beta = L*Omega */
179: BLASgemv_("N",&blasint,&blasint2,&one,&mfqP->L[(mfqP->n+1)*blasint],&blasint,mfqP->omega,&ione,&zero,mfqP->beta,&ione);
180:
181: /* solve M'*Alpha = RESk - N'*Beta */
182: BLASgemv_("T",&blasint,&blasnp,&negone,mfqP->N,&blasint,mfqP->beta,&ione,&one,&mfqP->RES[mfqP->npmax*k],&ione);
183: LAPACKgetrs_("T",&blasnplus1,&ione,mfqP->M,&blasnplus1,mfqP->npmaxiwork,&mfqP->RES[mfqP->npmax*k],&blasnplus1,&blasinfo);
184: if (blasinfo != 0) {
185: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine getrs returned with value %D\n",blasinfo);
186: }
188: /* Gdel(:,k) = Alpha(2:n+1) */
189: for (i=0;i<mfqP->n;i++) {
190: mfqP->Gdel[i + mfqP->n*k] = mfqP->RES[mfqP->npmax*k + i+1];
191: }
192:
193: /* Set Hdels */
194: num=0;
195: for (i=0;i<mfqP->n;i++) {
196: /* H[i,i,k] = Beta(num) */
197: mfqP->Hdel[(i*mfqP->n + i)*mfqP->m + k] = mfqP->beta[num];
198: num++;
199: for (j=i+1;j<mfqP->n;j++) {
200: /* H[i,j,k] = H[j,i,k] = Beta(num)/sqrt(2) */
201: mfqP->Hdel[(j*mfqP->n + i)*mfqP->m + k] = mfqP->beta[num]/sqrt2;
202: mfqP->Hdel[(i*mfqP->n + j)*mfqP->m + k] = mfqP->beta[num]/sqrt2;
203: num++;
204: }
205: }
206: }
207:
208: return(0);
209: }
213: PetscErrorCode morepoints(TAO_POUNDERS *mfqP) {
214: /* Assumes mfqP->model_indices[0] is minimum index
215: Finishes adding points to mfqP->model_indices (up to npmax)
216: Computes L,Z,M,N
217: np is actual number of points in model (should equal npmax?) */
218: PetscInt point,i,j,offset;
219: PetscInt reject;
220: PetscBLASInt blasn=mfqP->n,blasnpmax=mfqP->npmax,blasnplus1=mfqP->n+1,blasinfo,blasnpmax_x_5=mfqP->npmax*5,blasint,blasint2,blasnp,blasmaxmn;
221: PetscReal *x,normd;
225: for (i=0;i<mfqP->n+1;i++) {
226: VecGetArray(mfqP->Xhist[mfqP->model_indices[i]],&x);
227: mfqP->M[(mfqP->n+1)*i] = 1.0;
228: for (j=0;j<mfqP->n;j++) {
229: mfqP->M[j+1+((mfqP->n+1)*i)] = (x[j] - mfqP->xmin[j]) / mfqP->delta;
230: }
231: VecRestoreArray(mfqP->Xhist[mfqP->model_indices[i]],&x);
232: phi2eval(&mfqP->M[1+((mfqP->n+1)*i)],mfqP->n,&mfqP->N[mfqP->n*(mfqP->n+1)/2 * i]);
233:
236: }
238: /* Now we add points until we have npmax starting with the most recent ones */
239: point = mfqP->nHist-1;
240: mfqP->nmodelpoints = mfqP->n+1;
241:
242: while (mfqP->nmodelpoints < mfqP->npmax && point>=0) {
243: /* Reject any points already in the model */
244: reject = 0;
245: for (j=0;j<mfqP->n+1;j++) {
246: if (point == mfqP->model_indices[j]) {
247: reject = 1;
248: break;
249: }
250: }
252: /* Reject if norm(d) >c2 */
253: if (!reject) {
254: VecCopy(mfqP->Xhist[point],mfqP->workxvec);
255: VecAXPY(mfqP->workxvec,-1.0,mfqP->Xhist[mfqP->minindex]);
256: VecNorm(mfqP->workxvec,NORM_2,&normd);
257: normd /= mfqP->delta;
258: if (normd > mfqP->c2) {
259: reject =1;
260: }
261: }
262: if (reject)
263: {
264: point--;
265: continue;
266: }
268:
269: VecGetArray(mfqP->Xhist[point],&x);
270: mfqP->M[(mfqP->n+1)*mfqP->nmodelpoints] = 1.0;
271: for (j=0;j<mfqP->n;j++) {
272: mfqP->M[j+1+((mfqP->n+1)*mfqP->nmodelpoints)] = (x[j] - mfqP->xmin[j]) / mfqP->delta;
273: }
275: VecRestoreArray(mfqP->Xhist[point],&x);
276: phi2eval(&mfqP->M[1+(mfqP->n+1)*mfqP->nmodelpoints],mfqP->n,&mfqP->N[mfqP->n*(mfqP->n+1)/2 * (mfqP->nmodelpoints)]);
278: /* Update QR factorization */
280: /* Copy M' to Q_tmp */
281: for (i=0;i<mfqP->n+1;i++) {
282: for (j=0;j<mfqP->npmax;j++) {
283: mfqP->Q_tmp[j+mfqP->npmax*i] = mfqP->M[i+(mfqP->n+1)*j];
284: }
285: }
286: blasnp = mfqP->nmodelpoints+1;
287: /* Q_tmp,R = qr(M') */
288: blasmaxmn=PetscMax(mfqP->m,mfqP->n+1);
289: LAPACKgeqrf_(&blasnp,&blasnplus1,mfqP->Q_tmp,&blasnpmax,mfqP->tau_tmp,mfqP->mwork,&blasmaxmn,&blasinfo);
291: if (blasinfo != 0) {
292: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine geqrf returned with value %D\n",blasinfo);
293: }
294:
295: /* Reject if min(svd(N*Q(:,n+2:np+1)) <= theta2 */
296: /* L = N*Qtmp */
297: blasint2 = mfqP->n * (mfqP->n+1) / 2;
298: /* Copy N to L_tmp */
299: for (i=0;i<mfqP->n*(mfqP->n+1)/2 * mfqP->npmax;i++) {
300: mfqP->L_tmp[i]= mfqP->N[i];
301: }
302:
303: /* Copy L_save to L_tmp */
305: /* L_tmp = N*Qtmp' */
306: LAPACKormqr_("R","N",&blasint2,&blasnp,&blasnplus1,mfqP->Q_tmp,&blasnpmax,mfqP->tau_tmp,
307: mfqP->L_tmp,&blasint2,mfqP->npmaxwork,&blasnpmax_x_5,&blasinfo);
309: if (blasinfo != 0) {
310: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine ormqr returned with value %D\n",blasinfo);
311: }
313: /* Copy L_tmp to L_save */
314: for (i=0;i<mfqP->npmax * mfqP->n*(mfqP->n+1)/2;i++) {
315: mfqP->L_save[i] = mfqP->L_tmp[i];
316: }
317:
318: /* Get svd for L_tmp(:,n+2:np+1) (L_tmp is modified in process) */
319: blasint = mfqP->nmodelpoints - mfqP->n;
320: LAPACKgesvd_("N","N",&blasint2,&blasint,&mfqP->L_tmp[(mfqP->n+1)*blasint2],&blasint2,
321: mfqP->beta,mfqP->work,&blasn,mfqP->work,&blasn,mfqP->npmaxwork,&blasnpmax_x_5,
322: &blasinfo);
323: if (blasinfo != 0) {
324: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine gesvd returned with value %D\n",blasinfo);
325: }
327: if (mfqP->beta[PetscMin(blasint,blasint2)-1] > mfqP->theta2) {
328: /* accept point */
329: mfqP->model_indices[mfqP->nmodelpoints] = point;
330: /* Copy Q_tmp to Q */
331: for (i=0;i<mfqP->npmax* mfqP->npmax;i++) {
332: mfqP->Q[i] = mfqP->Q_tmp[i];
333: }
334: for (i=0;i<mfqP->npmax;i++){
335: mfqP->tau[i] = mfqP->tau_tmp[i];
336: }
337: mfqP->nmodelpoints++;
338: blasnp = mfqP->nmodelpoints+1;
340: /* Copy L_save to L */
341: for (i=0;i<mfqP->npmax * mfqP->n*(mfqP->n+1)/2;i++) {
342: mfqP->L[i] = mfqP->L_save[i];
343: }
344: }
345: point--;
347: }
350: blasnp = mfqP->nmodelpoints;
351: /* Copy Q(:,n+2:np) to Z */
352: /* First set Q_tmp to I */
353: for (i=0;i<mfqP->npmax*mfqP->npmax;i++) {
354: mfqP->Q_tmp[i] = 0.0;
355: }
356: for (i=0;i<mfqP->npmax;i++) {
357: mfqP->Q_tmp[i + mfqP->npmax*i] = 1.0;
358: }
360: /* Q_tmp = I * Q */
361: LAPACKormqr_("R","N",&blasnp,&blasnp,&blasnplus1,mfqP->Q,&blasnpmax,mfqP->tau,mfqP->Q_tmp,&blasnpmax,mfqP->npmaxwork,&blasnpmax_x_5,&blasinfo);
363: if (blasinfo != 0) {
364: SETERRQ1(PETSC_COMM_SELF,1,"LAPACK routine ormqr returned with value %D\n",blasinfo);
365: }
367: /* Copy Q_tmp(:,n+2:np) to Z) */
368: offset = mfqP->npmax * (mfqP->n+1);
369: for (i=offset;i<mfqP->npmax*mfqP->npmax;i++) {
370: mfqP->Z[i-offset] = mfqP->Q_tmp[i];
371: }
372: return(0);
373: }
377: /* Only call from modelimprove, addpoint() needs ->Q_tmp and ->work to be set */
378: PetscErrorCode addpoint(TaoSolver tao, TAO_POUNDERS *mfqP, PetscInt index)
379: {
382: /* Create new vector in history: X[newidx] = X[mfqP->index] + delta*X[index]*/
383:
384: VecDuplicate(mfqP->Xhist[0],&mfqP->Xhist[mfqP->nHist]);
385: VecSetValues(mfqP->Xhist[mfqP->nHist],mfqP->n,mfqP->indices,&mfqP->Q_tmp[index*mfqP->npmax],INSERT_VALUES);
386: VecAssemblyBegin(mfqP->Xhist[mfqP->nHist]);
387: VecAssemblyEnd(mfqP->Xhist[mfqP->nHist]);
388: VecAYPX(mfqP->Xhist[mfqP->nHist],mfqP->delta,mfqP->Xhist[mfqP->minindex]);
390: /* Project into feasible region */
391: if (tao->XU && tao->XL) {
392: VecMedian(mfqP->Xhist[mfqP->nHist], tao->XL, tao->XU,
393: mfqP->Xhist[mfqP->nHist]);
394: }
396: /* Compute value of new vector */
397: VecDuplicate(mfqP->Fhist[0],&mfqP->Fhist[mfqP->nHist]);
398: CHKMEMQ;
399: TaoComputeSeparableObjective(tao,mfqP->Xhist[mfqP->nHist],mfqP->Fhist[mfqP->nHist]);
400: VecNorm(mfqP->Fhist[mfqP->nHist],NORM_2,&mfqP->Fres[mfqP->nHist]);
401: mfqP->Fres[mfqP->nHist]*=mfqP->Fres[mfqP->nHist];
403: /* Add new vector to model */
404: mfqP->model_indices[mfqP->nmodelpoints] = mfqP->nHist;
405: mfqP->nmodelpoints++;
406: mfqP->nHist++;
408:
409: return(0);
410: }
414: PetscErrorCode modelimprove(TaoSolver tao, TAO_POUNDERS *mfqP, PetscInt addallpoints) {
415: /* modeld = Q(:,np+1:n)' */
417: PetscInt i,j,minindex=0;
418: PetscReal dp,half=0.5,one=1.0,minvalue=TAO_INFINITY;
419: PetscBLASInt blasn=mfqP->n, blasnpmax = mfqP->npmax, blask,info;
420: PetscBLASInt blas1=1,blasnpmax_x_5 = mfqP->npmax*5;
421: blask = mfqP->nmodelpoints;
423:
424: /* Qtmp = I(n x n) */
425: for (i=0;i<mfqP->n;i++) {
426: for (j=0;j<mfqP->n;j++) {
427: mfqP->Q_tmp[i + mfqP->npmax*j] = 0.0;
428: }
429: }
430: for (j=0;j<mfqP->n;j++) {
431: mfqP->Q_tmp[j + mfqP->npmax*j] = 1.0;
432: }
433:
435: /* Qtmp = Q * I */
436: LAPACKormqr_("R","N",&blasn,&blasn,&blask,mfqP->Q,&blasnpmax,
437: mfqP->tau, mfqP->Q_tmp, &blasnpmax, mfqP->npmaxwork,&blasnpmax_x_5, &info);
438:
439: for (i=mfqP->nmodelpoints;i<mfqP->n;i++) {
440: dp = BLASdot_(&blasn,&mfqP->Q_tmp[i*mfqP->npmax],&blas1,mfqP->Gres,&blas1);
441: if (dp>0.0) { /* Model says use the other direction! */
442: for (j=0;j<mfqP->n;j++) {
443: mfqP->Q_tmp[i*mfqP->npmax+j] *= -1;
444: }
445: }
446: /* mfqP->work[i] = Cres+Modeld(i,:)*(Gres+.5*Hres*Modeld(i,:)') */
447: for (j=0;j<mfqP->n;j++) {
448: mfqP->work2[j] = mfqP->Gres[j];
449: }
450: BLASgemv_("N",&blasn,&blasn,&half,mfqP->Hres,&blasn,
451: &mfqP->Q_tmp[i*mfqP->npmax], &blas1, &one, mfqP->work2,&blas1);
452: mfqP->work[i] = BLASdot_(&blasn,&mfqP->Q_tmp[i*mfqP->npmax],&blas1,
453: mfqP->work2,&blas1);
454: if (i==mfqP->nmodelpoints || mfqP->work[i] < minvalue) {
455: minindex=i;
456: minvalue = mfqP->work[i];
457: }
458:
459: if (addallpoints != 0) {
460: addpoint(tao,mfqP,i);
461: }
463: }
464:
465: if (!addallpoints) {
466: addpoint(tao,mfqP,minindex);
467: }
468: return(0);
469: }
474: PetscErrorCode affpoints(TAO_POUNDERS *mfqP, PetscReal *xmin,
475: PetscReal c) {
476: PetscInt i,j;
477: PetscBLASInt blasm=mfqP->m,blasj,blask,blasn=mfqP->n,ione=1,info;
478: PetscBLASInt blasnpmax = mfqP->npmax,blasmaxmn;
479: PetscReal proj,normd;
480: PetscReal *x;
484: for (i=mfqP->nHist-1;i>=0;i--) {
485: VecGetArray(mfqP->Xhist[i],&x);
486: for (j=0;j<mfqP->n;j++) {
487: mfqP->work[j] = (x[j] - xmin[j])/mfqP->delta;
488: }
489: VecRestoreArray(mfqP->Xhist[i],&x);
490: BLAScopy_(&blasn,mfqP->work,&ione,mfqP->work2,&ione);
491: normd = BLASnrm2_(&blasn,mfqP->work,&ione);
492: if (normd <= c*c) {
493: blasj=(mfqP->n - mfqP->nmodelpoints);
494: if (!mfqP->q_is_I) {
495: /* project D onto null */
496: blask=(mfqP->nmodelpoints);
497: LAPACKormqr_("R","N",&ione,&blasn,&blask,mfqP->Q,&blasnpmax,mfqP->tau,
498: mfqP->work2,&ione,mfqP->mwork,&blasm,&info);
499:
500: if (info < 0) {
501: SETERRQ1(PETSC_COMM_SELF,1,"ormqr returned value %D\n",info);
502: }
503: }
504: proj = BLASnrm2_(&blasj,&mfqP->work2[mfqP->nmodelpoints],&ione);
506: if (proj >= mfqP->theta1) { /* add this index to model */
507: mfqP->model_indices[mfqP->nmodelpoints]=i;
508: mfqP->nmodelpoints++;
509: BLAScopy_(&blasn,mfqP->work,&ione,&mfqP->Q_tmp[mfqP->npmax*(mfqP->nmodelpoints-1)],&ione);
510: blask=mfqP->npmax*(mfqP->nmodelpoints);
511: BLAScopy_(&blask,mfqP->Q_tmp,&ione,mfqP->Q,&ione);
512: blask = mfqP->nmodelpoints;
513: blasmaxmn = PetscMax(mfqP->m,mfqP->n);
514: LAPACKgeqrf_(&blasn,&blask,mfqP->Q,&blasnpmax,mfqP->tau,mfqP->mwork,
515: &blasmaxmn,&info);
516: mfqP->q_is_I = 0;
517: if (info < 0) {
518: SETERRQ1(PETSC_COMM_SELF,1,"geqrf returned value %D\n",info);
519: }
521:
522: }
523: if (mfqP->nmodelpoints == mfqP->n) {
524: break;
525: }
526: }
527: }
528:
529: return(0);
530: }
533: static PetscErrorCode TaoSolve_POUNDERS(TaoSolver tao)
534: {
535: TAO_POUNDERS *mfqP = (TAO_POUNDERS *)tao->data;
537: PetscInt i,ii,j,k,l,iter=0;
538: PetscReal step=1.0;
539: TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING;
541: PetscInt low,high;
542: PetscReal minnorm;
543: PetscReal *x,*f,*fmin,*xmint;
544: PetscReal cres,deltaold;
545: PetscReal gnorm;
546: PetscBLASInt info,ione=1,iblas;
547: PetscBool valid;
548: PetscReal mdec, rho, normxsp;
549: PetscReal one=1.0,zero=0.0,ratio;
550: PetscBLASInt blasm,blasn,blasnpmax,blasn2;
552:
553:
554: /* n = # of parameters
555: m = dimension (components) of function */
556:
558:
559: if (tao->XL && tao->XU) {
560: /* Check x0 + delta < XU */
561: PetscReal val;
562: VecSet(mfqP->Xhist[0],mfqP->delta);
563: VecAXPY(mfqP->Xhist[0],1.0,tao->solution);
564: VecAXPY(mfqP->Xhist[0],-1.0,tao->XU);
565: VecMin(mfqP->Xhist[0],PETSC_NULL,&val);
566: if (val > 1e-10) {
567: SETERRQ(PETSC_COMM_SELF,1,"X0 + delta > upper bound");
568: }
569: /* Check x0 > xl */
570: VecCopy(tao->solution,mfqP->Xhist[0]);
571: VecAXPY(mfqP->Xhist[0],-1.0,tao->XL);
572: VecMax(mfqP->Xhist[0],PETSC_NULL,&val);
573: if (val < -1e-10) {
574: SETERRQ(PETSC_COMM_SELF,1,"X0 + delta > upper bound");
575: }
576:
577: }
579: CHKMEMQ;
580: blasm = mfqP->m; blasn=mfqP->n; blasnpmax = mfqP->npmax;
581: for (i=0;i<mfqP->n*mfqP->n*mfqP->m;i++) mfqP->H[i]=0;
582:
583: VecCopy(tao->solution,mfqP->Xhist[0]);
584: CHKMEMQ;
585: TaoComputeSeparableObjective(tao,tao->solution,mfqP->Fhist[0]);
586: VecNorm(mfqP->Fhist[0],NORM_2,&mfqP->Fres[0]);
587: mfqP->Fres[0]*=mfqP->Fres[0];
588: mfqP->minindex = 0;
589: minnorm = mfqP->Fres[0];
590: VecGetOwnershipRange(mfqP->Xhist[0],&low,&high);
591: for (i=1;i<mfqP->n+1;i++) {
592: VecCopy(tao->solution,mfqP->Xhist[i]);
593: if (i-1 >= low && i-1 < high) {
594: VecGetArray(mfqP->Xhist[i],&x);
595: x[i-1-low] += mfqP->delta;
596: VecRestoreArray(mfqP->Xhist[i],&x);
597: }
598: CHKMEMQ;
599: TaoComputeSeparableObjective(tao,mfqP->Xhist[i],mfqP->Fhist[i]);
600: VecNorm(mfqP->Fhist[i],NORM_2,&mfqP->Fres[i]);
601: mfqP->Fres[i]*=mfqP->Fres[i];
602: if (mfqP->Fres[i] < minnorm) {
603: mfqP->minindex = i;
604: minnorm = mfqP->Fres[i];
605: }
606: }
608: VecCopy(mfqP->Xhist[mfqP->minindex],tao->solution);
609: VecCopy(mfqP->Fhist[mfqP->minindex],tao->sep_objective);
610: /* Gather mpi vecs to one big local vec */
612:
614: /* Begin serial code */
616: /* Disp[i] = Xi-xmin, i=1,..,mfqP->minindex-1,mfqP->minindex+1,..,n */
617: /* Fdiff[i] = (Fi-Fmin)', i=1,..,mfqP->minindex-1,mfqP->minindex+1,..,n */
618: /* (Column oriented for blas calls) */
619: ii=0;
621: if (mfqP->mpisize == 1) {
622: VecGetArray(mfqP->Xhist[mfqP->minindex],&xmint);
623: for (i=0;i<mfqP->n;i++) mfqP->xmin[i] = xmint[i];
624: VecRestoreArray(mfqP->Xhist[mfqP->minindex],&xmint);
625: VecGetArray(mfqP->Fhist[mfqP->minindex],&fmin);
626: for (i=0;i<mfqP->n+1;i++) {
627: if (i == mfqP->minindex) continue;
629: VecGetArray(mfqP->Xhist[i],&x);
630: for (j=0;j<mfqP->n;j++) {
631: mfqP->Disp[ii+mfqP->npmax*j] = (x[j] - mfqP->xmin[j])/mfqP->delta;
632: }
633: VecRestoreArray(mfqP->Xhist[i],&x);
635: VecGetArray(mfqP->Fhist[i],&f);
636: for (j=0;j<mfqP->m;j++) {
637: mfqP->Fdiff[ii+mfqP->n*j] = f[j] - fmin[j];
638: }
639: VecRestoreArray(mfqP->Fhist[i],&f);
640: mfqP->model_indices[ii++] = i;
642: }
643: for (j=0;j<mfqP->m;j++) {
644: mfqP->C[j] = fmin[j];
645: }
646: VecRestoreArray(mfqP->Fhist[mfqP->minindex],&fmin);
648: } else {
649: VecScatterBegin(mfqP->scatterx,mfqP->Xhist[mfqP->minindex],mfqP->localxmin,INSERT_VALUES,SCATTER_FORWARD);
650: VecScatterEnd(mfqP->scatterx,mfqP->Xhist[mfqP->minindex],mfqP->localxmin,INSERT_VALUES,SCATTER_FORWARD);
651: VecGetArray(mfqP->localxmin,&xmint);
652: for (i=0;i<mfqP->n;i++) mfqP->xmin[i] = xmint[i];
653: VecRestoreArray(mfqP->localxmin,&mfqP->xmin);
657: VecScatterBegin(mfqP->scatterf,mfqP->Fhist[mfqP->minindex],mfqP->localfmin,INSERT_VALUES,SCATTER_FORWARD);
658: VecScatterEnd(mfqP->scatterf,mfqP->Fhist[mfqP->minindex],mfqP->localfmin,INSERT_VALUES,SCATTER_FORWARD);
659: VecGetArray(mfqP->localfmin,&fmin);
660: for (i=0;i<mfqP->n+1;i++) {
661: if (i == mfqP->minindex) continue;
662:
663: mfqP->model_indices[ii++] = i;
664: VecScatterBegin(mfqP->scatterx,mfqP->Xhist[ii],mfqP->localx,
665: INSERT_VALUES, SCATTER_FORWARD);
666: VecScatterEnd(mfqP->scatterx,mfqP->Xhist[ii],mfqP->localx,
667: INSERT_VALUES, SCATTER_FORWARD);
668: VecGetArray(mfqP->localx,&x);
669: for (j=0;j<mfqP->n;j++) {
670: mfqP->Disp[i+mfqP->npmax*j] = (x[j] - mfqP->xmin[j])/mfqP->delta;
671: }
672: VecRestoreArray(mfqP->localx,&x);
674:
675: VecScatterBegin(mfqP->scatterf,mfqP->Fhist[ii],mfqP->localf,
676: INSERT_VALUES, SCATTER_FORWARD);
677: VecScatterEnd(mfqP->scatterf,mfqP->Fhist[ii],mfqP->localf,
678: INSERT_VALUES, SCATTER_FORWARD);
679: VecGetArray(mfqP->localf,&f);
680: for (j=0;j<mfqP->m;j++) {
681: mfqP->Fdiff[i*mfqP->n+j] = f[j] - fmin[j];
683: }
684: VecRestoreArray(mfqP->localf,&f);
685: }
686: for (j=0;j<mfqP->m;j++) {
687: mfqP->C[j] = fmin[j];
688: }
689:
690: VecRestoreArray(mfqP->localfmin,&fmin);
691:
692: }
694:
695: /* Determine the initial quadratic models */
696:
697: /* G = D(ModelIn,:) \ (F(ModelIn,1:m)-repmat(F(xkin,1:m),n,1)); */
698: /* D (nxn) Fdiff (nxm) => G (nxm) */
699: blasn2 = blasn;
700: LAPACKgesv_(&blasn,&blasm,mfqP->Disp,&blasnpmax,mfqP->iwork,mfqP->Fdiff,&blasn2,&info);
701: PetscInfo1(tao,"gesv returned %D\n",info);
703: cres = minnorm;
704: /* Gres = G*F(xkin,1:m)' G (nxm) Fk (m) */
705: BLASgemv_("N",&blasn,&blasm,&one,mfqP->Fdiff,&blasn,mfqP->C,&ione,&zero,mfqP->Gres,&ione);
707: /* Hres = G*G' */
708: BLASgemm_("N","T",&blasn,&blasn,&blasm,&one,mfqP->Fdiff, &blasn,mfqP->Fdiff,&blasn,&zero,mfqP->Hres,&blasn);
710: valid = PETSC_TRUE;
712: VecSetValues(tao->gradient,mfqP->n,mfqP->indices,mfqP->Gres,INSERT_VALUES);
713: VecAssemblyBegin(tao->gradient);
714: VecAssemblyEnd(tao->gradient);
715: VecNorm(tao->gradient,NORM_2,&gnorm);
716: gnorm *= mfqP->delta;
717: VecCopy(mfqP->Xhist[mfqP->minindex],tao->solution);
718: TaoMonitor(tao, iter, minnorm, gnorm, 0.0, step, &reason);
719: mfqP->nHist = mfqP->n+1;
720: mfqP->nmodelpoints = mfqP->n+1;
722: while (reason == TAO_CONTINUE_ITERATING) {
724: iter++;
725: PetscReal gnm = 10;
726: /* Solve the subproblem min{Q(s): ||s|| <= delta} */
727: gqtwrap(tao,&gnm,&mdec);
728: /* Evaluate the function at the new point */
730: for (i=0;i<mfqP->n;i++) {
731: mfqP->work[i] = mfqP->Xsubproblem[i]*mfqP->delta + mfqP->xmin[i];
732: }
733: VecDuplicate(tao->solution,&mfqP->Xhist[mfqP->nHist]);
734: VecDuplicate(tao->sep_objective,&mfqP->Fhist[mfqP->nHist]);
735: VecSetValues(mfqP->Xhist[mfqP->nHist],mfqP->n,mfqP->indices,mfqP->work,INSERT_VALUES);
736: VecAssemblyBegin(mfqP->Xhist[mfqP->nHist]);
737: VecAssemblyEnd(mfqP->Xhist[mfqP->nHist]);
739: TaoComputeSeparableObjective(tao,mfqP->Xhist[mfqP->nHist],mfqP->Fhist[mfqP->nHist]);
740: VecNorm(mfqP->Fhist[mfqP->nHist],NORM_2,&mfqP->Fres[mfqP->nHist]);
741: mfqP->Fres[mfqP->nHist]*=mfqP->Fres[mfqP->nHist];
742: rho = (mfqP->Fres[mfqP->minindex] - mfqP->Fres[mfqP->nHist]) / mdec;
743: mfqP->nHist++;
745: /* Update the center */
746: if ((rho >= mfqP->eta1) || (rho > mfqP->eta0 && valid==PETSC_TRUE)) {
747: /* Update model to reflect new base point */
748: for (i=0;i<mfqP->n;i++) {
749: mfqP->work[i] = (mfqP->work[i] - mfqP->xmin[i])/mfqP->delta;
750: }
751: for (j=0;j<mfqP->m;j++) {
752: /* C(j) = C(j) + work*G(:,j) + .5*work*H(:,:,j)*work';
753: G(:,j) = G(:,j) + H(:,:,j)*work' */
754: for (k=0;k<mfqP->n;k++) {
755: mfqP->work2[k]=0.0;
756: for (l=0;l<mfqP->n;l++) {
757: mfqP->work2[k]+=mfqP->H[j + mfqP->m*(k + l*mfqP->n)]*mfqP->work[l];
758: }
759: }
760: for (i=0;i<mfqP->n;i++) {
761: mfqP->C[j]+=mfqP->work[i]*(mfqP->Fdiff[i + mfqP->n* j] + 0.5*mfqP->work2[i]);
762: mfqP->Fdiff[i+mfqP->n*j] +=mfqP-> work2[i];
763: }
764: }
765: /* Cres += work*Gres + .5*work*Hres*work';
766: Gres += Hres*work'; */
768: BLASgemv_("N",&blasn,&blasn,&one,mfqP->Hres,&blasn,mfqP->work,&ione,
769: &zero,mfqP->work2,&ione);
770: for (i=0;j<mfqP->n;j++) {
771: cres += mfqP->work[i]*(mfqP->Gres[i] + 0.5*mfqP->work2[i]);
772: mfqP->Gres[i] += mfqP->work2[i];
773: }
774: mfqP->minindex = mfqP->nHist-1;
775: minnorm = mfqP->Fres[mfqP->minindex];
776: /* Change current center */
777: VecGetArray(mfqP->Xhist[mfqP->minindex],&xmint);
778: for (i=0;i<mfqP->n;i++) {
779: mfqP->xmin[i] = xmint[i];
780: }
781: VecRestoreArray(mfqP->Xhist[mfqP->minindex],&xmint);
784: }
786: /* Evaluate at a model-improving point if necessary */
787: if (valid == PETSC_FALSE) {
788: mfqP->q_is_I = 1;
789: affpoints(mfqP,mfqP->xmin,mfqP->c1);
790: if (mfqP->nmodelpoints < mfqP->n) {
791: PetscInfo(tao,"Model not valid -- model-improving");
792: modelimprove(tao,mfqP,1);
793: }
794: }
795:
797: /* Update the trust region radius */
798: deltaold = mfqP->delta;
799: normxsp = 0;
800: for (i=0;i<mfqP->n;i++) {
801: normxsp += mfqP->Xsubproblem[i]*mfqP->Xsubproblem[i];
802: }
803: normxsp = PetscSqrtReal(normxsp);
804: if (rho >= mfqP->eta1 && normxsp > 0.5*mfqP->delta) {
805: mfqP->delta = PetscMin(mfqP->delta*mfqP->gamma1,mfqP->deltamax);
806: } else if (valid == PETSC_TRUE) {
807: mfqP->delta = PetscMax(mfqP->delta*mfqP->gamma0,mfqP->deltamin);
808: }
810: /* Compute the next interpolation set */
811: mfqP->q_is_I = 1;
812: mfqP->nmodelpoints=0;
813: affpoints(mfqP,mfqP->xmin,mfqP->c1);
814: if (mfqP->nmodelpoints == mfqP->n) {
815: valid = PETSC_TRUE;
816: } else {
817: valid = PETSC_FALSE;
818: affpoints(mfqP,mfqP->xmin,mfqP->c2);
819: if (mfqP->n > mfqP->nmodelpoints) {
820: PetscInfo(tao,"Model not valid -- adding geometry points");
821: modelimprove(tao,mfqP,mfqP->n - mfqP->nmodelpoints);
822: }
823: }
824: for (i=mfqP->nmodelpoints;i>0;i--) {
825: mfqP->model_indices[i] = mfqP->model_indices[i-1];
826: }
827: mfqP->nmodelpoints++;
828: mfqP->model_indices[0] = mfqP->minindex;
829: morepoints(mfqP);
830: if (mfqP->nmodelpoints - mfqP->n - 1 == 0) {
831: reason = TAO_DIVERGED_USER;
832: tao->reason = TAO_DIVERGED_USER;
833: continue;
834: }
835: for (i=0;i<mfqP->nmodelpoints;i++) {
836: VecGetArray(mfqP->Xhist[mfqP->model_indices[i]],&x);
837: for (j=0;j<mfqP->n;j++) {
838: mfqP->Disp[i + mfqP->npmax*j] = (x[j] - mfqP->xmin[j]) / deltaold;
839: }
840: VecRestoreArray(mfqP->Xhist[mfqP->model_indices[i]],&x);
841: VecGetArray(mfqP->Fhist[mfqP->model_indices[i]],&f);
842: for (j=0;j<mfqP->m;j++) {
843: for (k=0;k<mfqP->n;k++) {
844: mfqP->work[k]=0.0;
845: for (l=0;l<mfqP->n;l++) {
846: mfqP->work[k] += mfqP->H[j + mfqP->m*(k + mfqP->n*l)] * mfqP->Disp[i + mfqP->npmax*l];
847: }
848: }
849: mfqP->RES[j*mfqP->npmax + i] = -mfqP->C[j] - BLASdot_(&blasn,&mfqP->Fdiff[j*mfqP->n],&ione,&mfqP->Disp[i],&blasnpmax) - 0.5*BLASdot_(&blasn,mfqP->work,&ione,&mfqP->Disp[i],&blasnpmax) + f[j];
850: }
851: VecRestoreArray(mfqP->Fhist[mfqP->model_indices[i]],&f);
852: }
855: /* Update the quadratic model */
856: getquadpounders(mfqP);
857: VecGetArray(mfqP->Fhist[mfqP->minindex],&fmin);
858: BLAScopy_(&blasm,fmin,&ione,mfqP->C,&ione);
859: /* G = G*(delta/deltaold) + Gdel */
860: ratio = mfqP->delta/deltaold;
861: iblas = blasm*blasn;
862: BLASscal_(&iblas,&ratio,mfqP->Fdiff,&ione);
863: BLASaxpy_(&iblas,&one,mfqP->Gdel,&ione,mfqP->Fdiff,&ione);
864: /* H = H*(delta/deltaold) + Hdel */
865: iblas = blasm*blasn*blasn;
866: ratio *= ratio;
867: BLASscal_(&iblas,&ratio,mfqP->H,&ione);
868: BLASaxpy_(&iblas,&one,mfqP->Hdel,&ione,mfqP->H,&ione);
871: /* Get residuals */
872: cres = mfqP->Fres[mfqP->minindex];
873: /* Gres = G*F(xkin,1:m)' */
874: BLASgemv_("N",&blasn,&blasm,&one,mfqP->Fdiff,&blasn,mfqP->C,&ione,&zero,mfqP->Gres,&ione);
875: /* Hres = sum i=1..m {F(xkin,i)*H(:,:,i)} + G*G' */
876: BLASgemm_("N","T",&blasn,&blasn,&blasm,&one,mfqP->Fdiff,&blasn,mfqP->Fdiff,&blasn,
877: &zero,mfqP->Hres,&blasn);
879: iblas = mfqP->n*mfqP->n;
881: for (j=0;j<mfqP->m;j++) {
882: BLASaxpy_(&iblas,&fmin[j],&mfqP->H[j],&blasm,mfqP->Hres,&ione);
883: }
884:
885: /* Export solution and gradient residual to TAO */
886: VecCopy(mfqP->Xhist[mfqP->minindex],tao->solution);
887: VecSetValues(tao->gradient,mfqP->n,mfqP->indices,mfqP->Gres,INSERT_VALUES);
888: VecAssemblyBegin(tao->gradient);
889: VecAssemblyEnd(tao->gradient);
890: VecNorm(tao->gradient,NORM_2,&gnorm);
891: gnorm *= mfqP->delta;
892: /* final criticality test */
893: TaoMonitor(tao, iter, minnorm, gnorm, 0.0, step, &reason);
894: }
895: return(0);
896: }
900: static PetscErrorCode TaoSetUp_POUNDERS(TaoSolver tao)
901: {
902: TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data;
903: PetscInt i;
904: IS isfloc,isfglob,isxloc,isxglob;
909: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
910: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
911: VecGetSize(tao->solution,&mfqP->n);
912: VecGetSize(tao->sep_objective,&mfqP->m);
913: mfqP->c1 = PetscSqrtReal((PetscReal)mfqP->n);
914: if (mfqP->npmax == PETSC_DEFAULT) {
915: mfqP->npmax = 2*mfqP->n + 1;
916: }
917: mfqP->npmax = PetscMin((mfqP->n+1)*(mfqP->n+2)/2,mfqP->npmax);
918: mfqP->npmax = PetscMax(mfqP->npmax, mfqP->n+2);
920: PetscMalloc((tao->max_funcs+10)*sizeof(Vec),&mfqP->Xhist);
921: PetscMalloc((tao->max_funcs+10)*sizeof(Vec),&mfqP->Fhist);
922: for (i=0;i<mfqP->n +1;i++) {
923: VecDuplicate(tao->solution,&mfqP->Xhist[i]);
924: VecDuplicate(tao->sep_objective,&mfqP->Fhist[i]);
925: }
926: VecDuplicate(tao->solution,&mfqP->workxvec);
927: mfqP->nHist = 0;
929: PetscMalloc((tao->max_funcs+10)*sizeof(PetscReal),&mfqP->Fres);
930: PetscMalloc(mfqP->npmax*mfqP->m*sizeof(PetscReal),&mfqP->RES);
931: PetscMalloc(mfqP->n*sizeof(PetscReal),&mfqP->work);
932: PetscMalloc(mfqP->n*sizeof(PetscReal),&mfqP->work2);
933: PetscMalloc(mfqP->n*sizeof(PetscReal),&mfqP->work3);
934: PetscMalloc(PetscMax(mfqP->m,mfqP->n+1)*sizeof(PetscReal),&mfqP->mwork);
935: PetscMalloc((mfqP->npmax - mfqP->n - 1)*sizeof(PetscReal),&mfqP->omega);
936: PetscMalloc((mfqP->n * (mfqP->n+1) / 2)*sizeof(PetscReal),&mfqP->beta);
937: PetscMalloc((mfqP->n + 1) *sizeof(PetscReal),&mfqP->alpha);
939: PetscMalloc(mfqP->n*mfqP->n*mfqP->m*sizeof(PetscReal),&mfqP->H);
940: PetscMalloc(mfqP->npmax*mfqP->npmax*sizeof(PetscReal),&mfqP->Q);
941: PetscMalloc(mfqP->npmax*mfqP->npmax*sizeof(PetscReal),&mfqP->Q_tmp);
942: PetscMalloc(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax)*sizeof(PetscReal),&mfqP->L);
943: PetscMalloc(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax)*sizeof(PetscReal),&mfqP->L_tmp);
944: PetscMalloc(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax)*sizeof(PetscReal),&mfqP->L_save);
945: PetscMalloc(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax)*sizeof(PetscReal),&mfqP->N);
946: PetscMalloc(mfqP->npmax * (mfqP->n+1) * sizeof(PetscReal),&mfqP->M);
947: PetscMalloc(mfqP->npmax * (mfqP->npmax - mfqP->n - 1) * sizeof(PetscReal), &mfqP->Z);
948: PetscMalloc(mfqP->npmax*sizeof(PetscReal),&mfqP->tau);
949: PetscMalloc(mfqP->npmax*sizeof(PetscReal),&mfqP->tau_tmp);
950: PetscMalloc(5*mfqP->npmax*sizeof(PetscReal),&mfqP->npmaxwork);
951: PetscMalloc(5*mfqP->npmax*sizeof(PetscBLASInt),&mfqP->npmaxiwork);
952: PetscMalloc(mfqP->n*sizeof(PetscReal),&mfqP->xmin);
953: PetscMalloc(mfqP->m*sizeof(PetscReal),&mfqP->C);
954: PetscMalloc(mfqP->m*mfqP->n*sizeof(PetscReal),&mfqP->Fdiff);
955: PetscMalloc(mfqP->npmax*mfqP->n*sizeof(PetscReal),&mfqP->Disp);
956: PetscMalloc(mfqP->n*sizeof(PetscReal),&mfqP->Gres);
957: PetscMalloc(mfqP->n*mfqP->n*sizeof(PetscReal),&mfqP->Hres);
958: PetscMalloc(mfqP->n*mfqP->n*sizeof(PetscReal),&mfqP->Gpoints);
959: PetscMalloc(mfqP->npmax*sizeof(PetscInt),&mfqP->model_indices);
960: PetscMalloc(mfqP->n*sizeof(PetscReal),&mfqP->Xsubproblem);
961: PetscMalloc(mfqP->m*mfqP->n*sizeof(PetscReal),&mfqP->Gdel);
962: PetscMalloc(mfqP->n*mfqP->n*mfqP->m*sizeof(PetscReal), &mfqP->Hdel);
963: PetscMalloc(PetscMax(mfqP->m,mfqP->n)*sizeof(PetscInt),&mfqP->indices);
964: PetscMalloc(mfqP->n*sizeof(PetscBLASInt),&mfqP->iwork);
965: for (i=0;i<PetscMax(mfqP->m,mfqP->n);i++) {
966: mfqP->indices[i] = i;
967: }
968: MPI_Comm_size(((PetscObject)tao)->comm,&mfqP->mpisize);
969: if (mfqP->mpisize > 1) {
970: VecCreateSeq(PETSC_COMM_SELF,mfqP->n,&mfqP->localx);
971: VecCreateSeq(PETSC_COMM_SELF,mfqP->n,&mfqP->localxmin);
972: VecCreateSeq(PETSC_COMM_SELF,mfqP->m,&mfqP->localf);
973: VecCreateSeq(PETSC_COMM_SELF,mfqP->m,&mfqP->localfmin);
974: ISCreateStride(MPI_COMM_SELF,mfqP->n,0,1,&isxloc);
975: ISCreateStride(MPI_COMM_SELF,mfqP->n,0,1,&isxglob);
976: ISCreateStride(MPI_COMM_SELF,mfqP->m,0,1,&isfloc);
977: ISCreateStride(MPI_COMM_SELF,mfqP->m,0,1,&isfglob);
978:
980: VecScatterCreate(tao->solution,isxglob,mfqP->localx,isxloc,&mfqP->scatterx);
981: VecScatterCreate(tao->sep_objective,isfglob,mfqP->localf,isfloc,&mfqP->scatterf);
983: ISDestroy(&isxloc);
984: ISDestroy(&isxglob);
985: ISDestroy(&isfloc);
986: ISDestroy(&isfglob);
988: }
990: if (!mfqP->usegqt) {
991: KSP ksp;
992: PC pc;
993: TaoCreate(PETSC_COMM_SELF,&mfqP->subtao);
994: TaoSetType(mfqP->subtao,"tao_bqpip");
995: TaoSetOptionsPrefix(mfqP->subtao,"pounders_subsolver_");
996: TaoSetObjectiveAndGradientRoutine(mfqP->subtao,pounders_fg,(void*)mfqP);
997: TaoSetMaximumIterations(mfqP->subtao,mfqP->gqt_maxits);
998: TaoGetKSP(mfqP->subtao,&ksp);
999: KSPGetPC(ksp,&pc);
1000: PCSetType(pc,PCNONE);
1001: TaoSetFromOptions(mfqP->subtao);
1002: }
1003: return(0);
1004: }
1008: static PetscErrorCode TaoDestroy_POUNDERS(TaoSolver tao)
1009: {
1010: TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data;
1011: PetscInt i;
1013:
1016: if (!mfqP->usegqt) {
1017: TaoDestroy(&mfqP->subtao);
1018: }
1019: PetscFree(mfqP->Fres);
1020: PetscFree(mfqP->RES);
1021: PetscFree(mfqP->work);
1022: PetscFree(mfqP->work2);
1023: PetscFree(mfqP->work3);
1024: PetscFree(mfqP->mwork);
1025: PetscFree(mfqP->omega);
1026: PetscFree(mfqP->beta);
1027: PetscFree(mfqP->alpha);
1028: PetscFree(mfqP->H);
1029: PetscFree(mfqP->Q);
1030: PetscFree(mfqP->Q_tmp);
1031: PetscFree(mfqP->L);
1032: PetscFree(mfqP->L_tmp);
1033: PetscFree(mfqP->L_save);
1034: PetscFree(mfqP->N);
1035: PetscFree(mfqP->M);
1036: PetscFree(mfqP->Z);
1037: PetscFree(mfqP->tau);
1038: PetscFree(mfqP->tau_tmp);
1039: PetscFree(mfqP->npmaxwork);
1040: PetscFree(mfqP->npmaxiwork);
1041: PetscFree(mfqP->xmin);
1042: PetscFree(mfqP->C);
1043: PetscFree(mfqP->Fdiff);
1044: PetscFree(mfqP->Disp);
1045: PetscFree(mfqP->Gres);
1046: PetscFree(mfqP->Hres);
1047: PetscFree(mfqP->Gpoints);
1048: PetscFree(mfqP->model_indices);
1049: PetscFree(mfqP->Xsubproblem);
1050: PetscFree(mfqP->Gdel);
1051: PetscFree(mfqP->Hdel);
1052: PetscFree(mfqP->indices);
1053: PetscFree(mfqP->iwork);
1054:
1055: for (i=0;i<mfqP->nHist;i++) {
1056: VecDestroy(&mfqP->Xhist[i]);
1057: VecDestroy(&mfqP->Fhist[i]);
1058: }
1059: if (mfqP->workxvec) {
1060: VecDestroy(&mfqP->workxvec);
1061: }
1062: PetscFree(mfqP->Xhist);
1063: PetscFree(mfqP->Fhist);
1065: if (mfqP->mpisize > 1) {
1066: VecDestroy(&mfqP->localx);
1067: VecDestroy(&mfqP->localxmin);
1068: VecDestroy(&mfqP->localf);
1069: VecDestroy(&mfqP->localfmin);
1070: }
1074: if (tao->data) {
1075: PetscFree(tao->data);
1076: }
1077: tao->data = PETSC_NULL;
1078: return(0);
1079: }
1083: static PetscErrorCode TaoSetFromOptions_POUNDERS(TaoSolver tao)
1084: {
1085: TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data;
1089: PetscOptionsHead("POUNDERS method for least-squares optimization");
1090: PetscOptionsReal("-tao_pounders_delta","initial delta","",mfqP->delta,&mfqP->delta,0);
1091: mfqP->npmax = PETSC_DEFAULT;
1092: PetscOptionsInt("-tao_pounders_npmax","max number of points in model","",mfqP->npmax,&mfqP->npmax,0);
1093: mfqP->usegqt = PETSC_FALSE;
1094: PetscOptionsBool("-tao_pounders_gqt","use gqt algorithm for subproblem","",mfqP->usegqt,&mfqP->usegqt,0);
1095: PetscOptionsTail();
1096: return(0);
1097: }
1101: static PetscErrorCode TaoView_POUNDERS(TaoSolver tao, PetscViewer viewer)
1102: {
1103: TAO_POUNDERS *mfqP = (TAO_POUNDERS *)tao->data;
1104: PetscBool isascii;
1109: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1110: if (isascii) {
1111: PetscViewerASCIIPushTab(viewer);
1112: if (mfqP->usegqt) {
1113: PetscViewerASCIIPrintf(viewer, "Subproblem solver: gqt\n");
1114: } else {
1115: PetscViewerASCIIPrintf(viewer, "Subproblem solver: bqpip\n");
1116: }
1117: PetscViewerASCIIPopTab(viewer);
1118: } else {
1119: SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NLS",((PetscObject)viewer)->type_name);
1120: }
1121: return(0);
1122: }
1127: PetscErrorCode TaoCreate_POUNDERS(TaoSolver tao)
1128: {
1129: TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data;
1131:
1134: tao->ops->setup = TaoSetUp_POUNDERS;
1135: tao->ops->solve = TaoSolve_POUNDERS;
1136: tao->ops->view = TaoView_POUNDERS;
1137: tao->ops->setfromoptions = TaoSetFromOptions_POUNDERS;
1138: tao->ops->destroy = TaoDestroy_POUNDERS;
1141: PetscNewLog(tao, TAO_POUNDERS, &mfqP);
1142: tao->data = (void*)mfqP;
1143: tao->max_it = 2000;
1144: tao->max_funcs = 4000;
1145: tao->fatol = 1e-8;
1146: tao->frtol = 1e-8;
1147: mfqP->delta = 0.1;
1148: mfqP->deltamax=1e3;
1149: mfqP->deltamin=1e-6;
1150: mfqP->c2 = 100.0;
1151: mfqP->theta1=1e-5;
1152: mfqP->theta2=1e-4;
1153: mfqP->gamma0=0.5;
1154: mfqP->gamma1=2.0;
1155: mfqP->eta0 = 0.0;
1156: mfqP->eta1 = 0.1;
1157: mfqP->gqt_rtol = 0.001;
1158: mfqP->gqt_maxits = 50;
1159: mfqP->workxvec = 0;
1160: return(0);
1161: }