[petsc-dev] [petsc-users] ublas sparse matrix bindings?

[Jed Brown]

On Sat, 3 Jul 2010 15:41:53 +0300, Aron Ahmadia <aja2111 at columbia.edu> wrote:
> I am not sure I  understand your argument here (I apologize if I talk past
> you, please feel free to clarify).  I am surmising from the rest of your
> paragraph that you feel it would be most beneficial to users to have the
> full flexibility of mixed-precision with matrices/vectors, i.e. A = A1 + A2,
> where A1 is single-precision and A2 is double-precision.  The layer at which
> you express your linear algebra algorithms should be unaware of the data
> that composes a given matrix A.

Consider a sparse matrix consisting of 12x12 blocks operating on vector
blocks with 3 single precision reals, 6 double precision reals, and 3
double precision complex.  Now write a templated factorization that
operates on these (user-defined) blocks.  Ideally this would merely
involve the user defining

  struct Field {
    float a[3];
    double b[6];
    complex double c[3];
  };
  Matrix<int64_t,Field,Field> A;

but I don't know how to make that work.  All I can think of is an
expression template solution operating on recursive types, which strikes
me as terribly confusing to write and use, as well as perhaps having
padding and ordering issues.

Jed