fix block petsc interface in complex arithmetic

Solver/linearSystemPETSc.cpp

 #include "GmshConfig.h"
-#if defined HAVE_PETSC
+
+#if defined(HAVE_PETSC)
 #include "linearSystemPETSc.h"
 #include "fullMatrix.h"
+
 template <>
-void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col, const fullMatrix<double> &val)
+void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
 {
-  fullMatrix<double> &modval = *const_cast<fullMatrix<double> *>(&val);
-  for (int ii = 0; ii<val.size1(); ii++)
-    for (int jj = 0; jj < ii; jj ++) {
-      double buff = modval(ii,jj);
-      modval(ii,jj) = modval (jj,ii);
-      modval(jj,ii) = buff;
+  fullMatrix<PetscScalar> &modval = *const_cast<fullMatrix<PetscScalar> *>(&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;
     }
   PetscInt i = row, j = col;
   _try(MatSetValuesBlocked(_a, 1, &i, 1, &j, &modval(0,0), ADD_VALUES));
   //transpose back so that the original matrix is not modified
-  for (int ii = 0; ii<val.size1(); ii++)
-    for (int jj = 0; jj < ii; jj ++) {
-      double buff = modval(ii,jj);
-      modval(ii,jj) = modval (jj,ii);
-      modval(jj,ii) = buff;
+  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;
     }
 }
+
 template<>
-void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row, const fullMatrix<double> &val)
+void linearSystemPETSc<fullMatrix<PetscScalar> >::addToRightHandSide(int row, const fullMatrix<PetscScalar> &val)
 {
-  for (int ii = 0; ii < _blockSize; ii ++) {
+  for (int ii = 0; ii < _blockSize; ii++) {
     PetscInt i = row * _blockSize + ii;
-     double v = val(ii,0);
+    PetscScalar v = val(ii, 0);
     _try(VecSetValues(_b, 1, &i, &v, ADD_VALUES));
   }
 }
 
 template<>
-void linearSystemPETSc<fullMatrix<double> >::getFromMatrix(int row, int col, fullMatrix<double> &val ) const
+void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromMatrix(int row, int col, fullMatrix<PetscScalar> &val ) const
 {
   Msg::Error("getFromMatrix not implemented for PETSc");
 }
 
 template<>
-void linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide(int row, fullMatrix<double> &val) const
+void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromRightHandSide(int row, fullMatrix<PetscScalar> &val) const
 {
   PetscScalar *tmp;
   _try(VecGetArray(_b, &tmp));
-  for (int i=0; i<_blockSize; i++) {
-    PetscScalar s = tmp[row*_blockSize+i];
+  for (int i = 0; i < _blockSize; i++) {
+    PetscScalar s = tmp[row * _blockSize + i];
 #if defined(PETSC_USE_COMPLEX)
     val(i,0) = s.real();
 #else
@@ -55,12 +58,12 @@ void linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide(int row, fullM
 }
 
 template<>
-void linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row, fullMatrix<double> &val) const
+void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullMatrix<PetscScalar> &val) const
 {
   PetscScalar *tmp;
   _try(VecGetArray(_x, &tmp));
-  for (int i=0; i<_blockSize; i++) {
-    PetscScalar s = tmp[row*_blockSize+i];
+  for (int i = 0; i < _blockSize; i++) {
+    PetscScalar s = tmp[row * _blockSize + i];
 #if defined(PETSC_USE_COMPLEX)
     val(i,0) = s.real();
 #else
@@ -71,10 +74,11 @@ void linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row, fullMatrix
 }
 
 template<>
-void linearSystemPETSc<fullMatrix<double> >::allocate(int nbRows) {
+void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows) 
+{
   clear();
   _try(MatCreate(PETSC_COMM_WORLD, &_a)); 
-  _try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows*_blockSize, nbRows*_blockSize));
+  _try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows * _blockSize, nbRows * _blockSize));
   _try(MatSetType(_a, MATSEQBAIJ));
   // override the default options with the ones from the option
   // database (if any)
@@ -82,7 +86,7 @@ void linearSystemPETSc<fullMatrix<double> >::allocate(int nbRows) {
   _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 4, NULL)); //todo preAllocate off-diagonal part
   //_try(MatMPIBAIJSetPreallocation(_a, _blockSize, 4, NULL, 0, NULL)); //todo preAllocate off-diagonal part
   _try(VecCreate(PETSC_COMM_WORLD, &_x));
-  _try(VecSetSizes(_x, PETSC_DECIDE, nbRows*_blockSize));
+  _try(VecSetSizes(_x, PETSC_DECIDE, nbRows * _blockSize));
   // override the default options with the ones from the option
   // database (if any)
   _try(VecSetFromOptions(_x));
@@ -91,24 +95,26 @@ void linearSystemPETSc<fullMatrix<double> >::allocate(int nbRows) {
 }
 
 #include "Bindings.h"
-void linearSystemPETScRegisterBindings(binding *b) {
+void linearSystemPETScRegisterBindings(binding *b) 
+{
   classBinding *cb;
   methodBinding *cm;
-  cb = b->addClass<linearSystemPETSc<double> >("linearSystemPETSc");
+  cb = b->addClass<linearSystemPETSc<PetscScalar> >("linearSystemPETSc");
   cb->setDescription("A linear system solver, based on PETSc");
-  cm = cb->setConstructor<linearSystemPETSc<double> >();
-  cm->setDescription ("A new PETSc<double> solver");
-  cb->setParentClass<linearSystem<double> >();
+  cm = cb->setConstructor<linearSystemPETSc<PetscScalar> >();
+  cm->setDescription ("A new PETSc<PetscScalar> solver");
+  cb->setParentClass<linearSystem<PetscScalar> >();
   cm->setArgNames(NULL);
-  cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<double> >::systemSolve);
+  cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<PetscScalar> >::systemSolve);
   cm->setDescription("compute x = A^{-1}b");
 
-  cb = b->addClass<linearSystemPETSc<fullMatrix<double> > >("linearSystemPETScBlock");
+  cb = b->addClass<linearSystemPETSc<fullMatrix<PetscScalar> > >("linearSystemPETScBlock");
   cb->setDescription("A linear system solver, based on PETSc");
-  cm = cb->setConstructor<linearSystemPETSc<fullMatrix<double> >, int>();
-  cm->setDescription ("A new PETScBlock<double> solver (we probably should get rid of the blockSize argument)");
+  cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> >, int>();
+  cm->setDescription ("A new PETScBlock<PetscScalar> solver (we probably should get rid of the blockSize argument)");
   cm->setArgNames("blockSize", NULL);
-  cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<double> >::systemSolve);
+  cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<PetscScalar> >::systemSolve);
   cm->setDescription("compute x = A^{-1}b");
 }
+
 #endif // HAVE_PETSC