corrections of amdis-use with petsc-seq, added globalmatrix insertion for petsc

AMDiS/AMDISUse.cmake

@@ -33,7 +33,7 @@ if(AMDIS_HAS_PARALLEL_DOMAIN)
 		endif(METIS_LIB)
 	endif()
 	
-elseif(AMDIS_NEED_PETSC)
+elseif(AMDIS_NEED_SEQ_PETSC)
     find_package(MPI REQUIRED)
     if(MPI_FOUND)
 	    list(APPEND AMDIS_LIBRARIES ${MPI_LIBRARIES})

AMDiS/src/solver/MatrixStreams.h

@@ -185,29 +185,96 @@ namespace AMDiS {
   /// create and fill \ref PetscMatrix
   inline void operator<<(PetscMatrix& mat, const SolverMatrix<Matrix<DOFMatrix*> >& Asolver)
   {
-    const Matrix< DOFMatrix* >& A = *(Asolver.getOriginalMat());
-    
-    std::vector<PetscInt> nnz;      
-    mat.nestMat.resize(A.getNumRows() * A.getNumCols());
-    
-    for (size_t i = 0; i < static_cast<size_t>(A.getNumRows()); i++) {
-      for (size_t j = 0; j < static_cast<size_t>(A.getNumCols()); j++) {
-	size_t idx = i * A.getNumCols() + j;
-	
-	if (A[i][j] == nullptr 
-	    || num_rows(A[i][j]->getBaseMatrix()) == 0 
-	    || num_cols(A[i][j]->getBaseMatrix()) == 0
-	    || A[i][j]->getBaseMatrix().nnz() == 0) {
-	  mat.nestMat[idx] = PETSC_NULL;
-	  continue;
-	}
-	
-	mat.nestMat[idx] << *(A[i][j]);
-      }
-    }
-      
-    // create nested matrix from a vector of block matrices
-    MatCreateNest(PETSC_COMM_SELF, A.getNumRows(), PETSC_NULL, A.getNumCols(), PETSC_NULL, &(mat.nestMat[0]), &mat.matrix);
+	  const Matrix< DOFMatrix* >& A = *(Asolver.getOriginalMat());
+	  if(mat.isNested) {
+
+		  std::vector<PetscInt> nnz;      
+		  mat.nestMat.resize(A.getNumRows() * A.getNumCols());
+
+		  for (size_t i = 0; i < static_cast<size_t>(A.getNumRows()); i++) {
+			  for (size_t j = 0; j < static_cast<size_t>(A.getNumCols()); j++) {
+				  size_t idx = i * A.getNumCols() + j;
+
+				  if (A[i][j] == nullptr 
+						  || num_rows(A[i][j]->getBaseMatrix()) == 0 
+						  || num_cols(A[i][j]->getBaseMatrix()) == 0
+						  || A[i][j]->getBaseMatrix().nnz() == 0) {
+					  mat.nestMat[idx] = PETSC_NULL;
+					  continue;
+				  }
+
+				  mat.nestMat[idx] << *(A[i][j]);
+			  }
+		  }
+
+		  // create nested matrix from a vector of block matrices
+		  MatCreateNest(PETSC_COMM_SELF, A.getNumRows(), PETSC_NULL, A.getNumCols(), PETSC_NULL, &(mat.nestMat[0]), &mat.matrix);
+	  }else {
+		  bool initMatrix = false;
+		  unsigned nRows(0);
+		  unsigned nEntries(0);
+		  std::vector< int > nRowsBlock(A.getNumRows());
+		  for(int i(0); i < A.getNumRows(); ++i) {
+			  int j(0);
+			  for( ; j<A.getNumCols() && A[i][j] == NULL; ++j) ;
+			  if( j == A.getNumCols() ) {
+				  std::stringstream ss; ss << "ERROR: solver matrix has empty row " << i << "\n";
+				  throw std::runtime_error(ss.str());
+			  }
+			  nRowsBlock[i]=num_rows(A[i][j]->getBaseMatrix());
+			  nRows+=nRowsBlock[i];
+		  }
+		  std::vector<PetscInt> nnz(nRows);
+		  //initialize the nnz (could be avoided, but gets complicated
+		  for (int i(0); i < nRows; ++i) nnz[i] = 0;
+
+		  //build a list of pairs (col, value) for each row
+		  std::vector< std::list< std::pair<int, double > > > rowData(nRows);
+		  //use the standard mapper to build global structure
+		  VectorialMapper mapper(Asolver);
+		  PetscInt maxNnz(0);
+		  for (int i(0); i < A.getNumRows(); ++i) {
+			  mapper.setRow(i);
+			  //compute the number of nonzeros per global row
+			  for (int j(0); j < A.getNumCols(); ++j) {
+				  if (A[i][j] != NULL) {
+					  mapper.setCol(j);
+					  const DOFMatrix::base_matrix_type& mtlMat(A[i][j]->getBaseMatrix());
+					  for (unsigned k(0);  k < nRowsBlock[i]; ++k) {
+						  unsigned curRow(mapper.row(k));
+						  nnz[curRow] += mtlMat.nnz_local(k);
+						  maxNnz = std::max(maxNnz, nnz[curRow]);
+						  unsigned curPos(mtlMat.ref_major()[k]);
+						  for(unsigned l(0); l < mtlMat.nnz_local(k); ++l, ++curPos) {
+							  rowData[curRow].push_back(std::make_pair(mapper.col(mtlMat.ref_minor()[curPos]), mtlMat.data[curPos]));
+						  }
+					  }
+				  }
+			  }
+
+		  }
+		  if (mat.matrix == PETSC_NULL) {
+			  MatCreateSeqAIJ(PETSC_COMM_SELF, nRows, nRows, 0, &(nnz[0]), &mat.matrix);
+			  initMatrix = true;
+		  }
+		  //reduce mallocs..
+		  std::vector< PetscInt > colIndices(maxNnz); 
+		  std::vector< PetscScalar > values(maxNnz);
+		  std::list< std::pair< int, double > >::iterator it;
+		  unsigned j;
+		  for (PetscInt i(0); i < nRows; ++i) {
+			  j=0; it = rowData[i].begin();
+			  for (; it != rowData[i].end(); ++it, ++j) {
+				  colIndices[j] = it->first;
+				  values[j] = it->second;
+			  }
+			  MatSetValues(mat.matrix, 1, &i, nnz[i], &colIndices[0], &values[0], INSERT_VALUES);
+		  }
+		  if (initMatrix) {      
+			  MatAssemblyBegin(mat.matrix, MAT_FINAL_ASSEMBLY);
+			  MatAssemblyEnd(mat.matrix, MAT_FINAL_ASSEMBLY);
+		  }
+	  }
   }
   
   
@@ -311,7 +378,7 @@ namespace AMDiS {
     size_t nComponents = rhs.getSize();
     
     // only create nested vector if requested, i.e. flag createNestedVector is set to true
-    if (petscVec.createNestedVector) {
+    if (petscVec.isNested) {
       petscVec.nestVec.resize(nComponents);
 
       for (size_t i = 0; i < nComponents; i++) {

AMDiS/src/solver/PetscSolver.h

@@ -88,7 +88,7 @@ namespace AMDiS {
 
     LinearSolver& oem;
     
-  private:    
+  //private:    
     /// PETSc solver object
     KSP ksp;
 
@@ -143,10 +143,19 @@ namespace AMDiS {
     
     void setNestedVectors(bool createNestedVector = true)
     {
-      vecSol.createNestedVector = createNestedVector;
-      vecRhs.createNestedVector = createNestedVector;
+      vecSol.isNested= createNestedVector;
+      vecRhs.isNested= createNestedVector;
+    }
+ 
+    void setNestedMatrix(bool nested=true)
+    {
+	    petscMat.isNested=nested;
+    } 
+
+    void setNested(bool n)
+    {
+	    setNestedVectors(n); setNestedMatrix(n); 
     }
-  
   protected:    
     Runner runner;
 

AMDiS/src/solver/PetscTypes.h

@@ -34,25 +34,29 @@ namespace AMDiS {
   /// pair of nested matrices and blockmatrix
   struct PetscMatrix
   {
-    PetscMatrix() : matrix(PETSC_NULL) {}
+    PetscMatrix() : 
+	 matrix(PETSC_NULL),
+	 isNested(false)
+	  {}
     
     std::vector<Mat> nestMat;
     Mat matrix;
     
     void destroy();
+    bool isNested;
   };
   
   
   /// pair of nested vectors and blockvector
   struct PetscVector
   {
-    PetscVector() : vector(PETSC_NULL), createNestedVector(false) {}
+    PetscVector() : vector(PETSC_NULL), isNested(false) {}
     
     std::vector<Vec> nestVec;
     Vec vector;
     
     void destroy();
-    bool createNestedVector;
+    bool isNested;
   };