Actual source code: lmvmmat.c

  1: #include "lmvmmat.h"
  2: #include "tao_general.h"
 3:  #include taovec.h
  4: #include "stdio.h"

  8: TaoLMVMMat::TaoLMVMMat( TaoVec *tv, int nlm){

 10:   TaoFunctionBegin;

 12:   lm = nlm;
 13:   lmnow = 0;
 14:   eps = 2.2e-16;
 15:   eps = 2.2e-11;
 16:   iter = 0;
 17:   rejects = 0;

 19:   tv->CloneVecs(nlm+1,&S);
 20:   tv->CloneVecs(nlm+1,&Y);

 22:   rho = new double[nlm+1];
 23:   beta = new double[nlm+1];

 25:   Gprev = 0;
 26:   Xprev = 0;
 27:   H0 = 0;
 28:   H0Vec = 0;

 30:   H0default=TAO_TRUE;

 32:   return;
 33: }

 37: TaoLMVMMat::~TaoLMVMMat(){
 38:   int info;

 40:   info = S[0]->DestroyVecs(lm+1,S);
 41:   info = Y[0]->DestroyVecs(lm+1,Y);
 42:   /*  if (H0) delete H0; */
 43:   if (H0Vec) delete H0Vec;

 45:   delete[] rho;
 46:   delete[] beta;
 47:   if (info) rho=0;
 48: }


 53: int TaoLMVMMat::InitialApproximation(TaoVec *HH0){
 54:   TaoFunctionBegin;
 55:   TaoFunctionReturn(0);
 56: }

 60: int TaoLMVMMat::SetH0(TaoMat* HH0){
 61:   int info;
 62:   TaoFunctionBegin;
 63:   if (this->H0){
 64:     /*    delete this->H0; Application deletes it */
 65:     this->H0=0;
 66:   };
 67:   if (this->H0Vec){ delete this->H0Vec;this->H0Vec=0;};
 68:   this->H0=HH0;
 69:   this->H0default=TAO_TRUE;
 70:   if (HH0){
 71:     info=this->S[0]->Clone(&this->H0Vec);CHKERRQ(info);
 72:     this->H0default=TAO_FALSE;
 73:   }
 74:   TaoFunctionReturn(0);
 75: }

 79: int TaoLMVMMat::reset(){
 80:   int i;
 81:   TaoFunctionBegin;
 82:   Gprev = Y[lm];
 83:   Xprev = S[lm];
 84:   for (i=0;i<lm;i++){
 85:     rho[i]=0.0;
 86:   }
 87:   rho[0]=1.0;
 88:   y0normsquared = 1.0;
 89:   iter=0;
 90:   rejects=0;
 91:   lmnow = 0;
 92:   TaoFunctionReturn(0);
 93: }


 98: int TaoLMVMMat::update(TaoVec* x, TaoVec* g){
 99:   int i,info;
100:   double rhotemp,rhotol;
101:   double y0temp;
102:   TaoFunctionBegin;

104:   if (this->iter==0){

106:     info=this->reset(); CHKERRQ(info);

108:   } else {
109: 
110:     info = Gprev->Aypx(-1.0,g); CHKERRQ(info);
111:     info = Xprev->Aypx(-1.0,x); CHKERRQ(info);
112:     info = Gprev->Dot(Xprev,&rhotemp); CHKERRQ(info);
113:     info = Gprev->Norm2squared(&y0temp); CHKERRQ(info);
114:     rhotol=eps*y0temp;
115:     if (rhotemp > rhotol){
116:       lmnow = TaoMin(lmnow+1,lm);
117:       for (i=lm-1;i>=0;i--){
118:         S[i+1]=S[i];
119:         Y[i+1]=Y[i];
120:         rho[i+1]=rho[i];
121:       }
122:       S[0]=Xprev;
123:       Y[0]=Gprev;
124:       rho[0]=1.0/rhotemp;
125:       y0normsquared=y0temp;
126:       Xprev=S[lm]; Gprev=Y[lm];

128:     } else {
129:       info = PetscLogInfo((0,"LMVM update rejected: iter= %d, %5.2e<=%5.2e \n",iter,rhotemp,rhotol)); CHKERRQ(info);
130:       rejects++;
131:     }
132: 
133:   }
134: 
135:   iter++;
136:   info = Xprev->CopyFrom(x); CHKERRQ(info);
137:   info = Gprev->CopyFrom(g); CHKERRQ(info);

139:   TaoFunctionReturn(0);
140: }

144: int TaoLMVMMat::Solve(TaoVec* tb, TaoVec* dx, TaoTruth *tt){

146:   int      ll,info;
147:   double   sq, yq,dd;

149:   TaoFunctionBegin;
150:   if (lmnow<1){
151:     rho[0]=1.0; y0normsquared = 1.0;
152:   }

154:   info = dx->CopyFrom(tb); CHKERRQ(info);
155: 
156:   for (ll = 0; ll<lmnow; ll++){
157: 
158:     info = dx->Dot(S[ll],&sq); CHKERRQ(info);
159:     beta[ll] = sq*rho[ll];
160:     info = dx->Axpy(-beta[ll], Y[ll]); CHKERRQ(info);

162:   }
163: 
164:   if (this->H0default==TAO_TRUE){
165:     info = dx->Scale( 1.0 / (rho[0]*y0normsquared) ); CHKERRQ(info);
166:   } else {
167:     info=this->H0Vec->CopyFrom(dx);CHKERRQ(info);
168:     info=dx->Scale( 1.0 / (rho[0]*y0normsquared) ); CHKERRQ(info);
169:     info=this->H0->Solve(H0Vec,dx,tt);CHKERRQ(info);
170:     info=dx->Dot(H0Vec,&dd);CHKERRQ(info);
171:     if (dd<=0){
172:       info=this->H0Vec->CopyFrom(dx);CHKERRQ(info);
173:       info=dx->Scale( 1.0 / (rho[0]*y0normsquared) ); CHKERRQ(info);
174:       info = PetscLogInfo((0,"REJECTING Application Hessian Solve: Inner product= %4.2e \n",dd)); CHKERRQ(info);
175:     } else {
176:       info = PetscLogInfo((0,"Using Application Hessian Solve: Inner product= %4.2e \n",dd)); CHKERRQ(info);
177:     }
178:   }

180:   for (ll=lmnow-1; ll>=0; ll--){
181: 
182:     info = dx->Dot(Y[ll],&yq); CHKERRQ(info);
183:     info = dx->Axpy( beta[ll] - yq * rho[ll] , S[ll]); CHKERRQ(info);
184: 
185:   }

187:   *tt = TAO_TRUE;
188:   TaoFunctionReturn(0);
189: }


194: int TaoLMVMMat::GetX0(TaoVec* x){
195:   int i,info;

197:   TaoFunctionBegin;
198:   info=x->CopyFrom(Xprev);CHKERRQ(info);
199:   for (i=0;i<this->lmnow;i++){
200:     info = x->Axpy(-1.0, S[i]); CHKERRQ(info);
201:   }
202:   TaoFunctionReturn(0);
203: }

207: int TaoLMVMMat::Clone(TaoMat** tm){
208:   TaoFunctionBegin;
209:   SETERRQ(56,"Operation not defined");
210:   /* TaoFunctionReturn(1); */
211: }

215: int TaoLMVMMat::Presolve(){
216:   TaoFunctionBegin;
217:   TaoFunctionReturn(0);
218: }

222: int TaoLMVMMat::CopyFrom(TaoMat* tm){
223:   TaoFunctionBegin;
224:   SETERRQ(56,"Operation not defined");
225:   /* TaoFunctionReturn(1); */
226: }

230: int TaoLMVMMat::GetDimensions( int* m, int* n ){
231:   TaoFunctionBegin;
232:   SETERRQ(56,"Operation not defined");
233:   /* TaoFunctionReturn(1); */
234: }

238: int TaoLMVMMat::Multiply(TaoVec* tx,TaoVec* ty){
239:   TaoFunctionBegin;
240:   SETERRQ(56,"Operation not defined");
241:   /* TaoFunctionReturn(1); */
242: }

246: int TaoLMVMMat::MultiplyAdd(TaoVec* tx,TaoVec* ty,TaoVec* tw){
247:   TaoFunctionBegin;
248:   SETERRQ(56,"Operation not defined");
249:   /* TaoFunctionReturn(1); */
250: }

254: int TaoLMVMMat::MultiplyTranspose(TaoVec* tx,TaoVec* ty){
255:   TaoFunctionBegin;
256:   SETERRQ(56,"Operation not defined");
257:   /* TaoFunctionReturn(1); */
258: }

262: int TaoLMVMMat::MultiplyTransposeAdd(TaoVec* tx,TaoVec* ty,TaoVec* yw){
263:   TaoFunctionBegin;
264:   SETERRQ(56,"Operation not defined");
265:   /* TaoFunctionReturn(1); */
266: }


271: int TaoLMVMMat::SetDiagonal(TaoVec* tv){
272:   TaoFunctionBegin;
273:   SETERRQ(56,"Operation not defined");
274:   /* TaoFunctionReturn(1); */
275: }

279: int TaoLMVMMat::AddDiagonal(TaoVec* tv){
280:   TaoFunctionBegin;
281:   SETERRQ(56,"Operation not defined");
282:   /* TaoFunctionReturn(1); */
283: }

287: int TaoLMVMMat::View(){
288:   int i,info;
289:   TaoFunctionBegin;
290:   printf("Xprev \n");
291:   info = Xprev->View(); CHKERRQ(info);
292:   printf("Gprev \n");
293:   info = Gprev->View(); CHKERRQ(info);
294:   for (i=0;i<lmnow; i++){
295:     printf("rho: %4.2e \n",this->rho[i]);
296:     printf("S %d \n",i);
297:     info = this->S[i]->View(); CHKERRQ(info);
298:     printf("Y %d \n",i);
299:     info = this->Y[i]->View(); CHKERRQ(info);
300:   }
301:   TaoFunctionReturn(0);
302: }

306: int TaoLMVMMat::GetDiagonal(TaoVec* tv){
307:   TaoFunctionBegin;
308:   SETERRQ(56,"Operation not defined");
309:   /* TaoFunctionReturn(1); */
310: }

314: int TaoLMVMMat::ShiftDiagonal(double c){
315:   TaoFunctionBegin;
316:   SETERRQ(56,"Operation not defined");
317:   /* TaoFunctionReturn(1); */
318: }

322: int TaoLMVMMat::extractSubMatrix(TaoMat**M, TaoIndexSet* row,TaoIndexSet* col){
323:   TaoFunctionBegin;
324:   SETERRQ(56,"Operation not defined");
325:   /* TaoFunctionReturn(1); */
326: };

330: int TaoLMVMMat::RowScale(TaoVec* tv){
331:   TaoFunctionBegin;
332:   SETERRQ(56,"Operation not defined");
333:   /* TaoFunctionReturn(1); */
334: }

338: int TaoLMVMMat::ColScale(TaoVec* tv){
339:   TaoFunctionBegin;
340:   SETERRQ(56,"Operation not defined");
341:   /* TaoFunctionReturn(1); */
342: }

346: int TaoLMVMMat::D_Fischer(TaoVec *tv1, TaoVec *tv2, TaoVec *tv3, TaoVec *tv4,
347:                           TaoVec *tv5, TaoVec *tv6, TaoVec *tv7, TaoVec *tv8)
348: {
349:   TaoFunctionBegin;
350:   SETERRQ(56,"Operation not defined");
351:   /* TaoFunctionReturn(1); */
352: }

356: int TaoLMVMMat::D_SFischer(TaoVec *tv1, TaoVec *tv2, TaoVec *tv3, TaoVec *tv4,
357:                            double, TaoVec *tv5, TaoVec *tv6,
358:                            TaoVec *tv7, TaoVec *tv8, TaoVec *tv9)
359: {
360:   TaoFunctionBegin;
361:   SETERRQ(56,"Operation not defined");
362:   /* TaoFunctionReturn(1); */
363: }

367: int TaoLMVMMat::minQuadraticTrustRegion(TaoVec*,TaoVec*,double, TaoTruth*){
368:   TaoFunctionBegin;
369:   SETERRQ(56,"Operation not defined");
370:   /* TaoFunctionReturn(1); */
371: }