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: }