Merged

sample/Hartmann_Matrix.txt

@@ -0,0 +1,7 @@
+6
+ 14.8253   -6.4243   7.8746  -1.2498  10.2733  10.2733 
+ -6.4243   15.1024  -1.1155  -0.2761  -8.2117  -8.2117 
+  7.8746   -1.1155  51.8519 -23.3482  12.5902  12.5902 
+ -1.2498   -0.2761 -23.3482  22.7962  -9.8958  -9.8958 
+ 10.2733   -8.2117  12.5902  -9.8958  21.0656  21.0656 
+ 10.2733   -8.2117  12.5902  -9.8958  21.0656  21.0656 

sample/testEskow.cpp

@@ -1,3 +1,4 @@
+#include <fstream>
 #include <cstring>
 #include <iostream>
 #include <iomanip>
@@ -26,25 +27,39 @@ struct MatrixOp
   }
 };
 
-int main()
+int main(int argc, char **argv)
 {
   MatrixOp M;
   double norm_E;
+  ifstream fi(argv[1]);
+  ofstream f("eskowed.txt");
+  CholeskyEskow<double,MatrixOp> chol;
 
-  M.N = 6;
+  fi >> M.N;
   M.M.resize(M.N*M.N);
 
-  memcpy(&M.M[0], &Hartmann_Matrix[0][0], sizeof(double)*36);
-
-  CholeskyEskow::cholesky_eskow(M, M.N, norm_E);
-
   for (int i = 0; i < M.N; i++)
     {
       for (int j = 0; j < M.N; j++)
 	{
-	  cout << setprecision(25) << M(i,j) << " ";
+	  fi >> M(i,j);
+	  if (j > i)
+	    M(i,j) =0;
 	}
-      cout << endl;
+    }
+
+  chol.cholesky(M, M.N, norm_E);
+
+  for (int i = 0; i < M.N; i++)
+    {
+      for (int j = 0; j < M.N; j++)
+	{
+	  if (j > i) 
+            f << "0 ";
+          else
+ 	    f << setprecision(25) << M(i,j) << " ";
+	}
+      f << endl;
     }
   return 0;
 }

src/eskow.hpp

@@ -7,24 +7,49 @@
 
 /* Implementation of Schnabel & Eskow, 1999, Vol. 9, No. 4, pp. 1135-148, SIAM J. OPTIM. */
 
-namespace CholeskyEskow
+template<typename T, typename A>
+class CholeskyEskow
+{
+private:
+  static const bool verbose_eskow = true;
+  T tau, tau_bar, mu;
+
+  void print_matrix(A& m, int N)
+  {
+    using std::cout;
+    using std::endl;
+    using std::setprecision;
+
+    if (verbose_eskow)
       {
 
-  template<typename T, typename A>
+	for (int i = 0; i < N; i++)
+	  {
+	    for (int j = 0; j < N; j++)
+	      {
+		cout.width(6);
+		cout << setprecision(5) << m(i,j) << " ";
+	      }
+	    cout << endl;
+	  }
+	cout << endl;
+      }
+  }
+
   T max_diag(A& m, int j, int N)
   {
-    T maxval = m(j,j);
+    T maxval = std::abs(m(j,j));
 
     for (int k = j+1; k < N; k++)
       {
-	maxval = std::max(maxval, m(k,k));
+	maxval = std::max(maxval, std::abs(m(k,k)));
       }
     return maxval;
   }
 
-  template<typename T, typename A>
   void minmax_diag(A& m, int j, int N, T& minval, T& maxval, int& i_min, int& i_max)
   {
+    i_min = i_max = j;
     minval = maxval = m(j,j);
 
     for (int k = j+1; k < N; k++)
@@ -35,34 +60,30 @@ namespace CholeskyEskow
   
     for (int k = j; k < N; k++)
       {
-	if (m(k,k) == minval)
+	if (m(k,k) == minval && i_min < 0)
 	  i_min = k;
-	if (m(k,k) == maxval)
+	if (m(k,k) == maxval && i_max < 0)
 	  i_max = k;
       }
   }
 
-  template<typename T, typename A>
   void swap_rows(A& m, int N, int i0, int i1)
   {
     for (int r = 0; r < N; r++)
       std::swap(m(r,i0), m(r,i1));
   }
 
-  template<typename T, typename A>
   void swap_cols(A& m, int N, int i0, int i1)
   {
     for (int c = 0; c < N; c++)
       std::swap(m(i0,c), m(i1,c));
   }
 
-  template<typename T>
   T square(T x)
   {
     return x*x;
   }
 
-  template<typename T, typename A>
   T min_row(A& m, int j, int N)
   {
     T a = 1/m(j,j);
@@ -70,13 +91,12 @@ namespace CholeskyEskow
     
     for (int i = j+2; i < N; i++)
       {
-	v = std::max(v, m(i, i) - square(m(i,j))*a);
+	v = std::min(v, m(i, i) - square(m(i,j))*a);
       }
 
     return v;
   }
 
-  template<typename T>
   int g_max(const std::vector<T>& g, int j, int N)
   {
     T a = g[j];
@@ -93,14 +113,18 @@ namespace CholeskyEskow
     return k;
   }
 
-  template<typename T, typename A>
+public:
-  void cholesky_eskow(A& m, int N, T& norm_E)
+  CholeskyEskow()
+  {
+    tau = std::pow(mach_epsilon<T>(), 1./3);
+    tau_bar = std::pow(mach_epsilon<T>(), 2./3);
+    mu=0.1;
+  }
+
+  void cholesky(A& m, int N, T& norm_E)
   {
-    T tau_bar = std::pow(mach_epsilon<T>(), 2./3);
-    T tau = std::pow(mach_epsilon<T>(), 1./3);
-    T mu = 0.1;
     bool phaseone = true;
-    T gamma = max_diag<T,A>(m, 0, N);
+    T gamma = max_diag(m, 0, N);
     int j;
 
     norm_E = 0;
@@ -110,7 +134,9 @@ namespace CholeskyEskow
 	T minval, maxval;
 	int i_min, i_max;
 	
-	minmax_diag<T,A>(m, j, N, minval, maxval, i_min, i_max);
+	print_matrix(m, N);
+     
+	minmax_diag(m, j, N, minval, maxval, i_min, i_max);
 	if (maxval < tau_bar*gamma || minval < -mu*maxval)
 	  {
 	    phaseone = false;
@@ -119,11 +145,13 @@ namespace CholeskyEskow
       
 	if (i_max != j)	  
 	  {
-	    swap_cols<T,A>(m, N, i_max, j);
+	    std::cout << "Have to swap i=" << i_max << " and j=" << j << std::endl;
-	    swap_rows<T,A>(m, N, i_max, j);
+	    swap_cols(m, N, i_max, j);
+	    swap_rows(m, N, i_max, j);
+
 	  }
       
-	if (min_row<T,A>(m, j, N) < -mu*gamma)
+	if (min_row(m, j, N) < -mu*gamma)
 	  {
 	    phaseone = false;
 	    break;
@@ -135,7 +163,7 @@ namespace CholeskyEskow
 	for (int i = j+1; i < N; i++)
 	  {
 	    m(i,j) /= L_jj;
-	    for (int k = j+1; k < i; k++)
+	    for (int k = j+1; k <= i; k++)
 	      m(i,k) -= m(i,j)*m(k,j);
 	  }
       }
@@ -151,9 +179,10 @@ namespace CholeskyEskow
 
     
 
-
     if (!phaseone && j < (N-1))
       {
+	std::cout << "Phase two ! (j=" << j << ")" << std::endl;
+
 	int k = j-1;
 	std::vector<T> g(N);
 
@@ -173,12 +202,15 @@ namespace CholeskyEskow
 	    int i = g_max(g, j, N);
 	    T norm_j;
 
+	    print_matrix(m, N);
+
 	    if (i != j)
 	      {
-		swap_cols<T,A>(m, N, i, j);
+		swap_cols(m, N, i, j);
-		swap_rows<T,A>(m, N, i, j);		
+		swap_rows(m, N, i, j);		
 	      }
 
+
 	    for (int i = j+1; j < N; j++)
 	      {
 		norm_j += std::abs(m(i,j));
@@ -208,7 +240,7 @@ namespace CholeskyEskow
 	    for (int i = j+1; i < N; i++)
 	      {
 		m(i,j) /= L_jj;
-		for (int k = j+1; k < i; k++)
+		for (int k = j+1; k <= i; k++)
 		  m(i,k) -= m(i,j)*m(k,j);
 	      }
 	  }

src/h5_readFlash.cpp

@@ -179,13 +179,13 @@ void h5_read_flash3_particles (H5File* file,
                         int* totalparticles,
                         int* localnp,
                         int* particle_offset,
-                        float pos1[],
+                        float *pos1,
-                        float pos2[],
+                        float *pos2,
-                        float pos3[],
+                        float *pos3,
-                        float vel1[],
+                        float *vel1,
-                        float vel2[],
+                        float *vel2,
-                        float vel3[],
+                        float *vel3,
-                        int    id[])
+                        int    *id)
 {
 
   herr_t   status;
@@ -339,15 +339,27 @@ void h5_read_flash3_particles (H5File* file,
 
     /* convert buffer into particle struct */
 
+    if (id) {
     for(p=0; p < (pcount); p++) {
       id[p+poffset] = (int)  *(partBuffer+iptag-1+p*numProps);
+    } }
+
+    if (pos1 && pos2 && pos3) {
+    for(p=0; p < (pcount); p++) {
       pos1[p+poffset]  = (float) *(partBuffer+ipx-1+p*numProps);
       pos2[p+poffset]  = (float) *(partBuffer+ipy-1+p*numProps);
       pos3[p+poffset]  = (float) *(partBuffer+ipz-1+p*numProps);
+    }
+    }
+   
+    
+    if (vel1 && vel2 && vel3) {
+    for(p=0; p < (pcount); p++) {
       vel1[p+poffset]  = (float) *(partBuffer+ipvx-1+p*numProps);
       vel2[p+poffset]  = (float) *(partBuffer+ipvy-1+p*numProps);
       vel3[p+poffset]  = (float) *(partBuffer+ipvz-1+p*numProps);
     }
+    }
 
     memspace.close();
     //status = H5Sclose(memspace);

src/h5_readFlash.hpp

@@ -23,13 +23,13 @@ void h5_read_flash3_particles (H5File* file,
                         int* totalparticles,
                         int* localnp,
                         int* particle_offset,
-                        float pos1[],
+                        float *pos1,
-                        float pos2[],
+                        float *pos2,
-                        float pos3[],
+                        float *pos3,
-                        float vel1[],
+                        float *vel1,
-                        float vel2[],
+                        float *vel2,
-                        float vel3[],
+                        float *vel3,
-                        int    id[]);
+                        int    *id);
 
 void h5_read_flash3_header_info(H5File* file,
 				double* time,                   /* simulation time */

src/loadFlash.cpp

@@ -51,29 +51,36 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags)
     data->Omega_Lambda = omegalambda;
 
     // particle data
+    if (loadflags& NEED_POSITION) {
       for (int i = 0; i < 3; i++) {
         data->Pos[i] = new float[data->NumPart];
         if (data->Pos[i] == 0) {
 	  delete data;
 	return 0;
       }
-    }
+    } }
 
+
+    if (loadflags &NEED_VELOCITY) {
     for (int i = 0; i < 3; i++) {
       data->Vel[i] = new float[data->NumPart];
       if (data->Vel[i] == 0) {
 	delete data;
 	return 0;
       }
-    }
+    } }
  
+    if (loadflags & NEED_GADGET_ID) {
     data->Id = new int[data->NumPart];
     if (data->Id == 0) {
       delete data;
       return 0;
     }
+    }
 
     int offset = 0;
+
+    if (loadflags &(NEED_GADGET_ID|NEED_POSITION|NEED_VELOCITY))
     h5_read_flash3_particles(&file, &npart, &npart, &offset, 
 			     data->Pos[0], data->Pos[1], data->Pos[2], 
 			     data->Vel[0], data->Vel[1], data->Vel[2],
@@ -81,8 +88,8 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags)
 
     for (int i = 0; i < 3; i++) {
       for (int n = 0; n < data->NumPart; n++) {
-	data->Pos[i][n] = data->Pos[i][n] * data->Hubble / kpc2cm;
+	if (loadflags& NEED_POSITION) data->Pos[i][n] = data->Pos[i][n] * data->Hubble / kpc2cm;
-	data->Vel[i][n] = data->Vel[i][n] * data->Hubble / km2cm;
+	if (loadflags&NEED_VELOCITY) data->Vel[i][n] = data->Vel[i][n] * data->Hubble / km2cm;
       }
     }
 

src/mach.hpp

@@ -1,5 +1,7 @@
-#ifndef __MACHINE_INFORMATION_HPP
+#ifndef __COSMO_MACHINE_TEST_HPP
-#define __MACHINE_INFORMATION_HPP
+#define __COSMO_MACHINE_TEST_HPP
+
+#include <iostream>
 
 template<typename T>
 T mach_epsilon()