# HG changeset patch
# User jwe
# Date 1106326941 0
# Node ID cda8c0a823c569d6a2c248cf81ba5ea52b8e091e
# Parent  cd9a8ae4e6d48efc0e55640b1e3c469d7b4565d1
[project @ 2005-01-21 17:02:21 by jwe]

diff -r cd9a8ae4e6d4 -r cda8c0a823c5 src/ChangeLog
--- a/src/ChangeLog	Fri Jan 21 03:32:09 2005 +0000
+++ b/src/ChangeLog	Fri Jan 21 17:02:21 2005 +0000
@@ -1,3 +1,7 @@
+2005-01-21  John W. Eaton  <jwe@octave.org>
+
+	* DLD-FUNCTIONS/filter.cc (filter): Avoid slow Marray indexing ops.
+
 2005-01-20  John W. Eaton  <jwe@octave.org>
 
 	* Makefile.in (EXTRAS): Move ov-base-mat.cc ov-base-scalar.cc here.
diff -r cd9a8ae4e6d4 -r cda8c0a823c5 src/DLD-FUNCTIONS/filter.cc
--- a/src/DLD-FUNCTIONS/filter.cc	Fri Jan 21 03:32:09 2005 +0000
+++ b/src/DLD-FUNCTIONS/filter.cc	Fri Jan 21 17:02:21 2005 +0000
@@ -146,10 +146,51 @@
 
       if (a_len > 1)
 	{
-	  for (int i = 0; i < x_len; i++)
+	  T *py = y.fortran_vec ();
+	  T *psi = si.fortran_vec ();
+
+	  const T *pa = a.data ();
+	  const T *pb = b.data ();
+	  const T *px = x.data ();
+
+	  psi += si_offset;
+
+	  for (int i = 0, idx = x_offset; i < x_len; i++, idx += x_stride)
 	    {
-	      int idx = i * x_stride + x_offset; 
-	      y (idx) = si (si_offset) + b (0) * x (idx);
+	      py[idx] = psi[0] + pb[0] * px[idx];
+
+	      if (si_len > 0)
+		{
+		  for (int j = 0; j < si_len - 1; j++)
+		    {
+		      OCTAVE_QUIT;
+
+		      psi[j] = psi[j+1] - pa[j+1] * py[idx] + pb[j+1] * px[idx];
+		    }
+
+		  psi[si_len-1] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
+		}
+	      else
+		{
+		  OCTAVE_QUIT;
+
+		  psi[0] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
+		}
+	    }
+	}
+      else if (si_len > 0)
+	{
+	  T *py = y.fortran_vec ();
+	  T *psi = si.fortran_vec ();
+
+	  const T *pb = b.data ();
+	  const T *px = x.data ();
+
+	  psi += si_offset;
+
+	  for (int i = 0, idx = x_offset; i < x_len; i++, idx += x_stride)
+	    {
+	      py[idx] = psi[0] + pb[0] * px[idx];
 
 	      if (si_len > 1)
 		{
@@ -157,39 +198,17 @@
 		    {
 		      OCTAVE_QUIT;
 
-		      si (j + si_offset) =  si (j + 1 + si_offset) - 
-			a (j+1) * y (idx) + b (j+1) * x (idx);
+		      psi[j] = psi[j+1] + pb[j+1] * px[idx];
 		    }
 
-		  si (si_len - 1 + si_offset) = b (si_len) * x (idx)
-		    - a (si_len) * y (idx);
+		  psi[si_len-1] = pb[si_len] * px[idx];
 		}
 	      else
-		si (si_offset) = b (si_len) * x (idx)
-		  - a (si_len) * y (idx);
-	    }
-	}
-      else if (si_len > 0)
-	{
-	  for (int i = 0; i < x_len; i++)
-	    {
-	      int idx = i * x_stride + x_offset; 
-	      y (idx) = si (si_offset) + b (0) * x (idx);
+		{
+		  OCTAVE_QUIT;
 
-	      if (si_len > 1)
-		{
-		  for (int j = 0; j < si_len - 1; j++)
-		    {
-		      OCTAVE_QUIT;
-
-		      si (j + si_offset) = si (j + 1 + si_offset) + 
-			b (j+1) * x (idx);
-		    }
-
-		  si (si_len - 1 + si_offset) = b (si_len) * x (idx);
+		  psi[0] = pb[1] * px[idx];
 		}
-	      else
-		si (si_offset) = b (1) * x (idx);
 	    }
 	}
     }