From bug-request at octave dot org Wed Jan 26 10:37:11 2005
Subject: Re: Bug concatenation of NDArray's
From: David Bateman <David dot Bateman at motorola dot com>
To: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
Cc: bug at octave dot org
Date: Wed, 26 Jan 2005 17:38:22 +0100

This is a multi-part message in MIME format.
--------------080703030608070604030603
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit

John W. Eaton wrote:

>Yes, it seems that it could be.  Even if this is incorrect usage, I
>think it should be detected so we can issue an error message instead
>of (potentially) crashing.
>
>  
>
Ok, then consider the following patch, that detects the typical case 
that is used for matrix concatenation, and uses the fast existing code, 
and if that can't be used falls back to the generic code that is about 
twice as slow for 2D arrays (probably slower for ND)..

Regards
David

-- 

David Bateman                                David dot Bateman at motorola dot com
Motorola CRM                                 +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary


--------------080703030608070604030603
Content-Type: text/plain;
 name="patch.insert-20050126"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="patch.insert-20050126"

*** Array.cc~	2005-01-23 00:35:41.000000000 +0100
--- Array.cc	2005-01-26 14:00:23.739371977 +0100
***************
*** 1041,1046 ****
--- 1041,1047 ----
        dim_vector dva = a.dims ();
        dim_vector dv = dims ();
        int len_a = dva.length ();
+       int non_full_dim = 0;
  
        for (int i = 0; i < n; i++)
  	{
***************
*** 1051,1087 ****
  		("Array<T>::insert: range error for insert");
  	      return *this;
  	    }
  	}
  
        if (dva.numel ())
          {
! 	  const T *a_data = a.data ();   
! 	  int numel_to_move = 1;
! 	  int skip = 0;
! 	  for (int i = 0; i < len_a; i++) 
! 	    if (ra_idx(i) == 0 && dva(i) == dv(i))
! 	      numel_to_move *= dva(i);
! 	    else
! 	      {
! 		skip = numel_to_move * (dv(i) - dva(i));
! 		numel_to_move *= dva(i);
! 		break;
! 	      }
! 
! 	  int jidx = ra_idx (n - 1);
! 	  for (int i = n-2; i >= 0; i--)
  	    {
! 	      jidx *= dv (i);
! 	      jidx += ra_idx (i);
! 	    }
  
! 	  int iidx = 0;
! 	  int moves = dva.numel () / numel_to_move;
! 	  for (int i = 0; i < moves; i++)
  	    {
! 	      for (int j = 0; j < numel_to_move; j++)
! 		elem (jidx++) = a_data[iidx++];
! 	      jidx += skip;
  	    }
  	}
      }
--- 1052,1117 ----
  		("Array<T>::insert: range error for insert");
  	      return *this;
  	    }
+ 
+ 	  if (dv(i) != (i < len_a ? dva(i) : 1))
+ 	    non_full_dim++;
  	}
  
        if (dva.numel ())
          {
! 	  if (non_full_dim < 2)
  	    {
! 	      // Special case for fast concatenation
  
! 	      const T *a_data = a.data ();   
! 	      int numel_to_move = 1;
! 	      int skip = 0;
! 	      for (int i = 0; i < len_a; i++) 
! 		if (ra_idx(i) == 0 && dva(i) == dv(i))
! 		  numel_to_move *= dva(i);
! 		else
! 		  {
! 		    skip = numel_to_move * (dv(i) - dva(i));
! 		    numel_to_move *= dva(i);
! 		    break;
! 		  }
! 
! 	      int jidx = ra_idx (n - 1);
! 	      for (int i = n-2; i >= 0; i--)
! 		{
! 		  jidx *= dv (i);
! 		  jidx += ra_idx (i);
! 		}
! 
! 	      int iidx = 0;
! 	      int moves = dva.numel () / numel_to_move;
! 	      for (int i = 0; i < moves; i++)
! 		{
! 		  for (int j = 0; j < numel_to_move; j++)
! 		    elem (jidx++) = a_data[iidx++];
! 		  jidx += skip;
! 		}
! 	    }
! 	  else
  	    {
! 	      // Generic code
! 	      const T *a_data = a.data ();   
! 	      int nel = a.numel ();
! 	      Array<int> a_idx (n, 0);
! 
! 	      for (int i = 0; i < nel; i++)
! 		{
! 		  int iidx = a_idx(n-1) + ra_idx(n-1);
! 		  for (int j = n-2; j >= 0; j--)
! 		    {
! 		      iidx *= dv(j);
! 		      iidx += a_idx(j) + ra_idx(j);
! 		    }
! 		  
! 		  elem (iidx) = a_data[i];
! 
! 		  increment_index (a_idx, dva); 
! 		}
  	    }
  	}
      }

--------------080703030608070604030603--



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------