Faster linearSystemPETScBlock by removing transpose

Solver/linearSystemPETSc.cpp

 // Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
 //
-// See the LICENSE.txt file for license information. Please report all
+// ee the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@geuz.org>.
 
 #include "GmshConfig.h"
@@ -35,35 +35,19 @@ void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col,
 {
   if (!_entriesPreAllocated)
     preAllocateEntries();
+  PetscInt i = row, j = col;
   #ifdef PETSC_USE_COMPLEX
   fullMatrix<std::complex<double> > modval(val.size1(), val.size2());
   for (int ii = 0; ii < val.size1(); ii++) {
     for (int jj = 0; jj < val.size1(); jj++) {
-      modval(ii, jj) = val (jj, ii);
-      modval(jj, ii) = val (ii, jj);
+      modval(ii, jj) = val (ii, jj);
     }
   }
+  MatSetValuesBlocked(_a, 1, &i, 1, &j, modval.getDataPtr(), ADD_VALUES);
   #else
-  fullMatrix<double> &modval = *const_cast<fullMatrix<double> *>(&val);
-  for (int ii = 0; ii < val.size1(); ii++) {
-    for (int jj = 0; jj < ii; jj++) {
-      PetscScalar buff = modval(ii, jj);
-      modval(ii, jj) = modval (jj, ii);
-      modval(jj, ii) = buff;
-    }
-  }
-  #endif
-  PetscInt i = row, j = col;
-  MatSetValuesBlocked(_a, 1, &i, 1, &j, &modval(0,0), ADD_VALUES);
-  //transpose back so that the original matrix is not modified
-  #ifndef PETSC_USE_COMPLEX
-  for (int ii = 0; ii < val.size1(); ii++)
-    for (int jj = 0; jj < ii; jj++) {
-      PetscScalar buff = modval(ii,jj);
-      modval(ii, jj) = modval (jj,ii);
-      modval(jj, ii) = buff;
-    }
+  MatSetValuesBlocked(_a, 1, &i, 1, &j, val.getDataPtr(), ADD_VALUES);
   #endif
+
   _valuesNotAssembled = true;
 }
 

Solver/linearSystemPETSc.h

@@ -86,7 +86,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   void insertInSparsityPattern (int i, int j);
   bool isAllocated() const { return _isAllocated; }
   void preAllocateEntries();
-  void allocate(int nbRows);
+  virtual void allocate(int nbRows);
   void print();
   void clear();
   void getFromMatrix(int row, int col, scalar &val) const;
@@ -100,7 +100,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   void zeroRightHandSide();
   void zeroSolution();
   void printMatlab(const char *filename) const;
-  int systemSolve();
+  virtual int systemSolve();
   Mat &getMatrix(){ return _a; }
   //std::vector<scalar> getData();
   //std::vector<int> getRowPointers();
@@ -133,7 +133,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
   void zeroRightHandSide() {}
   void zeroSolution() {}
   void printMatlab(const char *filename) const{};
-  int systemSolve() { return 0; }
+  virtual int systemSolve() { return 0; }
   double normInfRightHandSide() const{return 0;}
 };
 #endif

Solver/linearSystemPETSc.hpp

@@ -156,6 +156,7 @@ void linearSystemPETSc<scalar>::allocate(int nbRows)
     else {
       MatSetType(_a, MATSEQBAIJ);
     }
+    MatSetOption(_a, MAT_ROW_ORIENTED, PETSC_FALSE);
   }
   // override the default options with the ones from the option
   // database (if any)