# HG changeset patch
# User Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
# Date 1405875557 -7200
# Node ID b1dc980506344ca4a31cab455ad2e23a58e28c83
# Parent  8f309b85bb7eb9c0ba7e1ddfa7aa29e7c66aea86
Add support for PETSc algebra back-end

diff -r 8f309b85bb7e -r b1dc98050634 inst/linear_algebra_backend.m
--- a/inst/linear_algebra_backend.m	Sun Jul 13 19:25:03 2014 +0200
+++ b/inst/linear_algebra_backend.m	Sun Jul 20 18:59:17 2014 +0200
@@ -19,12 +19,15 @@
 ##
 ## Query or set the internal variable that specifies the linear algebra back-end
 ## to use when assembling systems, matrices or vectors. It defaults to "uBLAS".
+## uBLAS and PETSc are currently available.
 ##
 ## @end deftypefn
 
 function output = linear_algebra_backend (varargin)
+
   persistent backend = "uBLAS";
-  persistent opts = {"uBLAS"};
+  persistent opts = {"uBLAS",
+                     "PETSc"};
 
   if (nargin > 1)
     print_usage ();
diff -r 8f309b85bb7e -r b1dc98050634 src/Makefile.in
--- a/src/Makefile.in	Sun Jul 13 19:25:03 2014 +0200
+++ b/src/Makefile.in	Sun Jul 20 18:59:17 2014 +0200
@@ -140,7 +140,10 @@
 uBLAS_factory.o: uBLAS_factory.cc uBLAS_factory.h femfenics_base_factory.h
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@
 
-libfemfenics_factories.a: femfenics_factory.o uBLAS_factory.o
+PETSc_factory.o: PETSc_factory.cc PETSc_factory.h femfenics_base_factory.h
+	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@
+
+libfemfenics_factories.a: femfenics_factory.o uBLAS_factory.o PETSc_factory.o
 	$(AR) $(ARFLAGS) $@ $^
 
 clean:
@@ -148,7 +151,7 @@
 
 cleanall:
 	$(RM) *.o core octave-core *.oct *~ *.xml *.status *.log \
-          octave-workspace configure *.pvd *.vtu
+          octave-workspace configure *.pvd *.vtu *.a
 	$(RM) -r autom4te.cache
 	$(RM) ../inst/private/get_vars.m
 	$(RM) Makefile
diff -r 8f309b85bb7e -r b1dc98050634 src/PETSc_factory.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/PETSc_factory.cc	Sun Jul 20 18:59:17 2014 +0200
@@ -0,0 +1,110 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ This program is free software; you can redistribute it and/or modify it under
+ the terms of the GNU General Public License as published by the Free Software
+ Foundation; either version 3 of the License, or (at your option) any later
+ version.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along with
+ this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "PETSc_factory.h"
+
+femfenics_base_factory const&
+PETSc_factory::instance (void)
+  {
+    static PETSc_factory const theinstance;
+    return theinstance;
+  }
+
+octave_value
+PETSc_factory::matrix (dolfin::Matrix const& A) const
+  {
+    octave_value retval;
+
+    std::size_t nr = A.size (0), nc = A.size (1);
+
+    dim_vector dims (0, 1);
+    Array <octave_idx_type> 
+      ridx (dims, 0), 
+      cidx (dims, 0);
+    Array<double> data (dims, 0);
+
+    for (std::size_t i = 0; i < nr; ++i)
+      {
+        std::vector <double> data_tmp;
+        std::vector <std::size_t> cidx_tmp;
+        A.getrow (i, cidx_tmp, data_tmp);
+        std::vector <std::size_t> ridx_tmp (cidx_tmp.size (), i);
+        add_to_arrays (ridx, cidx, data, ridx_tmp, cidx_tmp, data_tmp);
+      }
+
+    SparseMatrix sm (data, ridx, cidx, nr, nc);
+    retval = sm;
+
+    return retval;
+  }
+
+octave_value
+PETSc_factory::vector (dolfin::Vector const& b) const
+  {
+    octave_value retval;
+
+    dim_vector dims;
+    dims.resize (2);
+    dims(0) = b.size ();
+    dims(1) = 1;
+    Array <double> myb (dims);
+
+    for (std::size_t i = 0; i < b.size (); ++i)
+      myb.xelem (i) = b[i];
+
+    retval = myb;
+
+    return retval;
+  }
+
+void
+PETSc_factory::add_to_arrays (Array <octave_idx_type> & ridx,
+                              Array <octave_idx_type> & cidx,
+                              Array <double> & values,
+                              std::vector <std::size_t> const& row,
+                              std::vector <std::size_t> const& col,
+                              std::vector <double> const& val)
+  {
+    octave_idx_type pos = cidx.numel ();
+
+    dim_vector dims = cidx.dims ();
+    if (dims(1) == 1)
+      dims(0) += col.size ();
+    else
+      dims(1) += col.size ();
+    ridx.resize (dims, 0);
+    cidx.resize (dims, 0);
+    values.resize (dims, 0);
+
+    std::vector <std::size_t>::const_iterator rit = row.begin (),
+                                                    cit = col.begin ();
+    std::vector <double>::const_iterator vit = val.begin ();
+
+    while (cit != col.end ())
+      {
+        ridx(pos) = *rit;
+        cidx(pos) = *cit;
+        values(pos) = *vit;
+
+        ++rit;
+        ++cit;
+        ++vit;
+        ++pos;
+      }
+
+    return;
+  }
\ No newline at end of file
diff -r 8f309b85bb7e -r b1dc98050634 src/PETSc_factory.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/PETSc_factory.h	Sun Jul 20 18:59:17 2014 +0200
@@ -0,0 +1,48 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ This program is free software; you can redistribute it and/or modify it under
+ the terms of the GNU General Public License as published by the Free Software
+ Foundation; either version 3 of the License, or (at your option) any later
+ version.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along with
+ this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __FEMFENICS_PETSC_FACTORY__
+#define __FEMFENICS_PETSC_FACTORY__
+
+#include "femfenics_base_factory.h"
+
+class PETSc_factory : public femfenics_base_factory
+{
+  public:
+
+  virtual ~PETSc_factory () {}
+
+  octave_value matrix (dolfin::Matrix const&) const;
+  octave_value vector (dolfin::Vector const&) const;
+
+  static femfenics_base_factory const& instance (void);
+
+  private:
+
+  PETSc_factory () {}
+  PETSc_factory (PETSc_factory const&);
+  PETSc_factory & operator = (PETSc_factory const&);
+
+  static void add_to_arrays (Array <octave_idx_type> &,
+                             Array <octave_idx_type> &,
+                             Array <double> &,
+                             std::vector <std::size_t> const&,
+                             std::vector <std::size_t> const&,
+                             std::vector <double> const&);
+};
+
+#endif
\ No newline at end of file
diff -r 8f309b85bb7e -r b1dc98050634 src/femfenics_factory.cc
--- a/src/femfenics_factory.cc	Sun Jul 13 19:25:03 2014 +0200
+++ b/src/femfenics_factory.cc	Sun Jul 20 18:59:17 2014 +0200
@@ -17,16 +17,16 @@
 
 #include "femfenics_factory.h"
 #include "uBLAS_factory.h"
+#include "PETSc_factory.h"
 #include <octave/parse.h>
 
 femfenics_base_factory const&
 femfenics_factory::factory (void) const
   {
     std::string backend = linear_algebra_backend ();
-    /* Here go the returns for other back-ends
+
     if (backend == "PETSc")
-      return PETSc_factory::instance (); ...
-    */
+      return PETSc_factory::instance ();
 
     // Default back-end
     return uBLAS_factory::instance ();