#======================================================================
#                    G A U S S J . P L 
#                    doc: Wed Feb 24 17:06:55 1999
#                    dlm: Fri Jan  6 10:23:44 2012
#                    (c) 1999 A.M. Thurnherr
#                    uE-Info: 46 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# GAUSSJ routine from Numerical Recipes adapted to ANTS

# Notes:
#	- both @A and @B passed by ref

# HISTORY:
#	Feb 24, 1999: - apparently created
#	Jul 19, 2001: - apparently fiddled
#	Jan  6, 2011: - added code to check for numericity of input

sub gaussj($$)
{
	my($AR,$BR) = @_;
	my($n) = $#{$AR};
	my($m) = $#{$BR->[1]};
	my(@indxc,@indxr,@ipiv);
	my($i,$icol,$irow,$j,$k,$l,$ll);
	my($big,$dum,$pivinv);
	my($temp);

#	print(STDERR "n = $n, m = $m\n");
#	for ($i=1; $i<=$n; $i++) {
#		for ($j=1; $j<=$n; $j++) {
#			print(STDERR "A[$i][$j] = $AR->[$i][$j]\n");
#		}
#	}

	&vector(\@indxc,1,$n);
	&vector(\@indxr,1,$n);
	&vector(\@ipiv, 1,$n);
	for ($j=1; $j<=$n; $j++) { $ipiv[$j] = 0; }
	for ($i=1; $i<=$n; $i++) {
		$big = 0.0;
		for ($j=1; $j<=$n; $j++) {
			if ($ipiv[$j] != 1) {
				for ($k=1; $k<=$n; $k++) {
					if ($ipiv[$k] == 0) {
						croak("GAUSSJ: non-numeric A[$j][$k]\n")
							unless numberp($AR->[$j][$k]);
						if (abs($AR->[$j][$k]) >= $big) {
							$big = abs($AR->[$j][$k]);
							$irow = $j;
							$icol = $k;
						}
					} elsif ($ipiv[$k] > 1) {
						croak("GAUSSJ: Singular Matrix-1");
					}
				}
			}
		}
		++($ipiv[$icol]);
		if ($irow != $icol) {
			for ($l=1; $l<=$n; $l++) {
				$temp = $AR->[$irow][$l];
				$AR->[$irow][$l] = $AR->[$icol][$l];
				$AR->[$icol][$l] = $temp;
			}
			for ($l=1; $l<=$m; $l++) {
				croak("GAUSSJ: non-numeric B[$irow][$l]\n")
						unless numberp($BR->[$irow][$l]);
				croak("GAUSSJ: non-numeric B[$icol][$l]\n")
						unless numberp($BR->[$icol][$l]);
				$temp = $BR->[$irow][$l];
				$BR->[$irow][$l] = $BR->[$icol][$l];
				$BR->[$icol][$l] = $temp;
			}
		}
		$indxr[$i] = $irow;
		$indxc[$i] = $icol;
		if ($AR->[$icol][$icol] == 0.0) {
			croak("GAUSSJ: Singular Matrix-2");
		}
		$pivinv = 1.0/$AR->[$icol][$icol];
		$AR->[$icol][$icol] = 1.0;
		for ($l=1; $l<=$n; $l++) {
			$AR->[$icol][$l] *= $pivinv;
		}
		for ($l=1; $l<=$m; $l++) {
			$BR->[$icol][$l] *= $pivinv;
		}
		for ($ll=1; $ll<=$n; $ll++) {
			if ($ll != $icol) {
				$dum = $AR->[$ll][$icol];
				$AR->[$ll][$icol] = 0.0;
				for ($l=1; $l<=$n; $l++) {
					$AR->[$ll][$l] -= $AR->[$icol][$l]*$dum;
				}
				for ($l=1; $l<=$m; $l++) {
					$BR->[$ll][$l] -= $BR->[$icol][$l]*$dum;
				}
			}
		}
	}
	for ($l=$n; $l>=1; $l--) {
		if ($indxr[$l] != $indxc[$l]) {
			for ($k=1; $k<=$n; $k++) {
				$temp = $AR->[$k][$indxr[$l]];
				$AR->[$k][$indxr[$l]] = $AR->[$k][$indxc[$l]];
				$AR->[$k][$indxc[$l]] = $temp;
			}
		}
	}
}

1;