Re: Project progress

[c.]

On 5 Oct 2013, at 21:46, marco Vassallo <address@hidden> wrote:

>>> where the method data is described here [1]
>> from this link I can't see how can you tell if the Matrix is in CSR or CSC?
>> How can you tell the total number of elements?
>> c.
> 
> To be honest I'm not sure of how we can use it.
> As I understand it by now, , the last information should be the the total 
> number of element,
> which is what we need.

I see, I had missed the last input argument, so yes that should give us the 
number of nonzeros
but there still is a proble, it seems from the docs that the returned arrays 
can be either in CSC
or CSR form, while Octave uses only CSC.

Is there a way of knowing which form is being used? 
If there is you could make sure that the matrix is CSC 
and then do (N.B. this an examlpe and is untested):

---------------------------------------------------------------
boost::tuples::tuple<const std::size_t*, const std::size_t*, const double*, 
int> t;
t = fenix_matrix.data ();

mat = SparseMatrix (fenics_matrix.size (0), fenics_matrix.size (1), t.get<3> 
());

const std::size_t* ir = t.get<0> ();
const std::size_t* ic = t.get<1> ();
const double* iv = t.get<2> ();

octave_idx_type* or = mat.ridx ();
octave_idx_type* oc = mat.cidx ();
double* ov = mat.data ();

for (octave_idx_type ii =0; ii < t.get<3> (); ++ii)
{
 or[ii] = ir[ii];
 oc[ii] = ic[ii];
 ov[ii] = iv[ii];
}
---------------------------------------------------------------

If the matrix is CSR, you have invert the assignments to or and oc
and transpose the matrix after copy.


If you don't know whether the matrix is CSC or CSR, yu can just use 
the number of nonzeros to preallocate 
the correct space for the Octave SparseMatrix, but to fill it
it is safer to still use dolfin high level methods.

For example you could do (N.B. this an examlpe and is untested):

----------------------------------------------------------------
dim_vector dims (nnz, 1);

Array<octave_idx_type> I(dims, 0);
octave_idx_type *Iptr;

Array<octave_idx_type> J(dims, 0);
octave_idx_type *Jptr;

Array<double> V(dims, 0);
double *Vptr;

octave_idx_type ival = 0;
std::vector<octave_idx_type> tmpcols;
std::vector<double> tmpvals;
for (octave_idx_value ir = 0; ir < fenics_matrix.size (0); ++ir)
{
 fenics_matrix.getrow (ir, tmpcols, tmpvals);
 for (octave_idx_type ii = 0; ii < tmpcols.length (); ++ii)
 {
   Iptr[ival] = ir;
   Jptr[ival] = tmpcols[ii];
   Vptr[ival++] = tmpvals[ii]; 
 }
}

retval(0) = octave_value (SparseMatrix (V, I, J, nrows, ncols, true));
----------------------------------------------------------------

> But, as Dolfin as well as Octave are using the BLAS,
> I was just wondering if there was any way to avoid the copy of the matrix.

I'm afraid you're a bit confused here,
Octave is using BLAS, not uBLAS.

BLAS only does algebra with full matrices,
for sparse matrices Octave is using SuiteSparse, 
while dolfin can use a set of different backends.

> Marco                                           

HTH,
c.

P.S. did you have time to fix the problem with __makeinfo__ ?