From 25947adc6e75e5c762e88aaf63d72d84f1980526 Mon Sep 17 00:00:00 2001
From: Thierry Moisan <thierry.moisan@gmail.com>
Date: Sun, 29 Jun 2008 22:49:07 -0400
Subject: [PATCH] Adding the is_similar_relative function and making the
 QR_decomp tests finally work

---
 lib/Math/GSL.pm             | 26 +++++++++++++++++++++++++-
 lib/Math/GSL/Linalg/Test.pm |  7 ++++---
 2 files changed, 29 insertions(+), 4 deletions(-)

diff --git a/lib/Math/GSL.pm b/lib/Math/GSL.pm
index c46ec5d..fdc1178 100644
--- a/lib/Math/GSL.pm
+++ b/lib/Math/GSL.pm
@@ -12,7 +12,7 @@ require DynaLoader;
 require Exporter;
 our @ISA = qw(Exporter DynaLoader);
 our @EXPORT = qw();
-our @EXPORT_OK = qw( ok_similar is_similar verify_results $GSL_MODE_DEFAULT $GSL_PREC_DOUBLE $GSL_PREC_SINGLE $GSL_PREC_APPROX);
+our @EXPORT_OK = qw( ok_similar is_similar is_similar_relative verify_results $GSL_MODE_DEFAULT $GSL_PREC_DOUBLE $GSL_PREC_SINGLE $GSL_PREC_APPROX);
 our %EXPORT_TAGS = ( all => [ @EXPORT_OK ] );
 
 our ($GSL_PREC_DOUBLE, $GSL_PREC_SINGLE, $GSL_PREC_APPROX ) = 0..2;
@@ -216,6 +216,30 @@ sub ok_similar {
     ok(is_similar($x,$y,$eps), $msg);
 }
 
+sub is_similar_relative {
+    my ($x,$y, $eps) = @_;
+    $eps ||= 1e-8;
+    if (ref $x eq 'ARRAY' && ref $y eq 'ARRAY') {
+        if ( $#$x != $#$y ){
+            warn "is_similar(): argument of different lengths, $#$x != $#$y !!!";
+            return 0;
+        } else {
+            map { 
+                    my $delta = abs($x->[$_] - $y->[$_]);
+                    if($delta > $eps){
+                        warn "\n\tElements start differing at index $_, delta = $delta\n";
+                        warn qq{\t\t\$x->[$_] = } . $x->[$_] . "\n";
+                        warn qq{\t\t\$y->[$_] = } . $y->[$_] . "\n";
+                        return 0;
+                    }
+                } (0..$#$x);
+            return 1;
+        }
+    } else {
+        (abs($x-$y)/abs($y)) <= $eps ? return 1 : return 0;
+    }
+}
+
 sub is_valid_double
 {
     my $x=shift;
diff --git a/lib/Math/GSL/Linalg/Test.pm b/lib/Math/GSL/Linalg/Test.pm
index fd5d757..315c093 100644
--- a/lib/Math/GSL/Linalg/Test.pm
+++ b/lib/Math/GSL/Linalg/Test.pm
@@ -241,7 +241,6 @@ sub GSL_LINALG_LU_LNDET : Tests {
 }
 
 sub GSL_LINALG_QR_DECOMP : Tests {
-# stolen from 
     my $matrix = gsl_matrix_alloc(3,5);
     my ($i, $j);
     for($i=0; $i<3; $i++) {
@@ -254,6 +253,8 @@ sub GSL_LINALG_QR_DECOMP : Tests {
     my $q = gsl_matrix_alloc(3,3);    
     my $r = gsl_matrix_alloc(3,5);    
     my $a = gsl_matrix_alloc(3,5);
+    my $save = gsl_matrix_alloc(3, 5);
+    gsl_matrix_memcpy($save, $matrix);
 
     is(gsl_linalg_QR_decomp($matrix, $tau),0);
     is(gsl_linalg_QR_unpack($matrix, $tau, $q, $r), 0);
@@ -264,8 +265,8 @@ sub GSL_LINALG_QR_DECOMP : Tests {
   for($i=0; $i<3; $i++) {
     for($j=0; $j<5; $j++) {
       $aij = gsl_matrix_get($a, $i, $j);
-      $mij = gsl_matrix_get($matrix, $i, $j);
-      ok_similar($aij, $mij, "QR decomposition", 2 * 8.0 * $GSL_DBL_EPSILON);      
+      $mij = gsl_matrix_get($save, $i, $j);
+      ok(is_similar_relative($aij, $mij, 2 * 8.0 * $GSL_DBL_EPSILON));      
     }
   }
 
-- 
2.11.4.GIT