From d8fe36f48cfb1c743ca0864f2b386d71fdb55d68 Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Fri, 2 Jun 2023 09:35:37 +0800
Subject: [PATCH 1/2] Make TVAL GPU executeable

---
 .../reconstruction/buildMatrix/buildMatrix.cpp           | 2 +-
 .../reconstruction/reconstruction.cpp                    | 8 ++++----
 .../reconstruction/solvingEquationSystem/TVAL/TVAL.cpp   | 9 +++++----
 .../reconstruction/solvingEquationSystem/solve.cpp       | 1 +
 test/Reconstruction_Test.cpp                             | 9 +++++++--
 5 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
index 3a72f51..d186d4b 100644
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
@@ -309,7 +309,7 @@ namespace Recon
                     result.hitmap[linearIndices[i]]+=1;
                 }
             }
-            printf("Progress: %f (%zu of %zu)\r\n",(double)rayCount*100/(double)nTotalRays,rayCount,nTotalRays);
+            // printf("Progress: %f (%zu of %zu)\r\n",(double)rayCount*100/(double)nTotalRays,rayCount,nTotalRays);
             cnt = cnt + pathLenDisc;
 
             if (cnt < safesize)
diff --git a/src/transmissionReconstruction/reconstruction/reconstruction.cpp b/src/transmissionReconstruction/reconstruction/reconstruction.cpp
index 3fbdb05..66da0f0 100644
--- a/src/transmissionReconstruction/reconstruction/reconstruction.cpp
+++ b/src/transmissionReconstruction/reconstruction/reconstruction.cpp
@@ -129,10 +129,10 @@ namespace Recon {
         for (size_t i = 0; i < transParams::gpuSelectionList.getDataSize(); i++)
         {
             std::string msg;
-            // if (!resetGPUDevice((int)transParams::gpuSelectionList[i],msg))
-            // {
-            //     std::cerr<<msg<<std::endl;
-            // }
+            if (!resetGPUDevice((int)transParams::gpuSelectionList[i],msg))
+            {
+                std::cerr<<msg<<std::endl;
+            }
         }
         std::vector<int> potentialMapDataDims = {1,1,1};
         for(size_t i=0; i<dims.getDataSize(); ++i)
diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
index 29597a6..70a9859 100644
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
@@ -1,6 +1,7 @@
 #include "TVAL.h"
 
 #include "Function2D.h"
+#include "Matrix.h"
 #include "tval3gpu3d.h"
 
 #include "mkl_spblas.h"
@@ -153,10 +154,10 @@ namespace Recon
     Aurora::Matrix callTval3(Aurora::Sparse& M, Aurora::Matrix& b,const Aurora::Matrix& dims,int device, TVALOptions& opt)
     {
         checkAndScale(M,b,(size_t)Aurora::prod(dims).getScalar());
-        int * xIdxs = new int[M.getColVector().getDataSize()];
-        std::copy(M.getColVector().getData(),M.getColVector().getData()+M.getColVector().getDataSize(),xIdxs);
-        int * yIdxs = new int[M.getRowVector().getDataSize()];
-        std::copy(M.getRowVector().getData(),M.getRowVector().getData()+M.getRowVector().getDataSize(),yIdxs);
+        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();
         int nz = std::max(M.getColVector().getDataSize(),std::max(M.getRowVector().getDataSize(),M.getValVector().getDataSize()));
         float* bData = new float[b.getDataSize()];
diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp
index 2d691b0..ea5dc00 100644
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp
@@ -29,6 +29,7 @@ namespace Recon
         {
             //callTval3
             TVALOptions opt;
+            opt.nonneg = solverOptions.nonNeg;
             opt.bent = false;
             opt.tol = 1E-10;
             opt.maxit = niter;
diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp
index 74c7fcf..ece1767 100644
--- a/test/Reconstruction_Test.cpp
+++ b/test/Reconstruction_Test.cpp
@@ -23,6 +23,9 @@ inline double fourDecimalRound(double src){
 #define EXPECT_DOUBLE_AE(valueA,valueB)\
     EXPECT_DOUBLE_EQ(fourDecimalRound(valueA),fourDecimalRound(valueB))
 
+    #define ASSERT_DOUBLE_AE(valueA,valueB)\
+    ASSERT_DOUBLE_EQ(fourDecimalRound(valueA),fourDecimalRound(valueB))
+
 class Reconstruction_Test : public ::testing::Test {
 protected:
     static void SetUpReconstructionTester() {
@@ -48,11 +51,13 @@ TEST_F(Reconstruction_Test, reconstructArt) {
     auto res = m.read("res");
     auto senderList  = m.read("senderList");
     auto SOS_IN_WATER  = m.read("SOS_IN_WATER");
-    auto f1 = m.read("f1");
+    
+    MatlabReader m2("/home/krad/TestData/gpuresult.mat");
+    auto f1 = m2.read("out");
     auto result = Recon::reconstructArt(data, dataAtt, dims, senderList, receiverList, res, SOS_IN_WATER[0]);
     for(size_t i=0; i<f1.getDataSize(); ++i)
     {
-        EXPECT_DOUBLE_AE(f1[i], result.outSOS[i]);
+        ASSERT_DOUBLE_AE(f1[i], result.outSOS[i])<<"index:"<<i;
     }
 }
 

From c1eeadf72a136c6982b535bf1aed24bf8242a70d Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Fri, 2 Jun 2023 11:34:14 +0800
Subject: [PATCH 2/2] Fix TVAL sparse input(convert form coo to csr)

---
 .../solvingEquationSystem/TVAL/TVAL.cpp       | 20 ++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
index 70a9859..152b2e5 100644
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
@@ -159,13 +159,31 @@ namespace Recon
         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 cols = M.getM(), rows = 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;
+        mkl_sparse_d_create_coo(&A, sparse_index_base_t::SPARSE_INDEX_BASE_ZERO,  rows, cols,  nz,  yIdxs, xIdxs,values.getData());
+        mkl_sparse_convert_csr(A, sparse_operation_t::SPARSE_OPERATION_NON_TRANSPOSE, &csrA);
+        int n_rows,n_cols;
+        int *rows_start,*rows_end,*col_indx;
+        double * csrValues;
+        sparse_index_base_t index;
+        mkl_sparse_d_export_csr(csrA, &index, &n_rows, &n_cols, &rows_start, &rows_end, &col_indx, &csrValues);
+        int *row_idx = new int[n_rows+1];
+        std::copy(rows_start,rows_start+n_rows,row_idx);
+        row_idx[n_rows] = rows_end[n_rows-1];
         float* bData = new float[b.getDataSize()];
         std::copy(b.getData(),b.getData()+b.getDataSize(),bData);
         size_t bDims[3]={(size_t)b.getDimSize(0),(size_t)b.getDimSize(1),(size_t)b.getDimSize(2)};
         size_t rdims[3] = {(size_t)dims[0], (size_t)dims[1], (size_t)dims[2]};
         bool pagelocked = false;
-        auto result = TVALGPU(xIdxs, yIdxs, values.getData(), M.getM(), M.getN(), nz, bData, bDims, rdims, opt, device, false);
+
+
+        auto result = TVALGPU(row_idx, col_indx,  csrValues, M.getM(), M.getN(), nz, bData, bDims, rdims, opt, device, false);
+        mkl_sparse_destroy(A);
+        mkl_sparse_destroy(csrA);
+        delete [] row_idx;
         delete [] xIdxs;
         delete [] yIdxs;
         delete [] bData;