octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

sparse is_symmetric functions


From: David Bateman
Subject: sparse is_symmetric functions
Date: Mon, 04 Dec 2006 23:54:33 +0100
User-agent: Thunderbird 1.5.0.7 (X11/20060921)

While starting to write svds I found that I need to check if a matrix
passed to eigs is symmetric, rather than symmetric with a positive
diagonal as the cheap test for positive definite matrices does. At the
moment eigs uses arpack's dseupd only for real positive definite
matrices, where in fact it could use and should use it for the case of
[sparse(size(A,1),size(A,1)),A;A',sparse(size(A,2),size(A,2))] as
created by svds. I therefore need the is_symmetric function of the
Matrix and SparseMatrix classes.

However, the current sparse is_symmetric is very inefficient. The
attached patch makes this method of the Sparse classes much better.

D.
Index: liboctave/dSparse.cc
===================================================================
RCS file: /cvs/octave/liboctave/dSparse.cc,v
retrieving revision 1.29
diff -c -r1.29 dSparse.cc
*** liboctave/dSparse.cc        22 Nov 2006 18:57:26 -0000      1.29
--- liboctave/dSparse.cc        4 Dec 2006 22:48:43 -0000
***************
*** 170,181 ****
  bool
  SparseMatrix::is_symmetric (void) const
  {
!   if (is_square () && rows () > 0)
      {
!       for (octave_idx_type i = 0; i < rows (); i++)
!       for (octave_idx_type j = i+1; j < cols (); j++)
!         if (elem (i, j) != elem (j, i))
!           return false;
  
        return true;
      }
--- 170,205 ----
  bool
  SparseMatrix::is_symmetric (void) const
  {
!   octave_idx_type nr = rows ();
!   octave_idx_type nc = cols ();
! 
!   if (nr == nc && nr > 0)
      {
!       for (octave_idx_type j = 0; j < nc; j++)
!       {
!         for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
!           {
!             octave_idx_type ri = ridx(i);
! 
!             if (ri != j)
!               {
!                 bool found = false;
! 
!                 for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++)
!                   {
!                     if (ridx(k) == j)
!                       {
!                         if (data(i) == data(k))
!                           found = true;
!                         break;
!                       }
!                   }
! 
!                 if (! found)
!                   return false;
!               }
!           }
!       }
  
        return true;
      }
Index: liboctave/CSparse.cc
===================================================================
RCS file: /cvs/octave/liboctave/CSparse.cc,v
retrieving revision 1.31
diff -c -r1.31 CSparse.cc
*** liboctave/CSparse.cc        22 Nov 2006 18:57:27 -0000      1.31
--- liboctave/CSparse.cc        4 Dec 2006 22:48:54 -0000
***************
*** 180,191 ****
    octave_idx_type nr = rows ();
    octave_idx_type nc = cols ();
  
!   if (is_square () && nr > 0)
      {
!       for (octave_idx_type i = 0; i < nr; i++)
!       for (octave_idx_type j = i; j < nc; j++)
!         if (elem (i, j) != conj (elem (j, i)))
!           return false;
  
        return true;
      }
--- 180,212 ----
    octave_idx_type nr = rows ();
    octave_idx_type nc = cols ();
  
!   if (nr == nc && nr > 0)
      {
!       for (octave_idx_type j = 0; j < nc; j++)
!       {
!         for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
!           {
!             octave_idx_type ri = ridx(i);
! 
!             if (ri != j)
!               {
!                 bool found = false;
! 
!                 for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++)
!                   {
!                     if (ridx(k) == j)
!                       {
!                         if (data(i) == conj(data(k)))
!                           found = true;
!                         break;
!                       }
!                   }
! 
!                 if (! found)
!                   return false;
!               }
!           }
!       }
  
        return true;
      }

reply via email to

[Prev in Thread] Current Thread [Next in Thread]