/MatrixSparse.pm

Large files files are truncated, but you can click here to view the full file

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;
    $_ = <HB>;
    chomp();
    my ($ident, $key) = $_ =~ m/^(.*)(........)\s*$/;
    $key =~ s/\s*$//;
    $this->{name} = $key;
    $_ = <HB>;
    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*$/;
    $_ = <HB>;
    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 = <HB>;
    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 = <HB>;
	chomp($rhsspec);
    }
    #now read the column pointer data...
    my @colpointers;
    for my $i  (1..$pline) {
	$_ = <HB>;
	s/^\s*//;
	s/\s*$//;
	push @colpointers, split(/\s+/);
    }
    #...and the row data...
    my @rowindex;
    for my $i (1..$rline) {
	$_ = <HB>;
	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) {
	    $_ = <HB>;
	    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";
    $_=<MM>;
    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 (<MM>) {
	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 ($…