--- trunk/OOPSE/libmdtools/NLModel0.cpp	2004/01/28 22:44:44	996
+++ trunk/OOPSE/libmdtools/NLModel0.cpp	2004/01/30 21:47:22	1000
@@ -1,6 +1,5 @@
 #include "NLModel.hpp"
 
-
 /**
  * calculate gradient using backward finite difference
  * df(Xk)/dXi = (f(Xk) - f(Xk - h*ei)) /h   
@@ -34,7 +33,6 @@ vector<double> NLModel0::BackwardGrad(const vector<dou
 
       hi = copysign(h[i], tempX[i]);
 
-      //tempX[i] = x[i] + hi;
       tempX[i] -= hi;
     }
 
@@ -173,15 +171,111 @@ vector<double> NLModel0::CentralGrad(const vector<doub
  * calculate hessian using finite difference
  * d2f(Xk)/dxidxj = (df(Xk+h*ej)/dXi + df(Xk - h*ej)/dXi) /2h   
  * where  h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
+ *
  */
-SymMatrix NLModel0::FDHessian(vector<double>& h){
+SymMatrix NLModel0::FiniteHessian(vector<double>& x, double fx, vector<double>& h){
   SymMatrix H(ndim);
+  vector<double> tempX;
+  vector<double> fi(ndim);
+  vector<double> hi;
+  double fij;
+  
+  tempX = x;
 
+  //calculate f(X+ h*ei)
   for(int i = 0; i < ndim; i++){
 
+#ifndef IS_MPI    
+
+    hi[i] = copysign(h[i],tempX[i]);
+
+    tempX[i] += hi[i];
+
+    fi[i] = calcF(tempX);
+    
+    tempX[i] -= hi;
+
+#else
+
+    if(procMappingArray[i] == myRank){
+
+      hi[i] = copysign(h[i],tempX[i]);
+
+      tempX[i] += hi[i];
+    }
+    
+    fi[i] = calcF(tempX);
+
+    if(procMappingArray[i] == myRank){    
+      tempX[i] -= hi[i];
+    }
+    
+#endif 
+  }
+
+#ifndef IS_MPI
+
+  for(int i = 0; i < ndim; i++){
+
+    tempX[i] += hi[i] * 2;
+
+    fii = calcF(tempX);
+    
+    H(i,j) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
+
+    tempX[i] -= hi[i];
+    
     for(int j = i + 1; j < ndim; j++){
 
+      tempX[j] += hi[j];
+
+      fij = calcF(tempX);
+
+      H(i,j) = ((f - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
+
+      tempX[j] -= hi[j];
     }
+    
+      tempX[i] -= hi[i];
   }
-  
+
+#else
+  for(int i = 0; i < ndim; i++){
+
+    if(procMappingArray[i] == myRank){    
+      tempX[i] += hi[i] * 2;
+    }
+    
+    fii = calcF(tempX);
+    
+    H(i,j) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
+
+    if(procMappingArray[i] == myRank){    
+      tempX[i] -= hi[i];
+    }
+    
+    for(int j = i + 1; j < ndim; j++){
+
+      if(procMappingArray[j] == myRank){    
+            tempX[j] += hi[j];
+      }
+
+      fij = calcF(tempX);
+      H(i,j) = ((f - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
+
+      if(procMappingArray[j] == myRank){    
+            tempX[j] -= hi[j];
+      }
+      
+    }
+    
+    if(procMappingArray[i] == myRank){    
+      tempX[i] -= hi[i];
+    }
+
+  }
+
+#endif
+
+  return H;  
 }