improved matrix_type check

[Jaroslav Hajek-2]

attached is the improved matrix_type check remade against newest
mainstream repo.
Compiles and make check's.
I've noticed that float sparse matrices are not supported. Will they?

--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[matrix_type.diff]

---

# HG changeset patch
# User Jaroslav Hajek <highegg@...>
# Date 1211382848 -7200
# Node ID e0fdf5deced75945bf75268540bafa6529357630
# Parent  39c1026191e93605b94bad19be3bae85c80e10ff
improved matrix_type check

diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,17 @@
[phone]  Jaroslav Hajek <highegg@...>
+
+ * MatrixType.cc (matrix_real_probe, matrix_complex_probe):
+ New template functions.
+ (MatrixType::MatrixType (const Matrix&),
+ MatrixType::MatrixType (const FloatMatrix&)):
+ just call matrix_real_probe.
+ (MatrixType::MatrixType (const ComplexMatrix&),
+ MatrixType::MatrixType (const FloatComplexMatrix&)):
+ just call matrix_complex_probe.
+
+ * MatrixType.cc (MatrixType::MatrixType (matrix_type, bool)):
+ add missing test for Unknown.
+
 2008-05-20  David Bateman  <dbateman@...>
 
  * Array.cc (Array<T> Array<T>::transpose () const): Modify for tiled
diff --git a/liboctave/MatrixType.cc b/liboctave/MatrixType.cc
--- a/liboctave/MatrixType.cc
+++ b/liboctave/MatrixType.cc
@@ -55,13 +55,15 @@
     }
 }
 
-MatrixType::MatrixType (const Matrix &a)
-  : typ (MatrixType::Unknown),
-    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
-    dense (false), full (true), nperm (0), perm (0)
+template<class T>
+MatrixType::matrix_type
+matrix_real_probe (const MArray2<T>& a)
 {
+  MatrixType::matrix_type typ;
   octave_idx_type nrows = a.rows ();
   octave_idx_type ncols = a.cols ();
+
+  const T zero = 0;
 
   if (ncols == nrows)
     {
@@ -69,36 +71,89 @@
       bool lower = true;
       bool hermitian = true;
 
-      for (octave_idx_type j = 0; j < ncols; j++)
+      // do the checks for lower/upper/hermitian all in one pass.
+      ColumnVector diag(ncols);
+
+      for (octave_idx_type j = 0;
+           j < ncols && upper; j++)
  {
-  if (j < nrows)
-    {
-      if (a.elem (j,j) == 0.)
- {
-  upper = false;
-  lower = false;
-  hermitian = false;
-  break;
- }
-      if (a.elem (j,j) < 0.)
- hermitian = false;
-    }      
-  for (octave_idx_type i = 0; i < j; i++)
-    if (lower && a.elem (i,j) != 0.)
-      {
- lower = false;
- break;
-      }
-  for (octave_idx_type i = j+1; i < nrows; i++)
-    {
-      if (hermitian && a.elem (i, j) != a.elem (j, i))
- hermitian = false;
-      if (upper && a.elem (i,j) != 0)
- upper = false;
-    }
-  if (!upper && !lower && !hermitian)
-    break;
+          T d = a.elem (j,j);
+          upper = upper && (d != zero);
+          lower = lower && (d != zero);
+          hermitian = hermitian && (d > zero);
+          diag(j) = d;
+        }
+
+      for (octave_idx_type j = 0;
+           j < ncols && (upper || lower || hermitian); j++)
+ {
+          for (octave_idx_type i = 0; i < j; i++)
+            {
+              double aij = a.elem (i,j), aji = a.elem (j,i);
+              lower = lower && (aij == zero);
+              upper = upper && (aji == zero);
+              hermitian = hermitian && (aij == aji
+                                        && aij*aij < diag(i)*diag(j));
+            }
  }
+
+      if (upper)
+ typ = MatrixType::Upper;
+      else if (lower)
+ typ = MatrixType::Lower;
+      else if (hermitian)
+ typ = MatrixType::Hermitian;
+      else
+ typ = MatrixType::Full;
+    }
+  else
+    typ = MatrixType::Rectangular;
+
+  return typ;
+}
+
+template<class T>
+MatrixType::matrix_type
+matrix_complex_probe (const MArray2<T>& a)
+{
+  MatrixType::matrix_type typ;
+  octave_idx_type nrows = a.rows ();
+  octave_idx_type ncols = a.cols ();
+
+  const typename T::value_type zero = 0;
+
+  if (ncols == nrows)
+    {
+      bool upper = true;
+      bool lower = true;
+      bool hermitian = true;
+
+      // do the checks for lower/upper/hermitian all in one pass.
+      ColumnVector diag(ncols);
+
+      for (octave_idx_type j = 0;
+           j < ncols && upper; j++)
+ {
+          T d = a.elem (j,j);
+          upper = upper && (d != zero);
+          lower = lower && (d != zero);
+          hermitian = hermitian && (d.real() > zero && d.imag() == zero);
+          diag (j) = d.real();
+        }
+
+      for (octave_idx_type j = 0;
+           j < ncols && (upper || lower || hermitian); j++)
+ {
+          for (octave_idx_type i = 0; i < j; i++)
+            {
+              T aij = a.elem (i,j), aji = a.elem (j,i);
+              lower = lower && (aij == zero);
+              upper = upper && (aji == zero);
+              hermitian = hermitian && (aij == std::conj (aji)
+                                        && std::norm (aij) < diag(i)*diag(j));
+            }
+ }
+
 
       if (upper)
  typ = MatrixType::Upper;
@@ -111,6 +166,16 @@
     }
   else
     typ = MatrixType::Rectangular;
+
+  return typ;
+}
+
+MatrixType::MatrixType (const Matrix &a)
+  : typ (MatrixType::Unknown),
+    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
+    dense (false), full (true), nperm (0), perm (0)
+{
+  typ = matrix_real_probe (a);
 }
 
 MatrixType::MatrixType (const ComplexMatrix &a)
@@ -118,61 +183,7 @@
     sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
     dense (false), full (true), nperm (0), perm (0)
 {
-  octave_idx_type nrows = a.rows ();
-  octave_idx_type ncols = a.cols ();
-
-  if (ncols == nrows)
-    {
-      bool upper = true;
-      bool lower = true;
-      bool hermitian = true;
-
-      for (octave_idx_type j = 0; j < ncols; j++)
- {
-  if (j < ncols)
-    {
-      if (imag(a.elem (j,j)) == 0. &&
-  real(a.elem (j,j)) == 0.)
- {
-  upper = false;
-  lower = false;
-  hermitian = false;
-  break;
- }
-
-      if (imag(a.elem (j,j)) != 0. ||
-  real(a.elem (j,j)) < 0.)
-    hermitian = false;
-    }
-  for (octave_idx_type i = 0; i < j; i++)
-    if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
-      {
- lower = false;
- break;
-      }
-  for (octave_idx_type i = j+1; i < nrows; i++)
-    {
-      if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
- hermitian = false;
-      if (upper && (real(a.elem (i,j)) != 0 ||
-    imag(a.elem (i,j)) != 0))
- upper = false;
-    }
-  if (!upper && !lower && !hermitian)
-    break;
- }
-
-      if (upper)
- typ = MatrixType::Upper;
-      else if (lower)
- typ = MatrixType::Lower;
-      else if (hermitian)
- typ = MatrixType::Hermitian;
-      else if (ncols == nrows)
- typ = MatrixType::Full;
-    }
-  else
-    typ = MatrixType::Rectangular;
+  typ = matrix_complex_probe (a);
 }
 
 
@@ -181,57 +192,7 @@
     sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
     dense (false), full (true), nperm (0), perm (0)
 {
-  octave_idx_type nrows = a.rows ();
-  octave_idx_type ncols = a.cols ();
-
-  if (ncols == nrows)
-    {
-      bool upper = true;
-      bool lower = true;
-      bool hermitian = true;
-
-      for (octave_idx_type j = 0; j < ncols; j++)
- {
-  if (j < nrows)
-    {
-      if (a.elem (j,j) == 0.)
- {
-  upper = false;
-  lower = false;
-  hermitian = false;
-  break;
- }
-      if (a.elem (j,j) < 0.)
- hermitian = false;
-    }      
-  for (octave_idx_type i = 0; i < j; i++)
-    if (lower && a.elem (i,j) != 0.)
-      {
- lower = false;
- break;
-      }
-  for (octave_idx_type i = j+1; i < nrows; i++)
-    {
-      if (hermitian && a.elem (i, j) != a.elem (j, i))
- hermitian = false;
-      if (upper && a.elem (i,j) != 0)
- upper = false;
-    }
-  if (!upper && !lower && !hermitian)
-    break;
- }
-
-      if (upper)
- typ = MatrixType::Upper;
-      else if (lower)
- typ = MatrixType::Lower;
-      else if (hermitian)
- typ = MatrixType::Hermitian;
-      else if (ncols == nrows)
- typ = MatrixType::Full;
-    }
-  else
-    typ = MatrixType::Rectangular;
+  typ = matrix_real_probe (a);
 }
 
 MatrixType::MatrixType (const FloatComplexMatrix &a)
@@ -239,61 +200,7 @@
     sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
     dense (false), full (true), nperm (0), perm (0)
 {
-  octave_idx_type nrows = a.rows ();
-  octave_idx_type ncols = a.cols ();
-
-  if (ncols == nrows)
-    {
-      bool upper = true;
-      bool lower = true;
-      bool hermitian = true;
-
-      for (octave_idx_type j = 0; j < ncols; j++)
- {
-  if (j < ncols)
-    {
-      if (imag(a.elem (j,j)) == 0. &&
-  real(a.elem (j,j)) == 0.)
- {
-  upper = false;
-  lower = false;
-  hermitian = false;
-  break;
- }
-
-      if (imag(a.elem (j,j)) != 0. ||
-  real(a.elem (j,j)) < 0.)
-    hermitian = false;
-    }
-  for (octave_idx_type i = 0; i < j; i++)
-    if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
-      {
- lower = false;
- break;
-      }
-  for (octave_idx_type i = j+1; i < nrows; i++)
-    {
-      if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
- hermitian = false;
-      if (upper && (real(a.elem (i,j)) != 0 ||
-    imag(a.elem (i,j)) != 0))
- upper = false;
-    }
-  if (!upper && !lower && !hermitian)
-    break;
- }
-
-      if (upper)
- typ = MatrixType::Upper;
-      else if (lower)
- typ = MatrixType::Lower;
-      else if (hermitian)
- typ = MatrixType::Hermitian;
-      else if (ncols == nrows)
- typ = MatrixType::Full;
-    }
-  else
-    typ = MatrixType::Rectangular;
+  typ = matrix_complex_probe (a);
 }
 
 MatrixType::MatrixType (const SparseMatrix &a)
@@ -560,54 +467,49 @@
       typ == MatrixType::Tridiagonal ||
       typ == MatrixType::Banded))
  {
-  // Check for symmetry, with positive real diagonal, which
-  // has a very good chance of being symmetric positive
-  // definite..
   bool is_herm = true;
 
-  for (octave_idx_type j = 0; j < ncols; j++)
-    {
-      bool diag_positive = false;
+          // first, check whether the diagonal is positive & extract it
+          ColumnVector diag (ncols);
 
-      for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
- {
-  octave_idx_type ri = a.ridx(i);
+  for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            {
+              is_herm = false;
+              for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
+                {
+                  if (a.ridx(i) == j)
+                    {
+                      double d = a.data(i);
+                      is_herm = d > 0.;
+                      diag(j) = d;
+                      break;
+                    }
+                }
+            }
 
-  if (ri == j)
-    {
-      if (a.data(i) == std::abs(a.data(i)))
- diag_positive = true;
-      else
- break;
-    }
-  else
-    {
-      bool found = false;
 
-      for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++)
- {
-  if (a.ridx(k) == j)
-    {
-      if (a.data(i) == a.data(k))
- found = true;
-      break;
-    }
- }
+          // next, check symmetry and 2x2 positiveness
 
-      if (! found)
- {
-  is_herm = false;
-  break;
- }
-    }
- }
-
-      if (! diag_positive || ! is_herm)
- {
-  is_herm = false;
-  break;
- }
-    }
+          for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            for (octave_idx_type i = a.cidx(j); is_herm && i < a.cidx(j+1); i++)
+              {
+                octave_idx_type k = a.ridx(i);
+                is_herm = k == j;
+                if (is_herm)
+                  continue;
+                double d = a.data(i);
+                if (d*d < diag(j)*diag(k))
+                  {
+                    for (octave_idx_type l = a.cidx(k); l < a.cidx(k+1); l++)
+                      {
+                        if (a.ridx(l) == j)
+                          {
+                            is_herm = a.data(l) == d;
+                            break;
+                          }
+                      }
+                  }
+              }
 
   if (is_herm)
     {
@@ -886,54 +788,50 @@
       typ == MatrixType::Tridiagonal ||
       typ == MatrixType::Banded))
  {
-  // Check for symmetry, with positive real diagonal, which
-  // has a very good chance of being symmetric positive
-  // definite..
   bool is_herm = true;
 
-  for (octave_idx_type j = 0; j < ncols; j++)
-    {
-      bool diag_positive = false;
+          // first, check whether the diagonal is positive & extract it
+          ColumnVector diag (ncols);
 
-      for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
- {
-  octave_idx_type ri = a.ridx(i);
+  for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            {
+              is_herm = false;
+              for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
+                {
+                  if (a.ridx(i) == j)
+                    {
+                      Complex d = a.data(i);
+                      is_herm = d.real() > 0. && d.imag() == 0.;
+                      diag(j) = d.real();
+                      break;
+                    }
+                }
+            }
 
-  if (ri == j)
-    {
-      if (a.data(i) == std::abs(a.data(i)))
- diag_positive = true;
-      else
- break;
-    }
-  else
-    {
-      bool found = false;
+          // next, check symmetry and 2x2 positiveness
 
-      for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++)
- {
-  if (a.ridx(k) == j)
-    {
-      if (a.data(i) == conj(a.data(k)))
- found = true;
-      break;
-    }
- }
+          for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            for (octave_idx_type i = a.cidx(j); is_herm && i < a.cidx(j+1); i++)
+              {
+                octave_idx_type k = a.ridx(i);
+                is_herm = k == j;
+                if (is_herm)
+                  continue;
+                Complex d = a.data(i);
+                if (std::norm (d) < diag(j)*diag(k))
+                  {
+                    d = std::conj (d);
+                    for (octave_idx_type l = a.cidx(k); l < a.cidx(k+1); l++)
+                      {
+                        if (a.ridx(l) == j)
+                          {
+                            is_herm = a.data(l) == d;
+                            break;
+                          }
+                      }
+                  }
+              }
 
-      if (! found)
- {
-  is_herm = false;
-  break;
- }
-    }
- }
-
-      if (! diag_positive || ! is_herm)
- {
-  is_herm = false;
-  break;
- }
-    }
 
   if (is_herm)
     {
@@ -953,10 +851,11 @@
     bandden (0), upper_band (0), lower_band (0),
     dense (false), full (_full), nperm (0), perm (0)
 {
-  if (t == MatrixType::Full || t == MatrixType::Diagonal ||
-      t == MatrixType::Permuted_Diagonal || t == MatrixType::Upper ||
-      t == MatrixType::Lower || t == MatrixType::Tridiagonal ||
-      t == MatrixType::Tridiagonal_Hermitian || t == MatrixType::Rectangular)
+  if (t == MatrixType::Unknown || t == MatrixType::Full
+      || t == MatrixType::Diagonal || t == MatrixType::Permuted_Diagonal
+      || t == MatrixType::Upper || t == MatrixType::Lower
+      || t == MatrixType::Tridiagonal || t == MatrixType::Tridiagonal_Hermitian
+      || t == MatrixType::Rectangular)
     typ = t;
   else
     (*current_liboctave_warning_handler) ("Invalid matrix type");
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
[phone]  Jaroslav Hajek <highegg@...>
+
+ * DLD-FUNCTIONS/matrix_type.cc: Fix tests relying on the
+ older more optimistic hermitian check.
+
 2008-05-21  John W. Eaton  <jwe@...>
 
  * pt-idx.h (tree_index_expression::tree_index_expression (int, int)):
diff --git a/src/DLD-FUNCTIONS/matrix_type.cc b/src/DLD-FUNCTIONS/matrix_type.cc
--- a/src/DLD-FUNCTIONS/matrix_type.cc
+++ b/src/DLD-FUNCTIONS/matrix_type.cc
@@ -512,9 +512,9 @@
 %!test
 %! bnd=spparms("bandden");
 %! spparms("bandden",0.5);
-%! a = spdiags(randn(10,3),[-1,0,1],10,10);
+%! a = spdiags(rand(10,3)-0.5,[-1,0,1],10,10);
 %! assert(matrix_type(a),"Tridiagonal");
-%! assert(matrix_type(abs(a')+abs(a)),"Tridiagonal Positive Definite");
+%! assert(matrix_type(a'+a+2*speye(10)),"Tridiagonal Positive Definite");
 %! spparms("bandden",bnd);
 %!test
 %! bnd=spparms("bandden");
@@ -551,14 +551,14 @@
 %! bnd=spparms("bandden");
 %! spparms("bandden",0.5);
 %! assert(matrix_type(spdiags(1i*randn(10,3),[-1,0,1],10,10)),"Tridiagonal");
-%! a = 1i*randn(9,1);a=[[a;0],ones(10,1),[0;-a]];
+%! a = 1i*(rand(9,1)-0.5);a=[[a;0],ones(10,1),[0;-a]];
 %! assert(matrix_type(spdiags(a,[-1,0,1],10,10)),"Tridiagonal Positive Definite");
 %! spparms("bandden",bnd);
 %!test
 %! bnd=spparms("bandden");
 %! spparms("bandden",0.5);
 %! assert(matrix_type(spdiags(1i*randn(10,4),[-2:1],10,10)),"Banded");
-%! a = 1i*randn(9,2);a=[[a;[0,0]],ones(10,1),[[0;-a(:,2)],[0;0;-a(1:8,1)]]];
+%! a = 1i*(rand(9,2)-0.5);a=[[a;[0,0]],ones(10,1),[[0;-a(:,2)],[0;0;-a(1:8,1)]]];
 %! assert(matrix_type(spdiags(a,[-2:2],10,10)),"Banded Positive Definite");
 %! spparms("bandden",bnd);
 %!test

---

[John W. Eaton]

On 21-May-2008, Jaroslav Hajek wrote:

| attached is the improved matrix_type check remade against newest
| mainstream repo.
| Compiles and make check's.

I applied this patch.

Thanks,

jwe

[dbateman3]

Jaroslav Hajek wrote:
> attached is the improved matrix_type check remade against newest
> mainstream repo.
> Compiles and make check's.
> I've noticed that float sparse matrices are not supported. Will they?

I hesitated to do that, but I don't see why not to add that except for
the work it implies. It seems to me that float sparse matrices are
probably more useful than float matrices.

D.





>
>

[Jaroslav Hajek-2]

On Wed, May 21, 2008 at 8:53 PM, David Bateman <adb014@...> wrote:
> Jaroslav Hajek wrote:
>> attached is the improved matrix_type check remade against newest
>> mainstream repo.
>> Compiles and make check's.
>> I've noticed that float sparse matrices are not supported. Will they?
>
> I hesitated to do that, but I don't see why not to add that except for
> the work it implies.

I didn't mean my question as asking for them. I guessed that part of
the reason was the fact that UMFPack et al. doesn't support single
precision, so the linear algebra stuff would have to be done via
conversion anyway. OTOH, single precision support available for BLAS,
LAPACK, FFTW and libraries in libcruft.

> It seems to me that float sparse matrices are
> probably more useful than float matrices.

Why?

>
> D.
>


--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[David Bateman]

Jaroslav Hajek wrote:
> I didn't mean my question as asking for them. I guessed that part of
> the reason was the fact that UMFPack et al. doesn't support single
> precision, so the linear algebra stuff would have to be done via
> conversion anyway. OTOH, single precision support available for BLAS,
> LAPACK, FFTW and libraries in libcruft.
>  
Yes that is one of the major pains of trying to port single precision
types to sparse matrices/

>> It seems to me that float sparse matrices are
>> probably more useful than float matrices.
>>    
>
> Why?
>  
If you are using sparse matrices, one of your concerns is that you don't
have enough memory to store the dense version of the matrix... So having
single precision would allow a matrix size roughly sqrt(2) large at the
same density without having to use specialized out-of-core techniques.

Regards
David


--
David Bateman                                David.Bateman@...
Motorola Labs - Paris                        [phone] (Ph)
Parc Les Algorithmes, Commune de St Aubin    [phone] (Mob)
91193 Gif-Sur-Yvette FRANCE                  [phone] (Fax)

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary

[Jaroslav Hajek-2]

On Thu, May 22, 2008 at 2:33 PM, David Bateman
<David.Bateman@...> wrote:

> Jaroslav Hajek wrote:
>> I didn't mean my question as asking for them. I guessed that part of
>> the reason was the fact that UMFPack et al. doesn't support single
>> precision, so the linear algebra stuff would have to be done via
>> conversion anyway. OTOH, single precision support available for BLAS,
>> LAPACK, FFTW and libraries in libcruft.
>>
> Yes that is one of the major pains of trying to port single precision
> types to sparse matrices/
>
>>> It seems to me that float sparse matrices are
>>> probably more useful than float matrices.
>>>
>>
>> Why?
>>
> If you are using sparse matrices, one of your concerns is that you don't
> have enough memory to store the dense version of the matrix... So having
> single precision would allow a matrix size roughly sqrt(2) large at the
> same density without having to use specialized out-of-core techniques.
>

But you can say the same for dense matrices. I don't see why
stretching your memory to its limits should be more common with sparse
matrices.

Another concern may be speed - especially if it's memory traffic that
slows the computations down.

> Regards
> David
>
>
> --
> David Bateman                                David.Bateman@...
> Motorola Labs - Paris                        [phone] (Ph)
> Parc Les Algorithmes, Commune de St Aubin    [phone] (Mob)
> 91193 Gif-Sur-Yvette FRANCE                  [phone] (Fax)
>
> The information contained in this communication has been classified as:
>
> [x] General Business Information
> [ ] Motorola Internal Use Only
> [ ] Motorola Confidential Proprietary
>
>

--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[dbateman]

Jaroslav Hajek-2 wrote:

But you can say the same for dense matrices. I don't see why
stretching your memory to its limits should be more common with sparse
matrices.

Ok given it depends on the person and their problem. But if that's the case single precision sparse matrices are as useful as full ones.

Another concern may be speed - especially if it's memory traffic that
slows the computations down.

True, but the same is true for full matrices as well. Ok I give up single precision full matrices are as useful as single precision sparse matrices :-)

Cheers
David