initial sd-vbs

initial sd-vbs

add sd-vbs

sd-vbs

31 files changed, 1586 insertions, 0 deletions

diff --git a/SD-VBS/common/matlab/cycle.h b/SD-VBS/common/matlab/cycle.h
new file mode 100644
index 0000000..2652a04
--- /dev/null
+++ b/SD-VBS/common/matlab/cycle.h
@@ -0,0 +1,514 @@
+/*
+ * Copyright (c) 2003, 2007-8 Matteo Frigo
+ * Copyright (c) 2003, 2007-8 Massachusetts Institute of Technology
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ *
+ */
+
+
+/* machine-dependent cycle counters code. Needs to be inlined. */
+
+/***************************************************************************/
+/* To use the cycle counters in your code, simply #include "cycle.h" (this
+   file), and then use the functions/macros:
+
+                 ticks getticks(void);
+
+   ticks is an opaque typedef defined below, representing the current time.
+   You extract the elapsed time between two calls to gettick() via:
+
+                 double elapsed(ticks t1, ticks t0);
+
+   which returns a double-precision variable in arbitrary units.  You
+   are not expected to convert this into human units like seconds; it
+   is intended only for *comparisons* of time intervals.
+
+   (In order to use some of the OS-dependent timer routines like
+   Solaris' gethrtime, you need to paste the autoconf snippet below
+   into your configure.ac file and #include "config.h" before cycle.h,
+   or define the relevant macros manually if you are not using autoconf.)
+*/
+
+/***************************************************************************/
+/* This file uses macros like HAVE_GETHRTIME that are assumed to be
+   defined according to whether the corresponding function/type/header
+   is available on your system.  The necessary macros are most
+   conveniently defined if you are using GNU autoconf, via the tests:
+   
+   dnl ---------------------------------------------------------------------
+
+   AC_C_INLINE
+   AC_HEADER_TIME
+   AC_CHECK_HEADERS([sys/time.h c_asm.h intrinsics.h mach/mach_time.h])
+
+   AC_CHECK_TYPE([hrtime_t],[AC_DEFINE(HAVE_HRTIME_T, 1, [Define to 1 if hrtime_t is defined in <sys/time.h>])],,[#if HAVE_SYS_TIME_H
+#include <sys/time.h>
+#endif])
+
+   AC_CHECK_FUNCS([gethrtime read_real_time time_base_to_time clock_gettime mach_absolute_time])
+
+   dnl Cray UNICOS _rtc() (real-time clock) intrinsic
+   AC_MSG_CHECKING([for _rtc intrinsic])
+   rtc_ok=yes
+   AC_TRY_LINK([#ifdef HAVE_INTRINSICS_H
+#include <intrinsics.h>
+#endif], [_rtc()], [AC_DEFINE(HAVE__RTC,1,[Define if you have the UNICOS _rtc() intrinsic.])], [rtc_ok=no])
+   AC_MSG_RESULT($rtc_ok)
+
+   dnl ---------------------------------------------------------------------
+*/
+
+/***************************************************************************/
+
+#if TIME_WITH_SYS_TIME
+# include <sys/time.h>
+# include <time.h>
+#else
+# if HAVE_SYS_TIME_H
+#  include <sys/time.h>
+# else
+#  include <time.h>
+# endif
+#endif
+
+#define INLINE_ELAPSED(INL) static INL double elapsed(ticks t1, ticks t0) \
+{                                                                         \
+     return (double)t1 - (double)t0;                                      \
+}
+
+/*----------------------------------------------------------------*/
+/* Solaris */
+#if defined(HAVE_GETHRTIME) && defined(HAVE_HRTIME_T) && !defined(HAVE_TICK_COUNTER)
+typedef hrtime_t ticks;
+
+#define getticks gethrtime
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* AIX v. 4+ routines to read the real-time clock or time-base register */
+#if defined(HAVE_READ_REAL_TIME) && defined(HAVE_TIME_BASE_TO_TIME) && !defined(HAVE_TICK_COUNTER)
+typedef timebasestruct_t ticks;
+
+static __inline ticks getticks(void)
+{
+     ticks t;
+     read_real_time(&t, TIMEBASE_SZ);
+     return t;
+}
+
+static __inline double elapsed(ticks t1, ticks t0) /* time in nanoseconds */
+{
+     time_base_to_time(&t1, TIMEBASE_SZ);
+     time_base_to_time(&t0, TIMEBASE_SZ);
+     return (((double)t1.tb_high - (double)t0.tb_high) * 1.0e9 + 
+             ((double)t1.tb_low - (double)t0.tb_low));
+}
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * PowerPC ``cycle'' counter using the time base register.
+ */
+#if ((((defined(__GNUC__) && (defined(__powerpc__) || defined(__ppc__))) || (defined(__MWERKS__) && defined(macintosh)))) || (defined(__IBM_GCC_ASM) && (defined(__powerpc__) || defined(__ppc__))))  && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     unsigned int tbl, tbu0, tbu1;
+
+     do {
+          __asm__ __volatile__ ("mftbu %0" : "=r"(tbu0));
+          __asm__ __volatile__ ("mftb %0" : "=r"(tbl));
+          __asm__ __volatile__ ("mftbu %0" : "=r"(tbu1));
+     } while (tbu0 != tbu1);
+
+     return (((unsigned long long)tbu0) << 32) | tbl;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* MacOS/Mach (Darwin) time-base register interface (unlike UpTime,
+   from Carbon, requires no additional libraries to be linked). */
+#if defined(HAVE_MACH_ABSOLUTE_TIME) && defined(HAVE_MACH_MACH_TIME_H) && !defined(HAVE_TICK_COUNTER)
+#include <mach/mach_time.h>
+typedef uint64_t ticks;
+#define getticks mach_absolute_time
+INLINE_ELAPSED(__inline__)
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * Pentium cycle counter 
+ */
+#if (defined(__GNUC__) || defined(__ICC)) && defined(__i386__)  && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+
+     __asm__ __volatile__("rdtsc": "=A" (ret));
+     /* no input, nothing else clobbered */
+     return ret;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#define TIME_MIN 5000.0   /* unreliable pentium IV cycle counter */
+#endif
+
+/* Visual C++ -- thanks to Morten Nissov for his help with this */
+#if _MSC_VER >= 1200 && _M_IX86 >= 500 && !defined(HAVE_TICK_COUNTER)
+#include <windows.h>
+typedef LARGE_INTEGER ticks;
+#define RDTSC __asm __emit 0fh __asm __emit 031h /* hack for VC++ 5.0 */
+
+static __inline ticks getticks(void)
+{
+     ticks retval;
+
+     __asm {
+          RDTSC
+          mov retval.HighPart, edx
+          mov retval.LowPart, eax
+     }
+     return retval;
+}
+
+static __inline double elapsed(ticks t1, ticks t0)
+{  
+     return (double)t1.QuadPart - (double)t0.QuadPart;
+}  
+
+#define HAVE_TICK_COUNTER
+#define TIME_MIN 5000.0   /* unreliable pentium IV cycle counter */
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * X86-64 cycle counter
+ */
+#if (defined(__GNUC__) || defined(__ICC) || defined(__SUNPRO_C)) && defined(__x86_64__)  && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     unsigned a, d; 
+     asm volatile("rdtsc" : "=a" (a), "=d" (d)); 
+     return ((ticks)a) | (((ticks)d) << 32); 
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* PGI compiler, courtesy Cristiano Calonaci, Andrea Tarsi, & Roberto Gori.
+   NOTE: this code will fail to link unless you use the -Masmkeyword compiler
+   option (grrr). */
+#if defined(__PGI) && defined(__x86_64__) && !defined(HAVE_TICK_COUNTER) 
+typedef unsigned long long ticks;
+static ticks getticks(void)
+{
+    asm(" rdtsc; shl    $0x20,%rdx; mov    %eax,%eax; or     %rdx,%rax;    ");
+}
+INLINE_ELAPSED(__inline__)
+#define HAVE_TICK_COUNTER
+#endif
+
+/* Visual C++, courtesy of Dirk Michaelis */
+#if _MSC_VER >= 1400 && (defined(_M_AMD64) || defined(_M_X64)) && !defined(HAVE_TICK_COUNTER)
+
+#include <intrin.h>
+#pragma intrinsic(__rdtsc)
+typedef unsigned __int64 ticks;
+#define getticks __rdtsc
+INLINE_ELAPSED(__inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * IA64 cycle counter
+ */
+
+/* intel's icc/ecc compiler */
+#if (defined(__EDG_VERSION) || defined(__ECC)) && defined(__ia64__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+#include <ia64intrin.h>
+
+static __inline__ ticks getticks(void)
+{
+     return __getReg(_IA64_REG_AR_ITC);
+}
+ 
+INLINE_ELAPSED(__inline__)
+ 
+#define HAVE_TICK_COUNTER
+#endif
+
+/* gcc */
+#if defined(__GNUC__) && defined(__ia64__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+
+     __asm__ __volatile__ ("mov %0=ar.itc" : "=r"(ret));
+     return ret;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* HP/UX IA64 compiler, courtesy Teresa L. Johnson: */
+#if defined(__hpux) && defined(__ia64) && !defined(HAVE_TICK_COUNTER)
+#include <machine/sys/inline.h>
+typedef unsigned long ticks;
+
+static inline ticks getticks(void)
+{
+     ticks ret;
+
+     ret = _Asm_mov_from_ar (_AREG_ITC);
+     return ret;
+}
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* Microsoft Visual C++ */
+#if defined(_MSC_VER) && defined(_M_IA64) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned __int64 ticks;
+
+#  ifdef __cplusplus
+extern "C"
+#  endif
+ticks __getReg(int whichReg);
+#pragma intrinsic(__getReg)
+
+static __inline ticks getticks(void)
+{
+     volatile ticks temp;
+     temp = __getReg(3116);
+     return temp;
+}
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * PA-RISC cycle counter 
+ */
+#if defined(__hppa__) || defined(__hppa) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+
+#  ifdef __GNUC__
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+
+     __asm__ __volatile__("mfctl 16, %0": "=r" (ret));
+     /* no input, nothing else clobbered */
+     return ret;
+}
+#  else
+#  include <machine/inline.h>
+static inline unsigned long getticks(void)
+{
+     register ticks ret;
+     _MFCTL(16, ret);
+     return ret;
+}
+#  endif
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* S390, courtesy of James Treacy */
+#if defined(__GNUC__) && defined(__s390__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks cycles;
+     __asm__("stck 0(%0)" : : "a" (&(cycles)) : "memory", "cc");
+     return cycles;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+/*----------------------------------------------------------------*/
+#if defined(__GNUC__) && defined(__alpha__) && !defined(HAVE_TICK_COUNTER)
+/*
+ * The 32-bit cycle counter on alpha overflows pretty quickly, 
+ * unfortunately.  A 1GHz machine overflows in 4 seconds.
+ */
+typedef unsigned int ticks;
+
+static __inline__ ticks getticks(void)
+{
+     unsigned long cc;
+     __asm__ __volatile__ ("rpcc %0" : "=r"(cc));
+     return (cc & 0xFFFFFFFF);
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+#if defined(__GNUC__) && defined(__sparc_v9__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+     __asm__ __volatile__("rd %%tick, %0" : "=r" (ret));
+     return ret;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+#if (defined(__DECC) || defined(__DECCXX)) && defined(__alpha) && defined(HAVE_C_ASM_H) && !defined(HAVE_TICK_COUNTER)
+#  include <c_asm.h>
+typedef unsigned int ticks;
+
+static __inline ticks getticks(void)
+{
+     unsigned long cc;
+     cc = asm("rpcc %v0");
+     return (cc & 0xFFFFFFFF);
+}
+
+INLINE_ELAPSED(__inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+/*----------------------------------------------------------------*/
+/* SGI/Irix */
+#if defined(HAVE_CLOCK_GETTIME) && defined(CLOCK_SGI_CYCLE) && !defined(HAVE_TICK_COUNTER)
+typedef struct timespec ticks;
+
+static inline ticks getticks(void)
+{
+     struct timespec t;
+     clock_gettime(CLOCK_SGI_CYCLE, &t);
+     return t;
+}
+
+static inline double elapsed(ticks t1, ticks t0)
+{
+     return ((double)t1.tv_sec - (double)t0.tv_sec) * 1.0E9 +
+          ((double)t1.tv_nsec - (double)t0.tv_nsec);
+}
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* Cray UNICOS _rtc() intrinsic function */
+#if defined(HAVE__RTC) && !defined(HAVE_TICK_COUNTER)
+#ifdef HAVE_INTRINSICS_H
+#  include <intrinsics.h>
+#endif
+
+typedef long long ticks;
+
+#define getticks _rtc
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* MIPS ZBus */
+#if HAVE_MIPS_ZBUS_TIMER
+#if defined(__mips__) && !defined(HAVE_TICK_COUNTER)
+#include <sys/mman.h>
+#include <unistd.h>
+#include <fcntl.h>
+
+typedef uint64_t ticks;
+
+static inline ticks getticks(void)
+{
+  static uint64_t* addr = 0;
+
+  if (addr == 0)
+  {
+    uint32_t rq_addr = 0x10030000;
+    int fd;
+    int pgsize;
+
+    pgsize = getpagesize();
+    fd = open ("/dev/mem", O_RDONLY | O_SYNC, 0);
+    if (fd < 0) {
+      perror("open");
+      return NULL;
+    }
+    addr = mmap(0, pgsize, PROT_READ, MAP_SHARED, fd, rq_addr);
+    close(fd);
+    if (addr == (uint64_t *)-1) {
+      perror("mmap");
+      return NULL;
+    }
+  }
+
+  return *addr;
+}
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+#endif /* HAVE_MIPS_ZBUS_TIMER */
+
diff --git a/SD-VBS/common/matlab/fSelfCheck.m b/SD-VBS/common/matlab/fSelfCheck.m
new file mode 100644
index 0000000..0288328
--- /dev/null
+++ b/SD-VBS/common/matlab/fSelfCheck.m
@@ -0,0 +1,27 @@
+function ret = fSelfCheck(in1, path, tol)
+
+r1 = size(in1, 1);
+c1 = size(in1, 2);
+
+ret = 1;
+
+file = [path, '/expected.m'];
+fd = fopen(file, 'r');
+
+[in2, count] = fscanf(fd, '%f');
+
+if(count ~= (r1*c1) )
+    fprintf(1, 'Dimensions mismatch: Expected %d\t Observed %d\n', count, (r1*c1));
+    ret = -1;
+else
+    ret = 1;
+    for i=1:(r1*c1)
+        if( (abs(in1(i)) - abs(in2(i)) > tol) || (abs(in2(i)) - abs(in1(i))) > tol)
+            fprintf(1, 'Checking Error: Index %d\tExpected %f\tObserved %f\n', i, in2(i), in1(i));
+            ret = -1;
+            break;
+        end
+    end
+
+end
+
diff --git a/SD-VBS/common/matlab/fWriteMatrix.m b/SD-VBS/common/matlab/fWriteMatrix.m
new file mode 100644
index 0000000..3995aa6
--- /dev/null
+++ b/SD-VBS/common/matlab/fWriteMatrix.m
@@ -0,0 +1,18 @@
+function fWriteMatrix(input, inpath)
+
+file = [inpath '/expected.m'];
+disp(file);
+fd = fopen(file, 'w');
+
+[rows, cols] = size(input);
+
+for j=1:rows
+    for i=1:cols
+        fprintf(fd, '%f\t', input(j,i));
+    end
+    fprintf(fd, '\n');
+end
+fclose(fd);
+end
+
+
diff --git a/SD-VBS/common/matlab/photonEndTiming.c b/SD-VBS/common/matlab/photonEndTiming.c
new file mode 100644
index 0000000..bbf5231
--- /dev/null
+++ b/SD-VBS/common/matlab/photonEndTiming.c
@@ -0,0 +1,15 @@
+#include"mex.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
+{
+    uint32_t* cycles;
+    plhs[0] = mxCreateNumericMatrix(1, 2, mxUINT32_CLASS, mxREAL);
+    cycles = (uint32_t*)mxGetPr(plhs[0]);
+    __asm__ __volatile__( "rdtsc": "=a" (cycles[0]), "=d" (cycles[1]));    
+ 
+    return;
+}
diff --git a/SD-VBS/common/matlab/photonEndTiming.mexa64 b/SD-VBS/common/matlab/photonEndTiming.mexa64
new file mode 100755
index 0000000..e41dd3e
--- /dev/null
+++ b/SD-VBS/common/matlab/photonEndTiming.mexa64
diff --git a/SD-VBS/common/matlab/photonEndTiming.mexglx b/SD-VBS/common/matlab/photonEndTiming.mexglx
new file mode 100755
index 0000000..fa9d4b6
--- /dev/null
+++ b/SD-VBS/common/matlab/photonEndTiming.mexglx
diff --git a/SD-VBS/common/matlab/photonPrintTiming.m b/SD-VBS/common/matlab/photonPrintTiming.m
new file mode 100644
index 0000000..b7c2abb
--- /dev/null
+++ b/SD-VBS/common/matlab/photonPrintTiming.m
@@ -0,0 +1,11 @@
+%! _photonPrintTiming_NA_i2
+
+function photonPrintTiming(elapsed)
+    if(elapsed(2) == 0)
+        fprintf(1,'Cycles elapsed\t\t- %u\n',elapsed(1));
+    else
+        fprintf(1,'Cycles elapsed\t\t- %u%\u\n',elapsed(2),elapsed(1));
+    end
+end
+
+
diff --git a/SD-VBS/common/matlab/photonReportTiming.m b/SD-VBS/common/matlab/photonReportTiming.m
new file mode 100644
index 0000000..e991175
--- /dev/null
+++ b/SD-VBS/common/matlab/photonReportTiming.m
@@ -0,0 +1,7 @@
+%! _photonReportTiming_i2_i2i2
+
+function elapsed = photonReportTiming(startCycles, endCycles)
+    elapsed = zeros(1,2);
+    elapsed(1) = endCycles(1) - startCycles(1);
+    elapsed(2) = endCycles(2) - startCycles(2);
+end
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/photonStartTiming.c b/SD-VBS/common/matlab/photonStartTiming.c
new file mode 100644
index 0000000..bbf5231
--- /dev/null
+++ b/SD-VBS/common/matlab/photonStartTiming.c
@@ -0,0 +1,15 @@
+#include"mex.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
+{
+    uint32_t* cycles;
+    plhs[0] = mxCreateNumericMatrix(1, 2, mxUINT32_CLASS, mxREAL);
+    cycles = (uint32_t*)mxGetPr(plhs[0]);
+    __asm__ __volatile__( "rdtsc": "=a" (cycles[0]), "=d" (cycles[1]));    
+ 
+    return;
+}
diff --git a/SD-VBS/common/matlab/photonStartTiming.mexa64 b/SD-VBS/common/matlab/photonStartTiming.mexa64
new file mode 100755
index 0000000..58e0d06
--- /dev/null
+++ b/SD-VBS/common/matlab/photonStartTiming.mexa64
diff --git a/SD-VBS/common/matlab/photonStartTiming.mexglx b/SD-VBS/common/matlab/photonStartTiming.mexglx
new file mode 100755
index 0000000..8b55b98
--- /dev/null
+++ b/SD-VBS/common/matlab/photonStartTiming.mexglx
diff --git a/SD-VBS/common/matlab/product.m b/SD-VBS/common/matlab/product.m
new file mode 100755
index 0000000..65548b4
--- /dev/null
+++ b/SD-VBS/common/matlab/product.m
@@ -0,0 +1,30 @@
+function [matrixOut] = product(matrixIn, direction)
+
+% Initialize matrices
+inputDim = size(matrixIn);
+
+% For input matrix (m,n), we need to multiply horizontally if output matrix
+%is (m). Else, we multiply vertically.
+
+%if direction is 1, we multiply vertically.
+if direction == 1
+    matrixOut = zeros(1, inputDim(2));     %initialize the output matrix
+    for cols = 1:inputDim(2)
+        val = 1;
+        for rows = 1:inputDim(1)
+            val = matrixIn(rows, cols) * val;
+        end
+        matrixOut(cols) = val;
+    end
+
+%  else multiply horizontally
+else
+    matrixOut = zeros(inputDim(1), 1);     %initialize the output matrix
+    for rows = 1:inputDim(1)
+        val = 1;
+        for cols = 1:inputDim(2)
+            val = matrixIn(rows, cols) * val;
+        end
+        matrixOut(rows) = val;
+    end
+end
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/randWrapper.m b/SD-VBS/common/matlab/randWrapper.m
new file mode 100644
index 0000000..dbef132
--- /dev/null
+++ b/SD-VBS/common/matlab/randWrapper.m
@@ -0,0 +1,17 @@
+function out = randWrapper(m,n)
+
+out = zeros(m,n);
+
+seed = 0.9;
+for i=1:m
+    for j=1:n
+        if(i<j)
+            out(i,j) = seed * (i/j);
+        else
+            out(i,j) = seed * (j/i);
+        end
+    end
+end
+
+end
+
diff --git a/SD-VBS/common/matlab/randn.m b/SD-VBS/common/matlab/randn.m
new file mode 100644
index 0000000..e8eb473
--- /dev/null
+++ b/SD-VBS/common/matlab/randn.m
@@ -0,0 +1,22 @@
+%! _randn_f2_i1i1

+

+function retRand = randn(m,n)

+

+retRand = ones(m,n);

+

+for i=1:m

+    for j=1:n

+        w = 10;

+        while(w>=1.0)

+            x1 = 2.0 * rand(1,1) - 1.0;

+            x2 = 2.0 * rand(1,1) - 1.0;

+            w = x1*x1 + x2*x2;

+        end

+        w = sqrt((-2.0*log(w))/w);

+        retRand(i,j) = x1*w;

+    end

+end

+

+

+

+

diff --git a/SD-VBS/common/matlab/randnWrapper.m b/SD-VBS/common/matlab/randnWrapper.m
new file mode 100644
index 0000000..1755957
--- /dev/null
+++ b/SD-VBS/common/matlab/randnWrapper.m
@@ -0,0 +1,17 @@
+function out = randnWrapper(m, n)
+
+out = ones(m,n);
+temp = randWrapper(m,n);
+
+for i=1:m
+    for j=1:n
+        w = temp(i,j);
+        x1 = 1;
+        %w = sqrt((-2.0*log(w))/w);
+        w = ((-2.0*log(w))/w);
+        out(i,j) = x1*w;
+    end
+end
+
+end
+
diff --git a/SD-VBS/common/matlab/readFile.m b/SD-VBS/common/matlab/readFile.m
new file mode 100644
index 0000000..4888103
--- /dev/null
+++ b/SD-VBS/common/matlab/readFile.m
@@ -0,0 +1,22 @@
+function fid = readFile(path)
+
+file = fopen(path, 'r');
+
+full = fscanf(file,'%f');
+elapsed = zeros(1,2);
+
+rows = full(1);
+cols = full(2);
+fid = zeros(rows, cols);
+
+k = 3;
+for i=1:rows
+    for j =1:cols
+        fid(i,j) = full(k);
+        k = k+1;
+    end
+end
+fclose(file);
+
+end
+
diff --git a/SD-VBS/common/matlab/readImage.m b/SD-VBS/common/matlab/readImage.m
new file mode 100644
index 0000000..c70222d
--- /dev/null
+++ b/SD-VBS/common/matlab/readImage.m
@@ -0,0 +1,60 @@
+function srcImage = readImage(pathName)
+
+    %Reading BMP image 
+    input = fopen(pathName,'r');
+    %start of header information
+    signature = fread(input, 2, 'uchar');
+    file_size = fread(input, 1, 'uint32');          
+    reserved1 = fread(input, 1, 'uint16');
+    reserved2 = fread(input, 1, 'uint16');
+    loc_of_bitmap = fread(input, 1, 'uint32');
+
+    size_of_infoheader = fread(input, 1, 'uint32');
+    width = fread(input, 1, 'uint32');
+    height = fread(input, 1, 'uint32');
+    number_of_planes = fread(input, 1, 'uint16');
+    bits_per_pixel = fread(input, 1, 'uint16');
+    compression_method = fread(input, 1, 'uint32');
+    bytes_of_bitmap = fread(input, 1, 'uint32');
+
+    hori_reso = fread(input, 1, 'uint32');
+    vert_reso = fread(input, 1, 'uint32');
+    no_of_colors = fread(input, 1, 'uint32');
+    no_of_imp_colors = fread(input, 1, 'uint32');        
+
+    %end of header information
+
+    srcImage = zeros(height, width);
+
+    % Conditions to check whether the BMP is interleaved and handling few exceptions
+    if (height <= 0 || width <= 0 || signature(1) ~= 'B' || signature(2) ~= 'M'  || ( bits_per_pixel ==16))
+        disp('Error in file format');
+        srcImage = 0;
+    end
+
+    status = fseek(input,loc_of_bitmap,-1);
+
+    nI = 0;
+    nJ = 0;
+    
+    if(bits_per_pixel == 24)
+    for nI=height:-1:1
+        for nJ=1:width
+            tempb = fread(input, 1,'uchar');
+            tempg = fread(input, 1,'uchar');
+            tempr = fread(input, 1,'uchar');
+            srcImage(nI,nJ) = uint8((tempb + 6*tempg + 3*tempr)/10);
+            srcImage(nI,nJ) = uint8(tempg);
+        end
+    end
+    else 
+    for nI=height:-1:1
+        for nJ=1:width
+            tempg = fread(input, 1,'uchar');
+            srcImage(nI,nJ) = tempg;
+        end
+    end
+    end
+    
+    fclose(input);
+end
diff --git a/SD-VBS/common/matlab/read_image8u_bmp.m b/SD-VBS/common/matlab/read_image8u_bmp.m
new file mode 100755
index 0000000..28e5834
--- /dev/null
+++ b/SD-VBS/common/matlab/read_image8u_bmp.m
@@ -0,0 +1,66 @@
+function image = read_image8u_bmp(pathName)

+

+                fid = fopen (pathName,'r'); %FILE *input;

+%               //check for the input FILE pointer

+                if(fid == 0)

+                

+                disp('File pointer error');

+        

+                else

+%                       start of header information

+            BYTES_PER_PIXEL=3;

+            signature = fread(fid,2,'char=>char');

+                        file_size = fread(fid,1,'int=>int');

+            reserved1 = fread(fid,1,'short=>short');

+            reserved2 = fread(fid,1,'short=>short');

+                        loc_of_bitmap = fread(fid, 1, 'int=>double');

+                    size_of_infoheader = fread(fid,1,'int=>int');

+            width = fread(fid,1,'int=>double');

+            height = fread(fid,1,'int=>double');

+            number_of_planes = fread(fid,1,'short=>short');

+            bits_per_pixel = fread(fid,1,'short=>short');

+            compression_method = fread(fid,1,'int=>int');

+            bytes_of_bitmap = fread(fid,1,'int=>int');

+            hori_reso = fread(fid,1,'int=>int');    

+            vert_reso = fread(fid,1,'int=>int');                        

+            no_of_colors = fread(fid,1,'int=>int');                     

+            no_of_imp_colors = fread(fid,1,'int=>int');                         

+                        

+                        nRows = height;

+                        nCols = width;                  

+                        nPitch = nCols;

+            pad = 4 - (nCols*BYTES_PER_PIXEL - 4*floor(nCols*BYTES_PER_PIXEL/4)); %bitmap multiple-of-4 requirement

+                        

+                        pixSize=nRows *nCols * 3;

+                        

+                        % Conditions to check whether the BMP is interleaved and handling few exceptions

+                        if  (nRows <= 0 || nCols <= 0 || signature(1) ~= 'B' || signature(2) ~= 'M'  || bits_per_pixel ~=24)

+                                disp ('Error');

+                return;

+            else

+        

+%read image

+%             fseek(fid,loc_of_bitmap,'bof');

+

+                        for nI = nRows:-1:1

+                                for nJ=1:nCols

+                                

+                                  image(nI,nJ,3) = fread(fid,1,'char=>uint8'); %B

+                                  image(nI,nJ,2) = fread(fid,1,'char=>uint8'); %G

+                                  image(nI,nJ,1) = fread(fid,1,'char=>uint8'); %R

+                  

+                end

+                if pad~=4

+                    for i=1:pad

+                    fread(fid,1);               

+                    end

+                end

+            end

+            imshow(image);

+            imwrite(image,'image.bmp');

+                        fclose(fid);    

+

+            end

+            

+        end

+        
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/reshapeMatrix.m b/SD-VBS/common/matlab/reshapeMatrix.m
new file mode 100755
index 0000000..f355828
--- /dev/null
+++ b/SD-VBS/common/matlab/reshapeMatrix.m
@@ -0,0 +1,13 @@
+function [outputMatrix] = reshapeMatrix(inputMatrix, outputSize)
+
+inputDim = size(inputMatrix);
+outputMatrix = zeros(outputSize(1), outputSize(2));
+
+k = 1;
+
+for i = 1:outputSize(2)
+    for j = 1:outputSize(1)
+        outputMatrix(j, i) = inputMatrix(k);
+        k = k+1;
+    end
+end
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/selfCheck.m b/SD-VBS/common/matlab/selfCheck.m
new file mode 100644
index 0000000..55c8de8
--- /dev/null
+++ b/SD-VBS/common/matlab/selfCheck.m
@@ -0,0 +1,27 @@
+function ret = selfCheck(in1, path, tol)
+
+r1 = size(in1, 1);
+c1 = size(in1, 2);
+
+ret = 1;
+
+file = [path, '/expected.m'];
+fd = fopen(file, 'r');
+
+[in2, count] = fscanf(fd, '%d');
+
+if(count ~= (r1*c1) )
+    fprintf(1, 'Dimensions mismatch: Expected %d\t Observed %d\n', count, (r1*c1));
+    ret = -1;
+else
+    ret = 1;
+    for i=1:(r1*c1)
+        if( (abs(in1(i)) - abs(in2(i)) > tol) || (abs(in2(i)) - abs(in1(i))) > tol)
+            fprintf(1, 'Checking Error: Index %d\tExpected %d\tObserved %d\n', i, in2(i), in1(i));
+            ret = -1;
+            break;
+        end
+    end
+
+end
+
diff --git a/SD-VBS/common/matlab/testProductFunction.m b/SD-VBS/common/matlab/testProductFunction.m
new file mode 100755
index 0000000..682a005
--- /dev/null
+++ b/SD-VBS/common/matlab/testProductFunction.m
@@ -0,0 +1,21 @@
+function testProductFunction
+
+matrixSize = zeros(1, 2);
+matrixSize(1) = 4;
+matrixSize(2) = 5;
+inputMatrix = zeros(matrixSize(1), matrixSize(2));
+
+for i=1:matrixSize(1)
+    for j=1:matrixSize(2)
+        inputMatrix(i,j) = i+j;
+    end
+end
+
+outputMatrix = product(inputMatrix, 1);
+actualOut = prod(inputMatrix, 1);
+
+if outputMatrix == actualOut
+    disp('SUCCESS');
+else
+    disp('ERROR');
+end
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/testReshapeFunction.m b/SD-VBS/common/matlab/testReshapeFunction.m
new file mode 100755
index 0000000..1dac432
--- /dev/null
+++ b/SD-VBS/common/matlab/testReshapeFunction.m
@@ -0,0 +1,21 @@
+function testReshapeFunction
+
+matrixSize = zeros(1, 2);
+matrixSize(1) = 4;
+matrixSize(2) = 5;
+inputMatrix = zeros(matrixSize(1), matrixSize(2));
+
+for i=1:matrixSize(1)
+    for j=1:matrixSize(2)
+        inputMatrix(i,j) = i+j;
+    end
+end
+
+outputMatrix = reshapeMatrix(inputMatrix, matrixSize);
+actualOut = reshape(inputMatrix, matrixSize);
+
+if outputMatrix == actualOut
+    disp('SUCCESS');
+else
+    disp('ERROR');
+end
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/timingFuncs/cycle.h b/SD-VBS/common/matlab/timingFuncs/cycle.h
new file mode 100644
index 0000000..2652a04
--- /dev/null
+++ b/SD-VBS/common/matlab/timingFuncs/cycle.h
@@ -0,0 +1,514 @@
+/*
+ * Copyright (c) 2003, 2007-8 Matteo Frigo
+ * Copyright (c) 2003, 2007-8 Massachusetts Institute of Technology
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ *
+ */
+
+
+/* machine-dependent cycle counters code. Needs to be inlined. */
+
+/***************************************************************************/
+/* To use the cycle counters in your code, simply #include "cycle.h" (this
+   file), and then use the functions/macros:
+
+                 ticks getticks(void);
+
+   ticks is an opaque typedef defined below, representing the current time.
+   You extract the elapsed time between two calls to gettick() via:
+
+                 double elapsed(ticks t1, ticks t0);
+
+   which returns a double-precision variable in arbitrary units.  You
+   are not expected to convert this into human units like seconds; it
+   is intended only for *comparisons* of time intervals.
+
+   (In order to use some of the OS-dependent timer routines like
+   Solaris' gethrtime, you need to paste the autoconf snippet below
+   into your configure.ac file and #include "config.h" before cycle.h,
+   or define the relevant macros manually if you are not using autoconf.)
+*/
+
+/***************************************************************************/
+/* This file uses macros like HAVE_GETHRTIME that are assumed to be
+   defined according to whether the corresponding function/type/header
+   is available on your system.  The necessary macros are most
+   conveniently defined if you are using GNU autoconf, via the tests:
+   
+   dnl ---------------------------------------------------------------------
+
+   AC_C_INLINE
+   AC_HEADER_TIME
+   AC_CHECK_HEADERS([sys/time.h c_asm.h intrinsics.h mach/mach_time.h])
+
+   AC_CHECK_TYPE([hrtime_t],[AC_DEFINE(HAVE_HRTIME_T, 1, [Define to 1 if hrtime_t is defined in <sys/time.h>])],,[#if HAVE_SYS_TIME_H
+#include <sys/time.h>
+#endif])
+
+   AC_CHECK_FUNCS([gethrtime read_real_time time_base_to_time clock_gettime mach_absolute_time])
+
+   dnl Cray UNICOS _rtc() (real-time clock) intrinsic
+   AC_MSG_CHECKING([for _rtc intrinsic])
+   rtc_ok=yes
+   AC_TRY_LINK([#ifdef HAVE_INTRINSICS_H
+#include <intrinsics.h>
+#endif], [_rtc()], [AC_DEFINE(HAVE__RTC,1,[Define if you have the UNICOS _rtc() intrinsic.])], [rtc_ok=no])
+   AC_MSG_RESULT($rtc_ok)
+
+   dnl ---------------------------------------------------------------------
+*/
+
+/***************************************************************************/
+
+#if TIME_WITH_SYS_TIME
+# include <sys/time.h>
+# include <time.h>
+#else
+# if HAVE_SYS_TIME_H
+#  include <sys/time.h>
+# else
+#  include <time.h>
+# endif
+#endif
+
+#define INLINE_ELAPSED(INL) static INL double elapsed(ticks t1, ticks t0) \
+{                                                                         \
+     return (double)t1 - (double)t0;                                      \
+}
+
+/*----------------------------------------------------------------*/
+/* Solaris */
+#if defined(HAVE_GETHRTIME) && defined(HAVE_HRTIME_T) && !defined(HAVE_TICK_COUNTER)
+typedef hrtime_t ticks;
+
+#define getticks gethrtime
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* AIX v. 4+ routines to read the real-time clock or time-base register */
+#if defined(HAVE_READ_REAL_TIME) && defined(HAVE_TIME_BASE_TO_TIME) && !defined(HAVE_TICK_COUNTER)
+typedef timebasestruct_t ticks;
+
+static __inline ticks getticks(void)
+{
+     ticks t;
+     read_real_time(&t, TIMEBASE_SZ);
+     return t;
+}
+
+static __inline double elapsed(ticks t1, ticks t0) /* time in nanoseconds */
+{
+     time_base_to_time(&t1, TIMEBASE_SZ);
+     time_base_to_time(&t0, TIMEBASE_SZ);
+     return (((double)t1.tb_high - (double)t0.tb_high) * 1.0e9 + 
+             ((double)t1.tb_low - (double)t0.tb_low));
+}
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * PowerPC ``cycle'' counter using the time base register.
+ */
+#if ((((defined(__GNUC__) && (defined(__powerpc__) || defined(__ppc__))) || (defined(__MWERKS__) && defined(macintosh)))) || (defined(__IBM_GCC_ASM) && (defined(__powerpc__) || defined(__ppc__))))  && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     unsigned int tbl, tbu0, tbu1;
+
+     do {
+          __asm__ __volatile__ ("mftbu %0" : "=r"(tbu0));
+          __asm__ __volatile__ ("mftb %0" : "=r"(tbl));
+          __asm__ __volatile__ ("mftbu %0" : "=r"(tbu1));
+     } while (tbu0 != tbu1);
+
+     return (((unsigned long long)tbu0) << 32) | tbl;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* MacOS/Mach (Darwin) time-base register interface (unlike UpTime,
+   from Carbon, requires no additional libraries to be linked). */
+#if defined(HAVE_MACH_ABSOLUTE_TIME) && defined(HAVE_MACH_MACH_TIME_H) && !defined(HAVE_TICK_COUNTER)
+#include <mach/mach_time.h>
+typedef uint64_t ticks;
+#define getticks mach_absolute_time
+INLINE_ELAPSED(__inline__)
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * Pentium cycle counter 
+ */
+#if (defined(__GNUC__) || defined(__ICC)) && defined(__i386__)  && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+
+     __asm__ __volatile__("rdtsc": "=A" (ret));
+     /* no input, nothing else clobbered */
+     return ret;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#define TIME_MIN 5000.0   /* unreliable pentium IV cycle counter */
+#endif
+
+/* Visual C++ -- thanks to Morten Nissov for his help with this */
+#if _MSC_VER >= 1200 && _M_IX86 >= 500 && !defined(HAVE_TICK_COUNTER)
+#include <windows.h>
+typedef LARGE_INTEGER ticks;
+#define RDTSC __asm __emit 0fh __asm __emit 031h /* hack for VC++ 5.0 */
+
+static __inline ticks getticks(void)
+{
+     ticks retval;
+
+     __asm {
+          RDTSC
+          mov retval.HighPart, edx
+          mov retval.LowPart, eax
+     }
+     return retval;
+}
+
+static __inline double elapsed(ticks t1, ticks t0)
+{  
+     return (double)t1.QuadPart - (double)t0.QuadPart;
+}  
+
+#define HAVE_TICK_COUNTER
+#define TIME_MIN 5000.0   /* unreliable pentium IV cycle counter */
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * X86-64 cycle counter
+ */
+#if (defined(__GNUC__) || defined(__ICC) || defined(__SUNPRO_C)) && defined(__x86_64__)  && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     unsigned a, d; 
+     asm volatile("rdtsc" : "=a" (a), "=d" (d)); 
+     return ((ticks)a) | (((ticks)d) << 32); 
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* PGI compiler, courtesy Cristiano Calonaci, Andrea Tarsi, & Roberto Gori.
+   NOTE: this code will fail to link unless you use the -Masmkeyword compiler
+   option (grrr). */
+#if defined(__PGI) && defined(__x86_64__) && !defined(HAVE_TICK_COUNTER) 
+typedef unsigned long long ticks;
+static ticks getticks(void)
+{
+    asm(" rdtsc; shl    $0x20,%rdx; mov    %eax,%eax; or     %rdx,%rax;    ");
+}
+INLINE_ELAPSED(__inline__)
+#define HAVE_TICK_COUNTER
+#endif
+
+/* Visual C++, courtesy of Dirk Michaelis */
+#if _MSC_VER >= 1400 && (defined(_M_AMD64) || defined(_M_X64)) && !defined(HAVE_TICK_COUNTER)
+
+#include <intrin.h>
+#pragma intrinsic(__rdtsc)
+typedef unsigned __int64 ticks;
+#define getticks __rdtsc
+INLINE_ELAPSED(__inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * IA64 cycle counter
+ */
+
+/* intel's icc/ecc compiler */
+#if (defined(__EDG_VERSION) || defined(__ECC)) && defined(__ia64__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+#include <ia64intrin.h>
+
+static __inline__ ticks getticks(void)
+{
+     return __getReg(_IA64_REG_AR_ITC);
+}
+ 
+INLINE_ELAPSED(__inline__)
+ 
+#define HAVE_TICK_COUNTER
+#endif
+
+/* gcc */
+#if defined(__GNUC__) && defined(__ia64__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+
+     __asm__ __volatile__ ("mov %0=ar.itc" : "=r"(ret));
+     return ret;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* HP/UX IA64 compiler, courtesy Teresa L. Johnson: */
+#if defined(__hpux) && defined(__ia64) && !defined(HAVE_TICK_COUNTER)
+#include <machine/sys/inline.h>
+typedef unsigned long ticks;
+
+static inline ticks getticks(void)
+{
+     ticks ret;
+
+     ret = _Asm_mov_from_ar (_AREG_ITC);
+     return ret;
+}
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/* Microsoft Visual C++ */
+#if defined(_MSC_VER) && defined(_M_IA64) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned __int64 ticks;
+
+#  ifdef __cplusplus
+extern "C"
+#  endif
+ticks __getReg(int whichReg);
+#pragma intrinsic(__getReg)
+
+static __inline ticks getticks(void)
+{
+     volatile ticks temp;
+     temp = __getReg(3116);
+     return temp;
+}
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/*
+ * PA-RISC cycle counter 
+ */
+#if defined(__hppa__) || defined(__hppa) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+
+#  ifdef __GNUC__
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+
+     __asm__ __volatile__("mfctl 16, %0": "=r" (ret));
+     /* no input, nothing else clobbered */
+     return ret;
+}
+#  else
+#  include <machine/inline.h>
+static inline unsigned long getticks(void)
+{
+     register ticks ret;
+     _MFCTL(16, ret);
+     return ret;
+}
+#  endif
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* S390, courtesy of James Treacy */
+#if defined(__GNUC__) && defined(__s390__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks cycles;
+     __asm__("stck 0(%0)" : : "a" (&(cycles)) : "memory", "cc");
+     return cycles;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+/*----------------------------------------------------------------*/
+#if defined(__GNUC__) && defined(__alpha__) && !defined(HAVE_TICK_COUNTER)
+/*
+ * The 32-bit cycle counter on alpha overflows pretty quickly, 
+ * unfortunately.  A 1GHz machine overflows in 4 seconds.
+ */
+typedef unsigned int ticks;
+
+static __inline__ ticks getticks(void)
+{
+     unsigned long cc;
+     __asm__ __volatile__ ("rpcc %0" : "=r"(cc));
+     return (cc & 0xFFFFFFFF);
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+#if defined(__GNUC__) && defined(__sparc_v9__) && !defined(HAVE_TICK_COUNTER)
+typedef unsigned long ticks;
+
+static __inline__ ticks getticks(void)
+{
+     ticks ret;
+     __asm__ __volatile__("rd %%tick, %0" : "=r" (ret));
+     return ret;
+}
+
+INLINE_ELAPSED(__inline__)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+#if (defined(__DECC) || defined(__DECCXX)) && defined(__alpha) && defined(HAVE_C_ASM_H) && !defined(HAVE_TICK_COUNTER)
+#  include <c_asm.h>
+typedef unsigned int ticks;
+
+static __inline ticks getticks(void)
+{
+     unsigned long cc;
+     cc = asm("rpcc %v0");
+     return (cc & 0xFFFFFFFF);
+}
+
+INLINE_ELAPSED(__inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+/*----------------------------------------------------------------*/
+/* SGI/Irix */
+#if defined(HAVE_CLOCK_GETTIME) && defined(CLOCK_SGI_CYCLE) && !defined(HAVE_TICK_COUNTER)
+typedef struct timespec ticks;
+
+static inline ticks getticks(void)
+{
+     struct timespec t;
+     clock_gettime(CLOCK_SGI_CYCLE, &t);
+     return t;
+}
+
+static inline double elapsed(ticks t1, ticks t0)
+{
+     return ((double)t1.tv_sec - (double)t0.tv_sec) * 1.0E9 +
+          ((double)t1.tv_nsec - (double)t0.tv_nsec);
+}
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* Cray UNICOS _rtc() intrinsic function */
+#if defined(HAVE__RTC) && !defined(HAVE_TICK_COUNTER)
+#ifdef HAVE_INTRINSICS_H
+#  include <intrinsics.h>
+#endif
+
+typedef long long ticks;
+
+#define getticks _rtc
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+
+/*----------------------------------------------------------------*/
+/* MIPS ZBus */
+#if HAVE_MIPS_ZBUS_TIMER
+#if defined(__mips__) && !defined(HAVE_TICK_COUNTER)
+#include <sys/mman.h>
+#include <unistd.h>
+#include <fcntl.h>
+
+typedef uint64_t ticks;
+
+static inline ticks getticks(void)
+{
+  static uint64_t* addr = 0;
+
+  if (addr == 0)
+  {
+    uint32_t rq_addr = 0x10030000;
+    int fd;
+    int pgsize;
+
+    pgsize = getpagesize();
+    fd = open ("/dev/mem", O_RDONLY | O_SYNC, 0);
+    if (fd < 0) {
+      perror("open");
+      return NULL;
+    }
+    addr = mmap(0, pgsize, PROT_READ, MAP_SHARED, fd, rq_addr);
+    close(fd);
+    if (addr == (uint64_t *)-1) {
+      perror("mmap");
+      return NULL;
+    }
+  }
+
+  return *addr;
+}
+
+INLINE_ELAPSED(inline)
+
+#define HAVE_TICK_COUNTER
+#endif
+#endif /* HAVE_MIPS_ZBUS_TIMER */
+
diff --git a/SD-VBS/common/matlab/timingFuncs/photonEndTiming.c b/SD-VBS/common/matlab/timingFuncs/photonEndTiming.c
new file mode 100644
index 0000000..bbf5231
--- /dev/null
+++ b/SD-VBS/common/matlab/timingFuncs/photonEndTiming.c
@@ -0,0 +1,15 @@
+#include"mex.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
+{
+    uint32_t* cycles;
+    plhs[0] = mxCreateNumericMatrix(1, 2, mxUINT32_CLASS, mxREAL);
+    cycles = (uint32_t*)mxGetPr(plhs[0]);
+    __asm__ __volatile__( "rdtsc": "=a" (cycles[0]), "=d" (cycles[1]));    
+ 
+    return;
+}
diff --git a/SD-VBS/common/matlab/timingFuncs/photonEndTiming.mexglx b/SD-VBS/common/matlab/timingFuncs/photonEndTiming.mexglx
new file mode 100755
index 0000000..fa9d4b6
--- /dev/null
+++ b/SD-VBS/common/matlab/timingFuncs/photonEndTiming.mexglx
diff --git a/SD-VBS/common/matlab/timingFuncs/photonPrintTiming.m b/SD-VBS/common/matlab/timingFuncs/photonPrintTiming.m
new file mode 100644
index 0000000..b9c7b22
--- /dev/null
+++ b/SD-VBS/common/matlab/timingFuncs/photonPrintTiming.m
@@ -0,0 +1,5 @@
+%! _photonPrintTiming_NA_i2
+
+function photonPrintTiming(elapsed)
+    disp(elapsed);
+end
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/timingFuncs/photonReportTiming.m b/SD-VBS/common/matlab/timingFuncs/photonReportTiming.m
new file mode 100644
index 0000000..e991175
--- /dev/null
+++ b/SD-VBS/common/matlab/timingFuncs/photonReportTiming.m
@@ -0,0 +1,7 @@
+%! _photonReportTiming_i2_i2i2
+
+function elapsed = photonReportTiming(startCycles, endCycles)
+    elapsed = zeros(1,2);
+    elapsed(1) = endCycles(1) - startCycles(1);
+    elapsed(2) = endCycles(2) - startCycles(2);
+end
\ No newline at end of file
diff --git a/SD-VBS/common/matlab/timingFuncs/photonStartTiming.c b/SD-VBS/common/matlab/timingFuncs/photonStartTiming.c
new file mode 100644
index 0000000..bbf5231
--- /dev/null
+++ b/SD-VBS/common/matlab/timingFuncs/photonStartTiming.c
@@ -0,0 +1,15 @@
+#include"mex.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
+{
+    uint32_t* cycles;
+    plhs[0] = mxCreateNumericMatrix(1, 2, mxUINT32_CLASS, mxREAL);
+    cycles = (uint32_t*)mxGetPr(plhs[0]);
+    __asm__ __volatile__( "rdtsc": "=a" (cycles[0]), "=d" (cycles[1]));    
+ 
+    return;
+}
diff --git a/SD-VBS/common/matlab/timingFuncs/photonStartTiming.mexglx b/SD-VBS/common/matlab/timingFuncs/photonStartTiming.mexglx
new file mode 100755
index 0000000..8b55b98
--- /dev/null
+++ b/SD-VBS/common/matlab/timingFuncs/photonStartTiming.mexglx
diff --git a/SD-VBS/common/matlab/writeMatrix.m b/SD-VBS/common/matlab/writeMatrix.m
new file mode 100644
index 0000000..ec2bbae
--- /dev/null
+++ b/SD-VBS/common/matlab/writeMatrix.m
@@ -0,0 +1,19 @@
+function writeMatrix(input, inpath)
+
+file = [inpath '/expected.m'];
+disp(file);
+fd = fopen(file, 'w');
+
+[rows, cols] = size(input);
+
+for j=1:rows
+    for i=1:cols
+        fprintf(fd, '%d\t', input(j,i)-1);
+    end
+    fprintf(fd, '\n');
+end
+fclose(fd);
+
+end
+
+
diff --git a/SD-VBS/common/matlab/write_image8u_bmp.m b/SD-VBS/common/matlab/write_image8u_bmp.m
new file mode 100755
index 0000000..71b48ef
--- /dev/null
+++ b/SD-VBS/common/matlab/write_image8u_bmp.m
@@ -0,0 +1,88 @@
+function image = write_image8u_bmp(inpImage, pathName) % this is matlab's sequence, removed width, height since they can be inferred from image

+

+BYTES_PER_PIXEL = 3;

+

+fid = fopen (pathName,'w'); %FILE *input;

+%               //check for the input FILE pointer

+if(fid == 0)            

+    disp('File pointer error');

+else

+    height = size(inpImage,1);

+    width=size(inpImage,2);

+    signature(1) = 'B';

+    signature(2) = 'M';

+    

+    reserved1 = 0;

+    reserved2 = 0;

+    file_size = (height*width*BYTES_PER_PIXEL)+54;

+    offset=54;

+        size_of_infoheader=40;

+        number_of_planes=1;

+        bits_per_pixel= 8*BYTES_PER_PIXEL ;

+        compression_method=0;

+        hori_reso=2835;

+        vert_reso=2835;

+        no_of_colors=0;

+        no_of_imp_colors=0;

+        

+%       start of header information

+    fwrite(fid,signature(1),'char');  

+    fwrite(fid,signature(2),'char');  

+    fwrite(fid,file_size,'int');

+    fwrite(fid,reserved1,'short');

+    fwrite(fid,reserved2,'short');

+    fwrite(fid,offset,'int');

+    fwrite(fid,size_of_infoheader,'int');

+    fwrite(fid,width,'int');

+    fwrite(fid,height,'int');

+    fwrite(fid,number_of_planes,'short');

+    fwrite(fid,bits_per_pixel,'short');

+    fwrite(fid,compression_method,'int');

+    

+    bytes_of_bitmap=width*height*BYTES_PER_PIXEL;

+    

+    fwrite(fid,bytes_of_bitmap,'int');

+    fwrite(fid,hori_reso,'int');    

+    fwrite(fid,vert_reso,'int');                        

+    fwrite(fid,no_of_colors,'int');                     

+    fwrite(fid,no_of_imp_colors,'int');                         

+

+        % Conditions to check whether the BMP is gray scaled and handling few exceptions

+

+    if (width <= 0 || height <= 0 || signature(1) ~= 'B' || signature(2) ~= 'M')

+        return;

+    end

+                        

+        % total size of pixels

+        % pixSize=srcImage->height * srcImage->width * NO_OF_BYTE_PER_PIXEL; //*3 is done as it is a 3 channel image

+    %reverse the image back

+

+    for nI=1:height

+        count = 0;

+        for nJ=1:width

+          temp = inpImage((height - nI+1),nJ,3);

+          fwrite(fid,temp,'char');

+          temp = inpImage((height - nI+1),nJ,2);

+          fwrite(fid,temp,'char');

+          temp = inpImage((height - nI+1),nJ,1);

+          fwrite(fid,temp,'char');

+          

+        end

+       pad = 4 - (width*BYTES_PER_PIXEL - (4*floor(width*BYTES_PER_PIXEL/4))); %bitmap multiple-of-4 requirement

+       if pad ~= 0

+           temp=0;

+           for cnt = 1:pad

+              fwrite(fid,temp,'char');

+%               fwrite(fid,temp,'char');

+%               fwrite(fid,temp,'char');              

+           end

+       end

+    end

+    

+    outImage = imread(pathName);

+    imshow(outImage);

+    fclose(fid);        

+    

+end

+            

+        
\ No newline at end of file