From 423a8b02602f011d58d8e1f78bc576ddd9d74e4d Mon Sep 17 00:00:00 2001
From: Sebastien Blaise <sebastien.blaise@uclouvain.be>
Date: Fri, 9 May 2014 13:13:30 +0000
Subject: [PATCH] Faster linearSystemPETScBlock by removing transpose
---
Solver/linearSystemPETSc.cpp | 28 ++++++----------------------
Solver/linearSystemPETSc.h | 6 +++---
Solver/linearSystemPETSc.hpp | 1 +
3 files changed, 10 insertions(+), 25 deletions(-)
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index fd44f363e9..2e460cdb0a 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -1,6 +1,6 @@
// 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;
}
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index ec8b666b07..23c839747d 100644
--- a/Solver/linearSystemPETSc.h
+++ b/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
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index f271188419..dcbe582da3 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/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)
--
GitLab