diff --git a/smop/Makefile b/smop/Makefile
index 4e1d49fd..bdcc5998 100644
--- a/smop/Makefile
+++ b/smop/Makefile
@@ -1,34 +1,21 @@
 # for py3: make PYTHON=python3 CYTHON="cython -3" V=3.4
 
-VPATH = $(SCRIPTS)/general
 OCTAVE = /home/lei/octave-3.8.2
 SCRIPTS = $(OCTAVE)/scripts
 
 CYTHON = cython
 PYTHON = python
 #FLAGS = -S stft.m,datenum.m -X orderfields.m,legend.m,pack.m,unpack.m,__unimplemented__.m,assert.m,optimset.m
-FLAGS = -X interp1.m,interpn.m,nargchk.m,nargoutchk.m,profexplore.m,profshow.m,quadgk.m,quadl.m
+#FLAGS = -X interp1.m,interpn.m,nargchk.m,nargoutchk.m,profexplore.m,profshow.m,quadgk.m,quadl.m
+FLAGS = -Xquadgk.m,quadl.m,lcm.m,adaptlobstp.m,spline.m,ppval.m
 V = 2.7
 ALL = \*.m
-GENERAL = accumarray.m accumdim.m bicubic.m bincoeff.m\
-    bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m cart2sph.m\
-    cell2mat.m celldisp.m chop.m circshift.m common_size.m cplxpair.m\
-    cumtrapz.m curl.m dblquad.m deal.m del2.m display.m divergence.m\
-    fieldnames.m flipdim.m fliplr.m flipud.m gradient.m idivide.m\
-    int2str.m interp1.m interp2.m interp3.m interpft.m interpn.m\
-    isa.m iscolumn.m isdir.m isequal.m isequaln.m isrow.m isscalar.m\
-    issquare.m isvector.m loadobj.m logspace.m methods.m nargchk.m\
-    narginchk.m nargoutchk.m nextpow2.m nthargout.m num2str.m pol2cart.m\
-    polyarea.m postpad.m prepad.m profexplore.m profile.m profshow.m\
-    quadgk.m quadl.m quadv.m randi.m rat.m repmat.m rot90.m rotdim.m\
-    saveobj.m shiftdim.m shift.m sortrows.m sph2cart.m structfun.m\
-    subsindex.m trapz.m triplequad.m
 
 all: mybench.py 
 	$(PYTHON) -c "import mybench ; mybench.mybench()"
 
-general.py: $(GENERAL)
-	$(PYTHON) main.py $^ -o $@ -v $(FLAGS)
+octave.py:
+	find $(SCRIPTS)/general $(SCRIPTS)/specfun -name \*.m | xargs smop -v -o $@ -r core $(FLAGS)
 
 go_so: solver.so
 	$(PYTHON) go.py	
diff --git a/smop/core.py b/smop/core.py
new file mode 100644
index 00000000..5ac7939f
--- /dev/null
+++ b/smop/core.py
@@ -0,0 +1,632 @@
+# SMOP compiler runtime support library
+# Copyright 2014 Victor Leikehman
+
+# MIT license
+
+import __builtin__
+from numpy import abs,ceil,floor,sqrt
+from numpy.fft import fft2
+from numpy.linalg import inv
+from numpy.linalg import qr  as _qr 
+try:
+    from scipy.linalg import schur as _schur
+except ImportError:
+    pass
+
+import numpy as np
+import os,sys,copy,time
+from sys import stdin,stdout,stderr
+try:
+    from scipy.io import loadmat
+except:
+    pass
+import unittest
+
+def isvector_or_scalar(a):
+    """
+    one-dimensional arrays having shape [N],
+    row and column matrices having shape [1 N] and
+    [N 1] correspondingly, and their generalizations
+    having shape [1 1 ... N ... 1 1 1].
+    Scalars have shape [1 1 ... 1].
+    Empty arrays dont count
+    """
+    try:
+        return a.size and a.ndim-a.shape.count(1) <= 1
+    except:
+        return False
+def isvector(a):
+    """
+    one-dimensional arrays having shape [N],
+    row and column matrices having shape [1 N] and
+    [N 1] correspondingly, and their generalizations
+    having shape [1 1 ... N ... 1 1 1]
+    """
+    try:
+        return a.ndim-a.shape.count(1) == 1
+    except:
+        return False
+
+class matlabarray(np.ndarray):
+    """
+    >>> matlabarray()
+    matlabarray([], shape=(0, 0), dtype=float64)
+    """
+
+    def __new__(cls,a=[],dtype=None):
+        obj = np.array(a,
+                       dtype=dtype,
+                       copy=False,
+                       order="F",
+                       ndmin=2).view(cls).copy(order="F")
+        if obj.size == 0:
+            obj.shape = (0,0)
+        return obj
+
+    #def __array_finalize__(self,obj):
+
+    def __copy__(self):
+        return np.ndarray.copy(self,order="F")
+
+    def __iter__(self):
+        """ must define iter or char won't work"""
+        return np.asarray(self).__iter__()
+
+    def compute_indices(self,index):
+        if not isinstance(index,tuple):
+           index = index,
+        if len(index) != 1 and len(index) != self.ndim:
+            raise IndexError
+        indices = []
+        for i,ix in enumerate(index):
+            if ix.__class__ is end:
+                indices.append(self.shape[i]-1+ix.n)
+            elif ix.__class__ is slice:
+                if self.size == 0 and ix.stop is None:
+                    raise IndexError
+                if len(index) == 1:
+                    n = self.size
+                else:
+                    n = self.shape[i]
+                indices.append(np.arange((ix.start or 1)-1,
+                                          ix.stop  or n,
+                                          ix.step  or 1,
+                                          dtype=int))
+            else:
+                try:
+                    indices.append(int(ix)-1)
+                except:
+                    indices.append(np.asarray(ix).astype("int32")-1)
+        if len(indices) == 2 and isvector(indices[0]) and isvector(indices[1]):
+            indices[0].shape = (-1,1)
+            indices[1].shape = (-1,)
+        return tuple(indices)
+
+    def __getslice__(self,i,j):
+        if i == 0 and j == sys.maxsize:
+            return self.reshape(-1,1,order="F")
+        return self.__getitem__(slice(i,j))
+
+    def __getitem__(self,index):
+        return matlabarray(self.get(index))
+
+    def get(self,index):
+        #import pdb; pdb.set_trace()
+        indices = self.compute_indices(index)
+        if len(indices) == 1:
+            return np.ndarray.__getitem__(self.reshape(-1,order="F"),indices)
+        else:
+            return np.ndarray.__getitem__(self,indices)
+
+    def __setslice__(self,i,j,value):
+        if i == 0 and j == sys.maxsize:
+            index = slice(None,None)
+        else:
+            index = slice(i,j)
+        self.__setitem__(index,value)
+        
+    def sizeof(self,ix):
+        if isinstance(ix,int):
+            n = ix+1
+        elif isinstance(ix,slice):
+            n = ix.stop
+        elif isinstance(ix,(list,np.ndarray)):
+            n = max(ix)+1
+        else:
+            assert 0,ix
+        if not isinstance(n,int):
+            raise IndexError
+        return n 
+
+    def __setitem__(self,index,value):
+        #import pdb; pdb.set_trace()
+        indices = self.compute_indices(index)
+        try:
+            if len(indices) == 1:
+                np.asarray(self).reshape(-1,order="F").__setitem__(indices,value)
+            else:
+                np.asarray(self).__setitem__(indices,value)
+        except (ValueError,IndexError):
+            #import pdb; pdb.set_trace()
+            if not self.size:
+                new_shape = [self.sizeof(s) for s in indices]
+                self.resize(new_shape,refcheck=0)
+                np.asarray(self).__setitem__(indices,value)
+            elif len(indices) == 1:
+                # One-dimensional resize is only implemented for
+                # two cases:
+                #
+                # a. empty matrices having shape [0 0]. These
+                #    matries may be resized to any shape.  A[B]=C
+                #    where A=[], and B is specific -- A[1:10]=C
+                #    rather than A[:]=C or A[1:end]=C
+                if self.size and not isvector_or_scalar(self):
+                    raise IndexError("One-dimensional resize "
+                                     "works only on vectors, and "
+                                     "row and column matrices")
+                # One dimensional resize of scalars creates row matrices
+                # ai = 3
+                # a(4) = 1
+                # 3 0 0 1
+                n = self.sizeof(indices[0]) # zero-based
+                if max(self.shape) == 1:
+                    new_shape = list(self.shape)
+                    new_shape[-1] = n
+                else:
+                    new_shape = [(1 if s==1 else n) for s in self.shape]
+                self.resize(new_shape,refcheck=0)
+                np.asarray(self).reshape(-1,order="F").__setitem__(indices,value)
+            else:
+                new_shape = list(self.shape)
+                if self.flags["C_CONTIGUOUS"]:
+                    new_shape[0] = self.sizeof(indices[0])
+                elif self.flags["F_CONTIGUOUS"]:
+                    new_shape[-1] = self.sizeof(indices[-1])
+                self.resize(new_shape,refcheck=0)
+                np.asarray(self).__setitem__(indices,value)
+
+    def __repr__(self):
+        return self.__class__.__name__ + repr(np.asarray(self))[5:]
+
+    def __str__(self):
+        return str(np.asarray(self))
+ 
+    def __add__(self,other):
+        return matlabarray(np.asarray(self)+np.asarray(other))
+
+    def __neg__(self):
+        return matlabarray(np.asarray(self).__neg__())
+
+class end(object):
+    def __add__(self,n):
+        self.n = n
+        return self
+    def __sub__(self,n):
+        self.n = -n
+        return self
+####
+class cellarray(matlabarray):
+    """
+    Cell array corresponds to matlab ``{}``
+
+
+    """
+
+    def __new__(cls, a=[]):
+        """
+        Create a cell array and initialize it with a.
+        Without arguments, create an empty cell array.
+
+        Parameters:
+        a : list, ndarray, matlabarray, etc.
+
+        >>> a=cellarray([123,"hello"])
+        >>> print a.shape
+        (1, 2)
+
+        >>> print a[1]
+        123
+
+        >>> print a[2]
+        hello
+        """
+        obj = np.array(a,
+                       dtype=object,
+                       order="F",
+                       ndmin=2).view(cls).copy(order="F")
+        if obj.size == 0:
+            obj.shape = (0,0)
+        return obj
+
+    def __getitem__(self,index): 
+        return self.get(index)
+
+#    def __str__(self):
+#        if self.ndim == 0:
+#            return ""
+#        if self.ndim == 1:
+#            return "".join(s for s in self)
+#        if self.ndim == 2:
+#            return "\n".join("".join(s) for s in self)
+#        raise NotImplementedError
+
+
+class cellstr(matlabarray):
+    """
+    >>> s=cellstr(char('helloworldkitty').reshape(3,5))
+    >>> s
+    cellstr([['hello', 'world', 'kitty']], dtype=object)
+    >>> print s
+    hello
+    world
+    kitty
+    >>> s.shape
+    (1, 3)
+    """
+
+    def __new__(cls, a):
+        """
+        Given a two-dimensional char object,
+        create a cell array where each cell contains
+        a line.
+        """
+        obj = np.array(["".join(s) for s in a], 
+                       dtype=object,
+                       copy=False,
+                       order="C",
+                       ndmin=2).view(cls).copy(order="F")
+        if obj.size == 0:
+            obj.shape = (0,0)
+        return obj
+
+    def __str__(self):
+        return "\n".join("".join(s) for s in self.reshape(-1))
+
+    def __getitem__(self,index): 
+        return self.get(index)
+
+
+class char(matlabarray):
+    """
+    class char is a rectangular string matrix, which
+    inherits from matlabarray all its features except
+    dtype.
+
+    >>> s=char()
+    >>> s.shape
+    (0, 0)
+
+    >>> s=char('helloworld').reshape(2,5)
+    >>> print s
+    hello
+    world
+
+    >>> s.shape
+    (2, 5)
+    """
+
+    def __new__(cls, a=""):
+        if not isinstance(a,str):
+            raise NotImplementedError
+        obj = np.array(list(a),
+                       dtype='|S1',
+                       copy=False,
+                       order="F",
+                       ndmin=2).view(cls).copy(order="F")
+        if obj.size == 0:
+            obj.shape = (0,0)
+        return obj
+
+    def __getitem__(self,index): 
+        return self.get(index)
+
+    def __str__(self):
+        if self.ndim == 0:
+            return ""
+        if self.ndim == 1:
+            return "".join(s for s in self)
+        if self.ndim == 2:
+            return "\n".join("".join(s) for s in self)
+        raise NotImplementedError
+
+class struct_(object):
+    def __init__(self,*args):
+        for i in range(0,len(args),2):
+            setattr(self,str(args[i]),args[i+1])
+
+def arange(start,stop,step=1,**kwargs):
+    """
+    >>> a=arange(1,10) # 1:10
+    >>> size(a)
+    matlabarray([[ 1, 10]])
+    """
+    return matlabarray(np.arange(start,
+                                 stop+1,
+                                 step,
+                                 **kwargs).reshape(1,-1),**kwargs)
+
+def cell(*args):
+    if len(args) == 1:
+        args += args
+    return cellarray(np.zeros(args,dtype=object,order="F"))
+
+def clc():
+    pass
+
+def copy(a):
+    return matlabarray(np.asanyarray(a).copy(order="F"))
+
+def disp(*args):
+    print (args)
+
+def eig(a):
+    u,v = numpy.linalg.eig(a)
+    return u.T
+
+def exist(a,b):
+    if str(b) == 'builtin':
+        return str(a) in globals()
+    if str(b) == 'file':
+        return os.path.exists(str(a))
+    raise NotImplementedError
+
+def false(*args):
+    if not args:
+        return False # or matlabarray(False) ???
+    if len(args) == 1: 
+        args += args
+    return np.zeros(args,dtype=bool,order="F")
+
+def find(a,n=None,d=None,nargout=1):
+    if d:
+        raise NotImplementedError
+
+    # there is no promise that nonzero or flatnonzero
+    # use or will use indexing of the argument without
+    # converting it to array first.  So we use asarray
+    # instead of asanyarray
+    if nargout == 1:
+        i = np.flatnonzero(np.asarray(a)).reshape(1,-1)+1
+        if n is not None:
+            i = i.take(n)
+        return matlabarray(i)
+    if nargout == 2:
+        i,j = np.nonzero(np.asarray(a))
+        if n is not None:
+            i = i.take(n)
+            j = j.take(n)
+        return (matlabarray((i+1).reshape(-1,1)),
+                matlabarray((j+1).reshape(-1,1)))
+    raise NotImplementedError
+
+def fopen(*args):
+    try:
+        fp = open(*args)
+        assert fp != -1
+        return fp
+    except:
+        return -1
+
+def fflush(fp):
+    fp.flush()
+
+def fprintf(fp,fmt,*args):
+    if not isinstance(fp,file):
+        fp = stdout
+    fp.write(str(fmt) % args)
+
+def fullfile(*args):
+    return os.path.join(*args)
+
+# implemented in "scripts/set/intersect.m"
+#def intersect(a,b,nargout=1):
+#    if nargout == 1:
+#        c = sorted(set(a) & set(b))
+#        if isinstance(a,str):
+#            return "".join(c)
+#        elif isinstance(a,list):
+#            return c
+#        else:
+#            # FIXME: the result is a column vector if
+#            # both args are column vectors; otherwise row vector
+#            return np.array(c)
+#    raise NotImplementedError
+#
+def iscellstr(a):
+    # TODO return isinstance(a,cellarray) and all(ischar(t) for t in a.flat)
+    return isinstance(a,cellarray) and all(isinstance(t,str) for t in a.flat)
+
+def ischar(a):
+    try:
+        return a.dtype == "|S1"
+    except AttributeError:
+        return False
+# ----------------------------------------------------
+def isempty(a):
+    try:
+        return 0 in np.asarray(a).shape
+    except AttributeError:
+        return False
+
+def isequal(a,b):
+    return np.array_equal(np.asanyarray(a),
+                          np.asanyarray(b))
+                          
+def isfield(a,b):
+    return str(b) in a.__dict__.keys()
+
+def isscalar(a):
+    """np.isscalar returns True if a.__class__ is a scalar
+    type (i.e., int, and also immutable containers str and
+    tuple, but not list.) Our requirements are different"""
+    try:
+        return a.size == 1
+    except AttributeError:
+        return np.isscalar(a)
+
+def length(a):
+    try:
+        return __builtin__.max(np.asarray(a).shape)
+    except ValueError:
+        return 1
+
+try:
+    def load(a):
+        return loadmat(a) # FIXME
+except:
+    pass
+
+def max(a, d=0, nargout=0):
+    if d or nargout:
+        raise NotImplementedError
+    return np.amax(a)
+
+def min(a, d=0, nargout=0):
+    if d or nargout:
+        raise NotImplementedError
+    return np.amin(a)
+
+def mod(a,b):
+    try:
+        return a % b
+    except ZeroDivisionError:
+        return a
+
+def ndims(a):
+    return np.asarray(a).ndim
+
+def numel(a):
+    return np.asarray(a).size
+
+def ones(*args,**kwargs):
+    if not args:
+        return 1.0
+    if len(args) == 1:
+        args += args
+    return matlabarray(np.ones(args,order="F",**kwargs))
+
+#def primes2(upto):
+#    primes=np.arange(2,upto+1)
+#    isprime=np.ones(upto-1,dtype=bool)
+#    for factor in primes[:int(math.sqrt(upto))]:
+#        if isprime[factor-2]: isprime[factor*2-2::factor]=0
+#    return primes[isprime]
+#
+#def primes(*args):
+#    return _primes.primes(*args)
+
+def qr(a):
+    return matlabarray(_qr(np.asarray(a)))
+
+def rand(*args,**kwargs):
+    if not args:
+        return np.random.rand()
+    if len(args) == 1:
+        args += args
+    try:
+        return np.random.rand(np.prod(args)).reshape(args,order="F")
+    except:
+        pass
+
+def rand(*args,**kwargs):
+    if not args:
+        return np.random.rand()
+    if len(args) == 1:
+        args += args
+    try:
+        return np.random.rand(np.prod(args)).reshape(args,order="F")
+    except:
+        pass
+
+def randn(*args,**kwargs):
+    if not args:
+        return np.random.randn()
+    if len(args) == 1:
+        args += args
+    try:
+        return np.random.randn(np.prod(args)).reshape(args,order="F")
+    except:
+        pass
+
+def ravel(a):
+    return np.asanyarray(a).reshape(-1,1)
+
+def roots(a):
+
+    return matlabarray(np.roots(np.asarray(a).ravel()))
+
+def round(a):
+    return np.round(np.asanyarray(a))
+
+def rows(a):
+    return np.asarray(a).shape[0]
+
+def schur(a):
+    return matlabarray(_schur(np.asarray(a)))
+
+def size(a, b=0, nargout=1):
+    """
+    >>> size(zeros(3,3)) + 1
+    matlabarray([[4, 4]])
+    """
+    s = np.asarray(a).shape
+    if s is ():
+        return 1 if b else (1,)*nargout
+    # a is not a scalar
+    try:
+        if b:
+            return s[b-1]
+        else:
+            return matlabarray(s) if nargout <= 1 else s
+    except IndexError:
+        return 1
+
+from numpy import sqrt
+sort = __builtin__.sorted
+
+def strread(s, format="", nargout=1):
+    if format == "":
+        a = [float(x) for x in s.split()]
+        return tuple(a) if nargout > 1 else np.asanyarray([a])
+    raise NotImplementedError
+
+def strrep(a,b,c):
+    return str(a).replace(str(b),str(c))
+
+def sum(a, dim=None):
+    if dim is None:
+        return np.asanyarray(a).sum()
+    else:
+        return np.asanyarray(a).sum(dim-1)
+
+def toupper(a):
+    return char(str(a.data).upper())
+
+true = True
+
+def tic():
+    return time.clock()
+
+def toc(t):
+    return time.clock()-t
+
+def true(*args):
+    if len(args) == 1:
+        args += args
+    return matlabarray(np.ones(args,dtype=bool,order="F"))
+
+def version():
+    return char('0.26')
+
+def zeros(*args,**kwargs):
+    if not args:
+        return 0.0
+    if len(args) == 1:
+        args += args
+    return matlabarray(np.zeros(args,**kwargs))
+
+if __name__ == "__main__":
+    import doctest
+    doctest.testmod()
+
+# vim:et:sw=4:si:tw=60
diff --git a/smop/mybench.py b/smop/mybench.py
new file mode 100644
index 00000000..ebd1a3e5
--- /dev/null
+++ b/smop/mybench.py
@@ -0,0 +1,248 @@
+# Autogenerated with SMOP version 0.26.3-15-gcf44923
+# /home/lei/.local/bin/smop -r core -r octave -o mybench.py mybench.m
+from __future__ import division
+from core import *
+from octave import *
+def mybench(*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 0-[].count(None)+len(args)
+
+    global IS_OCTAVE,IS_IMAGE_PROCESSING
+    conf=struct_(char('noOfRepeats'),3,char('normalize'),true,char('imagep'),true,char('onlyTests'),[])
+    conf=getargs(conf,varargin)
+    IS_OCTAVE=exist(char('OCTAVE_VERSION'),char('builtin')) > 0
+    IS_IMAGE_PROCESSING=copy(false)
+    NO_REPETITIONS=conf.noOfRepeats
+    NORMALIZE_TIMES=conf.normalize
+    if exist(char('imrotate'),char('file')) > 0 and exist(char('imresize'),char('file')) > 0 and exist(char('imerode'),char('file')) > 0 and conf.imagep == true:
+        disp(char('Image processing toolbox found'))
+        IS_IMAGE_PROCESSING=copy(true)
+    if conf.noOfRepeats < 1:
+        conf.noOfRepeats = copy(1)
+    clc
+    mytests=getBenchmarkTests()
+    noOftests=length(mytests)
+    moVersio=strrep(version(),char(' '),char('_'))
+    outFileName=char('foo.txt')
+    fid=fopen(outFileName,char('w'))
+    if NORMALIZE_TIMES:
+        fprintf(fid,char('%s\t%s\t%s\n'),[char('Name_'),moVersio],char('Time'),char('Norm_time'))
+    else:
+        fprintf(fid,char('%s\t%s\n'),[char('Name_'),moVersio],char('Time_[s]'))
+    avarage_time=0
+    times_vector=matlabarray([])
+    times_vector1=matlabarray([])
+    if isempty(conf.onlyTests):
+        doTheseTests=arange(1,noOftests)
+    else:
+        doTheseTests=conf.onlyTests
+        noOftests=length(conf.onlyTests)
+    for i in doTheseTests.reshape(-1):
+        fprintf(1,char('Execute test %d/%d -  %s\n'),i,noOftests,mytests[i].name)
+        if IS_OCTAVE:
+            fflush(stdout)
+        try:
+            x=mytests[i].input()
+            cumulative_time=0
+            cumulative_time1=0
+            goldResult=1
+            for ii in arange(1,NO_REPETITIONS).reshape(-1):
+                fprintf(1,char('%d '),ii)
+                if IS_OCTAVE:
+                    fflush(stdout)
+                t0=tic()
+                mytests[i].test(x)
+                t1=toc(t0)
+                if isfield(mytests[i],char('goldResult')) and NORMALIZE_TIMES == true:
+                    goldResult=mytests[i].goldResult
+                cumulative_time=cumulative_time + t1 / goldResult
+                cumulative_time1=cumulative_time1 + t1
+            avarage_time=cumulative_time / NO_REPETITIONS
+            avarage_time1=cumulative_time1 / NO_REPETITIONS
+            times_vector[end() + 1]=avarage_time
+            times_vector1[end() + 1]=avarage_time1
+        finally:
+            pass
+        if isfield(mytests[i],char('post')):
+            mytests[i].post()
+        fprintf(1,char('\n\tTime %.2f [s]\n'),avarage_time1)
+        if NORMALIZE_TIMES == true:
+            fprintf(1,char('\tNormalized time %.2f \n'),avarage_time)
+        fprintf(1,char('\n'))
+        if IS_OCTAVE:
+            fflush(stdout)
+        if NORMALIZE_TIMES:
+            fprintf(fid,char('%s\t%f\t%f\n'),mytests[i].name,avarage_time1,avarage_time)
+        else:
+            fprintf(fid,char('%s\t%f\n'),mytests[i].name,avarage_time1)
+    times_product=prod(times_vector)
+    times_mean=mean(times_vector)
+    times_geometric_mean=times_product ** (1 / length(times_vector))
+    times_product1=prod(times_vector1)
+    times_mean1=mean(times_vector1)
+    times_geometric_mean1=times_product1 ** (1 / length(times_vector1))
+    res.norm_geometric_mean = copy(times_geometric_mean)
+    res.norm_arithmetic_mean = copy(times_mean)
+    res.norm_min = copy(min(times_vector))
+    res.norm_max = copy(max(times_vector))
+    fprintf(1,char('\n\t --- SUMMARY ---\n'))
+    fprintf(1,char('\n\tMean: geometric %.3f [s], arithmetic %.3f [s]'),times_geometric_mean1,times_mean1)
+    fprintf(1,char('\n\tMin %.3f [s], Max %.3f [s]\n\n'),min(times_vector1),max(times_vector1))
+    if NORMALIZE_TIMES == true:
+        fprintf(1,char('\tNormalized Mean: geometric %.3f, arithmetic %.3f'),times_geometric_mean,times_mean)
+        fprintf(1,char('\n\tNormalized Min %.3f [s], Max %.3f [s]\n\n'),min(times_vector),max(times_vector))
+    if IS_OCTAVE:
+        fflush(stdout)
+    if NORMALIZE_TIMES:
+        fprintf(fid,char('%s\t%f\t%f\n'),char('Geom_mean'),times_geometric_mean1,times_geometric_mean)
+    else:
+        fprintf(fid,char('%s\t%f\t%f\n'),char('Geom_mean'),times_geometric_mean1)
+    fclose(fid)
+    disp(char(''))
+    disp([char('End of test. File '),outFileName,char(' was created.')])
+    if exist(char('out_gray.png')):
+        delete(char('out_gray.png'))
+    if exist(char('out_1.png')):
+        delete(char('out_1.png'))
+    if exist(char('out_mtx')):
+        delete(char('out_mtx'))
+    if exist(char('out_mtx.mat')):
+        delete(char('out_mtx.mat'))
+    if exist(char('dlm.txt')):
+        delete(char('dlm.txt'))
+    clear(char('IS_OCTAVE'),char('IS_IMAGE_PROCESSING'))
+    return res
+def getBenchmarkTests(*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 0-[].count(None)+len(args)
+
+    global IS_OCTAVE,IS_IMAGE_PROCESSING
+    s=cellarray([])
+    s[end() + 1]=struct_(char('name'),char('rand'),char('test'),lambda x: rand(x),char('input'),lambda : 4000,char('goldResult'),1.12)
+    s[end() + 1]=struct_(char('name'),char('randn'),char('test'),lambda x: randn(x),char('input'),lambda : 4000,char('goldResult'),0.58)
+    s[end() + 1]=struct_(char('name'),char('primes'),char('test'),lambda x: primes(x),char('input'),lambda : 10000000.0,char('goldResult'),0.74)
+    s[end() + 1]=struct_(char('name'),char('fft2'),char('test'),lambda x: fft2(x),char('input'),lambda : rand(3000),char('goldResult'),2.28)
+    s[end() + 1]=struct_(char('name'),char('square'),char('test'),lambda x: x ** 2,char('input'),lambda : rand(1000),char('goldResult'),1.35)
+    s[end() + 1]=struct_(char('name'),char('inv'),char('test'),lambda x: inv(x),char('input'),lambda : rand(1000),char('goldResult'),0.87)
+    s[end() + 1]=struct_(char('name'),char('eig'),char('test'),lambda x: eig(x),char('input'),lambda : rand(1000),char('goldResult'),9.45)
+    s[end() + 1]=struct_(char('name'),char('qr'),char('test'),lambda x: qr(x),char('input'),lambda : rand(1000),char('goldResult'),0.79)
+    s[end() + 1]=struct_(char('name'),char('schur'),char('test'),lambda x: schur(x),char('input'),lambda : rand(600),char('goldResult'),2.67)
+    s[end() + 1]=struct_(char('name'),char('roots'),char('test'),lambda x: roots(x),char('input'),lambda : rand(600,1),char('goldResult'),2.08)
+    s[end() + 1]=struct_(char('name'),char('interp2'),char('test'),lambda x: interp2(x,2,char('spline')),char('input'),lambda : rand(600),char('goldResult'),4.2)
+    s[end() + 1]=struct_(char('name'),char('binary'),char('test'),lambda x: eye(x) < 1,char('input'),lambda : 5000,char('goldResult'),0.51)
+    s[end() + 1]=struct_(char('name'),char('forLoop'),char('test'),lambda x: forLoop(x),char('input'),lambda : 200,char('goldResult'),0.06)
+    s[end() + 1]=struct_(char('name'),char('makeSparse'),char('test'),lambda x: sparse(x),char('input'),lambda : eye(5000),char('goldResult'),0.49)
+    s[end() + 1]=struct_(char('name'),char('multiplySparse'),char('test'),lambda x: sparse(x) * sparse(x),char('input'),lambda : eye(5000) * rand(1),char('goldResult'),0.98)
+    s[end() + 1]=struct_(char('name'),char('sLinearEq'),char('test'),lambda x: magic(x) / rand(1,x),char('input'),lambda : 2000,char('goldResult'),1.94)
+    s[end() + 1]=struct_(char('name'),char('saveLoadMtx'),char('test'),lambda x: saveLoadMtx(x),char('input'),lambda : rand(1000),char('goldResult'),0.93)
+    s[end() + 1]=struct_(char('name'),char('dlmwriteRead'),char('test'),lambda x: dlmwriteRead(x),char('input'),lambda : rand(500),char('goldResult'),5.03)
+    s[end() + 1]=struct_(char('name'),char('median'),char('test'),lambda x: median(x[:]),char('input'),lambda : rand(4000),char('goldResult'),3.32)
+    s[end() + 1]=struct_(char('name'),char('std'),char('test'),lambda x: std(x[:]),char('input'),lambda : rand(4000),char('goldResult'),0.84)
+    if IS_IMAGE_PROCESSING:
+        s[end() + 1]=struct_(char('name'),char('doImgAndSaveAsPNG'),char('test'),lambda x: doImgAndSaveAsPNG(x),char('input'),lambda : rand(1500,1500,3),char('post'),lambda : pause(2),char('goldResult'),2.0)
+        s[end() + 1]=struct_(char('name'),char('imageToGray'),char('test'),lambda I: imageToGray(I),char('input'),lambda : imread(char('out_1.png')),char('goldResult'),0.56)
+        s[end() + 1]=struct_(char('name'),char('imageRotate'),char('test'),lambda I: imageRotate(I),char('input'),lambda : imread(char('out_gray.png')),char('goldResult'),2.94)
+        s[end() + 1]=struct_(char('name'),char('imresize'),char('test'),lambda I: imresize(I,1.2),char('input'),lambda : imread(char('out_gray.png')),char('goldResult'),1.24)
+        s[end() + 1]=struct_(char('name'),char('imageFilter'),char('test'),lambda I: imageFilter(I),char('input'),lambda : imread(char('out_gray.png')),char('goldResult'),0.2)
+        s[end() + 1]=struct_(char('name'),char('imageErode'),char('test'),lambda I: imageErode(I),char('input'),lambda : imread(char('out_gray.png')),char('goldResult'),0.46)
+        s[end() + 1]=struct_(char('name'),char('medfilt2'),char('test'),lambda I: medfilt2(I),char('input'),lambda : magic(2000),char('goldResult'),1.03)
+    return s
+def saveLoadMtx(out_mtx=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[out_mtx].count(None)+len(args)
+
+    save(char('out_mtx'))
+    clear(char('out_mtx'))
+    load(char('out_mtx'))
+    return
+def dlmwriteRead(x=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[x].count(None)+len(args)
+
+    dlmwrite(char('dlm.txt'),x)
+    dlmread(char('dlm.txt'))
+    return
+def forLoop(x=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[x].count(None)+len(args)
+
+    for i in arange(1,x).reshape(-1):
+        for j in arange(1,x).reshape(-1):
+            for k in arange(1,x).reshape(-1):
+                i + j + k
+    return
+def doImgAndSaveAsPNG(x=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[x].count(None)+len(args)
+
+    imwrite(x,char('out_1.png'))
+    return
+def solveNonLinearEq(*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 0-[].count(None)+len(args)
+
+    x,fval,info=fsolve(equationsToSolve,[[1],[1],[1],[1]],nargout=3)
+    return
+def equationsToSolve(x=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[x].count(None)+len(args)
+
+    y[1]=- 2 * x[1] ** 2 + 3 * x[1] * x[2] + 4 * sin(x[2]) + log(x[3]) - 6
+    y[2]=3 * x[1] ** 2 - 2 * x[2] * x[2] ** 2 + 3 * cos(x[1]) + 4
+    y[3]=1 * x[1] ** 2 - 2 * x[1] * x[2] * x[4] ** 2 + 3 * cos(x[1]) + 4
+    y[4]=1 * x[1] ** 2 - 2 * x[1] * x[2] * x[3] ** 2 + 3 * cos(x[4]) + 4
+    return y
+def imageToGray(I=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[I].count(None)+len(args)
+
+    Igray=rgb2gray(I)
+    imwrite(Igray,char('out_gray.png'))
+    return
+def imageRotate(I=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[I].count(None)+len(args)
+
+    I2=imrotate(I,2)
+    return
+def imageFilter(I=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[I].count(None)+len(args)
+
+    h=fspecial(char('sobel'))
+    filteredI=imfilter(I,h)
+    return
+def imageErode(I=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 1-[I].count(None)+len(args)
+
+    SE=eye(5)
+    erodedI=imerode(I,SE)
+    return
+def getargs(defs=None,varglist=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 2-[defs,varglist].count(None)+len(args)
+
+    l=length(varglist)
+    if l == 0:
+        return defs
+    if mod(l,2) != 0:
+        disp(char(' !!! Odd number of parameters !!!'))
+        defs=cellarray([])
+        return defs
+    varnames=cellarray([varglist[1:2:l]])
+    varvalues=cellarray([varglist[2:2:l]])
+    given_vars=zeros(1,l / 2)
+    for i in arange(1,l / 2,1).reshape(-1):
+        existss=isfield(defs,varnames[i])
+        given_vars[i]=existss
+    if min(given_vars) == 0:
+        disp(char('!!! No such parameter(s):'))
+        disp(varnames[not given_vars])
+        defs=cellarray([])
+        return defs
+    for i in arange(1,l / 2,1).reshape(-1):
+        setattr(defs,varnames[i],varvalues[i])
+    return defs
+
+mybench()
diff --git a/smop/r8_random.m b/smop/r8_random.m
new file mode 100644
index 00000000..82599ae1
--- /dev/null
+++ b/smop/r8_random.m
@@ -0,0 +1,55 @@
+function [ r, s1, s2, s3 ] = r8_random ( s1, s2, s3 )
+
+%*****************************************************************************80
+%
+%% R8_RANDOM returns a pseudorandom number between 0 and 1.
+%
+%  Discussion:
+%
+%    This function returns a pseudo-random number rectangularly distributed
+%    between 0 and 1.   The cycle length is 6.95E+12.  (See page 123
+%    of Applied Statistics (1984) volume 33), not as claimed in the
+%    original article.
+%
+%  Licensing:
+%
+%    This code is distributed under the GNU LGPL license.
+%
+%  Modified:
+%
+%    08 July 2008
+%
+%  Author:
+%
+%    Original FORTRAN77 version by Brian Wichman, David Hill.
+%    MATLAB version by John Burkardt.
+%
+%  Reference:
+%
+%    Brian Wichman, David Hill,
+%    Algorithm AS 183: An Efficient and Portable Pseudo-Random
+%    Number Generator,
+%    Applied Statistics,
+%    Volume 31, Number 2, 1982, pages 188-190.
+%
+%  Parameters:
+%
+%    Input, integer S1, S2, S3, three values used as the
+%    seed for the sequence.  These values should be positive
+%    integers between 1 and 30,000.
+%
+%    Output, real R, the next value in the sequence.
+%
+%    Output, integer S1, S2, S3, updated seed values.
+%
+  s1 = mod ( 171 * s1, 30269 );
+  s2 = mod ( 172 * s2, 30307 );
+  s3 = mod ( 170 * s3, 30323 );
+
+  r = mod ( s1 / 30269.0 ...
+          + s2 / 30307.0 ...
+          + s3 / 30323.0, 1.0 );
+
+  return
+end
+