[ublas] [patch] Fixing the UMFPACK bindings

From: Andreas Klöckner (kloeckner_at_[hidden])
Date: 2005-03-10 11:18:25

Hi there,

When using UBlas from Boost CVS as of some time mid-February and the
bindings from the boost-sandbox CVS as of around the same time, the
UMFPACK bindings are broken. Over the course of this email, I'll explain
my fixes.

UMFPACK deals with sparse matrices in roughly the same format as
ublas::compressed_matrix<double, ublas::column_major> >, except for two
little differences:

- By default, UBlas uses size_t in its index arrays, while UMFPACK uses
int. I tried to convince UBlas to also use ints for its indices by
passing ublas::unbounded_array<int> as an index array to
ublas::compressed_matrix<>, but that only gives tons of "comparison
between signed and unsigned" warnings, as well as a few errors. Judging
UBlas not ready signed index types, I settled for a cast-based solution.
Since casting to (int *) by default in the bindings code is not a good
idea, I gave most of the UMFPACK bindings functions an optional extra
parameter that allows to specify a class that performs the necessary
casts if desired. The class "concept" of this parameter is called
IndexArrayCaster (or IAC, for short). Two models of it are provided:
index_array_caster_no_op (to which the parameter defaults) and
index_array_caster_to_int. See the attached bindings-umfpack-iac.patch
against the boost-sandbox CVS.

- For compressed_matrix<> types in UBlas, the column index array is not
always completely filled. Instead, the filled1_ data member says how
many fields of this index array are actually used. UMFPACK, on the other
hand, demands a fully filled column index array. The previous bindings
code solved this by manually filling up the index array, which is
arguably a hack. To perform this task more cleanly, I've added a method
compressed_matrix<>::complete_index1_data(). I've decided to declare the
method const, since the invariant

   compressed_matrix<double> A, B;
   fillSomehow(A);
   B = A;
   B.complete_index1_data();
   assert(B == A);

holds. See the attached ublas-compressed-complete-index.patch. (The
patch also fixes one comment's indentation and spelling.)

On a different note, UMFPACK Version 4.3 and below only allowed complex
numbers by specifying a matrix as its real and imaginary part,
necessitating a copy of UBlas's matrices, which use std::complex (and
thus packed complex storage). This restriction has been lifted in
UMFPACK 4.4. It might yield substantial speed benefits to take advantage
of this improvement. Unfortunately, UMFPACK's author did not include a
version #define in his headers, so doing this automatically might prove
hard.

I'd be grateful for a commit or additional discussion.

Andreas
--
--------------------------------------+-----------------------------
Dipl.-Math. techn. Andreas Klöckner   | FON [phone]
Institut für Angewandte Mathematik II | FAX [phone]
Universität Karlsruhe                 | EMAIL kloeckner_at_[hidden]
Englerstraße 2                        |
76128 Karlsruhe                       | RAUM 022
Germany                               | Kollegiengebäude Mathematik
--------------------------------------+-----------------------------