File Coverage

blib/lib/Math/MatrixSparse.pm

Criterion	Covered	Total	%
statement	648	1595	40.6
branch	203	782	25.9
condition	40	191	20.9
subroutine	63	136	46.3
pod	0	106	0.0
total	954	2810	33.9

package Math::MatrixSparse;

use 5.006;
use strict;
use warnings;

use Carp;

require Exporter;

our @ISA = qw(Exporter);

# Items to export into callers namespace by default. Note: do not export
# names by default without a very good reason. Use EXPORT_OK instead.
# Do not simply export all your public functions/methods/constants.

# This allows declaration	use Math::MatrixSparse ':all';
# If you do not need this, moving things directly into @EXPORT or @EXPORT_OK
# will save memory.
our %EXPORT_TAGS = ( 'all' => [ qw(
				   splitkey
				   makekey
) ] );

our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} });

our @EXPORT = qw(

);
our $VERSION = '0.01';

use overload
    "+" => "add",
    "*" => "quickmultiply",
    "x" => "kronecker",
    "-" => "subtract",
    "neg" => "negate",
    "~"  => "transpose",
    "**" => "exponentiate",
    "==" => "equals",
    '""' => "print",
    "&" => "matrixand",
    "|" => "matrixor",
    "fallback" => undef;


# Preloaded methods go here.

### CREATION METHODS
sub new {
    my ($proto,$name) = @_;
    my $this;
    $this->{data} = undef;
    $this->{name} = $name;
    $this->{rows} = 0;
    $this->{columns} =0;
    $this->{special}->{bandwidth} = 0;
    $this->{special}->{structure} = "";
    $this->{special}->{shape} = "";
    $this->{special}->{sign} = "";
    $this->{special}->{pattern} = 0;
    $this->{special}->{field} = "real";

    bless ($this);
    return $this;
}


sub newfromstring {
    my ($proto, $string,$name) = @_;
    my $this=new Math::MatrixSparse;
    my @entries = split(/\n/,$string);
    $this->{special}->{structure} = ".";
    foreach my $entry (@entries) {
	my ($i,$j,$value) = $entry=~m/^\s*(\d+)\s+(\d+)\s+(.+)\s*$/;
	$this->assign($i,$j,$value);
    }
    $this->{name} =$name;
    $this->{special}->{field} = "real";
    return $this;
}

sub newdiag {
    my ($proto, $diag, $name) = @_;
    my @data = @{$diag};
    my $this = new Math::MatrixSparse;
    $this->name($name);
    my $i=1;
    foreach my $entry (@data){
	$this->assign($i,$i,$entry);
	$i++;
    }
    $this->{special}->{structure} = "symmetric";
    $this->{special}->{shape} = "diagonal";
    $this->{special}->{field} = "real";
    $this->{special}->{square} = 1;
    return $this;
}



sub newdiagfromstring {
    my ($proto, $string, $name) = @_;
    my $this=new Math::MatrixSparse;
    my @entries = split(/\n/,$string);
    foreach my $entry (@entries) {
	my ($i,$value) = $entry=~m/^(\d+)\s+(.+)$/;
	$this->assign($i,$i,$value);
    }
    $this->{name} =$name;
    $this->{special}->{structure} = "symmetric";
    $this->{special}->{field} = "real";
    $this->{special}->{shape} = "diagonal";
    $this->{special}->{square} = 1;
    return $this;
}


sub newidentity {
    my ($proto, $n,$m) = @_;
    my $this = new Math::MatrixSparse;
    $m=$n unless ($m);
    $this->{name} = "I$n";
    my  $min = $m<$n ? $m :$n;
    for my $i (1..$min) {
	$this->assign($i,$i,1);
    }
    $this->{rows} = $n;
    $this->{columns} = $m;
    $this->{special}->{structure} = "symmetric";
    $this->{special}->{field} = "real";
    $this->{special}->{shape} = "diagonal";
    $this->{special}->{sign} = "nonnegative";
    $this->{special}->{square} = 1 if $m==$n;
    $this->{special}->{pattern} = 0;
    return $this;
}

sub newrandom {
    my ($proto,$maxrow,$maxcol,$max,$density,$name) = @_;
    $maxcol = $maxrow unless defined $maxcol;
    $density = 1 unless (defined $density)&&($density>=0) && ($density<=1);
    my $stoch = new Math::MatrixSparse($name);
    unless (defined $max) {
	for my $i (1..$maxrow) {
	    for my $j (1..$maxcol) {
		my $uniform = rand;
		next unless $uniform<=$density;
		$stoch->assign($i,$j,rand);
	    }
	}
	return $stoch;
    }
    $name = "" unless defined $name;
    $stoch->assign($maxrow,$maxcol,0);
    my $k = 1;
    while ($k<=$max) {
	my ($i,$j) = (int($maxrow* rand)+1,int($maxcol* rand)+1);
	next if $stoch->element($i,$j);
	my $uniform = rand;
	if ($uniform>$density) {
	    carp "Math::MatrixReal::newrandom ignoring element";
	}
	next unless $uniform <= $density;
	$stoch->assign($i,$j,rand);
	$k++;
    }
    $stoch->{special}->{sign}="nonnegative";
    return $stoch;
}

### INPUT-OUTPUT METHODS
sub newharwellboeing {
    my ($proto, $filename) = @_;
    open(HB,"<$filename" ) || croak "Math::MatrixSparse::newharwellboeing Can't open $filename\n";
    my $this = new Math::MatrixSparse;
    $_ = ;
    chomp();
    my ($ident, $key) = $_ =~ m/^(.*)(........)\s*$/;
    $key =~ s/\s*$//;
    $this->{name} = $key;
    $_ = ;
    chomp();
    my ($lines, $pline, $rline, $nvline, $rhsline) = $_ =~
     m/^\s*(\d{1,14})\s*(\d{1,14})\s*(\d{1,14})\s*(\d{1,14})\s*(\d{1,14})\s*$/;
    $_ = ;
    chomp();
    my ($mattype, $rows, $columns, $rvi, $eme) = $_ =~
         m/^\s*(...)\s*(\d{1,14})\s*(\d{1,14})\s*(\d{1,14})\s*(\d{1,14})\s*$/;
    if ($mattype =~ /^C/) {
	croak "Math::MatrixSparse::newharwellboeing Complex types not implemented, exiting\n";
    } else {
	$this->{special}->{field} = "real";
    }
    $this->{rows} = $rows;
    $this->{columns} = $columns;
    my $formatline = ;
    chomp($formatline);
    my ($colft, $rowft, $valft,$rhsft) = $formatline =~ m/(\([a-zA-Z0-9.+-+]+\))/g;
    unless (($colft =~ /I/i)&&($rowft =~ /I/i)) {
	carp "Math::MatrixSparse::newharwellboeing non-integer format for rows and columns";
    }
    $valft =~ s/\d+P//g;
    my ($valrep, $valsize) = $valft =~ m/(\d+)[a-zA-Z](\d+)/;
    my $valregex = '.' x $valsize;
    my $rhsspec;
    if ($rhsline) {
	$rhsspec = ;
	chomp($rhsspec);
    }
    #now read the column pointer data...
    my @colpointers;
    for my $i  (1..$pline) {
	$_ = ;
	s/^\s*//;
	s/\s*$//;
	push @colpointers, split(/\s+/);
    }
    #...and the row data...
    my @rowindex;
    for my $i (1..$rline) {
	$_ = ;
	s/^\s*//;
	s/\s*$//;
	push @rowindex, split(/\s+/);
    }
    #...and any value data.  If the matrix is a pattern type, fill
    # @values with ones.
    my @values;
    if ($mattype =~ m/^P..$/i) {
	@values = (1) x $rvi;
    } else {
	for my $i (1..$nvline) {
	    $_ = ;
	    s/D/e/g;
	    push @values, map {s/\s//g; $_+0} m/$valregex/g;
	}
    }

    my $curcol = 1;
    foreach my $i (0..($#colpointers-1)) {
	my $thiscol = $colpointers[$i];
	my $nextcol = $colpointers[$i+1];
	my @rowslice = @rowindex[$thiscol-1..$nextcol-2];
	my @valslice = @values[$thiscol-1..$nextcol-2];
	foreach my $j (0..$#rowslice) {
	    $this->assign($rowslice[$j],$curcol,$valslice[$j]);
	}
	$curcol++;
    }
    if ($mattype =~ /^.S.$/) {
	 $this->_symmetrify();
    } elsif ($mattype =~ /^.Z.$/) {
	$this->_skewsymmetrify();
    }

    return $this;
}

sub newmatrixmarket {
    my ($proto, $filename) = @_;
    my $this = new Math::MatrixSparse;
    confess "Math::MatrixSparse::newmatrixmarket undefined filename" unless defined $filename;
    open(MM,"<$filename" ) || croak "Math::MatrixSparse::newmatrixmarket Can't open $filename for reading";
    $_=;
    unless (/^\%\%MatrixMarket/) {
	confess "Math::MatrixSparse::newmatrixmarket Invalid start of file";
    }
    unless (/coordinate/i) {
	confess "Math::MatrixSparse::newmatrixmarket dense format not implemented";
    }
    if (/complex/i) {
	carp "Math::MatrixSparse::newmatrixmarket Complex matrices not implemented, ignoring imaginary part\n";
	$this->{special}->{field} = "real";
    } else {
	$this->{special}->{field} = "real";
    }
    my $specifications = $_;
    my $ispattern;
    my $issymmetric;
    my $isskewsymmetric;
    if ($specifications =~ /pattern/i) {
	$ispattern = 1;
	$this->{special}->{pattern} = 1;
    }
    $this->{name} = $filename;
    my $startdata=0;
    my $entries =0;
    my ($rows,$columns);

    while () {
	next if /^\%/;
	unless ($startdata) {
	    s/^\s*//;
	    s/\s*$//;
	    ($rows,$columns,$entries) = split(/\s+/);
	    $this->{rows} = $rows;
	    $this->{columns} = $columns;
	    $startdata = 1;
	} else {
	    s/^\s*//;
	    s/\s*$//;
	    ($rows,$columns,$entries) = split(/\s+/);
	    $entries = 1 if $ispattern;
	    $this->assign($rows,$columns,$entries);
	}
    }
    if ($specifications =~ /\Wsymmetric/i) {
	$this->{special}->{structure} = "symmetric";
	return $this->symmetrify();
    }
    if ($specifications =~ /skewsymmetric/i) {
	$this->{special}->{structure} = "skewsymmetric";
	return $this->skewsymmetrify();
    }
#    if ($specifications =~ /hermetian/i) {
#	$this->{special}->{structure} = "hermetian";
#	return $this->symmetrify();
#    }
    return $this;
}

sub writematrixmarket {
    my ($matrix, $filename) = @_;
    open(MM,">$filename" ) || croak "Math::MatrixSparse::newmatrixmarket Can't open $filename for writing";
    print MM '%%MatrixMarket matrix coordinate ';
    print MM $matrix->{special}->{pattern} ? "pattern " : "real ";
    if ($matrix->{special}->{structure} =~ m/^symmetric/i) {
	$matrix=$matrix->symmetrify();
	print MM "symmetric";
    } elsif ($matrix->{special}->{structure} =~ m/skewsymmetric/i) {
	$matrix=$matrix->skewsymmetrify();
	print MM "skewsymmetric";
    } else {
	print MM "general\n";
    }
    print MM "$matrix->{rows} $matrix->{columns} ";
    print MM $matrix->count() , "\n";
    if ($matrix->{special}->{pattern}) {
	foreach my $key ($matrix->sortbycolumn()) {
	    my ($i,$j) = &splitkey($key);
	    next unless $matrix->element($i,$j);
	    print MM "$i $j\n";
	}
    } else {
	foreach my $key ($matrix->sortbycolumn()) {
	    my ($i,$j) = &splitkey($key);
	    print MM "$i $j ", $matrix->{data}{$key}, "\n";
	}
    }
    return;
}


sub copy {
    my ($proto,$name) = @_;
    my $this = new Math::MatrixSparse;
    return $this unless defined $proto;
    if (defined $proto->{data}) {
	%{$this->{data}} = %{$proto->{data}};
    }
    %{$this->{special}} = %{$proto->{special}};
    $this->{name} = defined $name ? $name : $proto->{name};
    $this->{rows} = $proto->{rows};
    $this->{columns} = $proto->{columns};
    return $this;
}

sub name {
    my ($object,$name) = @_;
    $object->{name} = $name;
    return $name;
}

### INSERTION AND LOOKUP METHODS
sub assign {
    my ($object, $i,$j,$x)=@_;
    return undef unless ((defined $i) && (defined $j)&&(defined $object));
    $x = 1 if $object->{special}->{pattern};
    return undef unless defined $x;
    #update matrix's shape if necessary.
    if ((defined $object->{special}) &&
	 ($object->{special}->{shape} =~ m/diagonal/i) &&
	 ($i!= $j) ) {
 	if ($i<$j) {
 	    $object->{special}->{shape} = "upper"
 	} else  {
 	    $object->{special}->{shape} = "lower"
 	}
     } elsif (($object->{special}->{shape} =~ m/strictlower/i)
	      && ($i<= $j) ) {
 	if ($i==$j) {
  	    $object->{special}->{shape} = "lower";
 	} else {
 	    $object->{special}->{shape}="";
 	}
     } elsif (($object->{special}->{shape} =~ m/strictupper/i)
	      && ($i>= $j) ) {
 	if ($i==$j) {
 	    $object->{special}->{shape} = "upper";
 	} else {
 	    $object->{special}->{shape}="";
 	}
     } elsif (($object->{special}->{shape} =~ m/^lower/i) && ($i< $j) ) {
 	$object->{special}->{shape}="";
     } elsif (($object->{special}->{shape} =~ m/^upper/i) && ($i= $j) ) {
        $object->{special}->{shape}="";
     }
    $object->{special}->{pattern} = 0 unless ($x==1);
    #update bandwidth
    if (abs($i-$j) > $object->{special}->{bandwidth}) {
	$object->{special}->{bandwidth} = abs($i-$j);
    }
    #update symmetric and skew-symmetric structure if necessary
     if (
 	($i!=$j) && ( defined $object->{special}->{structure}) &&
  	(
 	 ($object->{special}->{structure} =~ /symmetric/i)
# 	 ||($object->{special}->{structure} =~ /hermetian/i)
 	 )
 	) {
 	$object->{special}->{structure} = "";
     } elsif (($i==$j)&&
 	     ($object->{special}->{structure} =~ /skewsymmetric/i)&&
	      ($x)) {
 	#skew-symmetric matrices must have zero diagonal
 	$object->{special}->{structure} = "";
     }
    #update sign if necessary
     if ( ($object->{special}->{sign} =~ /^positive/) && ($x<=0) ) {
 	if ($x<0) {
 	    $object->{special}->{sign}="";
 	} else {
 	    $object->{special}->{sign} = "nonnegative";
 	}
     } elsif ( ($object->{special}->{sign} =~ /^negative/) && ($x>=0) ) {
 	if ($x>0) {
 	     $object->{special}->{sign}="";
 	} else {
 	    $object->{special}->{sign} = "nonpositive";
 	}
     } elsif ( ($object->{special}->{sign} =~ /nonnegative/i) && ($x<0) ) {
 	 $object->{special}->{sign}="";
     } elsif ( ($object->{special}->{sign} =~ /nonpositive/i) && ($x>0) ) {
 	 $object->{special}->{sign}="";
     }
    my $key = &makekey($i,$j);
    delete $object->{sortedrows};
    delete $object->{sortedcolumns};
    delete $object->{data}->{$key};
    $object->{data}->{$key} = $x;
    #update size of matrix, and squareness, if necessary
    $object->{rows} = $i if ($i>$object->{rows});
    $object->{columns} = $j if ($j>$object->{columns});
    $object->{special}->{square} = ( $object->{columns}==$object->{rows});
    return $x;
}

sub assignspecial {
    #as assign, except that it respects special properties of
    #the matrix.  For example, symmetrical matrices are kept symmetric.
    my ($object, $i,$j,$x)=@_;
    return undef unless ((defined $i) && (defined $j)&&(defined $object));
    $x = 1 if $object->{special}->{pattern};
    return undef unless defined $x;
    my $key = &makekey($i,$j);
    my $revkey = &makekey($j,$i);
    if ($object->{special}->{structure} =~ m/^symmetric/i) {
	if ($i==$j)  {
	    $object->{data}{$key} = $x;
	} else {
	    $object->{data}{$key} = $x;
	    $object->{data}{$revkey} = $x;
	}
    } elsif ($object->{special}->{structure} =~ m/^symmetric/i) {
	if (($i==$j)&&($x)) {
	    croak "Math::MatrixSparse::assignspecial skewsymmetric matrices must have zero diagonal";
	} else {
	    $object->{data}{$key} = $x;
	    $object->{data}{$revkey} = -$x;
	}
    } else {
	$object->assign($i,$j,$x);
    }
    return $x;
}

sub assignkey {
    my ($object, $key,$x)=@_;
    my ($i,$j) = &splitkey($key);
    return unless ((defined $i) && (defined $j));
    $object->assign($i,$j,$x);
    return $x;
}

sub element {
    my ($object, $i,$j) = @_;
    my $key = &makekey($i,$j);
    if (defined $object->{data}{$key}) {
	return $object->{data}{$key};
    } else {
	return 0;
    }
}

sub elementkey {
    my ($object, $key) = @_;
    if (defined $object->{data}{$key}) {
	return $object->{data}{$key};
    } else {
	return 0;
    }
}

sub elements {
    my ($matrix) = @_;
    return keys %{$matrix->{data}};
}


#returns row number $row of matrix as a row matrix
sub row {
    my ($matrix, $row,$persist) = @_;
    my $matrow = new Math::MatrixSparse;
    $matrow->{columns} = $matrix->{columns};
    $matrow->{rows} = 1;
    $matrow->name($matrix->{name});
    reuturn $matrow unless $row;
    if ($persist) {
	@{$matrix->{sortedrows}} = $matrix->sortbyrow()
	    unless defined $matrix->{sortedrows};
    }
    if (defined $matrix->{sortedrows}) {
	#binary search for proper values
	my @rows = @{$matrix->{sortedrows}};
	my ($left,$right) = (0,$#rows);
	my $mid = int(($right+$left)/2.0);
	while (
	       ((&splitkey($rows[$mid]))[0] != $row )
	       && ($right-$left>0)
	       ) {
	    if ((&splitkey($rows[$mid]))[0] < $row) {
		$left = $mid;
	    } else {
		$right = $mid;
	    }
	    $mid = int(($right+$left)/2.0);
	}
	return $matrow unless (&splitkey($rows[$mid]))[0]==$row;
	$right = $mid;
	while (  ($right<=$#rows)&&
		 ((&splitkey($rows[$right]))[0]==$row)
		 )
	{
	    $matrow->assign(1,
			    (&splitkey($rows[$right]))[1],
			    $matrix->elementkey($rows[$right]));
	    $right++;
	}
	$left = $mid-1;
	while (  ($left) &&
		 ((&splitkey($rows[$left]))[0]==$row)
		 ){
	    $matrow->assign(1,
			    (&splitkey($rows[$left]))[1],
			    $matrix->elementkey($rows[$left]));
	    $left--;
	}
	return $matrow;
    }
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matrow->assignkey($key, $matrix->{data}{$key}) if ($i==$row);
    }
    return $matrow;
}

#returns column number $col of matrix as a column matrix
sub column {
    my ($matrix, $col,$persist) = @_;
    my $matcol = new Math::MatrixSparse;
    $matcol->{columns} = $matrix->{columns};
    $matcol->{cols} = 1;
    $matcol->name($matrix->{name});
    if ($persist) {
	@{$matrix->{sortedcols}} = $matrix->sortbycolumn() unless defined $matrix->{sortedcols};
    }
    if (defined $matrix->{sortedcols}) {
	#binary search for proper values
	my @cols = @{$matrix->{sortedcols}};
	my ($left,$right) = (0,$#cols);
	my $mid = int(($right+$left)/2.0);
	while (
	       ((&splitkey($cols[$mid]))[1] != $col )
	       && ($right-$left>0)
	       ) {
	    if ((&splitkey($cols[$mid]))[1] < $col) {
		$left = $mid;
	    } else {
		$right = $mid;
	    }
	    $mid = int(($right+$left)/2.0);
	}
	return $matcol unless (&splitkey($cols[$mid]))[1]==$col;
	$right = $mid;
	while (  ($right<=$#cols)&&
		 ((&splitkey($cols[$right]))[1]==$col)
		 )
	{
	    $matcol->assign((&splitkey($cols[$right]))[0],1,
			    $matrix->elementkey($cols[$right]));
	    $right++;
	}
	$left = $mid-1;
	while (  ($left) &&
		 ((&splitkey($cols[$left]))[1]==$col)
		 ){
	    $matcol->assign((&splitkey($cols[$left]))[0],1,
			    $matrix->elementkey($cols[$left]));
	    $left--;
	}
	return $matcol;
    }
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matcol->assignkey($key, $matrix->{data}{$key}) if ($j==$col);
    }
    return $matcol;
}


### PRINT

sub print {
    my ($object,$name,$rc) = @_;
    return unless defined $object;
    my $label = $name ? $name : $object->{name};
    my @order ;
    $rc = "n" unless defined($rc);
    if ($rc =~ /^r/i) {
	@order = $object->sortbyrow();
    } elsif ($rc =~ /^c/i) {
	@order = $object->sortbycolumn();

    } else {
	@order = keys %{$object->{data}};
    }
    foreach my $key (@order){
	my ($i,$j) = &splitkey($key);
	print "$label($i, $j) = ", $object->{data}{$key},"\n";
    }
    return "";
}


###ARITHMETIC METHODS

#left+$right, dimensions must be identical
sub add {
    my ($left,$right,$switch) = @_;
    if (($left->{rows} == $right->{rows})&&($left->{columns} == $right->{columns})) {

	my $sum= $left->addfree($right);

	$sum->{rows} = $left->{rows};
	$sum->{columns} = $left->{columns};
	return $sum;
    } else {
	return undef;
    }
}

#as add, but no restrictions on dimensions.
sub addfree {
    my ($left, $right,$switch) = @_;
    my $sum = new Math::MatrixSparse;
    $sum = $left->copy();
    if ((defined $left->{name})&&(defined $right->{name})) {
	$sum->{name} = $left->{name} . "+" . $right->{name};
    }
    foreach my $rightkey (keys %{$right->{data}}) {
	if (defined $sum->{data}{$rightkey}) {
	    $sum->{data}{$rightkey}+= $right->{data}{$rightkey};
	} else {
	    $sum->{data}{$rightkey} = $right->{data}{$rightkey};
	}
    }
    $sum->{rows} = $left->{rows} >$right->{rows} ?
	$left->{rows} : $right->{rows};
    $sum->{columns} = $left->{columns} >$right->{columns} ?
	$left->{columns} : $right->{columns};
    if ($left->{special}->{structure} eq $right->{special}->{structure}) {
	$sum->{special}->{structure} = $left->{special}->{structure};
    }
    if ($left->{special}->{shape} eq $right->{special}->{shape}) {
	$sum->{special}->{shape} = $left->{special}->{shape};
    }
    return $sum;
}


sub subtract {
    my ($left, $right, $switch) = @_;
    if ($switch) {
	($left,$right) = ($right, $left);
    }
    if (
	($left->{rows} == $right->{rows})&&
	($left->{columns} == $right->{columns})
	) {
	my $diff= $left->subtractfree($right);
	$diff->{rows} = $left->{rows};
	$diff->{columns} = $left->{columns};
	return $diff;
    } else {
	return undef;
    }

}

#as subtract, but no restrictions on dimensions.
sub subtractfree {
    my ($left, $right,$switch) = @_;
    my $diff = new Math::MatrixSparse;
    $diff = $left->copy();
    foreach my $rightkey (keys %{$right->{data}}) {
	if (defined $diff->{data}{$rightkey}) {
	    $diff->{data}{$rightkey}-= $right->{data}{$rightkey};
	} else {
	    $diff->{data}{$rightkey} = -1* $right->{data}{$rightkey};
	}
    }
    if ((defined $left->{name})&&(defined $right->{name})) {
	$diff->{name} = $left->{name} . "-" . $right->{name};
    }
    $diff->{rows} = $left->{rows} >$right->{rows} ?
	$left->{rows} : $right->{rows};
    $diff->{columns} = $left->{columns} >$right->{columns} ?
	$left->{columns} : $right->{columns};
    if ($left->{special}->{structure} eq $right->{special}->{structure}) {
	$diff->{special}->{structure} = $left->{special}->{structure};
    }
    if ($left->{special}->{shape} eq $right->{special}->{shape}) {
	$diff->{special}->{shape} = $left->{special}->{shape};
    }
    return $diff;
}

sub negate {
    my ($matrix) = @_;
    return $matrix->multiplyscalar(-1);
}

sub _negate {
    my ($matrix) = @_;
    return $matrix->_multiplyscalar(-1);
}

sub multiplyscalar {
    my ($matrix, $scalar) = @_;
    my $product = $matrix->copy();
    $scalar = 0 unless $scalar;
    foreach my $key (keys %{$product->{data}}) {
	$product->assignkey($key,$matrix->elementkey($key)*$scalar);
    }
    $product->{name} = "$scalar*".$product->{name} if defined $product->{name};
    if ($scalar <0 ) {
	if ($matrix->{special}->{sign} =~ m/positive/i) {
	    $product->{special}->{sign} = "negative";
	} elsif ($matrix->{special}->{sign} =~ m/negative/i) {
	    $product->{special}->{sign} = "positive";
	} elsif ($matrix->{special}->{sign} =~ m/nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($matrix->{special}->{sign} =~ m/nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    }  else {
	$product->{special}->{sign} =$matrix->{special}->{sign};
    }
    $product->{special}->{sign} = "zero" unless $scalar;
    return $product;
}

sub _multiplyscalar {
    my ($matrix, $scalar) = @_;
    $scalar = 0 unless $scalar;
    foreach my $key (keys %{$matrix->{data}}) {
	$matrix->{data}{$key} = $matrix->{data}{$key}*$scalar;
    }
    if ($scalar <0 ) {
	if ($matrix->{special}->{sign} =~ m/positive/i) {
	    $matrix->{special}->{sign} = "negative";
	} elsif ($matrix->{special}->{sign} =~ m/negative/i) {
	    $matrix->{special}->{sign} = "positive";
	} elsif ($matrix->{special}->{sign} =~ m/nonpositive/i) {
	    $matrix->{special}->{sign} = "nonnegative";
	} elsif ($matrix->{special}->{sign} =~ m/nonnegative/i) {
	    $matrix->{special}->{sign} = "nonpositive";
	}
    }
    $matrix->{special}->{sign} = "zero" unless $scalar;
    return $matrix;
}


#finds $left*$right, if compatible
sub multiply {
    my ($left,$right,$switch) = @_;
    unless (ref($right)) {
	return $left->multiplyscalar($right);
    }
    return undef if ($left->{columns} != $right->{rows});
    my $product = new Math::MatrixSparse;
    $product->{rows} = $left->{rows};
    $product->{columns} = $right->{columns};
    if ((defined $left->{name})&&(defined $right->{name})) {
	$product->{name} = $left->{name} . "*" . $right->{name};
    }
    foreach my $leftkey (keys %{$left->{data}}) {
	my ($li,$lj) = &splitkey($leftkey);
	foreach my $rightkey (keys %{$right->{data}}) {
	    my ($ri,$rj) = &splitkey($rightkey);
	    next unless ($lj==$ri);
	    my $thiskey = &makekey($li, $rj);
	    my $prod = $left->{data}{$leftkey}*$right->{data}{$rightkey};
	    if (defined $product->{data}{$thiskey}) {
		$product->{data}{$thiskey} += $prod;
	    } else {
		$product->{data}{$thiskey} = $prod;
	    }

	}

    }
    if (
	($left->{special}->{sign} =~ /zero/i)||
	($right->{special}->{sign} =~ /zero/i)) {
	$product->{special}->{sign} = "zero";
	return $product;
    }
    if ($left->{special}->{sign} =~ /^positive/i) {
	$product->{special}->{sign} = $right->{special}->{sign};
    } elsif ($left->{special}->{sign} =~ /nonpositive/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /^negative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "negative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "positive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /nonnegative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonnegative";
	}
    }
    return $product;

}


sub quickmultiply {
    my ($left,$right,$switch) = @_;
    unless (ref($right)) {
	return $left->multiplyscalar($right);
    }
    return undef if ($left->{columns} != $right->{rows});
    my $product = new Math::MatrixSparse;
    $product->{special}->{structure} = "";
    $product->{special}->{pattern}=0;
    $product->{special}->{shape} = "";
    $product->{rows} = $left->{rows};
    $product->{columns} = $right->{columns};
    my @leftcols = $left->sortbycolumn();
    my @rightrows = $right->sortbyrow();
    if ((defined $left->{name})&&(defined $right->{name})) {
	$product->{name} = $left->{name} . "*" . $right->{name};
    }
    my $lastrow = 0;
    foreach my $leftkey (@leftcols) {
	my ($li,$lj) = &splitkey($leftkey);
	my $i = 0;
	my $thiskey;
	if ($lj >$lastrow ) {
	    $lastrow = $lj;
	    #remove elements that won't be used again in multiplication
	    while (defined ($thiskey = $rightrows[0])){
		my ($ri,$rj) = &splitkey($thiskey);
		last if $ri>=$lj;
		shift @rightrows;
	    }
	}
	foreach my $rightkey (@rightrows) {
	    my ($ri,$rj) = &splitkey($rightkey);
	    last if ($ri>$lj);
	    next if ($ri<$lj);
	    my $thiskey = &makekey($li, $rj);
	    my $prod = $left->{data}{$leftkey}*$right->{data}{$rightkey};
	    if (defined $product->{data}{$thiskey}) {
		$product->{data}{$thiskey} += $prod;
	    } else {
		$product->{data}{$thiskey} = $prod;
	    }

	}

    }
    if (
	($left->{special}->{sign} =~ /zero/i)||
	($right->{special}->{sign} =~ /zero/i)) {
	$product->{special}->{sign} = "zero";
	return $product;
    }
    if ($left->{special}->{sign} =~ /^positive/i) {
	$product->{special}->{sign} = $right->{special}->{sign};
    } elsif ($left->{special}->{sign} =~ /nonpositive/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /^negative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "negative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "positive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /nonnegative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonnegative";
	}
    }
    return $product;

}

sub quickmultiplyfree {
    my ($left,$right,$switch) = @_;
    unless (ref($right)) {
	return $left->multiplyscalar($right);
    }
    my $product = new Math::MatrixSparse;
    $product->{rows} = $left->{rows};
    $product->{columns} = $right->{columns};
    my @leftcols = $left->sortbycolumn();
    my @rightrows = $right->sortbyrow();
    if ((defined $left->{name})&&(defined $right->{name})) {
	$product->{name} = $left->{name} . "*" . $right->{name};
    }
    my $lastrow = 0;
    foreach my $leftkey (@leftcols) {
	my ($li,$lj) = &splitkey($leftkey);
	my $i = 0;
	my $thiskey;
	if ($lj >$lastrow ) {
	    $lastrow = $lj;
	    #remove elements that won't be used again in multiplication
	    while (defined ($thiskey = $rightrows[0])){
		my ($ri,$rj) = &splitkey($thiskey);
		last if $ri>=$lj;
		shift @rightrows;
	    }
	}
	foreach my $rightkey (@rightrows) {
	    my ($ri,$rj) = &splitkey($rightkey);
	    last if ($ri>$lj);
	    next if ($ri<$lj);
	    my $thiskey = &makekey($li , $rj);
	    my $prod = $left->{data}{$leftkey}*$right->{data}{$rightkey};
	    if (defined $product->{data}{$thiskey}) {
		$product->{data}{$thiskey} += $prod;
	    } else {
		$product->{data}{$thiskey} = $prod;
	    }

	}

    }
    if (
	($left->{special}->{sign} =~ /zero/i)||
	($right->{special}->{sign} =~ /zero/i)) {
	$product->{special}->{sign} = "zero";
	return $product;
    }
    if ($left->{special}->{sign} =~ /^positive/i) {
	$product->{special}->{sign} = $right->{special}->{sign};
    } elsif ($left->{special}->{sign} =~ /nonpositive/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /^negative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "negative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "positive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /nonnegative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonnegative";
	}
    }
    return $product;

}


#as multiply, but no restrictions on dimensions.
sub multiplyfree {
    my ($left,$right,$switch) = @_;
    my $product = new Math::MatrixSparse;
    $product->{rows} = $left->{rows};
    $product->{columns} = $right->{columns};
    if ((defined $left->{name})&&(defined $right->{name})) {
	$product->{name} = $left->{name} . "*" . $right->{name};
    }
    foreach my $leftkey (keys %{$left->{data}}) {
	my ($li,$lj) = &splitkey($leftkey);
	foreach my $rightkey (keys %{$right->{data}}) {
	    my ($ri,$rj) = &splitkey($rightkey);
	    next unless ($lj==$ri);
	    my $thiskey = &makekey($li, $rj);
	    my $prod = $left->{data}{$leftkey}*$right->{data}{$rightkey};
	    if (defined $product->{data}{$thiskey}) {
		$product->{data}{$thiskey} += $prod;
	    } else {
		$product->{data}{$thiskey} = $prod;
	    }

	}
    }
    if (
	($left->{special}->{sign} =~ /zero/i)||
	($right->{special}->{sign} =~ /zero/i)) {
	$product->{special}->{sign} = "zero";
	return $product;
    }
    if ($left->{special}->{sign} =~ /^positive/i) {
	$product->{special}->{sign} = $right->{special}->{sign};
    } elsif ($left->{special}->{sign} =~ /nonpositive/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /^negative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "negative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "positive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonpositive";
	}
    } elsif ($left->{special}->{sign} =~ /nonnegative/i) {
	if ($right->{special}->{sign} =~ /^positive/i) {
	    $product->{special}->{sign} = "nonnegative";
	} elsif ($right->{special}->{sign} =~ /nonpositive/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /^negative/i) {
	    $product->{special}->{sign} = "nonpositive";
	} elsif ($right->{special}->{sign} =~ /nonnegative/i) {
	    $product->{special}->{sign} = "nonnegative";
	}
    }
    return $product;

}

sub kronecker {
    my ($left,$right);
    my ($rr,$rc) = $right->size();
    return undef unless ($rr&&$rc);
    my $kroprod = new Math::MatrixSparse;
    foreach my $key (keys %{$left->{data}}) {
	my ($i,$j) = &splitkey($key);
	$kroprod->_insert($i*$rr,$j*$rc,$right*$left->elementkey($key));
    }
    if ((defined $left->{name})&&(defined $right->{name})) {
	$kroprod->{name} = $left->{name} . "x" . $right->{name};
    }
    return $kroprod;
}

sub termmult {
    my ($left,$right) = @_;
    my $termprod = $left->copy();
    foreach my $key (keys %{$termprod->{data}}) {
	$termprod->assignkey($key,$termprod->elementkey($key)* $right->elementkey($key));
    }
    return $termprod;

}

sub exponentiate {
    my ($matrix,$power)  = @_;
    unless ($matrix->is_square()) {
	carp "Math::MatrixSparse::exponentiate matrix must be square";
	return undef ;
    }
    return Math::MatrixSparse->newidentity($matrix->{rows}) unless $power;
    unless ($power>0) {
	carp "Math::MatrixSparse::exponentiate exponent must be positive";
	return undef ;
    }
    unless ($power =~ /^[+]?\d+/) {
	carp "Math::MatrixSparse::exponentiate exponent must be an integer";
	return undef ;
    }
    my $product = $matrix->copy();
#    $product->clearspecials();
    for my $i (2..$power) {
	$product = $product->quickmultiply($matrix);
    }
    $product->{name} = $matrix->{name} ."^$power" if  $matrix->{name};
    return $product;
}

sub largeexponentiate {
    #find matrix^power using the square-and-multiply method
    my ($matrix,$power)  = @_;
    unless ($matrix->is_square()) {
	carp "Math::MatrixSparse::exponentiate matrix must be square";
	return undef ;
    }
    return Math::MatrixSparse->newidentity($matrix->{rows}) unless $power;
    unless ($power>0) {
	carp "Math::MatrixSparse::exponentiate exponent must be positive";
	return undef ;
    }
    unless ($power =~ /^[+]?\d+/) {
	carp "Math::MatrixSparse::exponentiate exponent must be an integer";
	return undef ;
    }
    #get a representation of $exponent in binary
    my $bitstr = unpack('B32',pack('N',$power));
    $bitstr =~s/^0*//;
    my @bitstr=split(//,$bitstr);
    my $z = Math::MatrixSparse->newidentity($matrix->{rows});
    foreach my $bit (@bitstr){
        $z = ($z*$z);
        if ($bit){
            $z = ($z*$matrix);
        }
    }
    $z->{name} = "";
    $z->{name} = $matrix->{name} . "^$power" if  $matrix->{name};
    return $z;
}


sub transpose {
    my ($matrix) = @_;
    my $this= new Math::MatrixSparse;
    $this->{rows} = $matrix->{columns};
    $this->{columns} = $matrix->{rows};
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$this->assign($j,$i,$matrix->{data}{$key});
    }
    $this->{name} = $matrix->{name} . "'" if $matrix->{name};
    $this->{special}->{structure} = $matrix->{special}->{structure};
    $this->{special}->{sign} = $matrix->{special}->{sign};
    $this->{special}->{pattern} = $matrix->{special}->{pattern};
    $this->{special}->{square} = $matrix->{special}->{square};
    if ($matrix->{special}->{shape} =~ /diagonal/i) {
	$this->{special}->{shape}  = "diagonal";
    } elsif ($matrix->{special}->{shape} =~ /^lower/i) {
	$this->{special}->{shape}  = "upper";
    } elsif  ($matrix->{special}->{shape} =~ /^upper/i) {
	$this->{special}->{shape}  = "lower";
    } elsif  ($matrix->{special}->{shape} =~ /strictupper/i) {
	$this->{special}->{shape}  = "strictlower";
    } elsif  ($matrix->{special}->{shape} =~ /strictlower/i) {
	$this->{special}->{shape}  = "strictlower";
    } else {
	$this->{special}->{shape} = "";
    }
    return $this;
}


sub terminvert {
    my ($matrix) = @_;
    my $this = $matrix->copy();
    foreach my $key (keys %{$this->{data}}) {
	next unless $this->{data}{$key};
	$this->{data}{$key} = 1.0/($this->{data}{$key});
    }
    return $this;
}

sub _terminvert {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}) {
	next unless $matrix->{data}{$key};
	$matrix->{data}{$key} = 1.0/($matrix->{data}{$key});
    }
    return $matrix;
}



### DISSECTION METHODS
sub diagpart {
    my ($matrix,$offset) = @_;
    my $diag = new Math::MatrixSparse;
    $offset = 0 unless defined $offset;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	next unless ($i == ($j+$offset));
	$diag->assign($i,$j,$matrix->{data}{$key});
    }
    $diag->{rows} = $matrix->{rows};
    $diag->{columns} = $matrix->{columns};
    $diag->{special}->{shape} = "diagonal";
    return $diag;
}

sub nondiagpart {
    my ($matrix,$offset) = @_;
    my $diag = new Math::MatrixSparse;
    $offset = 0 unless defined $offset;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	next if ($i == ($j+$offset));
	$diag->assign($i,$j,$matrix->{data}{$key});
    }
    $diag->{rows} = $matrix->{rows};
    $diag->{columns} = $matrix->{columns};
    return $diag;
}

sub lowerpart {
    my ($matrix) = @_;
    my $lower = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	next unless ($i > $j);
	$lower->assign($i,$j,$matrix->{data}{$key});
    }
    $lower->{rows} = $matrix->{rows};
    $lower->{columns} = $matrix->{columns};
    $lower->{special}->{shape} = "strictlower";

    return $lower;
}

sub nonlowerpart {
    my ($matrix) = @_;
    my $lower = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	next if ($i > $j);
	$lower->assign($i,$j,$matrix->{data}{$key});
    }
    $lower->{rows} = $matrix->{rows};
    $lower->{columns} = $matrix->{columns};
    $lower->{special}->{shape} = "upper";
    return $lower;
}

sub upperpart {
    my ($matrix) = @_;
    my $upper = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	next unless ($i < $j);
	$upper->assign($i,$j,$matrix->{data}{$key});
    }
    $upper->{rows} = $matrix->{rows};
    $upper->{columns} = $matrix->{columns};
    $upper->{special}->{shape} = "strictupper";
    return $upper;
}

sub nonupperpart {
    my ($matrix) = @_;
    my $upper = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	next if ($i < $j);
	$upper->assign($i,$j,$matrix->{data}{$key});
    }
    $upper->{rows} = $matrix->{rows};
    $upper->{columns} = $matrix->{columns};
    $upper->{special}->{shape} = "lower";

    return $upper;
}


sub _diagpart {
    my ($matrix,$offset) = @_;
    $offset = 0 unless defined $offset;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matrix->delete($i,$j) unless ($i == ($j+$offset));
    }
    $matrix->{special}->{shape} = "diagonal";
    return $matrix;
}

sub _nondiagpart {
    my ($matrix,$offset) = @_;
    $offset = 0 unless defined $offset;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matrix->delete($i,$j) if ($i == ($j+$offset));
    }
    $matrix->{special}->{shape}="" if
	$matrix->{special}->{shape} =~ m/diagonal/i;
    return $matrix;
}


sub _lowerpart {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matrix->delete($i,$j) unless ($i > $j);
    }
    $matrix->{special}->{shape} = "strictlower";
    return $matrix;
}

sub _nonlowerpart {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matrix->delete($i,$j) if ($i > $j);
    }
    $matrix->{special}->{shape} = "upper";
    return $matrix;
}

sub _upperpart {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matrix->delete($i,$j) unless ($i < $j);
    }
    $matrix->{special}->{shape} = "strictupper";
    return $matrix;
}

sub _nonupperpart {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}){
	my ($i,$j) = &splitkey($key);
	$matrix->delete($i,$j) if ($i < $j);
    }
    $matrix->{special}->{shape} = "lower";
    return $matrix;
}


sub symmetricpart {
    #.5*( A+A' )
    my ($matrix) = @_;
    my $sp =($matrix+$matrix->transpose());
    $sp = 0.5*$sp;
    $sp->{special}->{structure} = "symmetric";
    return $sp;

}
sub _symmetricpart {
    #.5*( A+A' )
    my ($matrix) = @_;
    my $sp =($matrix+$matrix->transpose());
    $sp = 0.5*$sp;
    $sp->{special}->{structure} = "symmetric";
    return $matrix=$sp->copy();
}

sub skewsymmetricpart {
    #.5*( A-A' )
    my ($matrix) = @_;
    my $ssp= 0.5*($matrix-$matrix->transpose());
    $ssp->{special}->{structure} = "skewsymmetric";
    return $ssp;
}

sub _skewsymmetricpart {
    #.5*( A-A' )
    my ($matrix) = @_;
    my $ssp= 0.5*($matrix-$matrix->transpose());
    $ssp->{special}->{structure} = "skewsymmetric";
    return $matrix=$ssp->copy();
}


sub positivepart {
    my ($matrix) = @_;
    my $pos = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}){
	next unless $matrix->elementkey($key) >0;
	$pos->assignkey($key,$matrix->{data}{$key});
    }
    $pos->{rows} = $matrix->{rows};
    $pos->{columns} = $matrix->{columns};
    return $pos;
}


sub negativepart {
    my ($matrix) = @_;
    my $neg = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}){
	next unless $matrix->elementkey($key) <0;
	$neg->assignkey($key,$matrix->{data}{$key});
    }
    $neg->{rows} = $matrix->{rows};
    $neg->{columns} = $matrix->{columns};
    return $neg;
}

sub _positivepart {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}){
	next if $matrix->elementkey($key) >0;
	$matrix->deletekey($key,$matrix->{data}{$key});
    }
    return $matrix;
}

sub _negativepart {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}){
	next if $matrix->elementkey($key) <0;
	$matrix->deletekey($key,$matrix->{data}{$key});
    }
    return $matrix;
}


sub mask{
    my ($matrix, $i1,$i2,$j1,$j2) = @_;
    return undef  unless (($i1<=$i2)&&($j1<=$j2));
    my $mask = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	next unless (($i>=$i1)&&($i<=$i2));
	next unless (($j>=$j1)&&($j<=$j2));
	$mask->assignkey($key,$matrix->{data}{$key});
    }
    return $mask;
}

sub _mask{
    my ($matrix, $i1,$i2,$j1,$j2) = @_;
    return undef  unless (($i1<=$i2)&&($j1<=$j2));
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	next if (($i>=$i1)&&($i<=$i2));
	next if (($j>=$j1)&&($j<=$j2));
	$matrix->assignkey($key,0);
    }
    return $matrix;
}

sub submatrix{
    my ($matrix, $i1,$i2,$j1,$j2) = @_;
    return undef  unless (($i1<=$i2)&&($j1<=$j2));
    my $subm = new Math::MatrixSparse;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	next unless (($i>=$i1)&&($i<=$i2));
	next unless (($j>=$j1)&&($j<=$j2));
	$subm->assign($i-$i1,$j-$j1,$matrix->{data}{$key});
    }
    return $subm;
}

sub insert {
    my ($big, $i0,$j0,$small) = @_;
    my $insert = $big->copy();
    foreach my $key (keys %{$small->{data}}){
	my ($i,$j) = &splitkey($key);
	$insert->assignkey($i+$i0,$j+$j0,$small->elementkey($key));
    }
    return $insert;
}

sub _insert {
    my ($big, $i0,$j0,$small) = @_;

    foreach my $key (keys %{$small->{data}}){
	my ($i,$j) = &splitkey($key);
	$big->assignkey($i+$i0,$j+$j0,$small->elementkey($key));
    }
    return $big;
}


sub shift {
    my ($matrix, $i1,$j1) = @_;
    my $this = new Math::MatrixSparse;
    $this->{name}=$matrix->{name};
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$this->assign($i+$i1,$j+$j1,$matrix->{data}{$key});
    }
    return $this;
}

sub each {
    my ($matrix,$coderef) = @_;
    my $this = $matrix->copy();
    foreach my $key (keys %{$this->{data}}) {
	my ($i,$j) = &splitkey($key);
	$this->assign($i,$j,&$coderef($this->{data}{$key},$i,$j));
    }
    return $this;
}

sub _each {
    my ($matrix,$coderef) = @_;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$matrix->assign($i,$j,&$coderef($matrix->{data}{$key},$i,$j));
    }
    return $matrix;
}

### INFORMATIONAL AND STATISTICAL METHODS


sub printstats {
    my ($matrix,$name) = @_;
    return unless defined $matrix;

    $name = $name || $matrix->{name} || "unnamed matrix";
    print "Statistics for $name :\n";
    print "rows: $matrix->{rows}\t";
    print "columns: $matrix->{columns}\tdefined elements: ";
    print $matrix->count(), "\n";
    my ($width) = $matrix->width();
    print "Bandwidth: $width\t";
    print "Symmetric " if $matrix->is_symmetric();
    print "Skew-Symmetric " if $matrix->is_skewsymmetric();
#    print "Hermetian " if $matrix->{special}->{structure} =~m/^hermetian/i;
    print "Real " if $matrix->{special}->{field} =~ m/real/i;
#    print "Complex " if $matrix->{special}->{field} =~ m/complex/i;
    print "Strictly " if $matrix->{special}->{shape} =~ m/strict/;
    print "Upper Triangular " if $matrix->{special}->{shape} =~ m/upper/;
    print "Lower Triangular " if $matrix->{special}->{shape} =~ m/lower/;
    print "Diagonal" if $matrix->is_diagonal();
    print "Pattern " if $matrix->is_pattern();
    print "Square " if $matrix->is_square();
    print "\n";
}


sub dim {
    my ($matrix) = @_;
    return ($matrix->{rows},$matrix->{columns});
}

sub exactdim {
    my ($matrix) = @_;
    my $rows=0;
    my $columns = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$rows = $i if $i>$rows;
	$columns = $j if $j > $rows;
    }
    return ($rows,$columns);
}


sub count {
    my ($matrix) = @_;
    return scalar keys %{$matrix->{data}};
}

sub width {
    my ($matrix) = @_;
    return $matrix->{special}->{bandwidth} if $matrix->{special}->{bandwidth};
    my $width = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$width = abs($i- $j) if abs($i-$j) > $width;
    }
    return ($width);
}

sub sum {
    my ($matrix) = @_;
    my $sum = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	$sum += $matrix->elementkey($key);
    }
    return $sum;
}

sub sumeach {
    my ($matrix,$coderef) = @_;
    my $sum = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey( $key);
	$sum += &$coderef($matrix->elementkey($key),$i,$j);
    }
    return $sum;
}

sub abssum {
    my ($matrix) = @_;
    my $sum = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	$sum += abs($matrix->elementkey($key));
    }
    return $sum;
}

sub rownorm {
    my ($matrix) = @_;
    my %rowsums;
    return  0 unless defined $matrix;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$rowsums{$i}+= abs($matrix->elementkey($key));
    }
    return (sort {$a <=> $b} values %rowsums)[0];
}


sub columnnorm {
    my ($matrix) = @_;
    my %columnsums;
    return 0 unless defined $matrix;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$columnsums{$i}+= $matrix->elementkey($key);
    }
    return (sort {$a <=> $b} values %columnsums)[0];
}

sub norm_max {
    return $_[0]->rownorm();
}

sub norm_one {
    return $_[0]->columnnorm();
}

sub norm_frobenius {
    return sqrt($_[0]->sumeach(sub {$_[0]*$_[0]}));
}


sub trace {
    my ($matrix) = @_;
    return $matrix->diagpart()->sum();
}

### BOOLEAN METHODS

sub equals {
    my ($left,$right) = @_;
    return 1 unless ( defined $left|| defined $right);
    return 0 unless  defined $left;
    return 0 unless  defined $right;
    my $truth = 1;
    foreach my $key (keys %{$left->{data}}, keys %{$right->{data}}) {
	$truth *= ($left->elementkey($key) == $right->elementkey($key));
    }
    return $truth;
}

sub is_square {
    my ($matrix)  = @_;
    return 0 unless  defined $matrix;
    return $matrix->{rows} == $matrix->{columns};
}
sub is_quadratic {
    my ($matrix)  = @_;
    return 0 unless  defined $matrix;
    return $matrix->{rows} == $matrix->{columns};
}

sub is_symmetric {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    return 1 if $matrix->{special}->{structure} =~ m/^symmetric/;
    return 0 if $matrix->{special}->{structure} =~ m/skewsymmetric/;
#    return 0 if $matrix->{special}->{structure} =~ m/hermetian/;

    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	my $reversekey = &makekey($j, $i);
	$truth *= ($matrix->elementkey($key)
		   == $matrix->elementkey($reversekey));
	return 0 unless $truth;
    }
    return $truth;
}

sub is_skewsymmetric {
    my ($matrix) = @_;
    return 0 unless  defined $matrix;
    return 0 if $matrix->{special}->{structure} =~ m/^symmetric/;
    return 1 if $matrix->{special}->{structure} =~ m/^skewsymmetric/;
    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	my $reversekey = &makekey($j , $i);
	$truth *= ($matrix->elementkey($key)
		   == -1*$matrix->elementkey($reversekey));
	return 0 unless $truth;
    }
    return $truth;
}

sub is_diagonal {
    my ($matrix) = @_;
    return 0 unless  defined $matrix;
    return 1 if $matrix->{special}->{shape}=~m/diagonal/i;
    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$truth *= ($i==$j);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_strictlowertriangular {
    my ($matrix) = @_;
    return 0 unless  defined $matrix;
    return 1 if $matrix->{special}->{shape}=~m/strictlower/i;
    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$truth *= ($i > $j);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_strictuppertriangular {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    return 1 if $matrix->{special}->{shape}=~m/strictupper/i;
    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$truth *= ($i < $j);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_lowertriangular {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    my $truth = 1;
    return 1 if $matrix->{special}->{shape}=~m/lower/i;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$truth *= ($i >= $j);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_uppertriangular {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    my $truth = 1;
    return 1 if $matrix->{special}->{shape}=~m/upper/i;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$truth *= ($i <= $j);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_positive {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    return 1 if $matrix->{special}->{sign} =~ m/^positive/i;
    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	$truth *= ($matrix->elementkey($key)>0);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_zero {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    my $truth = 1;
    return 1 if $matrix->{special}->{sign} =~ /zero/i;
    foreach my $key (keys %{$matrix->{data}}) {
	$truth *= ($matrix->elementkey($key)==0);
	return 0 unless $truth;
    }
    return $truth;
}


sub is_nonpositive {
    my ($matrix) = @_;
    my $truth = 1;
    return 0 unless defined $matrix;
    foreach my $key (keys %{$matrix->{data}}) {
	$truth *= ($matrix->elementkey($key)<=0);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_negative {
    my ($matrix) = @_;
    my $truth = 1;
    return 0 unless defined $matrix;
    foreach my $key (keys %{$matrix->{data}}) {
	$truth *= ($matrix->elementkey($key)<0);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_nonnegative {
    my ($matrix) = @_;
    my $truth = 1;
    return 0 unless defined $matrix;
    foreach my $key (keys %{$matrix->{data}}) {
	$truth *= ($matrix->elementkey($key)>=0);
	return 0 unless $truth;
    }
    return $truth;
}

sub is_boolean {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    return 1 if $matrix->{special}->{pattern};
    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	$truth *= (($matrix->elementkey($key) == 0)||($matrix->elementkey($key) == 1));
	return 0 unless $truth;
    }
    return $truth;
}

sub is_pattern {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    return 1 if $matrix->{special}->{pattern};
    my $truth = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	$truth *= (($matrix->elementkey($key) == 0)||($matrix->elementkey($key) == 1));
	return 0 unless $truth;
    }
    return $truth;
}


sub diagnose {
    my ($matrix) = @_;
    return 0 unless defined $matrix;
    my $boolean = 1;
    $matrix->_sizetofit();
    $matrix->{special}->{square}  = ($matrix->{rows} == $matrix->{columns});
    my $upper = 1;
    my $diagonal = 1;
    my $lower = 1;
    my $strictupper = 1;
    my $strictlower = 1;
    my $positive = 1;
    my $nonpositive = 1;
    my $negative = 1;
    my $nonnegative = 1;
    my $zero = 1;
    my $symmetric = 1;
    my $skewsymmetric = 1;
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	my $reversekey = &makekey($j,$i);
        $symmetric *= ($matrix->elementkey($key)
		       == $matrix->elementkey($reversekey));
        $skewsymmetric *= ($matrix->elementkey($key)
			   == -1*$matrix->elementkey($reversekey));
	$diagonal *= ($i==$j);
	$strictlower *= ($i > $j);
	$lower *= ($i >= $j);
	$strictupper *= ($i > $j);
	$upper *= ($i >= $j);
	$positive *= ($matrix->elementkey($key)>0);
	$nonpositive *= ($matrix->elementkey($key)<=0);
	$negative *= ($matrix->elementkey($key)<0);
	$nonnegative *= ($matrix->elementkey($key)>=0);
	$boolean *= (($matrix->elementkey($key) == 0)
		     ||($matrix->elementkey($key) == 1));
	$zero *= !($matrix->elementkey($key));
    }
    $matrix->{special}->{pattern} = $boolean;
    if ($diagonal) {
	$matrix->{special}->{shape} = "diagonal";
    } elsif ($lower) {
	if ($strictlower) {
	    $matrix->{special}->{shape} = "strictlower";
	} else {
	    $matrix->{special}->{shape} = "lower";
	}
    } elsif ($upper) {
	if ($strictlower) {
	    $matrix->{special}->{shape} = "strictupper";
	} else {
	    $matrix->{special}->{shape} = "upper";
	}
    } else {
	$matrix->{special}->{shape}="";
    }
    if ($symmetric) {
	$matrix->{special}->{structure}="symmetric";
    } elsif ($skewsymmetric) {
	$matrix->{special}->{structure}="skewsymmetric";
    } else {
	 $matrix->{special}->{structure}="";
    }
    if ($zero) {
	$matrix->{special}->{sign} = "zero";
    } elsif ($nonpositive) {
	if ($negative) {
	    $matrix->{special}->{sign}="negative";
	} else {
	    $matrix->{special}->{sign}="nonpositive";
	}
    } elsif ($nonnegative) {
	if ($positive) {
	    $matrix->{special}->{sign}="positive";
	} else {
	    $matrix->{special}->{sign}="nonnegative";
	}
    } else {
	$matrix->{special}->{sign}="";
    }
    return $matrix;
}


### BOOLEAN ARITHMETIC METHODS

sub matrixand {
    my ($matrix1, $matrix2) = @_;
    my $truth = new Math::MatrixSparse;
    $truth->{rows} = ($matrix1->{rows}<$matrix2->{rows}) ? $matrix1->{rows} : $matrix2->{rows};
    $truth->{columns} = ($matrix1->{columns}<$matrix2->{columns}) ? $matrix1->{columns} : $matrix2->{columns};
    foreach my $key (keys %{$matrix1->{data}}) {
	$truth->assignkey($key,$matrix1->elementkey($key)&&
			  $matrix2->elementkey($key));
    }
    if ((defined $matrix1->{name})&&(defined $matrix2->{name})) {
	$truth->{name} = $matrix1->{name} . "&" . $matrix2->{name};
    }
    #should be updated already.
    $truth->{special}->{pattern} = 1;
    return $truth;
}

sub matrixor {
    my ($matrix1, $matrix2) = @_;
    my $truth = new Math::MatrixSparse;
    $truth->{rows} = ($matrix1->{rows}<$matrix2->{rows}) ? $matrix1->{rows} : $matrix2->{rows};
    $truth->{columns} = ($matrix1->{columns}<$matrix2->{columns}) ? $matrix1->{columns} : $matrix2->{columns};
    foreach my $key (keys %{$matrix1->{data}}) {
	$truth->assignkey($key,$matrix1->elementkey($key)||
			  $matrix2->elementkey($key));
    }
    if ((defined $matrix1->{name})&&(defined $matrix2->{name})) {
	$truth->{name} = $matrix1->{name} . "|" . $matrix2->{name};
    }
    #should be updated already.
    $truth->{special}->{pattern} = 1;
    return $truth;
}


### DELETION FUNCTIONS
sub delete {
    my ($matrix,$i,$j) = @_;
    my $key = &makekey($i,$j);
    $matrix->deletekey($key);
    return;
}

sub deletekey {
    my ($matrix, $key) = @_;
    my $old = $matrix->elementkey($key);
    my ($i,$j) = &splitkey($key);
    delete $matrix->{data}{$key};
    #sign
    if ($matrix->{special}->{sign} =~ /^positive/i) {
	$matrix->{special}->{sign} = "nonnegative";
    } elsif ($matrix->{special}->{sign} =~ /^positive/i) {
	$matrix->{special}->{sign} = "nonpositive";
    }
    #structure
    if (
	($matrix->{special}->{structure} =~ m/symmetric/i)
	&&($i!=$j)
	) {
	$matrix->{special}->{structure} = "";
    }
    #band
    if (abs($i-$j) >= $matrix->{special}->{bandwidth}) {
	$matrix->{special}->{bandwidth} = 0;
    }
    #pattern--no change necessary
    #shape--no change necessary
    #persistent row and column data
    delete $matrix->{sortedrows};
    delete $matrix->{sortedcolumns};
    return undef;
}

sub cull {
    my ($matrix) = @_;
    my $this = $matrix->copy();
    foreach my $key (keys %{$this->{data}}) {
	next if $this->{data}{$key};
	$this->deletekey($key);
    }
    return $this;
}

sub _cull {
    my ($matrix) = @_;
    foreach my $key (keys %{$matrix->{data}}) {
	next if $matrix->{data}{$key};
	$matrix->deletekey($key);
    }
    return $matrix;
}

sub threshold {
    my ($matrix,$thresh) = @_;
    return undef if $thresh <0;
    my $this = $matrix->copy();
    foreach my $key (keys %{$this->{data}}) {
	next if abs($this->{data}{$key})>$thresh;
	$this->deletekey($key);
    }
    return $this;
}

sub _threshold {
    my ($matrix,$thresh) = @_;
    return undef if $thresh <0;
    foreach my $key (keys %{$matrix->{data}}) {
	next if abs($matrix->{data}{$key})>$thresh;
	$matrix->deletekey($key);
    }
    return $matrix;
}

sub sizetofit {
    my ($matrix) = @_;
    my $this = $matrix->copy();
    my ($maxrow,$maxcol) = $this->exactdim();
    $this->{rows} = $maxrow;
    $this->{columns} = $maxcol;
    return $this;
}


sub _sizetofit {
    my ($matrix) = @_;
    my ($maxrow,$maxcol) = $matrix->exactdim();
    $matrix->{rows} = $maxrow;
    $matrix->{columns} = $maxcol;
    return $matrix;
}


sub clearspecials {
    my ($matrix) = @_;
    $matrix->{special}->{pattern} = 0;
    $matrix->{special}->{sign} = "";
    $matrix->{special}->{structure} = "";
    $matrix->{special}->{shape} = "";
}

### SOLVER METHODS

sub gaussseidel {
    #solves Ax=b by Gauss-Seidel, or SOR with relaxation parameter 1.
    #b ($constant) and x0 ($guess) should be column vectors.
    my ($matrix, $constant, $guess,$tol, $steps) = @_;
    return $matrix->SOR($constant,$guess,1,$tol,$steps);
}

sub SOR {
    #solves Ax=b by Symmetric Over-Relaxation
    #b ($constant) and x0 ($guess) should be column vectors.
    my ($matrix, $constant, $guess, $omega,$tol, $steps) = @_;
    my $diag = $matrix->diagpart();
    my $lower = $matrix->lowerpart();
    my $upper = $matrix->upperpart();
    my $iterator;
    my $dinv = $diag->terminvert();
    my $soln = $guess->copy();
    foreach my $key (keys %{$constant->{data}}) {
	next  if defined $soln->{data}{$key};
	$soln->{data}{$key} = 0;
    }
    my @lowerkeys = $lower->sortbyrow();
    my @upperkeys = $upper->sortbyrow();
    my @solnkeys = $soln->sortbyrow();
    $steps = 100 unless $steps;
    for my $k (1 .. $steps) {
	my $oldsoln = $soln->copy();
	foreach my $thissolnkey (keys %{$soln->{data}}) {
	    my ($i,$j) = &splitkey($thissolnkey);
	    next unless $j == 1;
	    my $prev = $oldsoln->{data}{$thissolnkey}*(1.0-$omega);
	    my $scal = $omega * $dinv->element($i,$i);
	    my $lowersum=0;
	    my $uppersum = 0;
	    foreach my $lowerkey (@lowerkeys) {
		my ($li,$lj) = &splitkey($lowerkey);
		#unnecessary b/c of source of @lowerkeys
		next if $lj > ($i-1);
#		print "L $li $lj\n";
		$lowersum += $soln->element($lj,1)*$lower->{data}{$lowerkey};
	    }
	    foreach my $upperkey (@upperkeys) {
		my ($ui,$uj) = &splitkey($upperkey);
		#unnecessary b/c of source of @upperkeys
		next if $uj < ($i+1);
		$uppersum += $soln->element($uj,1)*$upper->{data}{$upperkey};
	    }
	    my $update = $prev+$scal*($constant->element($i,1)-$lowersum-$uppersum);
	    $soln->assign($i,1,$update);
	}
	my $err = &abssum($oldsoln-$soln);
	return $soln if $err<$tol;
    }
    return undef;
}


sub jacobi {
    #solves Ax=b by Jacobi iteration.
    #Note that b doesn't have to be a column vector
    #$guess (x0) should be the same size as $constant (b)
    my ($matrix, $constant, $guess,$tol, $steps) = @_;
    my $diag = $matrix->diagpart();
    my $lower = $matrix->lowerpart();
    my $upper = $matrix->upperpart();
    my $iterator;
    my $dinv = $diag->terminvert();
    my $nondiag = ($lower+$upper);
    my $soln = $guess->copy();
    my $oldsoln;
    $steps = 100 unless $steps;
    for my $i (1 .. $steps) {
	$oldsoln=$soln->copy();;
	$soln = $dinv*$constant - $dinv*$nondiag*$soln;
	my $err = &abssum($oldsoln-$soln);
	return $soln if $err<$tol;
    }
    return undef;
}


### SORTING METHODS

#note: column is a secondary sort criterion
# 1 2 3
# 4 5 6
# 7 8 9
sub sortbyrow {
    my ($matrix) = @_;
    my @sorted = map  { $_->[0] }
    sort { ($a->[1] <=> $b->[1])||($a->[2]<=>$b->[2]) }
    map  { [ $_, &splitkey($_) ] } keys %{$matrix->{data}};
    return @sorted;
}

#note: row is a secondary sort criterion
# 1 4 7
# 2 5 8
# 3 6 9
sub sortbycolumn {
    my ($matrix) = @_;
    my @sorted = map  { $_->[0] }
    sort { ($a->[2] <=> $b->[2])||($a->[1]<=>$b->[1]) }
    map  { [ $_, &splitkey($_) ] } keys %{$matrix->{data}};

    return @sorted;
}

#note: row is a secondary sort criterion
#lower diagonals are sorted in front of higher ones
#4 7 9
#2 5 8
#1 3 6
sub sortbydiagonal {
    my ($matrix) = @_;
    my @sorted = map  { $_->[0] }
    sort { (($a->[1]-$a->[2]) <=> ($b->[1]-$b->[2]))||($a->[1]<=>$b->[1]) }
    map  { [ $_, &splitkey($_) ] } keys %{$matrix->{data}};

    return @sorted;
}

#sorts the elements of $matrix by value, smallest first.
#row is a secondary criterion, and column is tertiaty.
sub sortbyvalue {
    my ($matrix) = @_;
    my @sorted = map  { $_->[0] }
    sort {
	($matrix->elementkey($a->[0]) <=> $matrix->elementkey($a->[0]))
	    || ($a->[1]<=>$b->[1])
		||($a->[2]<=>$b->[2])
	}
    map  { [ $_, &splitkey($_) ] } keys %{$matrix->{data}};

    return @sorted;
}


### SYMMETRIC METHODS
sub symmetrify {
    #takes a matrix, returns the matrix obtained by reflecting
    #non-lower part around main diagonal
    my ($matrix) = @_;
    return $matrix if $matrix->{special}->{shape} =~ /diagonal/i;
    my $this = $matrix->copy();
    $this->_symmetrify();
    return $this;
}


sub _symmetrify {
    #takes a matrix, reflects it about main diagonal
    my ($matrix) = @_;
    return $matrix if $matrix->{special}->{shape} =~ /diagonal/i;
    $matrix->_nonlowerpart();
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	unless ($i==$j) {
	    my $transkey = &makekey($j , $i);
	    $matrix->assign($j,$i,$matrix->element($i,$j));
	}
    }
    $matrix->_sizetofit();
    $matrix->{special}->{shape} = "";
    $matrix->{special}->{structure} = "symmetric";
    return $matrix;
}

sub skewsymmetrify {
    #takes a matrix, returns the matrix obtained by reflecting
    #upper part around main diagonal
    my ($matrix) = @_;
    my $this = $matrix->copy();
    return $this->_skewsymmetrify();
}

sub _skewsymmetrify {
    #takes a matrix, reflects upper part about main diagonal
    my ($matrix) = @_;
    $matrix->_upperpart();
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	my $transkey = &makekey($j , $i);
	$matrix->assign($j,$i,-1*$matrix->element($i,$j));
    }
    $matrix->{special}->{shape} = "";
    $matrix->{special}->{structure} = "skewsymmetric";
    return $matrix;
}


### PROBABILISTIC METHODS
sub normalize {
    my ($matrix) = @_;
    return undef unless defined $matrix;
    my $name = $matrix->{name};
    unless ($matrix->is_nonnegative()) {
	carp "Math::MatrixSparse::normalize matrix has negative elements";
	return undef;
    }
    my  $matsum = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	$matsum += $matrix->elementkey($key);
    }
    $matrix= $matrix->multiplyscalar(1.0/$matsum) if $matsum;
    $matrix->name($name . "/||" . $name . "||") if ($name);
    return $matrix;
}

sub _normalize {

    my ($matrix) = @_;
    return undef unless defined $matrix;
    my $name = $matrix->{name};
    unless ($matrix->is_nonnegative()) {
	carp "Math::MatrixSparse::normalize matrix has negative elements";
	return undef;
    }
    my  $matsum = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	$matsum += $matrix->elementkey($key);
    }
    $matrix= $matrix->_multiplyscalar(1.0/$matsum) if $matsum;
    $matrix->name($name . "/||" . $name . "||") if ($name);
    return $matrix;
}


sub normalizerows {
    my ($matrix) = @_;
    my %rowsums;
    return undef unless defined $matrix;
    unless ($matrix->is_nonnegative()) {
	carp "Math::MatrixSparse::normalizerows matrix has negative elements";
	return undef;
    }
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$rowsums{$i}+= $matrix->elementkey($key);
    }
    my $rownormed = $matrix->cull();
    foreach my $key (keys %{$rownormed->{data}}) {
	my ($i,$j) = &splitkey($key);
	next unless $rowsums{$i};
	$rownormed->assign($i,$j,($rownormed->elementkey($key))/$rowsums{$i});
    }
    my $name = $matrix->{name};
    $rownormed->name($name . "/||" . $name . "||") if ($name);
    return $rownormed;
}



sub _normalizerows {
    my ($matrix) = @_;
    my %rowsums;
    return undef unless defined $matrix;
    unless ($matrix->is_nonnegative()) {
	carp "Math::MatrixSparse::normalizerows matrix has negative elements";
	return undef;
    }
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$rowsums{$i}+= $matrix->{data}{$key};
    }
    $matrix->_cull();
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	next unless $rowsums{$i};
	$matrix->assign($i,$j,($matrix->elementkey($key))/$rowsums{$i});
    }
    my $name = $matrix->{name};
    $matrix->name($name . "/||" . $name . "||") if ($name);
    return $matrix;
}

sub normalizecolumns {
    my ($matrix) = @_;
    my %columnsums;
    return undef unless defined $matrix;
    unless ($matrix->is_nonnegative()) {
	carp "Math::MatrixSparse::normalizecolumns matrix has negative elements";
	return undef;
    }
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$columnsums{$j}+= $matrix->elementkey($key);
    }
    my $columnnormed = $matrix->copy();
    foreach my $key (keys %{$columnnormed->{data}}) {
	my ($i,$j) = &splitkey($key);
	$columnnormed->assign($i,$j,$columnnormed->elementkey($key)/$columnsums{$j});
    }
    my $name = $matrix->{name};
    $columnnormed->name($name . "/||" . $name . "||") if ($name);
    return $columnnormed;
}


sub _normalizecolumns {
    my ($matrix) = @_;
    my %columnsums;
    return undef unless defined $matrix;
    unless ($matrix->is_nonnegative()) {
	carp "Math::MatrixSparse::normalizecolumns matrix has negative elements";
	return undef;
    }
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$columnsums{$j}+= $matrix->elementkey($key);
    }
    $matrix->_cull();
    foreach my $key (keys %{$matrix->{data}}) {
	my ($i,$j) = &splitkey($key);
	$matrix->assign($i,$j,$matrix->elementkey($key)/$columnsums{$j});
    }
    my $name = $matrix->{name};
    $matrix->name($name . "/||" . $name . "||") if ($name);
    return $matrix;
}


sub discretepdf {
    my ($matrix) = @_;
    my $uniform = rand;

    my $current;
    my $curpos = 0;
    foreach my $key (keys %{$matrix->{data}}) {
	last if $curpos >=$uniform;
	$current = $key;
	$curpos += $matrix->{data}{$key};
    }
    return $current;

}

### KEY-INTERFACE METHODS

sub splitkey {
    my ($key) = @_;
    my ($i,$j) = split(/\0/,$key);
    return ($i,$j);
}

sub makekey {
    my ($i,$j) = @_;
    return ($i . "\0" . $j) if ((defined $i)&&(defined $j));

}


1;
__END__