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