diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp
index b130d47..999775f 100644
--- a/test/Reconstruction_Test.cpp
+++ b/test/Reconstruction_Test.cpp
@@ -2,6 +2,7 @@
 #include <gtest/gtest.h>
 #include <limits>
 
+#include "CudaMatrix.h"
 #include "Function1D.h"
 #include "Function3D.h"
 #include "MatlabReader.h"
@@ -14,6 +15,9 @@
 #include "transmissionReconstruction/reconstruction/buildMatrix/DGradient.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/FMM.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
+
+#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cuh"
+
 #include "transmissionReconstruction/reconstruction/reconstruction.h"
 #include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h"
 
@@ -23,6 +27,9 @@
 
 #include "reflectionReconstruction/preprocessData/determineOptimalPulse.h"
 
+#include <cusparse.h>
+#include <mkl_spblas.h>
+
 
 
 
@@ -52,6 +59,14 @@ protected:
     }
 };
 
+TEST_F(Reconstruction_Test, GPUMajor) {
+	int cudaDevice;
+    cudaGetDevice(&cudaDevice);
+    cudaDeviceProp prop;
+	cudaGetDeviceProperties(&prop, cudaDevice);
+    EXPECT_GT(prop.major,2)<<"GPU major:"<<prop.major;
+}
+
 TEST_F(Reconstruction_Test, determineOptimalPulse) {
     Recon::reflectParams::imageResolution = 8.925316499483377e-04;
     Recon::reflectParams::optPulseFactor = 48;
@@ -171,30 +186,85 @@ TEST_F(Reconstruction_Test, checkAndScale) {
 }
 
 TEST_F(Reconstruction_Test, buildMatrix) {
-    MatlabReader m("/home/sun/testData/buildMatrix.mat");
-    auto j1 = m.read("j1");
-    auto i1 = m.read("i1");
-    auto s1 = m.read("s1");
-    auto dims = m.read("dims");
+    Recon::transParams::bresenham = 1;
+    MatlabReader m("/home/krad/TestData/buildmatrixFloat.mat");
+    auto j1 = m.read("jCopy");
+    auto i1 = m.read("iCopy");
+    auto s1 = m.read("sCopy");
+    auto dims = Aurora::Matrix::fromRawData(new float[3]{128, 128, 73}, 1, 3);;
     auto receiverList  = m.read("receiverList");
     auto potentialMap = m.read("potentialMap");
     auto senderList  = m.read("senderList"); 
     auto res  = m.read("res");
     auto result = Recon::buildMatrix(senderList, receiverList, res, dims, false, potentialMap);
+    size_t count = 0;
+
     for(size_t i=0; i<i1.getDataSize(); ++i)
     {
-        EXPECT_FLOAT_AE(i1[i] - 1, result.M.getColVector()[i]);
+        // EXPECT_FLOAT_AE(i1[i] - 1, result.M.getRowVector()[i])<<"ROW error,"<<" index:"<<i;
+        if ((int)(i1[i] - 1)!= (int)result.M.getRowVector()[i]){
+            count++;
+        }
     }
-
+    EXPECT_EQ(count, 0)<<"ROW error percent,"<<((double)count/(double)i1.getDataSize())*100.0;
+    count=0;
     for(size_t i=0; i<j1.getDataSize(); ++i)
     {
-        EXPECT_FLOAT_AE(j1[i]- 1, result.M.getRowVector()[i]);
+        EXPECT_FLOAT_AE(j1[i]- 1, result.M.getColVector()[i])<<"COL error,"<<" index:"<<i;
+        if ((int)(j1[i] - 1)!= (int)result.M.getColVector()[i]){
+            count++;
+        }
     }
-
+    EXPECT_EQ(count, 0)<<"COL error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
+    count=0;
     for(size_t i=0; i<s1.getDataSize(); ++i)
     {
-        EXPECT_FLOAT_AE(s1[i], result.M.getValVector()[i]);
+        EXPECT_FLOAT_AE(s1[i], result.M.getValVector()[i])<<"Value error,"<<" index:"<<i;
+        if ((int)(s1[i] - 1)!= (int)result.M.getValVector()[i]){
+            count++;
+        }
     }
+    EXPECT_EQ(count, 0)<<"Value error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
+
+}
+
+TEST_F(Reconstruction_Test, buildMatrixCuda) {
+    Recon::transParams::bresenham = 1;
+
+    MatlabReader m("/home/krad/TestData/buildmatrixFloat.mat");
+    auto j1 = m.read("jCopy");
+    auto i1 = m.read("iCopy");
+    auto s1 = m.read("sCopy");
+    auto dims = Aurora::Matrix::fromRawData(new float[3]{128, 128, 73}, 1, 3).toDeviceMatrix();
+    auto receiverList  = m.read("receiverList").toDeviceMatrix();
+    auto potentialMap = m.read("potentialMap").toDeviceMatrix();
+    auto senderList  = m.read("senderList").toDeviceMatrix(); 
+    auto res  = m.read("res").toDeviceMatrix();
+    auto result = Recon::buildMatrix(senderList, receiverList, res, dims, false, potentialMap);
+    size_t count = 0;
+
+    for(size_t i=0; i<i1.getDataSize(); ++i)
+    {
+        EXPECT_FLOAT_AE(i1[i] - 1, result.M.getRowVector()[i])<<"ROW error,"<<" index:"<<i;
+        if ((int)(i1[i] - 1)!= (int)result.M.getRowVector()[i]){
+            count++;
+        }
+    }
+    EXPECT_EQ(count, 0)<<"ROW error percent,"<<((double)count/(double)i1.getDataSize())*100.0;
+    count=0;
+    for(size_t i=0; i<j1.getDataSize(); ++i)
+    {
+        EXPECT_FLOAT_AE(j1[i]- 1, result.M.getColVector()[i])<<"COL error,"<<" index:"<<i;
+        if ((int)(j1[i] - 1)!= (int)result.M.getColVector()[i]){
+            count++;
+        }
+    }
+    EXPECT_EQ(count, 0)<<"COL error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
+    for(size_t i=0; i<s1.getDataSize(); ++i)
+    {
+        EXPECT_FLOAT_EQ(s1[i], result.M.getValVector()[i])<<"Value error,"<<" index:"<<i;
+    }
+    EXPECT_EQ(count, 0)<<"Value error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
 
 }
 
@@ -371,17 +441,64 @@ TEST_F(Reconstruction_Test,callTval3){
     opt.mu = 100;
     opt.beta = 1;
     opt.beta0 = 1;
-    MatlabReader m("/home/sun/testData/buildMatrixTest.mat");
-    auto j1 = m.read("j1");
-    auto i1 = m.read("i1");
-    auto s1 = m.read("s1");
-    Aurora::Sparse M(i1-1,j1-1,s1,734989,1196032);
-    MatlabReader m2("/home/sun/testData/tval3gpu3d.mat");
-    auto b = m2.read("b");
-    auto dims = Aurora::Matrix::fromRawData(new float[3]{128,128,73}, 1, 3);
-    auto result = Recon::callTval3(M, b, dims, 0,opt);
-    auto outSOS = Recon::slownessToSOS(result, 1.498206569328594e+03) ;
-    MatlabWriter w2("/home/krad/transmissionSOS111.mat");
-    w2.setMatrix(outSOS, "SOS");
-    w2.write();
+    MatlabReader m("/home/krad/TestData/ijs.mat");
+    auto j = m.read("j");   
+    auto i = m.read("i");
+    auto s = m.read("s");
+    Aurora::Sparse M(i-1,j-1,s,1900092,1196032);
+    MatlabReader m2("/home/sun/TVALGPU.mat");
+    auto j1 = m2.read("colPtr");
+    auto i1 = m2.read("rowPtr");
+    auto s1 = m2.read("data");
+
+    
+    int * yIdxs = new int[M.getColVector().getDataSize()];
+    std::copy(M.getColVector().getData(),M.getColVector().getData()+M.getColVector().getDataSize(),yIdxs);
+    int * xIdxs = new int[M.getRowVector().getDataSize()];
+    std::copy(M.getRowVector().getData(),M.getRowVector().getData()+M.getRowVector().getDataSize(),xIdxs);
+    Aurora::Matrix values = M.getValVector();
+    size_t rows = M.getM(), cols = M.getN();        
+    int nz = std::max(M.getColVector().getDataSize(),std::max(M.getRowVector().getDataSize(),M.getValVector().getDataSize()));        
+     sparse_matrix_t A;
+    sparse_matrix_t csrA;
+    sparse_matrix_t csrA_t;
+    sparse_status_t status;
+    status = mkl_sparse_s_create_coo(&A, sparse_index_base_t::SPARSE_INDEX_BASE_ZERO,  rows,cols,  nz, xIdxs,yIdxs,values.getData());     
+  
+    status = mkl_sparse_convert_csr(A, sparse_operation_t::SPARSE_OPERATION_TRANSPOSE, &csrA); 
+    // status = mkl_sparse_order(csrA);  
+
+    int n_rows = 0;
+    int n_cols = 0;
+    int *col_start,*col_end,*row_idx;
+    float * csrValues;
+    sparse_index_base_t index;
+
+
+    status = mkl_sparse_s_export_csr(csrA, &index, &n_cols, &n_rows, &col_start, &col_end, &row_idx, &csrValues);      
+
+        
+    int *col_idx = new int[n_cols+1];
+    std::copy(col_start,col_start+ n_cols,col_idx);
+    col_idx[n_cols] = nz;
+    
+    printf("j1 count:%zu,i1 count:%zu,s1 count:%zu\r\n",j1.getDataSize(),i1.getDataSize(),s1.getDataSize());
+    printf("n_rows:%d,n_cols:%d\r\n", n_rows, n_cols);
+    for (size_t k = 0; k < j1.getDataSize(); k++)
+    {
+        ASSERT_FLOAT_EQ(j1.getData()[k],col_idx[k])<<"col error, index:"<<k;
+    }
+
+    for (size_t k = 0; k < i1.getDataSize(); k++)
+    {
+        ASSERT_FLOAT_EQ(i1.getData()[k],row_idx[k])<<"row error, index:"<<k;
+    }
+    for (size_t k = 0; k < s1.getDataSize(); k++)
+    {
+        ASSERT_FLOAT_EQ(s1.getData()[k],csrValues[k])<<"value error, index:"<<k;
+    }
+    mkl_sparse_destroy(A);
+    mkl_sparse_destroy(csrA);
+    delete [] xIdxs;
+    delete [] yIdxs;
 }