Re: rank(gf) doesn't really work?

[David Bateman]

Ok, I was trying to be way too smart before and was trying a rank 
revealing LU decomposition on the matrix. However, its makes much more 
sense to just find the linearly independent rows and columns to 
calculate the rank. I've written a m-file attached to do this and will 
commit it to the octave-forge CVS tonight, including the changes to 
remove the old version of the grank function..

I might convert it to an oct-file as the double for-loop that is needed 
in this function might be too slow in an m-file.

Regards
David

--
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

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

## Copyright (C) 2003 David Bateman
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

## -*- texinfo -*-
## @deftypefn {Loadable Function} address@hidden = } grank (@var{a})
## Compute the rank of the Galois array @var{a}, by extracting the independent
## rows and columns.
## @end deftypefn
## @seealso{rank}

## PKG_ADD: dispatch ("rank", "grank", "galois");
function rank = grank (A)
  rank = 0;
  if (nargin != 1)
    error ("grank: single galois matrix required as an argument");
    usage (gfrank);
  elseif (! isgalois (A))
    error ("grank: input must be a galois matrix");
  else
    [r, c] = size(A);
    gzero = gf (0, A.m, A.prim_poly);
    for i = 1:c
      ci = A(:,i);
      idx = find (ci != gzero);
      if (! isempty(idx))
        rank++;
        ci = ci ./ A(idx(1),i);
        for j = (i+1):c
          A(:,j) -= ci .* A(idx(1),j);
        endfor
      endif
    endfor
  endif
endfunction