Fixed PETSc nnz problem.

AMDiS/libtool

@@ -44,7 +44,7 @@ available_tags=" CXX F77"
 
 # ### BEGIN LIBTOOL CONFIG
 
-# Libtool was configured on host deimos104:
+# Libtool was configured on host p1d098:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -82,13 +82,13 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O3"
 
 # A language-specific compiler.
-CC="gcc"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
 
 # Is the compiler the GNU C compiler?
 with_gcc=yes
@@ -174,7 +174,7 @@ dlopen_self=unknown
 dlopen_self_static=unknown
 
 # Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
 
 # Compiler flag to turn off builtin functions.
 no_builtin_flag=" -fno-builtin"
@@ -6763,7 +6763,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
 # End:
 # ### BEGIN LIBTOOL TAG CONFIG: CXX
 
-# Libtool was configured on host deimos104:
+# Libtool was configured on host p1d098:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -6801,13 +6801,13 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O3"
 
 # A language-specific compiler.
-CC="g++"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicxx"
 
 # Is the compiler the GNU C compiler?
 with_gcc=yes
@@ -6893,7 +6893,7 @@ dlopen_self=unknown
 dlopen_self_static=unknown
 
 # Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
 
 # Compiler flag to turn off builtin functions.
 no_builtin_flag=" -fno-builtin"
@@ -6960,11 +6960,11 @@ predeps=""
 
 # Dependencies to place after the objects being linked to create a
 # shared library.
-postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
+postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -libverbs -lrt -lnuma -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
 
 # The library search path used internally by the compiler when linking
 # a shared library.
-compiler_lib_search_path=`echo "-L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/fastfs/wir/local/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+compiler_lib_search_path=`echo "-L/usr/lib64 -L/licsoft/libraries/openmpi/1.2.6/64bit/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/fastfs/wir/local/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
 
 # Method to check whether dependent libraries are shared objects.
 deplibs_check_method="pass_all"
@@ -7071,7 +7071,7 @@ include_expsyms=""
 
 # ### BEGIN LIBTOOL TAG CONFIG: F77
 
-# Libtool was configured on host deimos104:
+# Libtool was configured on host p1d098:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -7109,7 +7109,7 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O3"

AMDiS/src/parallel/PetscSolver.cc

@@ -24,10 +24,12 @@
 
 namespace AMDiS {
 
+  using namespace std;
+
   PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
   {    
     if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
-      std::cout << "[0]  Petsc-Iteration " << iter << ": " << rnorm << std::endl;
+      cout << "[0]  Petsc-Iteration " << iter << ": " << rnorm << endl;
 
     return 0;
   }
@@ -66,12 +68,12 @@ namespace AMDiS {
     typedef traits::range_generator<row, Matrix>::type cursor_type;
     typedef traits::range_generator<nz, cursor_type>::type icursor_type;
 
-    std::vector<int> cols;
-    std::vector<double> values;
+    vector<int> cols;
+    vector<double> values;
     cols.reserve(300);
     values.reserve(300);
     
-    std::vector<int> globalCols;
+    vector<int> globalCols;
 
 
     // === Traverse all rows of the dof matrix and insert row wise the values ===
@@ -86,7 +88,9 @@ namespace AMDiS {
       bool periodicRow = meshDistributor->isPeriodicDof(globalRowDof);
 
       if (!periodicRow) {
-	// Calculate petsc row index.
+	// === Row DOF index is not periodic. ===
+
+	// Calculate PETSc row index.
 	int rowIndex = globalRowDof * dispMult + dispAddRow;
 
 	cols.clear();
@@ -99,35 +103,54 @@ namespace AMDiS {
 	  int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
 	  // Test if the current col dof is a periodic dof.
 	  bool periodicCol = meshDistributor->isPeriodicDof(globalColDof);
+	  // Calculate PETSc col index.
+	  int colIndex = globalColDof * dispMult + dispAddCol;
 
-	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
+	  // Ignore all zero entries, expect it is a diagonal entry.
+	  //	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
+	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
 	    continue;
 
 	  if (!periodicCol) {
-	    // Calculate the exact position of the column index in the petsc matrix.
+	    // Calculate the exact position of the column index in the PETSc matrix.
 	    cols.push_back(globalColDof * dispMult + dispAddCol);
 	    values.push_back(value(*icursor));
 	  } else {
-	    std::set<int>& perColAsc = 
-	      meshDistributor->getPerDofAssociations(globalColDof);
+	    // === Row index is not periodic, but column index is. ===
 
+	    // Create set of all periodic associations of the column index.
 	    std::set<int> perAsc;
+	    std::set<int>& perColAsc = 
+	      meshDistributor->getPerDofAssociations(globalColDof);
 	    for (std::set<int>::iterator it = perColAsc.begin(); 
 		 it != perColAsc.end(); ++it)
 	      if (*it >= -3)
 		perAsc.insert(*it);
     
+	    // Scale value to the number of periodic associations of the column index.
 	    double scaledValue = 
-	      value(*icursor) * std::pow(0.5, static_cast<double>(perAsc.size()));
-	    std::vector<int> newCols;
+	      value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));
+
+	    
+	    // === Create set of all matrix column indices due to the periodic ===
+	    // === associations of the column DOF index.                       ===
+
+	    vector<int> newCols;
+
+	    // First, add the original matrix index.
 	    newCols.push_back(globalColDof);
 	    
+	    // And add all periodic matrix indices.
 	    for (std::set<int>::iterator it = perAsc.begin(); 
 		 it != perAsc.end(); ++it) {
-	      int nCols = static_cast<int>(newCols.size());
-	      for (int i = 0; i < nCols; i++) {
+	      for (unsigned int i = 0; i < newCols.size(); i++) {
  		TEST_EXIT(meshDistributor->isPeriodicDof(newCols[i], *it))
- 		  ("Should not happen: %d %d\n", *it, newCols[i]);
+ 		  ("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
+		   *it, newCols[i]);
+
+		MSG("MAP DOF %d -> %d\n",
+		    newCols[i],
+		    meshDistributor->getPeriodicMapping(newCols[i], *it));
 
 		newCols.push_back(meshDistributor->getPeriodicMapping(newCols[i], *it));
 	      }
@@ -143,18 +166,30 @@ namespace AMDiS {
 	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), 
 		     &(cols[0]), &(values[0]), ADD_VALUES);	
       } else {
-	std::map<int, std::vector<int> > colsMap;
-	std::map<int, std::vector<double> > valsMap;
+	// === Row DOF index is periodic. ===
 
+	// Because this row is periodic, we will have to add the entries of this 
+	// matrix row to multiple rows. The following maps store to each row an
+	// array of column indices and values of the entries that must be added to
+	// the PETSc matrix.
+	map<int, vector<int> > colsMap;
+	map<int, vector<double> > valsMap;
+
+	// Traverse all column entries.
 	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
 	     icursor != icend; ++icursor) {
 
 	  // Global index of the current column index.
 	  int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
 
+	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
 	    continue;
 
+
+	  // === Add all periodic associations of both, the row and the column ===
+	  // === indices to the set perAsc.                                    ===
+
 	  std::set<int> perAsc;
 
 	  if (meshDistributor->isPeriodicDof(globalColDof)) {
@@ -173,11 +208,21 @@ namespace AMDiS {
 	    if (*it >= -3)
 	      perAsc.insert(*it);
 
+	  // Scale the value with respect to the number of periodic associations.
 	  double scaledValue = 
-	    value(*icursor) * std::pow(0.5, static_cast<double>(perAsc.size()));
-	  std::vector<std::pair<int, int> > entry;
-	  entry.push_back(std::make_pair(globalRowDof, globalColDof));
+	    value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));
+
+
+	  // === Create all matrix entries with respect to the periodic  ===
+	  // === associations of the row and column indices.             ===
+
+	  vector<pair<int, int> > entry;
+	  
+	  // First, add the original entry.
+	  entry.push_back(make_pair(globalRowDof, globalColDof));
 
+	  // Then, traverse the periodic associations of the row and column indices
+	  // and create the corresponding entries.
 	  for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) {
 	    int nEntry = static_cast<int>(entry.size());
 	    for (int i = 0; i < nEntry; i++) {
@@ -194,23 +239,28 @@ namespace AMDiS {
 		perColDof = entry[i].second;	      	      
 	      
 
-	      entry.push_back(std::make_pair(perRowDof, perColDof));
+	      entry.push_back(make_pair(perRowDof, perColDof));
 	    }
 	  }
 
-	  for (std::vector<std::pair<int, int> >::iterator eIt = entry.begin(); 
+
+	  // === Translate the matrix entries to PETSc's matrix.
+
+	  for (vector<pair<int, int> >::iterator eIt = entry.begin(); 
 	       eIt != entry.end(); ++eIt) {
-	    // Calculate petsc row index.
+	    // Calculate row and column indices of the PETSc matrix.
 	    int rowIndex = eIt->first * dispMult + dispAddRow;
-
 	    int colIndex = eIt->second * dispMult + dispAddCol;
+
 	    colsMap[rowIndex].push_back(colIndex);
 	    valsMap[rowIndex].push_back(scaledValue);
 	  }
 	}
 
 
-	for (std::map<int, std::vector<int> >::iterator rowIt = colsMap.begin();
+	// === Finally, add all periodic rows to the PETSc matrix. ===
+
+	for (map<int, vector<int> >::iterator rowIt = colsMap.begin();
 	     rowIt != colsMap.end(); ++rowIt) {
 	  TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size())
 	    ("Should not happen!\n");
@@ -275,18 +325,18 @@ namespace AMDiS {
     using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
     namespace traits = mtl::traits;
     typedef DOFMatrix::base_matrix_type Matrix;
-    typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
-    typedef std::map<int, DofContainer> RankToDofContainer;
+    typedef vector<pair<int, int> > MatrixNnzEntry;
+    typedef map<int, DofContainer> RankToDofContainer;
 
     // Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
     // that this rank will send to. These nnz entries will be assembled on this rank,
     // but because the row DOFs are not DOFs of this rank they will be send to the
     // owner of the row DOFs.
-    std::map<int, MatrixNnzEntry> sendMatrixEntry;
+    map<int, MatrixNnzEntry> sendMatrixEntry;
 
     // Maps to each DOF that must be send to another rank the rank number of the
     // receiving rank.
-    std::map<DegreeOfFreedom, int> sendDofToRank;
+    map<DegreeOfFreedom, int> sendDofToRank;
 
 
     // First, create for all ranks we send data to MatrixNnzEntry object with 0 entries.
@@ -324,9 +374,9 @@ namespace AMDiS {
 	  
 	  int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor);
 
-	  std::vector<int> rows;
+	  vector<int> rows;
 	  rows.push_back(globalRowDof);
-	  std::vector<int> rowRank;
+	  vector<int> rowRank;
 	  if (meshDistributor->getIsRankDof(*cursor)) {
 	    rowRank.push_back(meshDistributor->getMpiRank());
 	  } else {
@@ -387,7 +437,7 @@ namespace AMDiS {
 		  meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j;
 		
 		sendMatrixEntry[sendToRank].
-		  push_back(std::make_pair(petscRowIdx, petscColIdx));
+		  push_back(make_pair(petscRowIdx, petscColIdx));
 	      }
 	    }
 	    
@@ -409,7 +459,7 @@ namespace AMDiS {
     // === Evaluate the nnz structure this rank got from other ranks and add it to ===
     // === the PETSc nnz data structure.                                           ===
 
-    for (std::map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
+    for (map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
 	 it != stdMpi.getRecvData().end(); ++it) {
       if (it->second.size() > 0) {
 	for (unsigned int i = 0; i < it->second.size(); i++) {