# HG changeset patch
# User Martin Albrecht <malb@informatik.uni-bremen.de>
# Date 1246293690 -3600
# Node ID 879c024ea778fdd968ac535cbc53424b4052f2a5
# Parent  52b259b6967803c1355cc5fb23ff491fa94d89d5
some experimental hackish heurstic to switch between M4RI and LQUP factorisation when computing the row echelon form.

diff -r 52b259b69678 -r 879c024ea778 src/brilliantrussian.c
--- a/src/brilliantrussian.c	Sat Jun 27 14:39:27 2009 +0100
+++ b/src/brilliantrussian.c	Mon Jun 29 17:41:30 2009 +0100
@@ -29,7 +29,7 @@
 
 #include "brilliantrussian.h"
 #include "grayflex.h"
-
+#include "lqup.h"
 
 /**
  * \brief Perform Gaussian reduction to reduced row echelon form on a
@@ -569,7 +569,7 @@
   }
 }
 
-size_t mzd_echelonize_m4ri(mzd_t *A, int full, int k) {
+size_t _mzd_echelonize_m4ri(mzd_t *A, int full, int k, int heuristic) {
   /**
    * \par General algorithm
    * \li Step 1.Denote the first column to be processed in a given
@@ -634,7 +634,29 @@
   size_t *L4 = (size_t *)m4ri_mm_calloc(TWOPOW(k), sizeof(size_t));
   size_t *L5 = (size_t *)m4ri_mm_calloc(TWOPOW(k), sizeof(size_t));
 
+  size_t last_check = 0;
+
+  if (heuristic) {
+    if (c<A->ncols && r< A->nrows && mzd_density(A,32)*A->ncols >= 10) {
+      mzd_t *Abar = mzd_init_window(A, r, (c/RADIX)*RADIX, A->nrows, A->ncols);
+      r += mzd_echelonize_pluq(Abar, full);
+      mzd_free(Abar);
+      c = ncols;
+    }
+  }
+
   while(c<ncols) {
+    if (heuristic && c > (last_check + 2048)) {
+      last_check = c;
+      if (c<A->ncols && r< A->nrows && mzd_density(A,32)*A->ncols >= 10) {
+        mzd_t *Abar = mzd_init_window(A, r, (c/RADIX)*RADIX, A->nrows, A->ncols);
+        r += mzd_echelonize_pluq(Abar, full);
+        c = ncols;
+        mzd_free(Abar);
+        break;
+      }
+    }
+
     if(c+kk > A->ncols) {
       kk = ncols - c;
     }
@@ -1265,3 +1287,11 @@
   return C;
 }
 
+
+size_t mzd_echelonize_m4ri(mzd_t *A, int full, int k) {
+  return _mzd_echelonize_m4ri(A, full, k, 0);
+}
+
+size_t mzd_echelonize(mzd_t *A, int full) {
+  return _mzd_echelonize_m4ri(A, full, 0, 1);
+}
diff -r 52b259b69678 -r 879c024ea778 src/brilliantrussian.h
--- a/src/brilliantrussian.h	Sat Jun 27 14:39:27 2009 +0100
+++ b/src/brilliantrussian.h	Mon Jun 29 17:41:30 2009 +0100
@@ -260,4 +260,6 @@
 
 #define M4RM_GRAY8
 
+size_t mzd_echelonize(mzd_t *A, int full);
+
 #endif //BRILLIANTRUSSIAN_H
diff -r 52b259b69678 -r 879c024ea778 src/lqup.c
--- a/src/lqup.c	Sat Jun 27 14:39:27 2009 +0100
+++ b/src/lqup.c	Mon Jun 29 17:41:30 2009 +0100
@@ -84,6 +84,7 @@
 
     size_t i, j;
     size_t n1 = (((ncols - 1) / RADIX + 1) >> 1) * RADIX;
+
     mzd_t *A0  = mzd_init_window(A,  0,  0, nrows,    n1);
     mzd_t *A1  = mzd_init_window(A,  0, n1, nrows, ncols);