# HG changeset patch
# User idamnjanovic
# Date 1300191659 0
# Node ID ad36f80e2ccfe627b34445040d8c6bcdb1e49d6a
# Parent  23f9dd7b9d785fdba6ffdde011c771d9e23620cf


diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/SolveFISTA.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/SolveFISTA.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,211 @@
+% Copyright ©2010. The Regents of the University of California (Regents). 
+% All Rights Reserved. Contact The Office of Technology Licensing, 
+% UC Berkeley, 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, 
+% (510) 643-7201, for commercial licensing opportunities.
+
+% Authors: Arvind Ganesh, Allen Y. Yang, Zihan Zhou.
+% Contact: Allen Y. Yang, Department of EECS, University of California,
+% Berkeley. <yang@eecs.berkeley.edu>
+
+% IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, 
+% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, 
+% ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF 
+% REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+% REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED
+% TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
+% PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY, 
+% PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO 
+% PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+
+%% This function is modified from Matlab code proximal_gradient_bp
+
+function [x_hat,nIter] = SolveFISTA(A,b, varargin)
+
+% b - m x 1 vector of observations/data (required input)
+% A - m x n measurement matrix (required input)
+%
+% tol - tolerance for stopping criterion.
+%     - DEFAULT 1e-7 if omitted or -1.
+% maxIter - maxilambdam number of iterations
+%         - DEFAULT 10000, if omitted or -1.
+% lineSearchFlag - 1 if line search is to be done every iteration
+%                - DEFAULT 0, if omitted or -1.
+% continuationFlag - 1 if a continuation is to be done on the parameter lambda
+%                  - DEFAULT 1, if omitted or -1.
+% eta - line search parameter, should be in (0,1)
+%     - ignored if lineSearchFlag is 0.
+%     - DEFAULT 0.9, if omitted or -1.
+% lambda - relaxation parameter
+%    - ignored if continuationFlag is 1.
+%    - DEFAULT 1e-3, if omitted or -1.
+% outputFileName - Details of each iteration are dumped here, if provided.
+%
+% x_hat - estimate of coeeficient vector
+% numIter - number of iterations until convergence
+%
+%
+% References
+% "Robust PCA: Exact Recovery of Corrupted Low-Rank Matrices via Convex Optimization", J. Wright et al., preprint 2009.
+% "An Accelerated Proximal Gradient Algorithm for Nuclear Norm Regularized Least Squares problems", K.-C. Toh and S. Yun, preprint 2009.
+%
+% Arvind Ganesh, Summer 2009. Questions? abalasu2@illinois.edu
+
+DEBUG = 0 ;
+
+STOPPING_GROUND_TRUTH = -1;
+STOPPING_DUALITY_GAP = 1;
+STOPPING_SPARSE_SUPPORT = 2;
+STOPPING_OBJECTIVE_VALUE = 3;
+STOPPING_SUBGRADIENT = 4;
+STOPPING_DEFAULT = STOPPING_SUBGRADIENT;
+
+stoppingCriterion = STOPPING_DEFAULT;
+maxIter = 1000 ;
+tolerance = 1e-3;
+[m,n] = size(A) ;
+x0 = zeros(n,1) ;
+xG = [];
+
+%% Initializing optimization variables
+t_k = 1 ; 
+t_km1 = 1 ;
+L0 = 1 ;
+G = A'*A ;
+nIter = 0 ;
+c = A'*b ;
+lambda0 = 0.99*L0*norm(c,inf) ;
+eta = 0.6 ;
+lambda_bar = 1e-4*lambda0 ;
+xk = zeros(n,1) ;
+lambda = lambda0 ;
+L = L0 ;
+beta = 1.5 ;
+
+% Parse the optional inputs.
+if (mod(length(varargin), 2) ~= 0 ),
+    error(['Extra Parameters passed to the function ''' mfilename ''' lambdast be passed in pairs.']);
+end
+parameterCount = length(varargin)/2;
+
+for parameterIndex = 1:parameterCount,
+    parameterName = varargin{parameterIndex*2 - 1};
+    parameterValue = varargin{parameterIndex*2};
+    switch lower(parameterName)
+        case 'stoppingcriterion'
+            stoppingCriterion = parameterValue;
+        case 'groundtruth'
+            xG = parameterValue;
+        case 'tolerance'
+            tolerance = parameterValue;
+        case 'linesearchflag'
+            lineSearchFlag = parameterValue;
+        case 'lambda'
+            lambda_bar = parameterValue;
+        case 'maxiteration'
+            maxIter = parameterValue;
+        case 'isnonnegative'
+            isNonnegative = parameterValue;
+        case 'continuationflag'
+            continuationFlag = parameterValue;
+        case 'initialization'
+            xk = parameterValue;
+            if ~all(size(xk)==[n,1])
+                error('The dimension of the initial xk does not match.');
+            end
+        case 'eta'
+            eta = parameterValue;
+            if ( eta <= 0 || eta >= 1 )
+                disp('Line search parameter out of bounds, switching to default 0.9') ;
+                eta = 0.9 ;
+            end
+        otherwise
+            error(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''.']);
+    end
+end
+clear varargin
+
+if stoppingCriterion==STOPPING_GROUND_TRUTH && isempty(xG)
+    error('The stopping criterion must provide the ground truth value of x.');
+end
+
+keep_going = 1 ;
+nz_x = (abs(xk)> eps*10);
+f = 0.5*norm(b-A*xk)^2 + lambda_bar * norm(xk,1);
+xkm1 = xk;
+while keep_going && (nIter < maxIter)
+    nIter = nIter + 1 ;
+    
+    yk = xk + ((t_km1-1)/t_k)*(xk-xkm1) ;
+    
+    stop_backtrack = 0 ;
+    
+    temp = G*yk - c ; % gradient of f at yk
+    
+    while ~stop_backtrack
+        
+        gk = yk - (1/L)*temp ;
+        
+        xkp1 = soft(gk,lambda/L) ;
+        
+        temp1 = 0.5*norm(b-A*xkp1)^2 ;
+        temp2 = 0.5*norm(b-A*yk)^2 + (xkp1-yk)'*temp + (L/2)*norm(xkp1-yk)^2 ;
+        
+        if temp1 <= temp2
+            stop_backtrack = 1 ;
+        else
+            L = L*beta ;
+        end
+        
+    end
+    
+    switch stoppingCriterion
+        case STOPPING_GROUND_TRUTH
+            keep_going = norm(xG-xkp1)>tolerance;
+        case STOPPING_SUBGRADIENT
+            sk = L*(yk-xkp1) + G*(xkp1-yk) ;
+            keep_going = norm(sk) > tolerance*L*max(1,norm(xkp1));
+        case STOPPING_SPARSE_SUPPORT
+            % compute the stopping criterion based on the change
+            % of the number of non-zero components of the estimate
+            nz_x_prev = nz_x;
+            nz_x = (abs(xkp1)>eps*10);
+            num_nz_x = sum(nz_x(:));
+            num_changes_active = (sum(nz_x(:)~=nz_x_prev(:)));
+            if num_nz_x >= 1
+                criterionActiveSet = num_changes_active / num_nz_x;
+                keep_going = (criterionActiveSet > tolerance);
+            end
+        case STOPPING_OBJECTIVE_VALUE
+            % compute the stopping criterion based on the relative
+            % variation of the objective function.
+            prev_f = f;
+            f = 0.5*norm(b-A*xkp1)^2 + lambda_bar * norm(xk,1);
+            criterionObjective = abs(f-prev_f)/(prev_f);
+            keep_going =  (criterionObjective > tolerance);
+        case STOPPING_DUALITY_GAP
+            error('Duality gap is not a valid stopping criterion for PGBP.');
+        otherwise
+            error('Undefined stopping criterion.');
+    end
+    
+    lambda = max(eta*lambda,lambda_bar) ;
+    
+    
+    t_kp1 = 0.5*(1+sqrt(1+4*t_k*t_k)) ;
+    
+    t_km1 = t_k ;
+    t_k = t_kp1 ;
+    xkm1 = xk ;
+    xk = xkp1 ;
+end
+
+x_hat = xk ;
+
+function y = soft(x,T)
+if sum(abs(T(:)))==0
+    y = x;
+else
+    y = max(abs(x) - T, 0);
+    y = sign(x).*y;
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/add_dc.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/add_dc.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,33 @@
+function x = add_dc(y,dc,columns)
+%ADD_DC Add DC channel to signals.
+%   X = ADD_DC(Y,DC) adds the specified DC value to the (possibly
+%   multi-dimensional) signal Y, returning the result as X. DC should be a
+%   scalar value.
+%
+%   X = ADD_DC(Y,DC,'columns') where Y is a 2D matrix and DC is an array of
+%   length size(Y,2), treats the columns of Y as individual 1D signals, 
+%   adding to each one the corresponding DC value from the DC array. X is
+%   the same size as Y and contains the resulting signals.
+%
+%   See also REMOVE_DC.
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
+
+
+if (nargin==3 && strcmpi(columns,'columns')), columns = 1;
+else columns = 0;
+end
+
+if (columns)
+  x = addtocols(y,dc);
+else
+  x = y + dc;
+end
+
+
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/addtocols.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/addtocols.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,85 @@
+/**************************************************************************
+ *
+ * File name: addtocols.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 19.4.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+
+
+/* Input Arguments */
+
+#define	X_IN	prhs[0]
+#define V_IN  prhs[1]
+
+
+/* Output Arguments */
+
+#define	Y_OUT	plhs[0]
+
+
+void mexFunction(int nlhs, mxArray *plhs[], 
+		             int nrhs, const mxArray*prhs[])
+     
+{ 
+    double *x, *y, *v, *xend;  
+    mwSize m,n,m1,n1;
+    mwIndex counter; 
+    
+    
+    /* Check for proper number of arguments */
+    
+    if (nrhs != 2) {
+      mexErrMsgTxt("Two input arguments required."); 
+    } else if (nlhs > 1) {
+      mexErrMsgTxt("Too many output arguments."); 
+    } 
+    
+    
+    /* Check the the input dimensions */ 
+    
+    m = mxGetM(X_IN);
+    n = mxGetN(X_IN);
+    if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || mxGetNumberOfDimensions(X_IN)>2) {
+      mexErrMsgTxt("ADDTOCOLS requires that X be a double matrix.");
+    }
+    m1 = mxGetM(V_IN);
+    n1 = mxGetN(V_IN);
+    if (!mxIsDouble(V_IN) || mxIsComplex(V_IN) || (m1!=1 && n1!=1)) {
+      mexErrMsgTxt("ADDTOCOLS requires that V be a double vector.");
+    }
+    if ((m1==1 && n1!=n) || (n1==1 && m1!=n)) {
+      mexErrMsgTxt("Error in ADDTOCOLS: dimensions of V and X must agree.");
+    }
+    
+    
+    /* Create a matrix for the return argument */ 
+    Y_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
+         
+    
+    /* Assign pointers to the various parameters */ 
+    x = mxGetPr(X_IN); 
+    v = mxGetPr(V_IN);
+    y = mxGetPr(Y_OUT);
+            
+    
+    /* Do the actual computation */
+    
+    xend = x+(m*n);
+    counter = 0;
+    while (x<xend) {
+      (*y) = (*x) + (*v);
+      y++; x++; counter++;
+      if (counter==m) {v++; counter=0;}
+    }
+    
+    return;
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/addtocols.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/addtocols.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,13 @@
+%ADDTOCOLS Add values to the columns of a matrix.
+%  Y=ADDTOCOLS(X,V) adds to each column of the MxN matrix X the
+%  corresponding value from the N-element vector V.
+%
+%  See also NORMCOLS.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  June 2005
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/addtocols.mexa64
Binary file DL/RLS-DLA/private/addtocols.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/col2imstep.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/col2imstep.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,146 @@
+/**************************************************************************
+ *
+ * File name: col2imstep.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 31.8.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+
+
+/* Input Arguments */
+
+#define	B_IN	 prhs[0]
+#define N_IN   prhs[1]
+#define SZ_IN  prhs[2]
+#define S_IN   prhs[3]
+
+
+/* Output Arguments */
+
+#define	X_OUT	plhs[0]
+
+
+void mexFunction(int nlhs, mxArray *plhs[], 
+		             int nrhs, const mxArray*prhs[])
+     
+{ 
+    double *x, *b, *s;
+    mwSize sz[3], stepsize[3], n[3], ndims;
+    mwIndex i, j, k, l, m, t, blocknum;
+    
+    
+    /* Check for proper number of arguments */
+    
+    if (nrhs < 3 || nrhs > 4) {
+      mexErrMsgTxt("Invalid number of input arguments."); 
+    } else if (nlhs > 1) {
+      mexErrMsgTxt("Too many output arguments."); 
+    } 
+    
+    
+    /* Check the the input dimensions */ 
+    
+    if (!mxIsDouble(B_IN) || mxIsComplex(B_IN) || mxGetNumberOfDimensions(B_IN)>2) {
+      mexErrMsgTxt("B should be a double matrix.");
+    }
+    if (!mxIsDouble(N_IN) || mxIsComplex(N_IN) || mxGetNumberOfDimensions(N_IN)>2) {
+      mexErrMsgTxt("Invalid output matrix size.");
+    }
+    ndims = mxGetM(N_IN)*mxGetN(N_IN);
+    if (ndims<2 || ndims>3) {
+      mexErrMsgTxt("Output matrix can only be 2-D or 3-D.");
+    }
+    if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    if (nrhs == 4) {
+      if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+    }
+    
+    
+    /* Get parameters */
+    
+    s = mxGetPr(N_IN);
+    if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+      mexErrMsgTxt("Invalid output matrix size.");
+    }
+    n[0] = (mwSize)(s[0] + 0.01);
+    n[1] = (mwSize)(s[1] + 0.01);
+    n[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+    
+    s = mxGetPr(SZ_IN);
+    if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    sz[0] = (mwSize)(s[0] + 0.01);
+    sz[1] = (mwSize)(s[1] + 0.01);
+    sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+    
+    if (nrhs == 4) {
+      s = mxGetPr(S_IN);
+      if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+      stepsize[0] = (mwSize)(s[0] + 0.01);
+      stepsize[1] = (mwSize)(s[1] + 0.01);
+      stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+    }
+    else {
+      stepsize[0] = stepsize[1] = stepsize[2] = 1;
+    }
+    
+    if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) {
+      mexErrMsgTxt("Block size too large.");
+    }
+
+    if (mxGetN(B_IN) != ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1)) {
+      mexErrMsgTxt("Invalid number of columns in B. Please use IM2COLSTEP to compute B.");
+    }
+    
+    
+    /* Create a matrix for the return argument */
+    
+    X_OUT = mxCreateNumericArray(ndims, n, mxDOUBLE_CLASS, mxREAL);
+    
+    
+    /* Assign pointers */
+    
+    b = mxGetPr(B_IN);
+    x = mxGetPr(X_OUT);
+            
+    
+    /* Do the actual computation */
+    
+    blocknum = 0;
+    
+    /* iterate over all blocks */
+    for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) {
+      for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) {
+        for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) {
+          
+          /* add single block */
+          for (m=0; m<sz[2]; m++) {
+            for (l=0; l<sz[1]; l++) {
+              for (t=0; t<sz[0]; t++) {
+                (x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i)[t] += (b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0])[t];
+              }
+            }
+          }
+          blocknum++;
+          
+        }
+      }
+    }
+    
+    return;
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/col2imstep.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/col2imstep.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,25 @@
+%COL2IMSTEP Rearrange matrix columns into blocks.
+%  A = COL2IMSTEP(B,[MM NN],[N1 N2]) rearranges the columns of B into
+%  sliding N1-by-N2 blocks producing the matrix A of size MM-by-NN. B is
+%  usually the result of calling IM2COLSTEP(...) with a stepsize of 1, or
+%  using Matlab's IM2COL(..,'sliding'). Overlapping blocks are summed in A.
+%
+%  A = COL2IMSTEP(B,[MM NN],[N1 N2],[S1 S2]) arranges the blocks in A with
+%  a step size of (S1,S2) between them. The first block is at A(1:N1,1:N2),
+%  and the rest are at A((1:N1)+i*S1,(1:N2)+j*S2). Overlapping blocks are
+%  summed in A. Note that B is usually the result of calling
+%  IM2COLSTEP(...) with a stepsize of [S1 S2].
+%
+%  A = IM2COLSTEP(B,[MM NN KK],[N1 N2 N3],[S1 S2 S3]) generates a 3-D
+%  output matrix A. The step size [S1 S2 S3] may be ommitted, and defaults
+%  to [1 1 1].
+%
+%  See also IM2COLSTEP, IM2COL, COUNTCOVER.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/col2imstep.mexa64
Binary file DL/RLS-DLA/private/col2imstep.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/collincomb.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/collincomb.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,165 @@
+/**************************************************************************
+ *
+ * File name: collincomb.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 21.5.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+
+
+/* Input Arguments */
+
+#define	A_IN	   prhs[0]
+#define ROWS_IN  prhs[1]
+#define COLS_IN1 prhs[1]
+#define COLS_IN2 prhs[2]
+#define X_IN1    prhs[2]
+#define X_IN2    prhs[3]
+
+
+/* Output Arguments */
+
+#define	Y_OUT	plhs[0]
+
+
+void mexFunction(int nlhs, mxArray *plhs[],
+int nrhs, const mxArray*prhs[])
+
+{
+  double *A, *x, *y, *rows, *cols;
+  mwSize m,n,m1,n1,m2,n2,rownum,colnum;
+  mwIndex col_id,*row_ids,i,j;
+  int rownumspecified=0;
+  
+  
+  /* Check for proper number of arguments */
+  
+  if (nrhs!=3 && nrhs!=4) {
+    mexErrMsgTxt("Invalid number of arguments.");
+  } else if (nlhs > 1) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+  
+  
+  /* Check the input dimensions */
+  
+  m = mxGetM(A_IN);
+  n = mxGetN(A_IN);
+  if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) {
+    mexErrMsgTxt("COLLINCOMB requires that A be a double matrix.");
+  }
+  
+  if (nrhs==3) {
+    
+    m1 = mxGetM(COLS_IN1);
+    n1 = mxGetN(COLS_IN1);
+    if (!mxIsDouble(COLS_IN1) || mxIsComplex(COLS_IN1) || (m1!=1 && n1!=1)) {
+      mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double.");
+    }
+    colnum = (m1 > n1) ? m1 : n1;   /* the number of columns in the linear combination */
+    
+    m2 = mxGetM(X_IN1);
+    n2 = mxGetN(X_IN1);
+    if (!mxIsDouble(X_IN1) || mxIsComplex(X_IN1) || (m2!=1 && n2!=1)) {
+      mexErrMsgTxt("COLLINCOMB requires that X be a double vector.");
+    }
+    
+    if (m2!=colnum && n2!=colnum) {
+      mexErrMsgTxt("The length of X does not match the number of columns in COLS.");
+    }
+    
+    rows = 0;
+    Y_OUT = mxCreateDoubleMatrix(m, 1, mxREAL);
+    cols = mxGetPr(COLS_IN1);
+    x = mxGetPr(X_IN1);
+  }
+  else {
+    
+    m1 = mxGetM(ROWS_IN);
+    n1 = mxGetN(ROWS_IN);
+    if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) {
+      mexErrMsgTxt("COLLINCOMB requires that ROWS be an index vector of type double.");
+    }
+    rownum = (m1 > n1) ? m1 : n1;   /* the number of rows in the linear combination */
+    rownumspecified = 1;
+    rows = mxGetPr(ROWS_IN);
+    
+    m1 = mxGetM(COLS_IN2);
+    n1 = mxGetN(COLS_IN2);
+    if (!mxIsDouble(COLS_IN2) || mxIsComplex(COLS_IN2) || (m1!=1 && n1!=1)) {
+      mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double.");
+    }
+    colnum = (m1 > n1) ? m1 : n1;   /* the number of columns in the linear combination */
+    
+    m2 = mxGetM(X_IN2);
+    n2 = mxGetN(X_IN2);
+    if (!mxIsDouble(X_IN2) || mxIsComplex(X_IN2) || (m2!=1 && n2!=1)) {
+      mexErrMsgTxt("COLLINCOMB requires that X be a double vector.");
+    }
+    
+    if (m2!=colnum && n2!=colnum) {
+      mexErrMsgTxt("The length of X does not match the number of columns in COLS.");
+    }
+    
+    Y_OUT = mxCreateDoubleMatrix(rownum, 1, mxREAL);
+    cols = mxGetPr(COLS_IN2);
+    x = mxGetPr(X_IN2);
+  }
+  
+  
+  /* Assign pointers to the various parameters */
+  A = mxGetPr(A_IN);
+  y = mxGetPr(Y_OUT);
+  
+  
+  if (rownumspecified) {
+    
+     /* check row indices */
+    
+    row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex));
+    
+    for (i=0; i<rownum; ++i) {
+      row_ids[i] = (mwIndex)(rows[i]+0.1)-1;
+      if (row_ids[i]<0 || row_ids[i]>=m) {
+        mexErrMsgTxt("Row index in ROWS is out of range.");
+      }
+    }
+    
+    /* Do the actual computation */
+    for (i=0; i<colnum; ++i) {
+      col_id = (mwIndex)(cols[i]+0.1)-1;
+      if (col_id<0 || col_id>=n) {
+        mexErrMsgTxt("Column index in COLS is out of range.");
+      }
+      for (j=0; j<rownum; ++j) {
+        y[j] += A[m*col_id+row_ids[j]]*x[i];
+      }
+    }
+    
+    mxFree(row_ids);
+  }
+  
+  else {
+    
+    /* Do the actual computation */
+    for (i=0; i<colnum; ++i) {
+      col_id = (mwIndex)(cols[i]+0.1)-1;
+      if (col_id<0 || col_id>=n) {
+        mexErrMsgTxt("Column index in COLS is out of range.");
+      }
+      for (j=0; j<m; ++j) {
+        y[j] += A[m*col_id+j]*x[i];
+      }
+    }
+  }
+  
+  return;
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/collincomb.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/collincomb.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,27 @@
+%COLLINCOMB Linear combination of matrix columns.
+%  Y = COLLINCOMB(A,COLS,X) computes a linear combination of the columns of
+%  the matrix A. The column indices are specified in the vector COLS, and
+%  the correspoinding coefficients are specified in the vector X. The
+%  vectors COLS and X must be of the same length. 
+%  
+%  The call Y = COLLINCOMB(A,COLS,X) is essentially equivalent to
+%
+%         Y = A(:,COLS)*X .
+%
+%  However, it is implemented much more efficiently.
+%
+%  Y = COLLINCOMB(A,ROWS,COLS,X) only works on the rows of A specified
+%  in ROWS, returning a vector of length equal to ROWS. This call is
+%  essentially equivalent to the command
+%
+%         Y = A(ROWS,COLS)*X .
+%
+%  See also ROWLINCOMB.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/collincomb.mexa64
Binary file DL/RLS-DLA/private/collincomb.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/completeOps.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/completeOps.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,78 @@
+function data = completeOps(data)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: completeOps.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+operators = {};
+flagM = 0; if isfield(data,'M'), flagM = 1; end;
+flagB = 0; if isfield(data,'B'), flagB = 1; end;
+
+if (~flagM) && (~flagB)
+   error('At least one of the operators M or B has be to given.');
+end
+
+% Define measurement matrix if needed
+if ~flagM
+  info   = data.B([],0);
+  data.M = opDirac(info{1});
+else
+  operators{end+1} = data.M;
+end
+
+% Define sparsity basis if needed
+if ~flagB
+  info   = data.M([],0);
+  data.B = opDirac(info{2});
+else
+  operators{end+1} = data.B;
+end
+
+% Define operator A if needed
+if ~isfield(data,'A')
+  if (length(operators) > 1)
+     data.A = opFoG(operators{:});
+  else
+     data.A = operators{1};
+  end
+end
+
+% Define empty solution if needed
+if ~isfield(data,'x0')
+  data.x0 = [];
+end
+
+% Define the operator size and string
+opInfo       = data.A([],0);
+data.sizeA   = [opInfo{1},opInfo{2}];
+opInfo       = data.B([],0);
+data.sizeB   = [opInfo{1},opInfo{2}];
+opInfo       = data.M([],0);
+data.sizeM   = [opInfo{1},opInfo{2}];
+data.op.strA = opToString(data.A);
+data.op.strB = opToString(data.B);
+data.op.strM = opToString(data.M);
+
+% Get the size of the desired signal
+if ~isfield(data,'signalSize')
+   if ~isfield(data,'signal')
+     error(['At least one of the fields signal ', ...
+            'or signalSize must be given.']);
+   end
+   data.signalSize = size(data.signal);
+end
+
+% Reconstruct signal from sparse coefficients
+if ~isfield(data,'reconstruct')
+   data.reconstruct = @(x) reshape(data.B(x,1),data.signalSize);
+end
+
+% Reorder the fields (sort alphabetically)
+m = fieldnames(data);
+m = sort(m);
+dataReorder = struct();
+for i=1:length(m)
+  eval(sprintf('dataReorder.%s = data.%s;',m{i},m{i}));
+end
+
+data = dataReorder;
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/countcover.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/countcover.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,33 @@
+function cnt = countcover(sz,blocksize,stepsize)
+%COUNTCOVER Covering of signal samples by blocks
+%  CNT = COUNTCOVER(SZ,BLOCKSIZE,STEPSIZE) assumes a p-dimensional signal
+%  of size SZ=[N1 N2 ... Np] covered by (possibly overlapping) blocks of
+%  size BLOCKSIZE=[M1 M2 ... Mp]. The blocks start at position (1,1,..,1)
+%  and are shifted between them by steps of size STEPSIZE=[S1 S2 ... Sp].
+%  COUNTCOVER returns a matrix the same size as the signal, containing in
+%  each entry the number of blocks covering that sample.
+%
+%  See also IM2COLSTEP, COL2IMSTEP, IM2COL.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2008
+
+
+cnt = ones(sz);
+for k = 1:length(sz)
+  
+  % this code is modified from function NDGRID, so it computes one
+  % output argument of NDGRID at a time (to conserve memory)
+  ids = (1:sz(k))';
+  s = sz; s(k) = [];
+  ids = reshape(ids(:,ones(1,prod(s))),[length(ids) s]);
+  ids = permute(ids,[2:k 1 k+1:length(sz)]);
+  
+  cnt = cnt .* max( min(floor((ids-1)/stepsize(k)),floor((sz(k)-blocksize(k))/stepsize(k))) - ...
+                    max(ceil((ids-blocksize(k))/stepsize(k)),0) + 1 , 0 );
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/diag_ids.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/diag_ids.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,30 @@
+function id = diag_ids(x,k)
+%DIAG_IDS Indices of matrix diagonal elements.
+%  ID = DIAG_IDS(X) returns the indices of the main diagonal of X.
+%
+%  ID = DIAG_IDS(X,K) returns the indices of the K-th diagonal. K=0
+%  represents the main diagonal, positive values are above the main
+%  diagonal and negative values are below the main diagonal.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  September 2006
+
+
+if (nargin==1), k=0; end
+
+[m,n] = size(x);
+l = max(m,n);
+
+if (k<=0)
+  id = (0:l-1)*m + (1:l) - k;
+else
+  id = (0:l-1)*m + (1:l) + k*m;
+end
+
+if (l-k>m), id = id(1:end-(l-k-m)); end
+if (l+k>n), id = id(1:end-(l+k-n)); end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/dictdist.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/dictdist.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,61 @@
+function [dist,ratio] = dictdist(approxD,D,epsilon)
+%DICTDIST Distance between dictionaries.
+%  [DIST,RATIO] = DICTDIST(APPROXD,D) computes the distance between the
+%  approximate dictionary APPROXD and the true dictionary D, where APPROXD
+%  is NxK and D is NxM.
+%
+%  The distance between the dictionary APPROXD and a single atom A of D is
+%  defined as:
+%
+%      DIST(APPROXD,A) = min  { 1-abs(APPROXD(:,i)' * A) }
+%                         i
+%
+%  The distance between the dictionaries APPROXD and D is defined as:
+%
+%      DIST(APPROXD,D) = sum { dist(APPROXD, D(:,k)) } / M
+%                         k
+%
+%  Note that 0 <= DIST(APPROXD,D) <= 1, where 0 implies that all atoms in D
+%  appear in APPROXD, and 1 implies that the atoms of D are orthogonal to
+%  APPROXD.
+%
+%  The similarity ratio between APPROXD and D is defined as:
+%
+%      RATIO(APPROXD,D) = #(atoms in D that appear in APPROXD) / M
+%
+%  where two atoms are considered identical when DIST(A1,A2) < EPSILON with
+%  EPSILON=0.01 by default. Note that 0 <= RATIO(APPROXD,D) <= 1, where 0
+%  means APPROXD and D have no identical atoms, and 1 means that all atoms
+%  of D appear in APPROXD.
+%
+%  [DIST,RATIO] = DICTDIST(DICT1,DICT2,EPSILON) specifies a different value
+%  for EPSILON.
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  October 2007
+
+
+if (nargin < 3), epsilon = 0.01; end
+
+[n,m] = size(D);
+
+approxD = normcols(approxD*spdiag(sign(approxD(1,:))));
+D = normcols(D*spdiag(sign(D(1,:))));
+
+identical_atoms = 0;
+dist = 0;
+
+for i = 1:m
+  atom = D(:,i);
+  distances = 1-abs(atom'*approxD);
+  mindist = min(distances);
+  dist = dist + mindist;
+  identical_atoms = identical_atoms + (mindist < epsilon);
+end
+
+dist = dist / m;
+ratio = identical_atoms / m;
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/genPDF.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/genPDF.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,113 @@
+function [pdf,val] = genPDF(imSize,p,pctg,distType,radius,disp)
+
+%[pdf,val] = genPDF(imSize,p,pctg [,distType,radius,disp])
+%
+%	generates a pdf for a 1d or 2d random sampling pattern
+%	with polynomial variable density sampling
+%
+%	Input:
+%		imSize - size of matrix or vector
+%		p - power of polynomial
+%		pctg - partial sampling factor e.g. 0.5 for half
+%		distType - 1 or 2 for L1 or L2 distance measure
+%		radius - radius of fully sampled center
+%		disp - display output
+%
+%	Output:
+%		pdf - the pdf
+%		val - min sampling density
+%
+% 
+%	Example:
+%	[pdf,val] = genPDF([128,128],2,0.5,2,0,1);
+%
+%	(c) Michael Lustig 2007
+
+%   This file is used with the kind permission of Michael Lustig
+%   (mlustig@stanford.edu), and originally appeared in the
+%   SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
+%
+% $Id: genPDF.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+
+if nargin < 4
+	distType = 2;
+end
+
+if nargin < 5
+	radius = 0;
+end
+
+if nargin < 6
+	disp = 0;
+end
+
+
+minval=0;
+maxval=1;
+val = 0.5;
+
+if length(imSize)==1
+	imSize = [imSize,1];
+end
+
+sx = imSize(1);
+sy = imSize(2);
+PCTG = floor(pctg*sx*sy);
+
+
+if sum(imSize==1)==0  % 2D
+	[x,y] = meshgrid(linspace(-1,1,sy),linspace(-1,1,sx));
+	switch distType
+		case 1
+			r = max(abs(x),abs(y));
+		otherwise
+			r = sqrt(x.^2+y.^2);
+			r = r/max(abs(r(:)));			
+	end
+
+else %1d
+	r = abs(linspace(-1,1,max(sx,sy)));
+end
+
+
+
+
+idx = find(r<radius);
+
+pdf = (1-r).^p; pdf(idx) = 1;
+if floor(sum(pdf(:))) > PCTG
+	error('infeasible without undersampling dc, increase p');
+end
+
+% begin bisection
+while(1)
+	val = minval/2 + maxval/2;
+	pdf = (1-r).^p + val; pdf(find(pdf>1)) = 1; pdf(idx)=1;
+	N = floor(sum(pdf(:)));
+	if N > PCTG % infeasible
+		maxval=val;
+	end
+	if N < PCTG % feasible, but not optimal
+		minval=val;
+	end
+	if N==PCTG % optimal
+		break;
+	end
+end
+
+if disp
+	figure, 
+	subplot(211), imshow(pdf)
+	if sum(imSize==1)==0
+		subplot(212), plot(pdf(end/2+1,:));
+	else
+		subplot(212), plot(pdf);
+	end
+end
+
+
+
+
+
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/genSampling.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/genSampling.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,53 @@
+function [minIntrVec,stat,actpctg] = genSampling(pdf,iter,tol)
+
+%[mask,stat,N] = genSampling(pdf,iter,tol)
+%
+% a monte-carlo algorithm to generate a sampling pattern with 
+% minimum peak interference. The number of samples will be
+% sum(pdf) +- tol
+%
+%	pdf - probability density function to choose samples from
+%	iter - number of tries
+%	tol  - the deviation from the desired number of samples in samples
+%
+% returns:
+%	mask - sampling pattern
+%	stat - vector of min interferences measured each try
+%	actpctg    - actual undersampling factor
+%
+%	(c) Michael Lustig 2007
+
+%   This file is used with the kind permission of Michael Lustig
+%   (mlustig@stanford.edu), and originally appeared in the
+%   SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
+%
+% $Id: genSampling.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+% h = waitbar(0);
+
+pdf(find(pdf>1)) = 1;
+K = sum(pdf(:));
+
+minIntr = 1e99;
+minIntrVec = zeros(size(pdf));
+
+for n=1:iter
+	tmp = zeros(size(pdf));
+	while abs(sum(tmp(:)) - K) > tol
+		tmp = rand(size(pdf))<pdf;
+	end
+	
+	TMP = ifft2(tmp./pdf);
+	if max(abs(TMP(2:end))) < minIntr
+		minIntr = max(abs(TMP(2:end)));
+		minIntrVec = tmp;
+	end
+	stat(n) = max(abs(TMP(2:end)));
+      % waitbar(n/iter,h);
+end
+
+actpctg = sum(minIntrVec(:))/prod(size(minIntrVec));
+
+% close(h);
+
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/im2colstep.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/im2colstep.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,131 @@
+/**************************************************************************
+ *
+ * File name: im2colstep.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 31.8.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+#include <string.h>
+
+
+/* Input Arguments */
+
+#define	X_IN	 prhs[0]
+#define SZ_IN  prhs[1]
+#define S_IN   prhs[2]
+
+
+/* Output Arguments */
+
+#define	B_OUT	plhs[0]
+
+
+void mexFunction(int nlhs, mxArray *plhs[], 
+		             int nrhs, const mxArray*prhs[])
+     
+{ 
+    double *x, *b, *s;
+    mwSize sz[3], stepsize[3], n[3], ndims;
+    mwIndex i, j, k, l, m, blocknum;
+    
+    
+    /* Check for proper number of arguments */
+    
+    if (nrhs < 2 || nrhs > 3) {
+      mexErrMsgTxt("Invalid number of input arguments."); 
+    } else if (nlhs > 1) {
+      mexErrMsgTxt("Too many output arguments."); 
+    } 
+    
+    
+    /* Check the the input dimensions */ 
+    
+    ndims = mxGetNumberOfDimensions(X_IN);
+    
+    if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ndims>3) {
+      mexErrMsgTxt("X should be a 2-D or 3-D double matrix.");
+    }
+    if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    if (nrhs == 3) {
+      if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+    }
+    
+    
+    /* Get parameters */
+    
+    s = mxGetPr(SZ_IN);
+    if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    sz[0] = (mwSize)(s[0] + 0.01);
+    sz[1] = (mwSize)(s[1] + 0.01);
+    sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+    
+    if (nrhs == 3) {
+      s = mxGetPr(S_IN);
+      if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+      stepsize[0] = (mwSize)(s[0] + 0.01);
+      stepsize[1] = (mwSize)(s[1] + 0.01);
+      stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+    }
+    else {
+      stepsize[0] = stepsize[1] = stepsize[2] = 1;
+    }
+    
+    n[0] = (mxGetDimensions(X_IN))[0];
+    n[1] = (mxGetDimensions(X_IN))[1];
+    n[2] = ndims==3 ? (mxGetDimensions(X_IN))[2] : 1;
+    
+    if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) {
+      mexErrMsgTxt("Block size too large.");
+    }
+    
+    
+    /* Create a matrix for the return argument */
+    
+    B_OUT = mxCreateDoubleMatrix(sz[0]*sz[1]*sz[2], ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1), mxREAL);
+    
+    
+    /* Assign pointers */
+    
+    x = mxGetPr(X_IN);
+    b = mxGetPr(B_OUT);
+            
+    
+    /* Do the actual computation */
+    
+    blocknum = 0;
+    
+    /* iterate over all blocks */
+    for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) {
+      for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) {
+        for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) {
+          
+          /* copy single block */
+          for (m=0; m<sz[2]; m++) {
+            for (l=0; l<sz[1]; l++) {
+              memcpy(b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0], x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i, sz[0]*sizeof(double));
+            }
+          }
+          blocknum++;
+          
+        }
+      }
+    }
+    
+    return;
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/im2colstep.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/im2colstep.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,31 @@
+%IM2COLSTEP Rearrange matrix blocks into columns.
+%  B = IM2COLSTEP(A,[N1 N2]) converts each sliding N1-by-N2 block of the
+%  2-D matrix A into a column of B, with no zero padding. B has N1*N2 rows
+%  and will contain as many columns as there are N1-by-N2 neighborhoods in
+%  A. Each column of B contains a neighborhood of A reshaped as NHOOD(:),
+%  where NHOOD is a matrix containing an N1-by-N2 neighborhood of A.
+%
+%  B = IM2COLSTEP(A,[N1 N2],[S1 S2]) extracts neighborhoods of A with a
+%  step size of (S1,S2) between them. The first extracted neighborhood is
+%  A(1:N1,1:N2), and the rest are of the form A((1:N1)+i*S1,(1:N2)+j*S2).
+%  Note that to ensure coverage of all A by neighborhoods,
+%  (size(A,i)-Ni)/Si must be whole for i=1,2. The default function behavior
+%  corresponds to [S1 S2] = [1 1]. Setting S1>=N1 and S2>=N2 results in no
+%  overlap between the neighborhoods.
+%
+%  B = IM2COLSTEP(A,[N1 N2 N3],[S1 S2 S3]) operates on a 3-D matrix A. The
+%  step size [S1 S2 S3] may be ommitted, and defaults to [1 1 1].
+%
+%  Note: the call IM2COLSTEP(A,[N1 N2]) produces the same output as
+%  Matlab's IM2COL(A,[N1 N2],'sliding'). However, it is significantly
+%  faster.
+%
+%  See also COL2IMSTEP, IM2COL, COUNTCOVER.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/im2colstep.mexa64
Binary file DL/RLS-DLA/private/im2colstep.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/imnormalize.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/imnormalize.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,19 @@
+function y = imnormalize(x)
+%IMNORMALIZE Normalize image values.
+%  Y = IMNORMALIZE(X) linearly transforms the intensity values of the image
+%  X to tightly cover the range [0,1]. If X has more than one channel, the
+%  channels are handled as one and normalized using the same transform.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  May 2004
+
+
+maxval = max(x(:));
+minval = min(x(:));
+
+y = (x-minval) / (maxval-minval);
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/iswhole.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/iswhole.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,17 @@
+function z = iswhole(x,epsilon)
+%ISWHOLE Determine whole numbers with finite precision.
+%  Z = ISWHOLE(X,EPSILON) returns a matrix the same size as X, containing
+%  1's where X is whole up to precision EPSILON, and 0's elsewhere. 
+%
+%  Z = ISWHOLE(X) uses the default value of EPSILON=1e-6.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  May 2005
+
+if (nargin<2), epsilon = 1e-6; end
+z = abs(round(x)-x) < epsilon;
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/make.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/make.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,40 @@
+function make
+%MAKE Build the KSVDBox MEX support files.
+%  MAKE compiles the KSVDBox supporting MEX functions, using Matlab's
+%  default MEX compiler. If the MEX compiler has not been set-up before,
+%  please run 
+%
+%    mex -setup
+%
+%  before using this MAKE file.
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
+
+
+% detect platform 
+
+compstr = computer;
+is64bit = strcmp(compstr(end-1:end),'64');
+
+
+% compilation parameters
+
+compile_params = cell(0);
+if (is64bit)
+  compile_params{1} = '-largeArrayDims';
+end
+
+
+% Compile files %
+
+sourcefiles = {{'addtocols.c'}, {'collincomb.c'}, {'rowlincomb.c'}, {'sprow.c','mexutils.c'}, {'im2colstep.c'}, {'col2imstep.c'}};
+
+for i = 1:length(sourcefiles)
+  printf('Compiling %s...', sourcefiles{i}{1});
+  mex(sourcefiles{i}{:},compile_params{:});
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/mexutils.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/mexutils.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,79 @@
+/**************************************************************************
+ *
+ * File name: mexutils.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 15.8.2009
+ *
+ *************************************************************************/
+
+#include "mexutils.h"
+#include <math.h>
+
+
+
+/* verify that the mxArray contains a double matrix */
+
+void checkmatrix(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-D double vector */
+
+void checkvector(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a double scalar */
+
+void checkscalar(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || 
+      mxGetM(param)!=1 || mxGetN(param)!=1) 
+  {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a sparse matrix */
+
+void checksparse(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname);
+  if (!mxIsSparse(param)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-dimensional cell array */
+
+void checkcell_1d(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname);
+  if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/mexutils.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/mexutils.h	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,103 @@
+/**************************************************************************
+ *
+ * File name: mexutils.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Utility functions for MEX files.
+ *
+ *************************************************************************/
+
+
+#ifndef __MEX_UTILS_H__
+#define __MEX_UTILS_H__
+
+#include "mex.h"
+
+
+
+/**************************************************************************
+ * Function checkmatrix:
+ *
+ * Verify that the specified mxArray is real, of type double, and has 
+ * no more than two dimensions. If not, an error message is printed
+ * and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkmatrix(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkvector:
+ *
+ * Verify that the specified mxArray is 1-D, real, and of type double. The
+ * vector may be a column or row vector. Otherwise, an error message is
+ * printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkvector(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkscalar:
+ *
+ * Verify that the specified mxArray represents a real double scalar value. 
+ * If not, an error message is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkscalar(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checksparse:
+ *
+ * Verify that the specified mxArray contains a sparse matrix. If not,
+ * an error message is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checksparse(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkcell_1d:
+ *
+ * Verify that the specified mxArray is a 1-D cell array. The cell array 
+ * may be arranged as either a column or a row. If not, an error message 
+ * is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkcell_1d(const mxArray *param, char *fname, char *pname);
+
+
+#endif
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/myblas.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/myblas.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,663 @@
+/**************************************************************************
+ *
+ * File name: myblas.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Version: 1.1
+ * Last updated: 13.8.2009
+ *
+ *************************************************************************/
+
+
+#include "myblas.h"
+#include <ctype.h>
+
+
+/* find maximum of absolute values */
+
+mwIndex maxabs(double c[], mwSize m)
+{
+  mwIndex maxid=0, k;
+  double absval, maxval = SQR(*c);   /* use square which is quicker than absolute value */
+
+  for (k=1; k<m; ++k) {
+    absval = SQR(c[k]);
+    if (absval > maxval) {
+      maxval = absval;
+      maxid = k;
+    }
+  }
+  return maxid;
+}
+
+
+/* compute y := alpha*x + y */
+
+void vec_sum(double alpha, double x[], double y[], mwSize n)
+{
+  mwIndex i;
+
+  for (i=0; i<n; ++i) {
+    y[i] += alpha*x[i];
+  }
+}
+
+
+/* compute y := alpha*A*x */
+
+void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m)
+{
+  mwIndex i, j, i_n;
+  double *Ax;
+
+  Ax = mxCalloc(n,sizeof(double));
+
+  for (i=0; i<m; ++i) {
+    i_n = i*n;
+    for (j=0; j<n; ++j) {
+      Ax[j] += A[i_n+j] * x[i];
+    }
+  }
+
+  for (j=0; j<n; ++j) {
+    y[j] = alpha*Ax[j];
+  }
+
+  mxFree(Ax);
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m)
+{
+  mwIndex i, j, n_i;
+  double sum0, sum1, sum2, sum3;
+
+  for (j=0; j<m; ++j) {
+    y[j] = 0;
+  }
+
+  /* use loop unrolling to accelerate computation */
+
+  for (i=0; i<m; ++i) {
+    n_i = n*i;
+    sum0 = sum1 = sum2 = sum3 = 0;
+    for (j=0; j+4<n; j+=4) {
+      sum0 += A[n_i+j]*x[j];
+      sum1 += A[n_i+j+1]*x[j+1];
+      sum2 += A[n_i+j+2]*x[j+2];
+      sum3 += A[n_i+j+3]*x[j+3];
+    }
+    y[i] += alpha * ((sum0 + sum1) + (sum2 + sum3));
+    while (j<n) {
+      y[i] += alpha*A[n_i+j]*x[j];
+      j++;
+    }
+  }
+}
+
+
+/* compute y := alpha*A*x */
+
+void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j1, j2;
+
+  for (i=0; i<n; ++i) {
+    y[i] = 0;
+  }
+  
+  j2 = jc[0];
+  for (i=0; i<m; ++i) {
+    j1 = j2; j2 = jc[i+1];
+    for (j=j1; j<j2; ++j) {
+      y[ir[j]] += alpha * pr[j] * x[i];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j1, j2;
+  
+  for (i=0; i<m; ++i) {
+    y[i] = 0;
+  }
+  
+  j2 = jc[0];
+  for (i=0; i<m; ++i) {
+    j1 = j2; j2 = jc[i+1];
+    for (j=j1; j<j2; ++j) {
+      y[i] += alpha * pr[j] * x[ir[j]];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A*x */
+
+void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j_n, k, kend;
+  
+  for (i=0; i<n; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jc[1];
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (k=0; k<kend; ++k) {
+    j = ir[k];
+    j_n = j*n;
+    for (i=0; i<n; ++i) {
+      y[i] += alpha * A[i+j_n] * pr[k];
+    }
+  }
+
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j_n, k, kend;
+  
+  for (i=0; i<m; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jc[1];
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (j=0; j<m; ++j) {
+    j_n = j*n;
+    for (k=0; k<kend; ++k) {
+      i = ir[k];
+      y[j] += alpha * A[i+j_n] * pr[k];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A*x */
+
+void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, k, kend, j1, j2;
+
+  for (i=0; i<n; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jcx[1]; 
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (k=0; k<kend; ++k) {
+    i = irx[k];
+    j1 = jc[i]; j2 = jc[i+1];
+    for (j=j1; j<j2; ++j) {
+      y[ir[j]] += alpha * pr[j] * prx[k];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, k, jend, kend, jadd, kadd, delta;
+  
+  for (i=0; i<m; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jcx[1];
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (i=0; i<m; ++i) {
+    j = jc[i]; 
+    jend = jc[i+1];
+    k = 0;
+    while (j<jend && k<kend) {
+      
+      delta = ir[j] - irx[k];
+      
+      if (delta) { /* if indices differ - increment the smaller one */
+        jadd = delta<0;
+        kadd = 1-jadd;
+        j += jadd;
+        k += kadd;
+      }
+      
+      else {    /* indices are equal - add to result and increment both */
+        y[i] += alpha * pr[j] * prx[k];
+        j++; k++;
+      }
+    }
+  }
+  
+}
+
+
+/* matrix-matrix multiplication */
+
+void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
+{
+  mwIndex i1, i2, i3, iX, iA, i2_n;
+  double b;
+  
+  for (i1=0; i1<n*k; i1++) {
+    X[i1] = 0;
+  }
+
+  for (i2=0; i2<m; ++i2) {
+    i2_n = i2*n;
+    iX = 0;
+    for (i3=0; i3<k; ++i3) {
+      iA = i2_n;
+      b = B[i2+i3*m];
+      for (i1=0; i1<n; ++i1) {
+        X[iX++] += A[iA++]*b;
+      }
+    }
+  }
+  
+  for (i1=0; i1<n*k; i1++) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* matrix-transpose-matrix multiplication */
+
+void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
+{
+  mwIndex i1, i2, i3, iX, iA, i2_n;
+  double *x, sum0, sum1, sum2, sum3;
+
+  for (i2=0; i2<m; ++i2) {
+    for (i3=0; i3<k; ++i3) {
+      sum0 = sum1 = sum2 = sum3 = 0;
+      for (i1=0; i1+4<n; i1+=4) {
+        sum0 += A[i1+0+i2*n]*B[i1+0+i3*n];
+        sum1 += A[i1+1+i2*n]*B[i1+1+i3*n];
+        sum2 += A[i1+2+i2*n]*B[i1+2+i3*n];
+        sum3 += A[i1+3+i2*n]*B[i1+3+i3*n];
+      }
+      X[i2+i3*m] = (sum0+sum1) + (sum2+sum3);
+      while(i1<n) {
+        X[i2+i3*m] += A[i1+i2*n]*B[i1+i3*n];
+        i1++;
+      }
+    }
+  }
+  
+  for (i1=0; i1<m*k; i1++) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* tensor-matrix product */
+
+void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l)
+{
+  mwIndex i1, i2, i3, i4, i2_n, nml;
+  double b;
+  
+  nml = n*m*l;
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] = 0;
+  }
+
+  for (i2=0; i2<m; ++i2) {
+    i2_n = i2*n;
+    for (i3=0; i3<k; ++i3) {
+      for (i4=0; i4<l; ++i4) {
+        b = B[i4+i3*l];
+        for (i1=0; i1<n; ++i1) {
+          X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b;
+        }
+      }
+    }
+  }
+  
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* tensor-matrix-transpose product */
+
+void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l)
+{
+  mwIndex i1, i2, i3, i4, i2_n, nml;
+  double b;
+  
+  nml = n*m*l;
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] = 0;
+  }
+
+  for (i2=0; i2<m; ++i2) {
+    i2_n = i2*n;
+    for (i4=0; i4<l; ++i4) {
+      for (i3=0; i3<k; ++i3) {
+        b = B[i3+i4*k];
+        for (i1=0; i1<n; ++i1) {
+          X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b;
+        }
+      }
+    }
+  }
+  
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* dot product */
+
+double dotprod(double a[], double b[], mwSize n)
+{
+  double sum = 0;
+  mwIndex i;
+  for (i=0; i<n; ++i)
+    sum += a[i]*b[i];
+  return sum;
+}
+
+
+/* find maximum of vector */
+
+mwIndex maxpos(double c[], mwSize m)
+{
+  mwIndex maxid=0, k;
+  double val, maxval = *c;
+
+  for (k=1; k<m; ++k) {
+    val = c[k];
+    if (val > maxval) {
+      maxval = val;
+      maxid = k;
+    }
+  }
+  return maxid;
+}
+
+
+/* solve L*x = b */
+
+void backsubst_L(double L[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=0; i<k; ++i) {
+    rhs = b[i];
+    for (j=0; j<i; ++j) {
+      rhs -= L[j*n+i]*x[j];
+    }
+    x[i] = rhs/L[i*n+i];
+  }
+}
+
+
+/* solve L'*x = b */
+
+void backsubst_Lt(double L[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=k; i>=1; --i) {
+    rhs = b[i-1];
+    for (j=i; j<k; ++j) {
+      rhs -= L[(i-1)*n+j]*x[j];
+    }
+    x[i-1] = rhs/L[(i-1)*n+i-1];
+  }
+}
+
+
+/* solve U*x = b */
+
+void backsubst_U(double U[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=k; i>=1; --i) {
+    rhs = b[i-1];
+    for (j=i; j<k; ++j) {
+      rhs -= U[j*n+i-1]*x[j];
+    }
+    x[i-1] = rhs/U[(i-1)*n+i-1];
+  }
+}
+
+
+/* solve U'*x = b */
+
+void backsubst_Ut(double U[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=0; i<k; ++i) {
+    rhs = b[i];
+    for (j=0; j<i; ++j) {
+      rhs -= U[i*n+j]*x[j];
+    }
+    x[i] = rhs/U[i*n+i];
+  }
+}
+
+
+/* back substitution solver */
+
+void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k)
+{
+  if (tolower(ul) == 'u') {
+    backsubst_U(A, b, x, n, k);
+  }
+  else if (tolower(ul) == 'l') {
+    backsubst_L(A, b, x, n, k);
+  }
+  else {
+    mexErrMsgTxt("Invalid triangular matrix type: must be ''U'' or ''L''");
+  }
+}
+
+
+/* solve equation set using cholesky decomposition */
+
+void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k)
+{
+  double *tmp;
+
+  tmp = mxMalloc(k*sizeof(double));
+
+  if (tolower(ul) == 'l') {
+    backsubst_L(A, b, tmp, n, k);
+    backsubst_Lt(A, tmp, x, n, k);
+  }
+  else if (tolower(ul) == 'u') {
+    backsubst_Ut(A, b, tmp, n, k);
+    backsubst_U(A, tmp, x, n, k);
+  }
+  else {
+    mexErrMsgTxt("Invalid triangular matrix type: must be either ''U'' or ''L''");
+  }
+
+  mxFree(tmp);
+}
+
+
+/* perform a permutation assignment y := x(ind(1:k)) */
+
+void vec_assign(double y[], double x[], mwIndex ind[], mwSize k)
+{
+  mwIndex i;
+
+  for (i=0; i<k; ++i)
+    y[i] = x[ind[i]];
+}
+
+
+/* matrix transpose */
+
+void transpose(double X[], double Y[], mwSize n, mwSize m)
+{
+  mwIndex i, j, i_m, j_n;
+  
+  if (n<m) {
+    for (j=0; j<m; ++j) {
+      j_n = j*n;
+      for (i=0; i<n; ++i) {
+        Y[j+i*m] = X[i+j_n];
+      }
+    }
+  }
+  else {
+    for (i=0; i<n; ++i) {
+      i_m = i*m;
+      for (j=0; j<m; ++j) {
+        Y[j+i_m] = X[i+j*n];
+      }
+    }
+  }
+}
+
+
+/* print contents of matrix */
+
+void printmat(double A[], int n, int m, char* matname)
+{
+  int i, j;
+  mexPrintf("\n%s = \n\n", matname);
+
+  if (n*m==0) {
+    mexPrintf("   Empty matrix: %d-by-%d\n\n", n, m);
+    return;
+  }
+
+  for (i=0; i<n; ++i) {
+    for (j=0; j<m; ++j)
+      mexPrintf("   %lf", A[j*n+i]);
+    mexPrintf("\n");
+  }
+  mexPrintf("\n");
+}
+
+
+/* print contents of sparse matrix */
+
+void printspmat(mxArray *a, char* matname)
+{
+  mwIndex *aJc = mxGetJc(a);
+  mwIndex *aIr = mxGetIr(a);
+  double *aPr = mxGetPr(a);
+
+  int i;
+
+  mexPrintf("\n%s = \n\n", matname);
+
+  for (i=0; i<aJc[1]; ++i)
+    printf("   (%d,1) = %lf\n", aIr[i]+1,aPr[i]);
+
+  mexPrintf("\n");
+}
+
+
+
+/* matrix multiplication using Winograd's algorithm */
+
+/*
+void mat_mat2(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
+{
+  
+  mwIndex i1, i2, i3, iX, iA, i2_n;
+  double b, *AA, *BB;
+  
+  AA = mxCalloc(n,sizeof(double));
+  BB = mxCalloc(k,sizeof(double));
+  
+  for (i1=0; i1<n*k; i1++) {
+    X[i1] = 0;
+  }
+  
+  for (i1=0; i1<n; ++i1) {
+    for (i2=0; i2<m/2; ++i2) {
+      AA[i1] += A[i1+2*i2*n]*A[i1+(2*i2+1)*n];
+    }
+  }
+
+  for (i2=0; i2<k; ++i2) {
+    for (i1=0; i1<m/2; ++i1) {
+      BB[i2] += B[2*i1+i2*m]*B[2*i1+1+i2*m];
+    }
+  }
+
+  for (i2=0; i2<k; ++i2) {
+    for (i3=0; i3<m/2; ++i3) {
+      for (i1=0; i1<n; ++i1) {
+        X[i1+i2*n] += (A[i1+(2*i3)*n]+B[2*i3+1+i2*m])*(A[i1+(2*i3+1)*n]+B[2*i3+i2*m]);
+      }
+    }
+  }
+  
+  if (m%2) {
+    for (i2=0; i2<k; ++i2) {
+      for (i1=0; i1<n; ++i1) {
+        X[i1+i2*n] += A[i1+(m-1)*n]*B[m-1+i2*m];
+      }
+    }
+  }
+  
+  for (i2=0; i2<k; ++i2) {
+    for (i1=0; i1<n; ++i1) {
+      X[i1+i2*n] -= (AA[i1] + BB[i2]);
+      X[i1+i2*n] *= alpha;
+    }
+  }
+  
+  mxFree(AA);
+  mxFree(BB);
+}
+*/
+
+
+
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/myblas.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/myblas.h	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,492 @@
+/**************************************************************************
+ *
+ * File name: myblas.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Version: 1.1
+ * Last updated: 17.8.2009
+ *
+ * A collection of basic linear algebra functions, in the spirit of the
+ * BLAS/LAPACK libraries.
+ *
+ *************************************************************************/
+
+
+
+#ifndef __MY_BLAS_H__
+#define __MY_BLAS_H__
+
+
+#include "mex.h"
+#include <math.h>
+
+
+
+/**************************************************************************
+ * Squared value.
+ **************************************************************************/
+#define SQR(X) ((X)*(X))
+
+
+
+/**************************************************************************
+ * Matrix-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   x - vector of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-transpose-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   x - vector of length n
+ *   y - output vector of length m
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * where A is a sparse matrix.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   x - vector of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-transpose-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * where A is a sparse matrix.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   x - vector of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * where A is a matrix and x is a sparse vector.
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   pr,ir,jc - sparse representation of the vector x, of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-transpose-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * where A is a matrix and x is a sparse vector.
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   pr,ir,jc - sparse representation of the vector x, of length n
+ *   y - output vector of length m
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * where A is a sparse matrix and x is a sparse vector.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   prx,irx,jcx - sparse representation of the vector x (of length m)
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-transpose-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * where A is a sparse matrix and x is a sparse vector.
+ *
+ * Importnant note: this function is provided for completeness, but is NOT efficient.
+ * If possible, convert x to non-sparse representation and use matT_vec_sp instead.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   prx,irx,jcx - sparse representation of the vector x (of length n)
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-matrix multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   X := alpha*A*B
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   B - matrix of size m X k
+ *   X - output matrix of size n X k
+ *   alpha - real constant
+ *   n, m, k - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k);
+
+
+
+/**************************************************************************
+ * Matrix-transpose-matrix multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   X := alpha*A*B
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   B - matrix of size m X k
+ *   X - output matrix of size n X k
+ *   alpha - real constant
+ *   n, m, k - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k);
+
+
+
+/**************************************************************************
+ * Tensor-matrix multiplication. 
+ *
+ * This function accepts a 3-D tensor A of size n X m X k
+ * and a 2-D matrix B of size l X k.
+ * The function computes the 3-D tensor X of size n X m X l, where
+ *
+ *   X(i,j,:) = B*A(i,j,:)
+ *
+ * for all i,j.
+ *
+ * Parameters:
+ *   A - tensor of size n X m X k
+ *   B - matrix of size l X k
+ *   X - output tensor of size n X m X l
+ *   alpha - real constant
+ *   n, m, k, l - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l);
+
+
+
+/**************************************************************************
+ * Tensor-matrix-transpose multiplication. 
+ *
+ * This function accepts a 3-D tensor A of size n X m X k
+ * and a 2-D matrix B of size k X l.
+ * The function computes the 3-D tensor X of size n X m X l, where
+ *
+ *   X(i,j,:) = B'*A(i,j,:)
+ *
+ * for all i,j.
+ *
+ * Parameters:
+ *   A - tensor of size n X m X k
+ *   B - matrix of size k X l
+ *   X - output tensor of size n X m X l
+ *   alpha - real constant
+ *   n, m, k, l - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l);
+
+
+
+/**************************************************************************
+ * Vector-vector sum.
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*x + y
+ *
+ * Parameters:
+ *   x - vector of length n
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n - length of x,y
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void vec_sum(double alpha, double x[], double y[], mwSize n);
+
+
+
+/**************************************************************************
+ * Triangular back substitution.
+ *
+ * Solve the set of linear equations
+ *
+ *   T*x = b
+ *
+ * where T is lower or upper triangular.
+ *
+ * Parameters:
+ *   ul - 'U' for upper triangular, 'L' for lower triangular
+ *   A  - matrix of size n x m containing T
+ *   b  - vector of length k
+ *   x  - output vector of length k
+ *   n  - size of first dimension of A
+ *   k  - the size of the equation set, k<=n,m
+ *
+ * Note:
+ *   The matrix A can be of any size n X m, as long as n,m >= k. 
+ *   Only the lower/upper triangle of the submatrix A(1:k,1:k) defines the
+ *   matrix T (depending on the parameter ul).
+ *
+ **************************************************************************/
+void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k);
+
+
+
+/**************************************************************************
+ * Solve a set of equations using a Cholesky decomposition.
+ *
+ * Solve the set of linear equations
+ *
+ *   M*x = b
+ *
+ * where M is positive definite with a known Cholesky decomposition:
+ * either M=L*L' (L lower triangular) or M=U'*U (U upper triangular).
+ *
+ * Parameters:
+ *   ul - 'U' for upper triangular, 'L' for lower triangular decomposition
+ *   A  - matrix of size n x m with the Cholesky decomposition of M
+ *   b  - vector of length k
+ *   x  - output vector of length k
+ *   n  - size of first dimension of A
+ *   k  - the size of the equation set, k<=n,m
+ *
+ * Note:
+ *   The matrix A can be of any size n X m, as long as n,m >= k. 
+ *   Only the lower/upper triangle of the submatrix A(1:k,1:k) is used as
+ *   the Cholesky decomposition of M (depending on the parameter ul).
+ *
+ **************************************************************************/
+void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k);
+
+
+
+/**************************************************************************
+ * Maximum absolute value.
+ *
+ * Returns the index of the coefficient with maximal absolute value in a vector.
+ *
+ * Parameters:
+ *   x - vector of length n
+ *   n - length of x
+ *
+ **************************************************************************/
+mwIndex maxabs(double x[], mwSize n);
+
+
+
+/**************************************************************************
+ * Maximum vector element.
+ *
+ * Returns the index of the maximal coefficient in a vector.
+ *
+ * Parameters:
+ *   x - vector of length n
+ *   n - length of x
+ *
+ **************************************************************************/
+mwIndex maxpos(double x[], mwSize n);
+
+
+
+/**************************************************************************
+ * Vector-vector dot product.
+ *
+ * Computes an operation of the form:
+ *
+ *   c = a'*b
+ *
+ * Parameters:
+ *   a, b - vectors of length n
+ *   n - length of a,b
+ *
+ * Returns: The dot product c.
+ *
+ **************************************************************************/
+double dotprod(double a[], double b[], mwSize n);
+
+
+
+/**************************************************************************
+ * Indexed vector assignment.
+ *
+ * Perform a permutation assignment of the form
+ *
+ *   y = x(ind)
+ *
+ * where ind is an array of indices to x.
+ *
+ * Parameters:
+ *   y - output vector of length k
+ *   x - input vector of arbitrary length
+ *   ind - array of indices into x (indices begin at 0)
+ *   k - length of the array ind
+ *
+ **************************************************************************/
+void vec_assign(double y[], double x[], mwIndex ind[], mwSize k);
+
+
+
+/**************************************************************************
+ * Matrix transpose.
+ *
+ * Computes Y := X'
+ *
+ * Parameters:
+ *   X - input matrix of size n X m
+ *   Y - output matrix of size m X n
+ *   n, m - dimensions of X
+ *
+ **************************************************************************/
+void transpose(double X[], double Y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Print a matrix.
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   n, m - dimensions of A
+ *   matname - name of matrix to display
+ *
+ **************************************************************************/
+void printmat(double A[], int n, int m, char* matname);
+
+
+
+/**************************************************************************
+ * Print a sparse matrix.
+ *
+ * Parameters:
+ *   A - sparse matrix of type double
+ *   matname - name of matrix to display
+ *
+ **************************************************************************/
+void printspmat(mxArray *A, char* matname);
+
+
+#endif
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/normcols.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/normcols.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,17 @@
+function y = normcols(x)
+%NORMCOLS Normalize matrix columns.
+%  Y = NORMCOLS(X) normalizes the columns of X to unit length, returning
+%  the result as Y.
+%
+%  See also ADDTOCOLS.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
+
+
+y = x*spdiag(1./sqrt(sum(x.*x)));
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/omp2mex.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/omp2mex.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,156 @@
+/**************************************************************************
+ *
+ * File name: omp2mex.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "ompcore.h"
+#include "omputils.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define	IN_D	        prhs[0]
+#define IN_X          prhs[1]
+#define IN_DtX        prhs[2]
+#define IN_XtX        prhs[3]
+#define IN_G          prhs[4]
+#define IN_EPS        prhs[5]
+#define IN_SPARSE_G   prhs[6]
+#define IN_MSGDELTA   prhs[7]
+#define IN_MAXATOMS   prhs[8]
+#define IN_PROFILE    prhs[9]
+
+
+/* Output Arguments */
+
+#define	GAMMA_OUT     plhs[0]
+
+
+/***************************************************************************************/
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
+
+{
+  double *D, *x, *DtX, *XtX, *G, eps, msgdelta;
+  int gmode, maxatoms, profile;
+  mwSize m, n, L;    /* D is n x m , X is n x L, DtX is m x L */
+
+  
+  /* check parameters */
+  
+  checkmatrix(IN_D, "OMP2", "D");
+  checkmatrix(IN_X, "OMP2", "X");
+  checkmatrix(IN_DtX, "OMP2", "DtX");
+  checkmatrix(IN_XtX, "OMP2", "XtX");
+  checkmatrix(IN_G, "OMP2", "G");
+  
+  checkscalar(IN_EPS, "OMP2", "EPSILON");
+  checkscalar(IN_SPARSE_G, "OMP2", "sparse_g");
+  checkscalar(IN_MSGDELTA, "OMP2", "msgdelta");
+  checkscalar(IN_MAXATOMS, "OMP2", "maxatoms");
+  checkscalar(IN_PROFILE, "OMP2", "profile");
+  
+  
+  /* get parameters */
+  
+  x = D = DtX = XtX = G = 0;
+  
+  if (!mxIsEmpty(IN_D))
+    D = mxGetPr(IN_D);
+  
+  if (!mxIsEmpty(IN_X))
+    x = mxGetPr(IN_X);
+  
+  if (!mxIsEmpty(IN_DtX))
+    DtX = mxGetPr(IN_DtX);
+  
+  if (!mxIsEmpty(IN_XtX))
+    XtX = mxGetPr(IN_XtX);
+  
+  if (!mxIsEmpty(IN_G))
+    G = mxGetPr(IN_G);
+  
+  eps = mxGetScalar(IN_EPS);
+  if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) {
+    gmode = SPARSE_GAMMA;
+  }
+  else {
+    gmode = FULL_GAMMA;
+  }
+  msgdelta = mxGetScalar(IN_MSGDELTA);
+  if (mxGetScalar(IN_MAXATOMS) < -1e-5) {
+    maxatoms = -1;
+  }
+  else {
+    maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2);
+  }
+  profile = (int)(mxGetScalar(IN_PROFILE)+1e-2);
+  
+  
+  /* check sizes */
+  
+  if (D && x) {
+    n = mxGetM(IN_D);
+    m = mxGetN(IN_D);
+    L = mxGetN(IN_X);
+    
+    if (mxGetM(IN_X) != n) {
+      mexErrMsgTxt("D and X have incompatible sizes.");
+    }
+    
+    if (G) {
+      if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+        mexErrMsgTxt("G must be a square matrix.");
+      }
+      if (mxGetN(IN_G) != m) {
+        mexErrMsgTxt("D and G have incompatible sizes.");
+      }
+    }
+  }
+  
+  else if (DtX && XtX) {
+    m = mxGetM(IN_DtX);
+    L = mxGetN(IN_DtX);
+    
+    /* set n to an arbitrary value that is at least the max possible number of selected atoms */
+    
+    if (maxatoms>0) {
+      n = maxatoms;
+    }
+    else {
+      n = m;
+    }
+    
+    if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) {
+      mexErrMsgTxt("DtX and XtX have incompatible sizes.");
+    }
+    
+    if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+      mexErrMsgTxt("G must be a square matrix.");
+    }
+    if (mxGetN(IN_G) != m) {
+      mexErrMsgTxt("DtX and G have incompatible sizes.");
+    }
+  }
+  
+  else {
+    mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified.");
+  }
+  
+  
+  /* Do OMP! */
+  
+  GAMMA_OUT = ompcore(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1);
+  
+  return;
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/omp2mex.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/omp2mex.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,23 @@
+%This is the Matlab interface to the OMP2 MEX implementation.
+%The function is not for independent use, only through omp2.m.
+
+
+%OMP2MEX Matlab interface to the OMP2 MEX implementation.
+%  GAMMA = OMP2MEX(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE)
+%  invokes the OMP2 MEX function according to the specified parameters. Not
+%  all the parameters are required. Those among D, X, DtX, XtX and G which
+%  are not specified should be passed as [].
+%
+%  EPSILON - the target error.
+%  SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero.
+%  MSGDELTA - the delay in secs between messages. Zero means no messages.
+%  MAXATOMS - the max number of atoms per signal, negative for no max.
+%  PROFILE - nonzero means that profiling information should be printed.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/omp2mex.mexa64
Binary file DL/RLS-DLA/private/omp2mex.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/ompcore.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/ompcore.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,409 @@
+/**************************************************************************
+ *
+ * File name: ompcore.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 25.8.2009
+ *
+ *************************************************************************/
+
+
+#include "ompcore.h"
+#include "omputils.h"
+#include "ompprof.h"
+#include "myblas.h"
+#include <math.h>
+#include <string.h>
+
+
+
+/******************************************************************************
+ *                                                                            *
+ *                           Batch-OMP Implementation                         *
+ *                                                                            *
+ ******************************************************************************/  
+
+mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L,
+                 int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp)
+{
+  
+  profdata pd;
+  mxArray *Gamma;
+  mwIndex i, j, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count;
+  mwSize allocated_coefs, allocated_cols;
+  int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms;
+  double *alpha, *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; 
+  double eps2, resnorm, delta, deltaprev, secs_remain;
+  int mins_remain, hrs_remain;
+  clock_t lastprint_time, starttime;
+ 
+  
+  
+  /*** status flags ***/
+  
+  DtX_specified = (DtX!=0);   /* indicates whether D'*x was provided */
+  XtX_specified = (XtX!=0);   /* indicates whether sum(x.*x) was provided */
+  
+  standardomp = (G==0);       /* batch-omp or standard omp are selected depending on availability of G */
+  batchomp = !standardomp;
+  
+  
+  
+  /*** allocate output matrix ***/
+  
+  
+  if (gamma_mode == FULL_GAMMA) {
+    
+    /* allocate full matrix of size m X L */
+    
+    Gamma = mxCreateDoubleMatrix(m, L, mxREAL);
+    gammaPr = mxGetPr(Gamma);
+    gammaIr = 0;
+    gammaJc = 0;
+  }
+  else {
+    
+    /* allocate sparse matrix with room for allocated_coefs nonzeros */
+    
+    /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */
+    allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T;
+    Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL);
+    gammaPr = mxGetPr(Gamma);
+    gammaIr = mxGetIr(Gamma);
+    gammaJc = mxGetJc(Gamma);
+    gamma_count = 0;
+    gammaJc[0] = 0;
+  }
+  
+  
+  /*** helper arrays ***/
+  
+  alpha = (double*)mxMalloc(m*sizeof(double));        /* contains D'*residual */
+  ind = (mwIndex*)mxMalloc(n*sizeof(mwIndex));        /* indices of selected atoms */
+  selected_atoms = (int*)mxMalloc(m*sizeof(int));     /* binary array with 1's for selected atoms */
+  c = (double*)mxMalloc(n*sizeof(double));            /* orthogonal projection result */
+  
+  /* current number of columns in Dsub / Gsub / Lchol */
+  allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T;
+  
+  /* Cholesky decomposition of D_I'*D_I */
+  Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double));
+
+  /* temporary vectors for various computations */
+  tempvec1 = (double*)mxMalloc(m*sizeof(double));
+  tempvec2 = (double*)mxMalloc(m*sizeof(double));
+  
+  if (batchomp) {
+    /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */
+    Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double));
+  }
+  else {
+    /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */
+    Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double));
+    
+    /* stores the residual */
+    r = (double*)mxMalloc(n*sizeof(double));        
+  }
+  
+  if (!DtX_specified) {
+    /* contains D'*x for the current signal */
+    DtX = (double*)mxMalloc(m*sizeof(double));  
+  }
+  
+  
+  
+  /*** initializations for error omp ***/
+  
+  if (erroromp) {
+    eps2 = eps*eps;        /* compute eps^2 */
+    if (T<0 || T>n) {      /* unspecified max atom num - set max atoms to n */
+      T = n;
+    }
+  }
+  
+  
+  
+  /*** initialize timers ***/
+  
+  initprofdata(&pd);             /* initialize profiling counters */
+  starttime = clock();           /* record starting time for eta computations */
+  lastprint_time = starttime;    /* time of last status display */
+  
+  
+  
+  /**********************   perform omp for each signal   **********************/
+  
+  
+  
+  for (signum=0; signum<L; ++signum) {
+    
+    
+    /* initialize residual norm and deltaprev for error-omp */
+    
+    if (erroromp) {
+      if (XtX_specified) {
+        resnorm = XtX[signum];
+      }
+      else {
+        resnorm = dotprod(x+n*signum, x+n*signum, n);
+        addproftime(&pd, XtX_TIME);
+      }
+      deltaprev = 0;     /* delta tracks the value of gamma'*G*gamma */
+    }
+    else {
+      /* ignore residual norm stopping criterion */
+      eps2 = 0;
+      resnorm = 1;
+    }
+    
+    
+    if (resnorm>eps2 && T>0) {
+      
+      /* compute DtX */
+      
+      if (!DtX_specified) {
+        matT_vec(1, D, x+n*signum, DtX, n, m);
+        addproftime(&pd, DtX_TIME);
+      }
+      
+      
+      /* initialize alpha := DtX */
+      
+      memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double));
+      
+      
+      /* mark all atoms as unselected */
+      
+      for (i=0; i<m; ++i) {
+        selected_atoms[i] = 0;
+      }
+      
+    }
+    
+
+    /* main loop */
+    
+    i=0;
+    while (resnorm>eps2 && i<T) {
+
+      /* index of next atom */
+      
+      pos = maxabs(alpha, m);
+      addproftime(&pd, MAXABS_TIME);
+      
+      
+      /* stop criterion: selected same atom twice, or inner product too small */
+      
+      if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) {
+        break;
+      }
+      
+      
+      /* mark selected atom */
+      
+      ind[i] = pos;
+      selected_atoms[pos] = 1;
+      
+      
+      /* matrix reallocation */
+      
+      if (erroromp && i>=allocated_cols) {
+        
+        allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01);
+        
+        Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double));
+        
+        batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) :
+                   (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ;
+      }
+      
+      
+      /* append column to Gsub or Dsub */
+      
+      if (batchomp) {
+        memcpy(Gsub+i*m, G+pos*m, m*sizeof(double));
+      }
+      else {
+        memcpy(Dsub+i*n, D+pos*n, n*sizeof(double));
+      }
+      
+      
+      /*** Cholesky update ***/
+      
+      if (i==0) {
+        *Lchol = 1;
+      }
+      else {
+        
+        /* incremental Cholesky decomposition: compute next row of Lchol */
+        
+        if (standardomp) {
+          matT_vec(1, Dsub, D+n*pos, tempvec1, n, i);      /* compute tempvec1 := Dsub'*d where d is new atom */
+          addproftime(&pd, DtD_TIME);
+        }
+        else {
+          vec_assign(tempvec1, Gsub+i*m, ind, i);          /* extract tempvec1 := Gsub(ind,i) */
+        }
+        backsubst('L', Lchol, tempvec1, tempvec2, n, i);   /* compute tempvec2 = Lchol \ tempvec1 */
+        for (j=0; j<i; ++j) {                              /* write tempvec2 to end of Lchol */
+          Lchol[j*n+i] = tempvec2[j];
+        }
+        
+        /* compute Lchol(i,i) */
+        sum = 0;
+        for (j=0; j<i; ++j) {         /* compute sum of squares of last row without Lchol(i,i) */
+          sum += SQR(Lchol[j*n+i]);
+        }
+        if ( (1-sum) <= 1e-14 ) {     /* Lchol(i,i) is zero => selected atoms are dependent */
+          break;
+        }
+        Lchol[i*n+i] = sqrt(1-sum);
+      }
+      
+      addproftime(&pd, LCHOL_TIME);
+
+      i++;
+      
+      
+      /* perform orthogonal projection and compute sparse coefficients */
+      
+      vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i);   /* extract tempvec1 = DtX(ind) */
+      cholsolve('L', Lchol, tempvec1, c, n, i);                     /* solve LL'c = tempvec1 for c */
+      addproftime(&pd, COMPCOEF_TIME);
+      
+
+      /* update alpha = D'*residual */
+      
+      if (standardomp) {
+        mat_vec(-1, Dsub, c, r, n, i);             /* compute r := -Dsub*c */
+        vec_sum(1, x+n*signum, r, n);              /* compute r := x+r */
+        
+        
+        /*memcpy(r, x+n*signum, n*sizeof(double));   /* assign r := x */
+        /*mat_vec1(-1, Dsub, c, 1, r, n, i);         /* compute r := r-Dsub*c */
+        
+        addproftime(&pd, COMPRES_TIME);
+        matT_vec(1, D, r, alpha, n, m);            /* compute alpha := D'*r */
+        addproftime(&pd, DtR_TIME);
+        
+        /* update residual norm */
+        if (erroromp) {
+          resnorm = dotprod(r, r, n);
+          addproftime(&pd, UPDATE_RESNORM_TIME);
+        }
+      }
+      else {
+        mat_vec(1, Gsub, c, tempvec1, m, i);                              /* compute tempvec1 := Gsub*c */
+        memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double));    /* set alpha = D'*x */
+        vec_sum(-1, tempvec1, alpha, m);                                  /* compute alpha := alpha - tempvec1 */
+        addproftime(&pd, UPDATE_DtR_TIME);
+        
+        /* update residual norm */
+        if (erroromp) {
+          vec_assign(tempvec2, tempvec1, ind, i);      /* assign tempvec2 := tempvec1(ind) */
+          delta = dotprod(c,tempvec2,i);               /* compute c'*tempvec2 */
+          resnorm = resnorm - delta + deltaprev;       /* residual norm update */
+          deltaprev = delta;
+          addproftime(&pd, UPDATE_RESNORM_TIME);
+        }
+      }
+    }
+    
+    
+    /*** generate output vector gamma ***/
+
+    if (gamma_mode == FULL_GAMMA) {    /* write the coefs in c to their correct positions in gamma */
+      for (j=0; j<i; ++j) {
+        gammaPr[m*signum + ind[j]] = c[j];
+      }
+    }
+    else {
+      /* sort the coefs by index before writing them to gamma */
+      quicksort(ind,c,i);
+      addproftime(&pd, INDEXSORT_TIME);
+      
+      /* gamma is full - reallocate */
+      if (gamma_count+i >= allocated_coefs) {
+        
+        while(gamma_count+i >= allocated_coefs) {
+          allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01);
+        }
+        
+        mxSetNzmax(Gamma, allocated_coefs);
+        mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double)));
+        mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex)));
+        
+        gammaPr = mxGetPr(Gamma);
+        gammaIr = mxGetIr(Gamma);
+      }
+      
+      /* append coefs to gamma and update the indices */
+      for (j=0; j<i; ++j) {
+        gammaPr[gamma_count] = c[j];
+        gammaIr[gamma_count] = ind[j];
+        gamma_count++;
+      }
+      gammaJc[signum+1] = gammaJc[signum] + i;
+    }
+    
+    
+    
+    /*** display status messages ***/
+    
+    if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta)
+    {
+      lastprint_time = clock();
+      
+      /* estimated remainig time */
+      secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) ,
+        &hrs_remain, &mins_remain, &secs_remain);
+      
+      mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n",        
+        signum+1, L, hrs_remain, mins_remain, secs_remain);
+      mexEvalString("drawnow;");
+    }
+    
+  }
+  
+  /* end omp */
+  
+  
+  
+  /*** print final messages ***/
+  
+  if (msg_delta>0) {
+    mexPrintf("omp: signal %d / %d\n", signum, L);
+  }
+  
+  if (profile) {
+    printprofinfo(&pd, erroromp, batchomp, L);
+  }
+  
+  
+  
+  /* free memory */
+  
+  if (!DtX_specified) {
+    mxFree(DtX);
+  }
+  if (standardomp) {
+    mxFree(r);
+    mxFree(Dsub);
+  }
+  else {
+    mxFree(Gsub);
+  }  
+  mxFree(tempvec2);
+  mxFree(tempvec1);
+  mxFree(Lchol);
+  mxFree(c);
+  mxFree(selected_atoms);
+  mxFree(ind);
+  mxFree(alpha);
+  
+  return Gamma;
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/ompcore.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/ompcore.h	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,80 @@
+/**************************************************************************
+ *
+ * File name: ompcore.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Contains the core implementation of Batch-OMP / OMP-Cholesky.
+ *
+ *************************************************************************/
+
+
+#ifndef __OMP_CORE_H__
+#define __OMP_CORE_H__
+
+
+#include "mex.h"
+
+
+
+/**************************************************************************
+ * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using
+ * either a fixed number of atoms or an error bound.
+ *
+ * Parameters (not all required):
+ *
+ *   D - the dictionary, of size n X m
+ *   x - the signals, of size n X L
+ *   DtX - D'*x, of size m X L
+ *   XtX - squared norms of the signals in x, sum(x.*x), of length L
+ *   G - D'*D, of size m X m
+ *   T - target sparsity, or maximal number of atoms for error-based OMP
+ *   eps - target residual norm for error-based OMP
+ *   gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA
+ *   profile - if non-zero, profiling info is printed
+ *   msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed
+ *   erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP
+ *
+ * Usage:
+ *
+ *   The function can be called using different parameters, and will have
+ *   different complexity depending on the parameters specified. Arrays which
+ *   are not specified should be passed as null (0). When G is specified, 
+ *   Batch-OMP is performed. Otherwise, OMP-Cholesky is performed.
+ *
+ *   Fixed-sparsity usage:
+ *   ---------------------
+ *   Either DtX, or D and x, must be specified. Specifying DtX is more efficient.
+ *   XtX does not need to be specified.
+ *   When D and x are specified, G is not required. However, not providing G
+ *   will significantly degrade efficiency.
+ *   The number of atoms must be specified in T. The value of eps is ignored.
+ *   Finally, set erroromp to 0.
+ *
+ *   Error-OMP usage:
+ *   ----------------
+ *   Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX
+ *   is more efficient.
+ *   When D and x are specified, G is not required. However, not providing G
+ *   will significantly degrade efficiency.
+ *   The target error must be specified in eps. A hard limit on the number
+ *   of atoms can also be specified via the parameter T. Otherwise, T should 
+ *   be negative. Finally, set erroromp to nonzero.
+ *
+ *
+ * Returns: 
+ *   An mxArray containing the sparse representations of the signals in x
+ *   (allocated using the appropriate mxCreateXXX() function).
+ *   The array is either full or sparse, depending on gamma_mode.
+ *
+ **************************************************************************/
+mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L,
+                 int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp);
+
+
+#endif
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/ompmex.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/ompmex.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,133 @@
+/**************************************************************************
+ *
+ * File name: ompmex.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "ompcore.h"
+#include "omputils.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define IN_D          prhs[0]
+#define IN_X          prhs[1]
+#define IN_DtX        prhs[2]
+#define IN_G          prhs[3]
+#define IN_T          prhs[4]
+#define IN_SPARSE_G   prhs[5]
+#define IN_MSGDELTA   prhs[6]
+#define IN_PROFILE    prhs[7]
+
+
+/* Output Arguments */
+
+#define	GAMMA_OUT     plhs[0]
+
+
+/***************************************************************************************/
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
+
+{
+  double *D, *x, *DtX, *G, msgdelta;
+  int gmode, profile, T;
+  mwSize m, n, L;   /* D is n x m , X is n x L, DtX is m x L */
+  
+  
+  /* check parameters */
+  
+  checkmatrix(IN_D, "OMP", "D");
+  checkmatrix(IN_X, "OMP", "X");
+  checkmatrix(IN_DtX, "OMP", "DtX");
+  checkmatrix(IN_G, "OMP", "G");
+  
+  checkscalar(IN_T, "OMP", "T");
+  checkscalar(IN_SPARSE_G, "OMP", "sparse_g");
+  checkscalar(IN_MSGDELTA, "OMP", "msgdelta");
+  checkscalar(IN_PROFILE, "OMP", "profile");
+
+  
+  /* get parameters */
+  
+  x = D = DtX = G = 0;
+  
+  if (!mxIsEmpty(IN_D))
+    D = mxGetPr(IN_D);
+  
+  if (!mxIsEmpty(IN_X))
+    x = mxGetPr(IN_X);
+  
+  if (!mxIsEmpty(IN_DtX))
+    DtX = mxGetPr(IN_DtX);
+  
+  if (!mxIsEmpty(IN_G))
+    G = mxGetPr(IN_G);
+  
+  T = (int)(mxGetScalar(IN_T)+1e-2);
+  if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) {
+    gmode = SPARSE_GAMMA;
+  }
+  else {
+    gmode = FULL_GAMMA;
+  }
+  msgdelta = mxGetScalar(IN_MSGDELTA);
+  profile = (int)(mxGetScalar(IN_PROFILE)+1e-2);
+  
+  
+  /* check sizes */
+  
+  if (D && x) {
+    n = mxGetM(IN_D);
+    m = mxGetN(IN_D);
+    L = mxGetN(IN_X);
+    
+    if (mxGetM(IN_X) != n) {
+      mexErrMsgTxt("D and X have incompatible sizes.");
+    }
+    
+    if (G) {
+      if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+        mexErrMsgTxt("G must be a square matrix.");
+      }
+      if (mxGetN(IN_G) != m) {
+        mexErrMsgTxt("D and G have incompatible sizes.");
+      }
+    }
+  }
+  
+  else if (DtX) {
+    m = mxGetM(IN_DtX);
+    L = mxGetN(IN_DtX);
+    
+    n = T;  /* arbitrary - it is enough to assume signal length is T */
+    
+    if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+      mexErrMsgTxt("G must be a square matrix.");
+    }
+    if (mxGetN(IN_G) != m) {
+      mexErrMsgTxt("DtX and G have incompatible sizes.");
+    }
+  }
+  
+  else {
+    mexErrMsgTxt("Either D and X, or DtX, must be specified.");
+  }
+  
+  
+  /* Do OMP! */
+  
+  GAMMA_OUT = ompcore(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0);
+  
+  return;
+}
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/ompmex.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/ompmex.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,22 @@
+%This is the Matlab interface to the OMP MEX implementation.
+%The function is not for independent use, only through omp.m.
+
+
+%OMPMEX Matlab interface to the OMP MEX implementation.
+%  GAMMA = OMPMEX(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP
+%  MEX function according to the specified parameters. Not all the
+%  parameters are required. Those among D, X, DtX and G which are not
+%  specified should be passed as [].
+%
+%  L - the target sparsity.
+%  SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero.
+%  MSGDELTA - the delay in secs between messages. Zero means no messages.
+%  PROFILE - nonzero means that profiling information should be printed.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/ompmex.mexa64
Binary file DL/RLS-DLA/private/ompmex.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/ompprof.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/ompprof.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,113 @@
+/**************************************************************************
+ *
+ * File name: ompprof.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 11.4.2009
+ *
+ *************************************************************************/
+
+
+#include "ompprof.h"
+
+
+/* initialize profiling information */
+
+void initprofdata(profdata *pd)
+{
+  pd->DtX_time = 0;
+  pd->XtX_time = 0;
+  pd->DtR_time = 0;
+  pd->maxabs_time = 0;
+  pd->DtD_time = 0;
+  pd->Lchol_time = 0;
+  pd->compcoef_time = 0;
+  pd->update_DtR_time = 0;
+  pd->update_resnorm_time = 0;
+  pd->compres_time = 0;
+  pd->indexsort_time = 0;
+  
+  pd->DtX_time_counted = 0;
+  pd->XtX_time_counted = 0;
+  pd->DtR_time_counted = 0;
+  pd->DtD_time_counted = 0;
+  pd->update_DtR_time_counted = 0;
+  pd->resnorm_time_counted = 0;
+  pd->compres_time_counted = 0;
+  pd->indexsort_time_counted = 0;
+  
+  pd->prevtime = clock();
+}
+
+
+/* add elapsed time to profiling data according to specified computation */
+
+void addproftime(profdata *pd, int comptype)
+{
+  switch(comptype) {
+    case DtX_TIME:            pd->DtX_time            += clock()-pd->prevtime; pd->DtX_time_counted = 1; break;
+    case XtX_TIME:            pd->XtX_time            += clock()-pd->prevtime; pd->XtX_time_counted = 1; break;
+    case DtR_TIME:            pd->DtR_time            += clock()-pd->prevtime; pd->DtR_time_counted = 1; break;
+    case DtD_TIME:            pd->DtD_time            += clock()-pd->prevtime; pd->DtD_time_counted = 1; break;
+    case COMPRES_TIME:        pd->compres_time        += clock()-pd->prevtime; pd->compres_time_counted = 1; break;
+    case UPDATE_DtR_TIME:     pd->update_DtR_time     += clock()-pd->prevtime; pd->update_DtR_time_counted = 1; break;
+    case UPDATE_RESNORM_TIME: pd->update_resnorm_time += clock()-pd->prevtime; pd->resnorm_time_counted = 1; break;
+    case INDEXSORT_TIME:      pd->indexsort_time      += clock()-pd->prevtime; pd->indexsort_time_counted = 1; break;
+    case MAXABS_TIME:         pd->maxabs_time         += clock()-pd->prevtime; break;
+    case LCHOL_TIME:          pd->Lchol_time          += clock()-pd->prevtime; break;
+    case COMPCOEF_TIME:       pd->compcoef_time       += clock()-pd->prevtime; break;
+  }
+  pd->prevtime = clock();
+}
+
+
+/* print profiling info */
+
+void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum)
+{
+  clock_t tottime;
+  
+  tottime = pd->DtX_time + pd->XtX_time + pd->DtR_time + pd->DtD_time + pd->compres_time + pd->maxabs_time + 
+            pd->Lchol_time + pd->compcoef_time + pd->update_DtR_time + pd->update_resnorm_time + pd->indexsort_time;
+  
+  mexPrintf("\n\n*****  Profiling information for %s  *****\n\n", erroromp? "OMP2" : "OMP");
+  
+  mexPrintf("OMP mode: %s\n\n", batchomp? "Batch-OMP" : "OMP-Cholesky");
+  
+  mexPrintf("Total signals processed: %d\n\n", signum);
+  
+  if (pd->DtX_time_counted) {
+    mexPrintf("Compute DtX time:      %7.3lf seconds\n", pd->DtX_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->XtX_time_counted) {
+    mexPrintf("Compute XtX time:      %7.3lf seconds\n", pd->XtX_time/(double)CLOCKS_PER_SEC);
+  }
+  mexPrintf("Max abs time:          %7.3lf seconds\n", pd->maxabs_time/(double)CLOCKS_PER_SEC);
+  if (pd->DtD_time_counted) {
+    mexPrintf("Compute DtD time:      %7.3lf seconds\n", pd->DtD_time/(double)CLOCKS_PER_SEC);
+  }
+  mexPrintf("Lchol update time:     %7.3lf seconds\n", pd->Lchol_time/(double)CLOCKS_PER_SEC);
+  mexPrintf("Compute coef time:     %7.3lf seconds\n", pd->compcoef_time/(double)CLOCKS_PER_SEC);
+  if (pd->compres_time_counted) {
+    mexPrintf("Compute R time:        %7.3lf seconds\n", pd->compres_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->DtR_time_counted) {
+    mexPrintf("Compute DtR time:      %7.3lf seconds\n", pd->DtR_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->update_DtR_time_counted) {
+    mexPrintf("Update DtR time:       %7.3lf seconds\n", pd->update_DtR_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->resnorm_time_counted) {
+    mexPrintf("Update resnorm time:   %7.3lf seconds\n", pd->update_resnorm_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->indexsort_time_counted) {
+    mexPrintf("Index sort time:       %7.3lf seconds\n", pd->indexsort_time/(double)CLOCKS_PER_SEC);
+  }
+  mexPrintf("---------------------------------------\n");
+  mexPrintf("Total time:            %7.3lf seconds\n\n", tottime/(double)CLOCKS_PER_SEC);
+}
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/ompprof.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/ompprof.h	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,106 @@
+/**************************************************************************
+ *
+ * File name: ompprof.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Collection of definitions and functions for profiling the OMP method.
+ *
+ *************************************************************************/
+
+
+#ifndef __OMP_PROF_H__
+#define __OMP_PROF_H__
+
+#include "mex.h"
+#include <time.h>
+
+
+
+/**************************************************************************
+ *
+ * Constants and data types.
+ *
+ **************************************************************************/
+
+
+/* constants denoting the various parts of the algorithm */
+
+enum { DtX_TIME, XtX_TIME, DtR_TIME, MAXABS_TIME, DtD_TIME, LCHOL_TIME, COMPCOEF_TIME, 
+       UPDATE_DtR_TIME, UPDATE_RESNORM_TIME, COMPRES_TIME, INDEXSORT_TIME };
+
+       
+       
+/* profiling data container with counters for each part of the algorithm */
+       
+typedef struct profdata 
+{
+  clock_t prevtime;  /* the time when last initialization/call to addproftime() was performed */
+  
+  clock_t DtX_time;
+  clock_t XtX_time;
+  clock_t DtR_time;
+  clock_t maxabs_time;
+  clock_t DtD_time;
+  clock_t Lchol_time;
+  clock_t compcoef_time;
+  clock_t update_DtR_time;
+  clock_t update_resnorm_time;
+  clock_t compres_time;
+  clock_t indexsort_time;
+  
+  /* flags indicating whether profiling data was gathered */
+  int DtX_time_counted;
+  int XtX_time_counted;
+  int DtR_time_counted;
+  int DtD_time_counted;
+  int update_DtR_time_counted;
+  int resnorm_time_counted;
+  int compres_time_counted;
+  int indexsort_time_counted;
+  
+} profdata;
+
+
+
+/**************************************************************************
+ *
+ * Initialize a profdata structure, zero all counters, and start its timer.
+ *
+ **************************************************************************/
+void initprofdata(profdata *pd);
+
+
+/**************************************************************************
+ *
+ * Add elapsed time from last call to addproftime(), or from initialization
+ * of profdata, to the counter specified by comptype. comptype must be one
+ * of the constants in the enumeration above.
+ *
+ **************************************************************************/
+void addproftime(profdata *pd, int comptype);
+
+
+/**************************************************************************
+ *
+ * Print the current contents of the counters in profdata.
+ *
+ * Parameters:
+ *   pd - the profdata to print
+ *   erroromp - indicates whether error-based (nonzero) or sparsity-based (zero)
+ *              omp was performed.
+ *   batchomp - indicates whether batch-omp (nonzero) or omp-cholesky (zero)
+ *              omp was performed.
+ *   signum   - number of signals processed by omp
+ *
+ **************************************************************************/
+void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum);
+
+
+#endif
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/omputils.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/omputils.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,89 @@
+/**************************************************************************
+ *
+ * File name: omputils.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "omputils.h"
+#include <math.h>
+
+
+const char FULL_GAMMA_STR[] = "full";
+const char SPARSE_GAMMA_STR[] = "sparse";
+
+
+/* convert seconds to hours, minutes and seconds */
+
+void secs2hms(double sectot, int *hrs, int *mins, double *secs)
+{
+  *hrs = (int)(floor(sectot/3600)+1e-2);
+  sectot = sectot - 3600*(*hrs);
+  *mins = (int)(floor(sectot/60)+1e-2);
+  *secs = sectot - 60*(*mins);
+}
+
+
+/* quicksort, public-domain C implementation by Darel Rex Finley. */
+/* modification: sorts the array data[] as well, according to the values in the array vals[] */
+
+#define  MAX_LEVELS  300
+
+void quicksort(mwIndex vals[], double data[], mwIndex n) {
+  
+  long piv, beg[MAX_LEVELS], end[MAX_LEVELS], i=0, L, R, swap ;
+  double datapiv;
+  
+  beg[0]=0;
+  end[0]=n;
+  
+  while (i>=0) {
+    
+    L=beg[i]; 
+    R=end[i]-1;
+    
+    if (L<R) {
+      
+      piv=vals[L];
+      datapiv=data[L];
+      
+      while (L<R) 
+      {
+        while (vals[R]>=piv && L<R) 
+          R--;
+        if (L<R) {
+          vals[L]=vals[R];
+          data[L++]=data[R];
+        }
+        
+        while (vals[L]<=piv && L<R) 
+          L++;
+        if (L<R) { 
+          vals[R]=vals[L];
+          data[R--]=data[L];
+        }
+      }
+      
+      vals[L]=piv;
+      data[L]=datapiv;
+      
+      beg[i+1]=L+1;
+      end[i+1]=end[i];
+      end[i++]=L;
+      
+      if (end[i]-beg[i] > end[i-1]-beg[i-1]) {
+        swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap;
+        swap=end[i]; end[i]=end[i-1]; end[i-1]=swap;
+      }
+    }
+    else {
+      i--;
+    }
+  }
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/omputils.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/omputils.h	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,77 @@
+/**************************************************************************
+ *
+ * File name: omputils.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Utility definitions and functions for the OMP library.
+ *
+ *************************************************************************/
+
+
+#ifndef __OMP_UTILS_H__
+#define __OMP_UTILS_H__
+
+#include "mex.h"
+
+
+/* constants for the representation mode of gamma */
+
+extern const char FULL_GAMMA_STR[];      /* "full" */
+extern const char SPARSE_GAMMA_STR[];    /* "sparse" */
+
+
+#define FULL_GAMMA 1
+#define SPARSE_GAMMA 2
+#define INVALID_MODE 3
+
+
+
+/**************************************************************************
+ * Memory management for OMP2.
+ *
+ * GAMMA_INC_FACTOR:
+ * The matrix GAMMA is allocated with sqrt(n)/2 coefficients per signal,
+ * for a total of nzmax = L*sqrt(n)/2 nonzeros. Whenever GAMMA needs to be
+ * increased, it is increased by a factor of GAMMA_INC_FACTOR.
+ *
+ * MAT_INC_FACTOR:
+ * The matrices Lchol, Gsub and Dsub are allocated with sqrt(n)/2
+ * columns each. If additional columns are needed, this number is 
+ * increased by a factor of MAT_INC_FACTOR.
+ **************************************************************************/
+
+#define GAMMA_INC_FACTOR (1.4)
+#define MAT_INC_FACTOR (1.6)
+
+
+
+/**************************************************************************
+ * Convert number of seconds to hour, minute and second representation.
+ *
+ * Parameters:
+ *   sectot - total number of seconds
+ *   hrs, mins, secs - output hours (whole) and minutes (whole) and seconds
+ *
+ **************************************************************************/
+void secs2hms(double sectot, int *hrs, int *mins, double *secs);
+
+
+
+/**************************************************************************
+ * QuickSort - public-domain C implementation by Darel Rex Finley.
+ *
+ * Modified to sort both the array vals[] and the array data[] according 
+ * to the values in the array vals[].
+ *
+ **************************************************************************/
+void quicksort(mwIndex vals[], double data[], mwIndex n);
+
+
+#endif
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/printf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/printf.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,26 @@
+function str = printf(varargin)
+%PRINTF Print formatted text to screen.
+%  PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to
+%  the format string FMT, and prints the result to the screen.
+%
+%  The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call
+%  DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the
+%  format string options see function SPRINTF.
+%
+%  STR = PRINTF(...) also returns the formatted string.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2008
+
+
+if (nargout>0)
+  str = sprintf(varargin{:});
+  disp(str);
+else
+  disp(sprintf(varargin{:}));
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/private/make.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/private/make.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,41 @@
+function make
+%MAKE Build the OMPBox package.
+%  MAKE compiles all OMPBox MEX functions, using Matlab's default MEX
+%  compiler. If the MEX compiler has not been set-up before, please run
+%
+%    mex -setup
+%
+%  before using this MAKE file.
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
+
+
+% detect platform 
+
+compstr = computer;
+is64bit = strcmp(compstr(end-1:end),'64');
+
+
+% compilation parameters
+
+compile_params = cell(0);
+if (is64bit)
+  compile_params{1} = '-largeArrayDims';
+end
+
+
+% Compile files %
+
+ompsources = {'mexutils.c','ompcore.c','omputils.c','myblas.c','ompprof.c'};
+
+disp('Compiling ompmex...');
+mex('ompmex.c', ompsources{:},compile_params{:});
+
+disp('Compiling omp2mex...');
+mex('omp2mex.c',ompsources{:},compile_params{:});
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/private/mexutils.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/private/mexutils.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,79 @@
+/**************************************************************************
+ *
+ * File name: mexutils.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 15.8.2009
+ *
+ *************************************************************************/
+
+#include "mexutils.h"
+#include <math.h>
+
+
+
+/* verify that the mxArray contains a double matrix */
+
+void checkmatrix(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-D double vector */
+
+void checkvector(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a double scalar */
+
+void checkscalar(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || 
+      mxGetM(param)!=1 || mxGetN(param)!=1) 
+  {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a sparse matrix */
+
+void checksparse(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname);
+  if (!mxIsSparse(param)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-dimensional cell array */
+
+void checkcell_1d(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname);
+  if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/private/mexutils.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/private/mexutils.h	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,103 @@
+/**************************************************************************
+ *
+ * File name: mexutils.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Utility functions for MEX files.
+ *
+ *************************************************************************/
+
+
+#ifndef __MEX_UTILS_H__
+#define __MEX_UTILS_H__
+
+#include "mex.h"
+
+
+
+/**************************************************************************
+ * Function checkmatrix:
+ *
+ * Verify that the specified mxArray is real, of type double, and has 
+ * no more than two dimensions. If not, an error message is printed
+ * and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkmatrix(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkvector:
+ *
+ * Verify that the specified mxArray is 1-D, real, and of type double. The
+ * vector may be a column or row vector. Otherwise, an error message is
+ * printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkvector(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkscalar:
+ *
+ * Verify that the specified mxArray represents a real double scalar value. 
+ * If not, an error message is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkscalar(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checksparse:
+ *
+ * Verify that the specified mxArray contains a sparse matrix. If not,
+ * an error message is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checksparse(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkcell_1d:
+ *
+ * Verify that the specified mxArray is a 1-D cell array. The cell array 
+ * may be arranged as either a column or a row. If not, an error message 
+ * is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkcell_1d(const mxArray *param, char *fname, char *pname);
+
+
+#endif
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/private/printf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/private/printf.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,26 @@
+function str = printf(varargin)
+%PRINTF Print formatted text to screen.
+%  PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to
+%  the format string FMT, and prints the result to the screen.
+%
+%  The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call
+%  DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the
+%  format string options see function SPRINTF.
+%
+%  STR = PRINTF(...) also returns the formatted string.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2008
+
+
+if (nargout>0)
+  str = sprintf(varargin{:});
+  disp(str);
+else
+  disp(sprintf(varargin{:}));
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/pwsmoothfield.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/pwsmoothfield.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,32 @@
+function f = pwsmoothfield(n,var,alpha)
+% PWSMOOTHFIELD(N,VAR,ALPHA)
+%   Generate an image of piecewise smooth filtered white noise
+%
+%   N = sidelength of the field in pixels
+%   VAR = variance of original white nose
+%   ALPHA = fraction of FFT coefficents to keep
+%
+%   Returns an N-by-N array
+%
+%   This file is used with the kind permission of Stephen J. Wright
+%   (swright@cs.wisc.edu), and was originally included in the GPSR
+%   v1.0 distribution: http://www.lx.it.pt/~mtf/GPSR .
+
+% $Id: pwsmoothfield.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+f = sqrt(var)*randn(n);
+F = fft2(f);
+a = ceil(n*alpha/2);
+b = fix(n*(1-alpha));
+F(a+1:a+b,:) = 0;
+F(:,a+1:a+b) = 0;
+f = real(ifft2(F));
+
+for i = 1:n
+    for j = 1:n
+        if (j/n >= 15*(i/n - 0.5)^3 + 0.4)
+            f(i,j) = f(i,j) + sqrt(var);
+        end
+    end
+end
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/reggrid.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/reggrid.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,136 @@
+function [varargout] = reggrid(sz,num,mode)
+%REGGRID Regular sampling grid.
+%  [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM) returns the indices
+%  of a regular uniform sampling grid over a p-dimensional matrix with
+%  dimensions N1xN2x...xNp. NUM is the minimal number of required samples,
+%  and it is ensured that the actual number of samples, given by
+%  length(I1)xlength(I2)x...xlength(Ip), is at least as large as NUM.
+%
+%  [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM,'MODE') specifies the
+%  method for distributing the samples along each dimension. Valid modes
+%  include 'eqdist' (the default mode) and 'eqnum'. 'eqdist' indicates an
+%  equal distance between the samples in each dimension, while 'eqnum'
+%  indicates an equal number of samples in each dimension.
+%
+%  Notes about MODE:
+%
+%    1. The 'eqnum' mode will generally fail when the p-th root of NUM
+%    (i.e. NUM^(1/p)) is larger than min([N1 N2 ... Np]). Thus 'eqdist' is
+%    the more useful choice for sampling an arbitrary number of samples
+%    from the matrix (up to the total number of matrix entries).
+%  
+%    2. In both modes, the equality (of the distance between samples, or
+%    the number of samples in each dimension) is only approximate. This is
+%    because REGGRID attempts to maintain the appropriate equality while at
+%    the same time find a sampling pattern where the total number of
+%    samples is as close as possible to NUM. In general, the larger {Ni}
+%    and NUM are, the tighter the equality.
+%
+%  Example: Sample a set of blocks uniformly from a 2D image.
+%
+%    n = 512; blocknum = 20000; blocksize = [8 8];
+%    im = rand(n,n);
+%    [i1,i2] = reggrid(size(im)-blocksize+1, blocknum);
+%    blocks = sampgrid(im, blocksize, i1, i2);
+%
+%  See also SAMPGRID.
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  November 2007
+
+dim = length(sz);
+
+if (nargin<3)
+  mode = 'eqdist';
+end
+
+if (any(sz<1))
+  error(['Invalid matrix size : [' num2str(sz) ']']);
+end
+
+if (num > prod(sz))
+  warning(['Invalid number of samples, returning maximum number of samples.']);
+elseif (num <= 0)
+  if (num < 0)
+    warning('Invalid number of samples, assuming 0 samples.');
+  end
+  for i = 1:length(sz)
+    varargout{i} = [];
+  end
+  return;
+end
+
+
+if (strcmp(mode,'eqdist'))
+  
+  % approximate distance between samples: total volume divided by number of
+  % samples gives the average volume per sample. then, taking the p-th root
+  % gives the average distance between samples
+  d = (prod(sz)/num)^(1/dim);
+  
+  % compute the initial guess for number of samples in each dimension.
+  % then, while total number of samples is too large, decrese the number of
+  % samples by one in the dimension where the samples are the most crowded.
+  % finally, do the opposite process until just passing num, so the final
+  % number of samples is the closest to num from above.
+  
+  n = min(max(round(sz/d),1),sz);   % set n so that it saturates at 1 and sz
+  
+  active_dims = find(n>1);    % dimensions where the sample num can be reduced
+  while(prod(n)>num && ~isempty(active_dims))
+    [y,id] = min((sz(active_dims)-1)./n(active_dims));
+    n(active_dims(id)) = n(active_dims(id))-1;
+    if (n(active_dims(id)) < 2)
+      active_dims = find(n>1);
+    end
+  end
+
+  active_dims = find(n<sz);    % dimensions where the sample num can be increased
+  while(prod(n)<num && ~isempty(active_dims))
+    [y,id] = max((sz(active_dims)-1)./n(active_dims));
+    n(active_dims(id)) = n(active_dims(id))+1;
+    if (n(active_dims(id)) >= sz(active_dims(id)))
+      active_dims = find(n<sz);
+    end
+  end
+
+  for i = 1:dim
+    varargout{i} = round((1:n(i))/n(i)*sz(i));
+    varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2);
+  end
+  
+elseif (strcmp(mode,'eqnum'))
+  
+  % same idea as above
+  n = min(max( ones(size(sz)) * round(num^(1/dim)) ,1),sz);
+
+  active_dims = find(n>1);
+  while(prod(n)>num && ~isempty(active_dims))
+    [y,id] = min((sz(active_dims)-1)./n(active_dims));
+    n(active_dims(id)) = n(active_dims(id))-1;
+    if (n(active_dims(id)) < 2)
+      active_dims = find(n>1);
+    end
+  end
+  
+  active_dims = find(n<sz);
+  while(prod(n)<num && ~isempty(active_dims))
+    [y,id] = max((sz(active_dims)-1)./n(active_dims));
+    n(active_dims(id)) = n(active_dims(id))+1;
+    if (n(active_dims(id)) >= sz(active_dims(id)))
+      active_dims = find(n<sz);
+    end
+  end
+  
+  for i = 1:dim
+    varargout{i} = round((1:n(i))/n(i)*sz(i));
+    varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2);
+  end
+else
+  error('Invalid sampling mode');
+end
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/remove_dc.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/remove_dc.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,31 @@
+function [y,dc] = remove_dc(x,columns)
+%REMOVE_DC Remove DC channel from signals.
+%   [Y,DC] = REMOVE_DC(X) removes the DC channel (i.e. the mean) from the
+%   specified (possibly multi-dimensional) signal X. Y is the DC-free
+%   signal and is the same size as X. DC is a scalar containing the mean of
+%   the signal X.
+%
+%   [Y,DC] = REMOVE_DC(X,'columns') where X is a 2D matrix, treats the
+%   columns of X as a set of 1D signals, removing the DC channel from each
+%   one individually. Y is the same size as X and contains the DC-free
+%   signals. DC is a row vector of length size(X,2) containing the means of
+%   the signals in X.
+%
+%   See also ADD_DC.
+
+
+if (nargin==2 && strcmpi(columns,'columns')), columns = 1;
+else columns = 0;
+end
+
+if (columns)
+  dc = mean(x);
+  y = addtocols(x,-dc);
+else
+  if (ndims(x)==2)  % temporary, will remove in future
+    warning('Treating 2D matrix X as a single signal and not each column individually');
+  end
+  dc = mean(x(:));
+  y = x-dc;
+end
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/rowlincomb.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/rowlincomb.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,148 @@
+/**************************************************************************
+ *
+ * File name: rowlincomb.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 21.5.2009
+ *
+ *************************************************************************/
+
+#include "mex.h"
+
+
+/* Input Arguments */
+
+#define	X_IN	   prhs[0]
+#define A_IN     prhs[1]
+#define ROWS_IN  prhs[2]
+#define COLS_IN  prhs[3]
+
+
+/* Output Arguments */
+
+#define	Y_OUT	plhs[0]
+
+
+void mexFunction(int nlhs, mxArray *plhs[],
+int nrhs, const mxArray*prhs[])
+
+{
+  double *A, *x, *y, *rows, *cols;
+  mwSize m,n,m1,n1,m2,n2,rownum,colnum;
+  mwIndex *row_ids,*col_ids,i,j;
+  int colnumspecified=0;
+  
+  
+  /* Check for proper number of arguments */
+  
+  if (nrhs!=3 && nrhs!=4) {
+    mexErrMsgTxt("Invalid number of input arguments.");
+  } else if (nlhs > 1) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+  
+  
+  /* Check the input dimensions */
+  
+  m = mxGetM(A_IN);
+  n = mxGetN(A_IN);
+  if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) {
+    mexErrMsgTxt("ROWLINCOMB requires that A be a double matrix.");
+  }
+  
+  m1 = mxGetM(ROWS_IN);
+  n1 = mxGetN(ROWS_IN);
+  if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) {
+    mexErrMsgTxt("ROWLINCOMB requires that ROWS be an index vector of type double.");
+  }
+  rownum = (m1 > n1) ? m1 : n1;   /* the number of rows in the linear combination */
+  
+  m2 = mxGetM(X_IN);
+  n2 = mxGetN(X_IN);
+  if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ((m2!=1) && (n2!=1))) {
+    mexErrMsgTxt("ROWLINCOMB requires that X be a double vector.");
+  }
+  
+  if (m2 != rownum && n2 != rownum) {
+    mexErrMsgTxt("The length of X does not match the number of rows in ROWS.");
+  }
+  
+  if (nrhs==4) {
+    m1 = mxGetM(COLS_IN);
+    n1 = mxGetN(COLS_IN);
+    if (!mxIsDouble(COLS_IN) || mxIsComplex(COLS_IN) || (m1!=1 && n1!=1)) {
+      mexErrMsgTxt("ROWLINCOMB requires that COLS be an index vector of type double.");
+    }
+    colnum = (m1 > n1) ? m1 : n1;   /* the number of columns */
+    colnumspecified = 1;
+    cols = mxGetPr(COLS_IN);
+    
+    Y_OUT = mxCreateDoubleMatrix(1, colnum, mxREAL);
+  }
+  else {
+    cols = 0;
+    Y_OUT = mxCreateDoubleMatrix(1, n, mxREAL);
+  }
+  
+  
+  /* Assign pointers to the various parameters */
+  A = mxGetPr(A_IN);
+  rows = mxGetPr(ROWS_IN);
+  x = mxGetPr(X_IN);
+  y = mxGetPr(Y_OUT);
+  
+  
+  /* check row indices */
+  
+  row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex));
+  
+  for (i=0; i<rownum; ++i) {
+    row_ids[i] = (mwIndex)(rows[i]+0.1)-1;
+    if (row_ids[i]<0 || row_ids[i]>=m) {
+      mexErrMsgTxt("Row index in ROWS is out of range.");
+    }
+  }
+  
+  
+  
+  if (colnumspecified) {
+    
+    /* check column indices */
+    col_ids = (mwIndex*)mxMalloc(colnum*sizeof(mwIndex));
+    
+    for (i=0; i<colnum; ++i) {
+      col_ids[i] = (mwIndex)(cols[i]+0.1)-1;
+      if (col_ids[i]<0 || col_ids[i]>=n) {
+        mexErrMsgTxt("Column index in COLS is out of range.");
+      }
+    }
+    
+    /* Do the actual computation */
+    for (j=0; j<colnum; ++j) {
+      for (i=0; i<rownum; ++i) {
+        y[j] += A[m*col_ids[j]+row_ids[i]]*x[i];
+      }
+    }
+    
+    mxFree(col_ids);
+  }
+  
+  else {
+    
+    /* Do the actual computation */
+    for (j=0; j<n; ++j) {
+      for (i=0; i<rownum; ++i) {
+        y[j] += A[m*j+row_ids[i]]*x[i];
+      }
+    }
+  }
+  
+  
+  mxFree(row_ids);
+  
+  return;
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/rowlincomb.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/rowlincomb.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,26 @@
+%ROWLINCOMB Linear combination of matrix rows.
+%  Y = ROWLINCOMB(X,A,ROWS) computes a linear combination of the rows of
+%  the matrix A. The row indices are specified in the vector ROWS, and the
+%  correspoinding coefficients are specified in the vector X. The vectors
+%  ROWS and X must be of the same length. The call Y = ROWLINCOMB(X,A,ROWS)
+%  is essentially equivalent to the command
+%
+%         Y = X'*A(ROWS,:) .
+%
+%  However, it is implemented much more efficiently.
+%
+%  Y = ROWLINCOMB(X,A,ROWS,COLS) only works on the columns of A specified
+%  in COLS, returning a vector of length equal to COLS. This call is
+%  essentially equivalent to the command
+%
+%         Y = X'*A(ROWS,COLS) .
+%
+%  See also COLLINCOMB.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/rowlincomb.mexa64
Binary file DL/RLS-DLA/private/rowlincomb.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/sampgrid.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/sampgrid.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,72 @@
+function y = sampgrid(x,blocksize,varargin)
+%SAMPGRID Sample a multi-dimensional matrix on a regular grid.
+%  Y = SAMPGRID(X,BLOCKSIZE,I1,I2,...,Ip) extracts block samples of size
+%  BLOCKSIZE from the p-dimensional matrix X, arranging the samples as the
+%  column vectors of the matrix Y. The locations of the (1,1,..,1)-th
+%  elements of each block are given in the index vectors I1,I2,..Ip. The
+%  total number of samples taken is length(I1)xlength(I2)x...xlength(Ip).
+%  BLOCKSIZE should either be a p-element vector of the form [N1,N2,...Np],
+%  or a scalar N which is shorthand for the square block size [N N ... N].
+%
+%  Example: Sample a set of blocks uniformly from a 2D image.
+%
+%    n = 512; blocknum = 20000; blocksize = [8 8];
+%    im = rand(n,n);
+%    [i1,i2] = reggrid(size(im)-blocksize+1, blocknum);
+%    blocks = sampgrid(im, blocksize, i1, i2);
+%
+%  See also REGGRID.
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  November 2007
+
+
+p = ndims(x);
+if (p==2 && any(size(x)==1) && length(blocksize)==1)
+  p = 1;
+end
+
+if (numel(blocksize)==1)
+  blocksize = ones(1,p)*blocksize;
+end
+
+n = zeros(1,p);
+for i = 1:p
+  n(i) = length(varargin{i});
+end
+
+nsamps = prod(n);
+
+% create y of the same class as x
+y = zeros(prod(blocksize),nsamps,class(x));
+
+% ids() contains the index of the current block in I1..Ip
+ids = ones(p,1);
+
+% block_ids contains the indices of the current block in X
+block_ids = cell(p,1);
+for j = 1:p
+  block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1;
+end
+
+for k = 1:nsamps
+  block = x(block_ids{:});
+  y(:,k) = block(:);
+  
+  % increment ids() and block_ids{}
+  if (k<nsamps)
+    j = 1;
+    while (ids(j) == n(j))
+      ids(j) = 1;
+      block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1;
+      j = j+1;
+    end
+    ids(j) = ids(j)+1;
+    block_ids{j} = varargin{j}(ids(j)) : varargin{j}(ids(j))+blocksize(j)-1;
+  end
+end
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/scalarToRGB.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/scalarToRGB.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,18 @@
+function s = scalarToRGB(x,colors)
+% input values are assumed to lie between 0 and 1
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: scalarToRGB.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+l  = size(colors,1);
+m  = size(x,1);
+n  = size(x,2);
+s  = zeros(m,n,3);
+
+for i=1:m
+  for j=1:n
+     idx = max(1,min(l,1+floor((l-1) * x(i,j))));
+     s(i,j,:) = colors(idx,:);
+  end
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/secs2hms.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/secs2hms.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,28 @@
+function [h,m,s] = secs2hms(t)
+%SECS2HMS Convert seconds to hours, minutes and seconds.
+%  [H,M,S] = SECS2HMS(T) converts the specified number of seconds T to
+%  hours, minutes and seconds. H and M are whole numbers, and S is real.
+%
+%  Example: Estimate the remaining time of a loop
+%
+%    n = 10; tic;
+%    for i = 1:n
+%      pause(1);
+%      [h,m,s] = secs2hms( (n-i)*toc/i );
+%      printf('estimated remaining time: %02d:%02d:%05.2f',h,m,s);
+%    end
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2008
+
+
+s = t;
+h = fix(s/3600);
+s = rem(s,3600);
+m = fix(s/60);
+s = rem(s,60);
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/spdiag.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/spdiag.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,38 @@
+function Y = spdiag(V,K)
+%SPDIAG Sparse diagonal matrices.
+%   SPDIAG(V,K) when V is a vector with N components is a sparse square
+%   matrix of order N+ABS(K) with the elements of V on the K-th diagonal. 
+%   K = 0 is the main diagonal, K > 0 is above the main diagonal and K < 0
+%   is below the main diagonal. 
+%
+%   SPDIAG(V) is the same as SPDIAG(V,0) and puts V on the main diagonal.
+%
+%   See also DIAG, SPDIAGS.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  June 2008
+
+
+if (nargin<2)
+  K = 0;
+end
+
+n = length(V) + abs(K);
+
+if (K>0)
+  i = 1:length(V);
+  j = K+1:n;
+elseif (K<0)
+  i = -K+1:n;
+  j = 1:length(V);
+else
+  i = 1:n;
+  j = 1:n;
+end
+
+Y = sparse(i,j,V(:),n,n);
\ No newline at end of file
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/sprow.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/sprow.c	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,94 @@
+/**************************************************************************
+ *
+ * File name: sprow.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 24.8.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define	A_IN	prhs[0]
+#define J_IN  prhs[1]
+
+
+/* Output Arguments */
+
+#define	X_OUT	  plhs[0]
+#define	ID_OUT	plhs[1]
+
+
+void mexFunction(int nlhs, mxArray *plhs[], 
+		             int nrhs, const mxArray*prhs[])
+     
+{ 
+  double *pr, *x, *id, rowid;
+  mwIndex *ir, *jc;
+  mwSize m, n;
+  mwIndex i, j, k, l, rowlen;
+  
+  if (nrhs != 2) {
+    mexErrMsgTxt("GETSPROW requires two input arguments.");
+  } else if (nlhs > 2) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+  
+  checkmatrix(A_IN, "GETSPROW", "A");
+  checksparse(A_IN, "GETSPROW", "A");
+  checkscalar(J_IN, "GETSPROW", "J");
+  
+  m = mxGetM(A_IN);
+  n = mxGetN(A_IN);
+  
+  rowid = mxGetScalar(J_IN);
+  if (rowid < 0) {
+    mexErrMsgTxt("Invalid row index.");
+  }
+  j = (mwIndex)(rowid + 1e-2);
+  if (j<1 || j>m) {
+    mexErrMsgTxt("Row index is out of range.");
+  }
+  j--;
+  
+  pr = mxGetPr(A_IN);
+  ir = mxGetIr(A_IN);
+  jc = mxGetJc(A_IN);
+  
+  /* Determine length of row */
+  rowlen = 0;
+  for (i=0; i<jc[n]; ++i) {
+    rowlen += (ir[i]==j);
+  }
+  
+  /* Allocate output parameters */
+  X_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL);
+  ID_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL);
+  
+  x = mxGetPr(X_OUT);
+  id = mxGetPr(ID_OUT);
+  
+  /* Compute j-th row */
+  k=0;
+  for (l=1; l<=n; ++l) {
+    i = jc[l-1];
+    while (i<jc[l] && ir[i]<j) {
+      i++;
+    }
+    if (i<jc[l] && ir[i]==j) {
+      x[k] = pr[i];
+      id[k] = l;
+      k++;
+    }
+  }
+  
+}
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/sprow.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/sprow.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,19 @@
+%SPROW Extract row of sparse matrix.
+%  X = SPROW(A,J) where A is a sparse matrix, returns the nonzero values in
+%  the row A(J,:).
+%
+%  [X,ID] = SPROW(A,J) also returns the column indices of the nonzeros.
+%
+%  Note that the call [X,ID] = SPROW(A,J) is equivalent (but more efficient
+%  than) the Matlab code
+%
+%    IDS = find(A(J,:)); 
+%    X = A(J,IDS);
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/sprow.mexa64
Binary file DL/RLS-DLA/private/sprow.mexa64 has changed
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/thumbFromOp.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/thumbFromOp.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,42 @@
+function P = thumbFromOp(op,m,n,sm,sn,grayscale)
+% Output matrix P, of size m x n, sampled at
+% sm x sn upper top
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbFromOp.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+if nargin < 6, grayscale = 0; end
+    
+info = op([],0);
+
+sm = min(info{1},sm);
+sn = min(info{2},sn);
+M  = zeros(sm,sn);
+
+for i=1:sn
+   v = zeros(info{2},1); v(i) = 1;
+   w = real(op(v,1));
+
+   M(:,i) = w(1:sm);
+end
+
+mn = min(min(M));
+mx = max(max(M));
+M  = (M - mn) / (mx-mn);
+
+idxm = floor(linspace(1,sm+1,m+1)); idxm = idxm(1:end-1);
+idxn = floor(linspace(1,sn+1,n+1)); idxn = idxn(1:end-1);
+
+if grayscale
+   P = 1-M(idxm,idxn);
+else
+   clrmap = hsv;
+   M = 1 + round(M * (length(clrmap)-1));
+   P = zeros(m,n,3);
+   for j1=1:m
+       for j2=1:n
+           P(j1,j2,:) = clrmap(M(idxm(j1),idxn(j2)),:);
+       end
+   end
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/thumbPlot.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/thumbPlot.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,38 @@
+function P = thumbPlot(P,x,y,color)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbPlot.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+m = size(P,1);
+n = size(P,2);
+if (size(P,3) == 0) & (length(color) == 3)
+  % Convert to gray-scale
+  color = 0.30*color(1) + 0.59*color(2) + 0.11*color(3);
+end
+
+mnx = min(x);  % Minimum x
+mxx = max(x);  % Maximum x
+mny = min(y);  % Minimum y
+mxy = max(y);  % Maximum y
+dy  = (mxy - mny) * 0.1;   % Offset on vertical axis
+sx  = (mxx - mnx) * 1.0;   % Scale of horizontal axis
+sy  = (mxy - mny) * 1.2;   % Scale of vertical axis
+
+if (sx < 1e-6), sx = 1; end
+if (sy < 1e-6), sy = 1; end
+
+for i=1:length(x)-1
+   x0 = floor(1 + (n-1) * (x(i  ) - mnx) / sx);
+   x1 = floor(1 + (n-1) * (x(i+1) - mnx) / sx);
+   y0 = floor(    (n-1) * (y(i  ) - mny + dy) / sy);
+   y1 = floor(    (n-1) * (y(i+1) - mny + dy) / sy);
+   
+   samples = 1+2*max(abs(x1-x0)+1,abs(y1-y0)+1);
+   c       = linspace(0,1,samples);
+   idx     = round((1-c)*x0 + c*x1);
+   idy     = n - round((1-c)*y0 + c*y1);
+   for j=1:samples
+      P(idy(j),idx(j),:) = color;
+   end
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/thumbwrite.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/thumbwrite.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,18 @@
+function thumbwrite(data,name,opts)
+
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbwrite.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+%
+% data       Thumbnail data in range 0-1
+% name       Name of file (no extension)
+% opts
+%   .thumbtype  Type of image (png, eps, ps, ...)
+%   .thumbdir   Output directory
+%
+
+[type,ext] = getFigureExt(opts.thumbtype);
+data = round(data * 255) / 255;
+imwrite(data,[opts.thumbpath,name,'.',ext],type);
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/timerclear.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/timerclear.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,37 @@
+function timerclear()
+%TIMERCLEAR Clear all timers.
+%   TIMERCLEAR clears all currenly registered timers, invalidating all
+%   timer ids.
+%
+%   Note: since registered timers do not consume CPU power except for when
+%   the TIMER<*> functions are called, this function is only useful in
+%   situations where a large number of timers have been initialized, and
+%   there is a need to reclaim memory.
+%
+%   See also TIMERINIT, TIMERETA.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  June 2008
+
+
+global utiltbx_timer_start_times    % start times
+global utiltbx_time_lastdisp        % last display times
+global utiltbx_timer_iternums       % iteration numbers
+global utiltbx_timer_lastiter       % last queried iteration numbers
+global utiltbx_timer_name           % timer names
+global utiltbx_timer_callfun        % timer calling functions
+
+
+% clear all timers %
+
+utiltbx_timer_start_times = [];
+utiltbx_time_lastdisp = [];
+utiltbx_timer_iternums = [];
+utiltbx_timer_lastiter = [];
+utiltbx_timer_name = [];
+utiltbx_timer_callfun = [];
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/timereta.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/timereta.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,98 @@
+function varargout = timereta(tid,iter,delay)
+%TIMERETA Estimated remaining time.
+%   S = TIMERETA(TID,ITER) returns the estimated remaining time (in
+%   seconds) for the process associated with timer TID, assuming the
+%   process has completed ITER iterations. Note: the function will exit
+%   with an error if the timer TID does not have an associated number of
+%   iterations (see function TIMERINIT).
+%
+%   [H,M,S] = TIMERETA(TID,ITER) returns the estimated remaining time in
+%   hours, minutes and seconds.
+%
+%   TIMERETA(TID,ITER), with no output arguments, prints the estimated
+%   remaining time to the screen. The time is displayed in the format
+%
+%     TIMERNAME: iteration ITER / ITERNUM, estimated remaining time: HH:MM:SS.SS
+%
+%   If the timer has no assigned name, the display format changes to
+%
+%     Iteration ITER / ITERNUM, estimated remaining time: HH:MM:SS.SS
+%
+%   TIMERETA(TID,ITER,DELAY) only displays the remaining time if the
+%   time elapsed since the previous printout is at least DELAY seconds.
+%
+%   See also TIMERINIT, TIMERCLEAR.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  June 2008
+
+
+global utiltbx_timer_start_times 
+global utiltbx_timer_iternums
+global utiltbx_timer_lastiter
+global utiltbx_time_lastdisp
+global utiltbx_timer_name
+
+
+if (tid<1 || tid>length(utiltbx_timer_iternums))
+  error('Unknown timer id');
+end
+
+if (utiltbx_timer_iternums(tid) < 0)
+  error('Specified timer does not have an associated number of iterations');
+end
+
+% update last reported iteration number
+utiltbx_timer_lastiter(tid) = iter;
+
+% compute elapsed time
+starttime = utiltbx_timer_start_times(tid,:);
+currtime = clock;
+timediff = etime(currtime, starttime);
+
+% total iteration number
+iternum = utiltbx_timer_iternums(tid);
+
+% compute eta
+timeremain = (iternum-iter)*timediff/iter;
+
+% return eta in seconds
+if (nargout==1)
+  varargout{1} = timeremain;
+  
+% return eta in hms
+elseif (nargout==3)
+  [varargout{1}, varargout{2}, varargout{3}] = secs2hms(timeremain);
+  
+  
+% print eta
+elseif (nargout==0)
+  
+  % check last display time
+  lastdisptime = utiltbx_time_lastdisp(tid,:);
+  if (nargin>2 && etime(currtime,lastdisptime) < delay)
+    return;
+  end
+  
+  % update last display time
+  utiltbx_time_lastdisp(tid,:) = currtime;
+
+  % display timer
+  [hrs,mins,secs] = secs2hms(timeremain);
+  if (isempty(utiltbx_timer_name{tid}))
+    printf('Iteration %d / %d, estimated remaining time: %02d:%02d:%05.2f', iter, iternum, hrs, mins, secs);
+  else
+    timername = utiltbx_timer_name{tid};
+    printf('%s: iteration %d / %d, estimated remaining time: %02d:%02d:%05.2f', timername, iter, iternum, hrs, mins, secs);
+  end
+   
+% invalid number of outputs
+else
+  error('Invalid number of output arguments');
+end
+
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/timerinit.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/timerinit.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,110 @@
+function tid = timerinit(par1,par2)
+%TIMERINIT Initialize a new timer.
+%   TID = TIMERINIT() initializes a new timer for counting elapsed time,
+%   and returns its id.
+%
+%   TID = TIMERINIT('TIMERNAME') sets the timer name to the specified
+%   string for display purposes.
+%
+%   TID = TIMERINIT(ITERNUM) initializes a new ETA timer for a process with
+%   ITERNUM iterations. An ETA timer can be used for both counting elapsed
+%   time and estimating remaining time.
+%
+%   TID = TIMERINIT('TIMERNAME',ITERNUM) sets the ETA timer name to the
+%   specified string for display purposes.
+%
+%   Example:
+%
+%     tid = timerinit(100); 
+%     for i = 1:100
+%       pause(0.07);
+%       timereta(tid,i,1);
+%     end
+%     timereta(tid,i);
+%
+%   See also TIMERETA, TIMERCLEAR.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  June 2008
+
+
+global utiltbx_timer_start_times    % start times
+global utiltbx_time_lastdisp        % last display times
+global utiltbx_timer_iternums       % iteration numbers
+global utiltbx_timer_lastiter       % last queried iteration numbers
+global utiltbx_timer_name           % timer names
+global utiltbx_timer_callfun        % timer calling functions
+
+
+% parse function arguments %
+
+if (nargin==0)
+
+  iternum = -1;
+  timername = '';
+
+elseif (nargin==1)
+
+  if (ischar(par1))
+    iternum = -1;
+    timername = par1;
+
+  elseif (isnumeric(par1) && numel(par1)==1 && par1>0)
+    iternum = par1;
+    timername = '';
+
+  else
+    error('Invalid number of iterations');
+  end
+
+elseif (nargin==2)
+
+  if (ischar(par1) && isnumeric(par2))
+    if (numel(par2)==1 && par2>0)
+      timername = par1;
+      iternum = par2;
+    else
+      error('Invalid number of iterations');
+    end
+  else
+    error('Invalid function syntax');
+  end
+
+else
+  error('Too many function parameters');
+end
+
+
+% register the timer %
+
+if (isempty(utiltbx_timer_start_times))
+  utiltbx_timer_start_times = clock;
+  utiltbx_time_lastdisp = utiltbx_timer_start_times;
+  utiltbx_timer_iternums = double(iternum);
+  utiltbx_timer_lastiter = 0;
+  utiltbx_timer_name = { timername };
+  utiltbx_timer_callfun = {};
+  tid = 1;
+else
+  utiltbx_timer_start_times(end+1,:) = clock;
+  utiltbx_time_lastdisp(end+1,:) = utiltbx_timer_start_times(end,:);
+  utiltbx_timer_iternums(end+1) = double(iternum);
+  utiltbx_timer_lastiter(end+1) = 0;
+  utiltbx_timer_name{end+1} = timername;
+  tid = size(utiltbx_timer_start_times,1);
+end
+
+
+% detect timer calling function %
+
+st = dbstack;
+if (length(dbstack) >= 2)
+  utiltbx_timer_callfun{end+1} = st(2).name;
+else
+  utiltbx_timer_callfun{end+1} = '';
+end
diff -r 23f9dd7b9d78 -r ad36f80e2ccf DL/RLS-DLA/private/updateFigure.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/RLS-DLA/private/updateFigure.m	Tue Mar 15 12:20:59 2011 +0000
@@ -0,0 +1,93 @@
+function updateFigure(opts, figTitle, figFilename)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: updateFigure.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+% Ensure default values are available
+opts.linewidth  = getOption(opts,'linewidth', []);
+opts.fontsize   = getOption(opts,'fontsize',  []);
+opts.markersize = getOption(opts,'markersize',[]);
+
+% Output the plots
+if opts.update
+  % Set the line width, font size and marker size
+  chld = [gca; get(gca,'Children')];
+  lnwd = ones(length(chld),1) * NaN;
+  fnts = ones(length(chld),1) * NaN;
+  mrks = ones(length(chld),1) * NaN;
+  for i=1:length(chld)
+    conf = get(chld(i));
+    if ~isempty(opts.linewidth) && isfield(conf,'LineWidth')
+      lnwd(i) = get(chld(i),'LineWidth');
+      if (lnwd(i) == 0.5) % Default
+        set(chld(i),'Linewidth',opts.linewidth);
+      end
+    end
+    if ~isempty(opts.fontsize) && isfield(conf,'FontSize')
+      fnts(i) = get(chld(i),'FontSize');
+      if (fnts(i) == 10) % Default
+        set(chld(i),'FontSize',opts.fontsize);
+      end
+    end
+    if ~isempty(opts.markersize) && isfield(conf,'MarkerSize')
+      mrks(i) = get(chld(i),'MarkerSize');
+      if (mrks(i) == 6) % Default
+        set(chld(i),'MarkerSize',opts.markersize);
+      end
+    end
+  end
+  
+  for i=1:length(opts.figtype)
+     updateFigureType(opts.update, 0, opts.figtype{i}, ...
+                      opts.figpath, figTitle, figFilename);
+  end
+
+  % Restore the line-widths, font size
+  for i=1:length(chld)
+    if ~isnan(lnwd(i))
+      set(chld(i),'LineWidth',lnwd(i));
+    end
+    if ~isnan(fnts(i))
+      set(chld(i),'FontSize',fnts(i));
+    end
+    if ~isnan(mrks(i))
+      set(chld(i),'MarkerSize',mrks(i));
+    end
+  end
+    
+end
+
+% Show the plot
+if opts.show
+  updateFigureType(0,opts.show,'','',figTitle,'');
+end
+
+
+
+function updateFigureType(update,show,figtype,figpath,figTitle,figFilename)
+filename = [figpath,figFilename];
+
+switch lower(figtype)
+ case {'pdf'}
+   cmdPostprocess = sprintf('!pdfcrop %s.pdf %s.pdf >& /dev/null', ...
+                            filename, filename);
+ otherwise
+   cmdPostprocess = [];
+end
+
+[figtype,figext] = getFigureExt(figtype);
+
+% Print the figure for output (without title)
+if update
+  evalc(sprintf('print -d%s %s.%s;', figtype, filename, figext));
+  
+  if ~isempty(cmdPostprocess)
+    eval(cmdPostprocess);
+  end
+end
+
+% Add title if needed
+if show
+  title(figTitle);
+end