diff src/uBLAS_factory.cc @ 247:8ca45824938e

Add factories hierarchy for matrices' and vectors' assembly
author Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
date Sat, 12 Jul 2014 12:30:13 +0200
parents
children 61830a4f9ab9
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/uBLAS_factory.cc	Sat Jul 12 12:30:13 2014 +0200
@@ -0,0 +1,98 @@
+/*
+ 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 "uBLAS_factory.h"
+
+femfenics_base_factory const&
+uBLAS_factory::instance (void)
+  {
+    static uBLAS_factory const theinstance;
+    return theinstance;
+  }
+
+octave_value
+uBLAS_factory::matrix (dolfin::Matrix const& A) const
+  {
+    octave_value retval;
+
+    // Get capacity of the dolfin sparse matrix
+    boost::tuples::tuple<const std::size_t*, 
+                         const std::size_t*, 
+                         const double*, int> 
+      aa = A.data ();
+
+    int nnz = aa.get<3> ();
+    std::size_t nr = A.size (0), nc = A.size (1);
+    std::vector<double> data_tmp;
+    std::vector<std::size_t> cidx_tmp;
+
+    dim_vector dims (nnz, 1);
+    octave_idx_type nz = 0, ii = 0;
+    Array<octave_idx_type> 
+      ridx (dims, 0), 
+      cidx (dims, 0);
+    Array<double> data (dims, 0);
+
+    octave_idx_type* orow = ridx.fortran_vec ();
+    octave_idx_type* oc = cidx.fortran_vec ();
+    double* ov = data.fortran_vec ();
+
+    for (std::size_t i = 0; i < nr; ++i)
+      {
+        A.getrow (i, cidx_tmp, data_tmp);
+        nz += cidx_tmp.size ();
+
+        for (octave_idx_type j = 0; 
+             j < cidx_tmp.size (); ++j)
+          {
+            orow [ii + j] = i;
+            oc [ii + j] = cidx_tmp [j];
+            ov [ii + j] = data_tmp [j];
+          }
+
+        ii = nz;
+      }
+
+    dims(0) = ii;
+    ridx.resize (dims);
+    cidx.resize (dims);
+    data.resize (dims);
+
+    SparseMatrix sm (data, ridx, cidx, nr, nc);
+    retval = sm;
+
+    return retval;
+  }
+
+octave_value
+uBLAS_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;
+  }
\ No newline at end of file