ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NLModel0.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NLModel0.cpp (file contents):
Revision 996 by tim, Wed Jan 28 22:44:44 2004 UTC vs.
Revision 1023 by tim, Wed Feb 4 22:26:00 2004 UTC

# Line 1 | Line 1
1   #include "NLModel.hpp"
2 <
3 <
2 > #include <math.h>
3 > #include "Utility.hpp"
4   /**
5   * calculate gradient using backward finite difference
6   * df(Xk)/dXi = (f(Xk) - f(Xk - h*ei)) /h  
# Line 22 | Line 22 | vector<double> NLModel0::BackwardGrad(const vector<dou
22      //tempX[i] = x[i] + hi;
23      tempX[i] -= hi;
24  
25 <    fminus = (*objfunc)(tempX);
25 >    fminus = calcF(tempX);
26  
27      partialGrad[i] =  (fx - fminus) / hi;
28  
# Line 34 | Line 34 | vector<double> NLModel0::BackwardGrad(const vector<dou
34  
35        hi = copysign(h[i], tempX[i]);
36  
37 –      //tempX[i] = x[i] + hi;
37        tempX[i] -= hi;
38      }
39  
40 <    fminus = (*objfunc)(tempX);
40 >    fminus = calcF(tempX);
41  
42      if(procMappingArray[i] == myRank){
43        partialGrad[i] =  (fx - fminus) / hi;
# Line 74 | Line 73 | vector<double> NLModel0::ForwardGrad(const vector<doub
73      //tempX[i] = x[i] + hi;
74      tempX[i] += hi;
75  
76 <    fplus = (*objfunc)(tempX);
76 >    fplus = calcF(tempX);
77  
78      partialGrad[i] =  (fplus - fx) / hi;
79  
# Line 90 | Line 89 | vector<double> NLModel0::ForwardGrad(const vector<doub
89        tempX[i] += hi;
90      }
91  
92 <    fminus = (*objfunc)(tempX);
92 >    fplus = calcF(tempX);
93  
94      if(procMappingArray[i] == myRank){
95 <      partialGrad[i] =  (fx - fminus) / hi;
95 >      partialGrad[i] =  (fplus - fx) / hi;
96  
97       //restore tempX to its original value
98        tempX[i] -= hi;
# Line 126 | Line 125 | vector<double> NLModel0::CentralGrad(const vector<doub
125      //tempX[i] = x[i] + hi
126      tempX[i] += hi;
127      
128 <    fplus = (*objfunc)(tempX);
128 >    fplus = calcF(tempX);
129  
130      //tempX[i] = x[i] -hi
131      tempX[i] -= 2*hi;
132 <    fminus = (*objfunc)(tempX);
132 >    fminus = calcF(tempX);
133      
134      partialGrad[i] =  (fplus + fminus) / (2*hi);
135  
# Line 146 | Line 145 | vector<double> NLModel0::CentralGrad(const vector<doub
145        tempX[i] += hi;
146      }
147  
148 <    fplus = (*objfunc)(tempX);
148 >    fplus = calcF(tempX);
149  
150      if(procMappingArray[i] == myRank){
151        partialGrad[i] =  (fx - fminus) / hi;
# Line 155 | Line 154 | vector<double> NLModel0::CentralGrad(const vector<doub
154        tempX[i] -= 2*hi;
155      }
156  
157 <    fminus = (*objfunc)(tempX);
157 >    fminus = calcF(tempX);
158      
159      if(procMappingArray[i] == myRank){
160        partialGrad[i] =  (fx - fminus) / (2*hi);
# Line 173 | Line 172 | vector<double> NLModel0::CentralGrad(const vector<doub
172   * calculate hessian using finite difference
173   * d2f(Xk)/dxidxj = (df(Xk+h*ej)/dXi + df(Xk - h*ej)/dXi) /2h  
174   * where  h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
175 + *
176   */
177 < SymMatrix NLModel0::FDHessian(vector<double>& h){
177 > SymMatrix NLModel0::FiniteHessian(vector<double>& x, double fx, vector<double>& h){
178    SymMatrix H(ndim);
179 +  vector<double> tempX;
180 +  vector<double> fi(ndim);
181 +  vector<double> hi;
182 +  double fii, fij;
183 +  
184 +  tempX = x;
185  
186 +  //calculate f(X+ h*ei)
187    for(int i = 0; i < ndim; i++){
188  
189 + #ifndef IS_MPI    
190 +
191 +    hi[i] = copysign(h[i],tempX[i]);
192 +
193 +    tempX[i] += hi[i];
194 +
195 +    fi[i] = calcF(tempX);
196 +    
197 +    tempX[i] -=  hi[i];
198 +
199 + #else
200 +
201 +    if(procMappingArray[i] == myRank){
202 +
203 +      hi[i] = copysign(h[i],tempX[i]);
204 +
205 +      tempX[i] += hi[i];
206 +    }
207 +    
208 +    fi[i] = calcF(tempX);
209 +
210 +    if(procMappingArray[i] == myRank){    
211 +      tempX[i] -= hi[i];
212 +    }
213 +    
214 + #endif
215 +  }
216 +
217 + #ifndef IS_MPI
218 +
219 +  for(int i = 0; i < ndim; i++){
220 +
221 +    tempX[i] += hi[i] * 2;
222 +
223 +    fii = calcF(tempX);
224 +    
225 +    H(i, i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
226 +
227 +    tempX[i] -= hi[i];
228 +    
229      for(int j = i + 1; j < ndim; j++){
230  
231 +      tempX[j] += hi[j];
232 +
233 +      fij = calcF(tempX);
234 +
235 +      H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
236 +
237 +      tempX[j] -= hi[j];
238      }
239 +    
240 +      tempX[i] -= hi[i];
241    }
242 <  
242 >
243 > #else
244 >  for(int i = 0; i < ndim; i++){
245 >
246 >    if(procMappingArray[i] == myRank){    
247 >      tempX[i] += hi[i] * 2;
248 >    }
249 >    
250 >    fii = calcF(tempX);
251 >    
252 >    H(i,i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
253 >
254 >    if(procMappingArray[i] == myRank){    
255 >      tempX[i] -= hi[i];
256 >    }
257 >    
258 >    for(int j = i + 1; j < ndim; j++){
259 >
260 >      if(procMappingArray[j] == myRank){    
261 >            tempX[j] += hi[j];
262 >      }
263 >
264 >      fij = calcF(tempX);
265 >      H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
266 >
267 >      if(procMappingArray[j] == myRank){    
268 >            tempX[j] -= hi[j];
269 >      }
270 >      
271 >    }
272 >    
273 >    if(procMappingArray[i] == myRank){    
274 >      tempX[i] -= hi[i];
275 >    }
276 >
277 >  }
278 >
279 > #endif
280 >
281 >  return H;  
282   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines