Changes in ProblemVec matrix debugging functions.

AMDiS/src/ProblemVec.cc

@@ -660,8 +660,6 @@ namespace AMDiS {
 	// will be set to false).
 	DOFMatrix *matrix = (*systemMatrix)[i][j];
 
-	std::cout << "   i = " << i << "   j = " << j << "   matrix " << matrix << std::endl;
-
 	if (matrix) 
 	  matrix->calculateNnz();
 	
@@ -765,29 +763,6 @@ namespace AMDiS {
 	
 	if (matrix)
 	  nnz += matrix->getBaseMatrix().nnz();
-
-	std::cout << "   i = " << i << "   j = " << j << "   matrix " << matrix << std::endl;
-
-	if (i == 0 && matrix) {
-	  using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
-	  namespace traits= mtl::traits;
-	  typedef DOFMatrix::base_matrix_type Matrix;
-	  
-	  traits::row<Matrix>::type                                 row(matrix->getBaseMatrix());
-	  traits::col<Matrix>::type                                 col(matrix->getBaseMatrix());
-	  traits::const_value<Matrix>::type                         value(matrix->getBaseMatrix());
-	  
-	  typedef traits::range_generator<major, Matrix>::type      cursor_type;
-	  typedef traits::range_generator<nz, cursor_type>::type    icursor_type;
-	  
-	  std::cout.precision(10);
-	  for (cursor_type cursor = begin<major>(matrix->getBaseMatrix()), cend = end<major>(matrix->getBaseMatrix()); cursor != cend; ++cursor)
-	    for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
-	      if (row(*icursor) == 0)
-		std::cout << "C = " << col(*icursor) << " : " << value(*icursor) << std::endl;    	  
-	  std::cout << "============================" << std::endl;
-	}
-
       }
 
       // And now assemble boundary conditions on the vectors
@@ -1464,19 +1439,27 @@ namespace AMDiS {
     typedef traits::range_generator<major, Matrix>::type cursor_type;
     typedef traits::range_generator<nz, cursor_type>::type icursor_type;
 
+    // Create map with all world coordinates of all DOFs.
     DofToCoord dofToCoords;
     createDofToCoordMap(dofToCoords);
 
+    // Start writing output file.
     std::ofstream out(filename.c_str());
+    out.precision(15);
+
+    // First, we write the number of DOFs the file contains.
     out << dofToCoords.size() << std::endl;
 
     int dimOfWorld = Global::getGeo(WORLD);
+    // And now write all the DOF's coords. The first line contains the dof index, all
+    // the other lines contain one component of the coordinates.
     for (DofToCoord::iterator it = dofToCoords.begin(); it != dofToCoords.end(); ++it) {
       out << it->first << std::endl;
       for (int j = 0; j < dimOfWorld; j++)
 	out << (it->second)[j] << std::endl;
     }
 
+    // Write the matrices.
     for (int i = 0; i < nComponents; i++) {
       for (int j = 0; j < nComponents; j++) {
 	DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
@@ -1488,6 +1471,7 @@ namespace AMDiS {
 	int nnz = matrix.nnz();
 	int testNnz = 0;
 	
+	// Write to file, how many entries the matrix conatins.
 	out << nnz << std::endl;
 
 	traits::row<Matrix>::type row(matrix);
@@ -1503,6 +1487,7 @@ namespace AMDiS {
 	  }
 	}
 
+	// Just to check, if nnz() is correct.
 	TEST_EXIT(testNnz == nnz)("NNZ does not fit together!\n");
       }
     }
@@ -1519,6 +1504,7 @@ namespace AMDiS {
     typedef traits::range_generator<major, Matrix>::type cursor_type;
     typedef traits::range_generator<nz, cursor_type>::type icursor_type;
 
+    // Create a map from coords of all DOFs, to the DOF indices in this problem.
     CoordToDof coordToDof;
     createCoordToDofMap(coordToDof);
 
@@ -1527,9 +1513,12 @@ namespace AMDiS {
     DegreeOfFreedom dof;
     int dimOfWorld = Global::getGeo(WORLD);
 
+    // Open file and read the number of DOFs the file contains.
     std::ifstream in(filename.c_str());
     int nReadDofs;
     in >> nReadDofs;
+
+    // Read all DOF indices and their world coordinates.
     for (int i = 0; i < nReadDofs; i++) {
       in >> dof;
       for (int j = 0; j < dimOfWorld; j++)
@@ -1537,6 +1526,7 @@ namespace AMDiS {
       dofToCoord[dof] = coords;
     }
 
+    // Now we traverse all matrices and check them.
     for (int i = 0; i < nComponents; i++) {
       for (int j = 0; j < nComponents; j++) {
 	DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
@@ -1547,6 +1537,8 @@ namespace AMDiS {
 	int readNnz;
 	in >> readNnz;
 
+	// Read each entry in file and check it with the corresponding matrix value 
+	// in this problem.
 	for (int k = 0; k < readNnz; k++) {
 	  DegreeOfFreedom row, col;
 	  double value;
@@ -1558,12 +1550,34 @@ namespace AMDiS {
 	  WorldVector<double> rowCoords = dofToCoord[row];
 	  WorldVector<double> colCoords = dofToCoord[col];
 
-	  if (coordToDof.find(rowCoords) == coordToDof.end())
-	    std::cout << "CANNOT FIND" << std::endl;
+	  if (coordToDof.find(rowCoords) == coordToDof.end()) {
+	    std::cout << "CANNOT ROW FIND" << std::endl;
+	    rowCoords.print();
+	    exit(0);
+	  }
+
+	  if (coordToDof.find(colCoords) == coordToDof.end()) {
+	    std::cout << "CANNOT COL FIND" << std::endl;
+	    rowCoords.print();
+	    exit(0);
+	  }
+
+	  DegreeOfFreedom rowHere = coordToDof[rowCoords];
+	  DegreeOfFreedom colHere = coordToDof[colCoords];
+	  double valueHere = (dofmatrix->getBaseMatrix())[rowHere][colHere];
+	  double diff = fabs(valueHere - value);
+
+	  // And so we check the difference between the value read from file and the
+	  // corresponding value in the matrix.
+	  if (diff > 1e-10) {
+	    std::cout << "WRONG value in matrix[" << i << "][" << j << "]!" << std::endl
+		      << "  DOF in other matrix [" << row << "][" << col << "] = " << value << std::endl
+		      << "  DOF in this matrix  [" << rowHere << "][" << colHere << "] = " << valueHere << std::endl;
+	      
+
+	    exit(0);
+	  }
 
-	  if (coordToDof.find(colCoords) == coordToDof.end())
-	    std::cout << "CANNOT FIND" << std::endl;
-	  
 	}
       }
     }