From 367176f86c88a767c99c780792c3f4c459f0a5db Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Thu, 1 Jun 2023 13:24:56 +0800
Subject: [PATCH] Add Check and Scale(unfinished).

---
 .../solvingEquationSystem/TVAL/TVAL.cpp       | 126 +++++++++++++++++-
 .../solvingEquationSystem/TVAL/TVAL.h         |   4 +-
 test/Reconstruction_Test.cpp                  |  14 ++
 3 files changed, 140 insertions(+), 4 deletions(-)

diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
index 817b3fc..29597a6 100644
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
@@ -1,25 +1,144 @@
 #include "TVAL.h"
 
 #include "Function2D.h"
-#include "Matrix.h"
 #include "tval3gpu3d.h"
+
+#include "mkl_spblas.h"
+#include "mkl_solvers_ee.h"
 #include <algorithm>
 #include <cstddef>
+#include <iostream>
+#include <vector>
 
 namespace Recon
 {
     void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n){
+        return;
         //TODO:暂时啥都没做
         
         //无用的代码？
         //opts.scale_A = true;
         //opts.consist_mu = false;
-        
+        sparse_matrix_t A;
+        size_t length = M.getRowVector().getDataSize();
+        int *row_indx = new int[length];
+        int *col_indx = new int[length];
+        size_t rows = M.getM(), cols = M.getN();
+        std::copy(M.getRowVector().getData(),M.getRowVector().getData()+length,row_indx);
+        std::copy(M.getColVector().getData(),M.getColVector().getData()+length,col_indx);
+        auto status = mkl_sparse_d_create_coo(&A,sparse_index_base_t::SPARSE_INDEX_BASE_ZERO,rows,cols, M.getValVector().getDataSize(),
+             row_indx, col_indx, M.getValVector().getData());
+        if(status != SPARSE_STATUS_SUCCESS)
+        {
+            std::cerr<<"create mkl sparse error!checkAndScale fail!"<<std::endl;
+            return;
+        }
+        int *row_indx2 = new int[n];
+        int *col_indx2 = new int[n];
+        for (size_t i = 0; i < n; i++)
+        {
+            row_indx2[i] = i;
+            col_indx2[i] = i;
+        }
+        double * value2 = new double[n]{1};
+        sparse_matrix_t eye;
+        status = mkl_sparse_d_create_coo(&eye,sparse_index_base_t::SPARSE_INDEX_BASE_ZERO,n, n, n,
+             row_indx2, col_indx2, value2);
+        if(status != SPARSE_STATUS_SUCCESS)
+        {
+            std::cerr<<"create eye mkl sparse error!checkAndScale fail!"<<std::endl;
+            return;
+        }
+        std::vector<int> vCols;
+        {
+            int *col_indx2 = new int[length];
+            std::copy(M.getColVector().getData(),M.getColVector().getData()+length,col_indx2);
+            std::sort(col_indx2,col_indx2+length);
+            vCols.push_back(col_indx2[0]);
+            double lastV = col_indx2[0];
+            for (size_t i = 1; i < length; i++)
+            {
+                if(col_indx2[i]==lastV)continue;
+                vCols.push_back(col_indx2[i]);
+                lastV = col_indx2[i];
+            }
+                    
+        }
+        double * vec = new double[M.getN()]{0};
+        double * vecOut1 = new double[M.getM()]{0};
+        double * vecOut2 = new double[M.getN()]{0};
+        int colIdxIdx = 0;
+        for (size_t i = 0; i < 6000; i++)
+        {
+            if ( vCols[colIdxIdx]>i)continue;
+            else colIdxIdx++;
+            vec[i]=1;
+            if (i-1>0)vec[i-1] = 0;
+            matrix_descr descr;
+            descr.type = SPARSE_MATRIX_TYPE_GENERAL;
+            status = mkl_sparse_d_mv(sparse_operation_t::SPARSE_OPERATION_NON_TRANSPOSE, 1.0, A, descr, vec, 0.0, vecOut1);
+            if(status != SPARSE_STATUS_SUCCESS)
+            {
+                std::cerr<<"sparse * vector error!checkAndScale fail!"<<std::endl;
+                return;
+            }
+            status = mkl_sparse_d_mv(sparse_operation_t::SPARSE_OPERATION_TRANSPOSE, 1.0, A, descr, vecOut1, 0.0, vecOut2);
+            if(status != SPARSE_STATUS_SUCCESS)
+            {
+                std::cerr<<"sparse * vector error!checkAndScale fail!"<<std::endl;
+                return;
+            }
+            for (size_t j = 0; j < M.getN(); j++)
+            {
+                if (vecOut2[j]!=0.0){
+                    mkl_sparse_d_set_value(eye, j, i, vecOut2[j]);
+                }   
+            }
+        }
+        delete [] row_indx;
+        delete [] col_indx;
+        delete [] vec;
+        delete [] vecOut1;
+        delete [] vecOut2;
+        matrix_descr descr;
+        descr.type = SPARSE_MATRIX_TYPE_GENERAL;
+        int k = 0;
+        int pm[128]{0};
+        mkl_sparse_ee_init(pm);
+        pm[1]=3;
+        double s2 = 0.0, res = 0.0;
+        double *X = new double[n];
+        char which = 'L';
+        status = mkl_sparse_d_ev(&which, pm, eye, descr, 1, &k, &s2, X, &res);
+        if(status != SPARSE_STATUS_SUCCESS)
+        {
+            std::cerr<<"sparse eigen error!checkAndScale fail!"<<std::endl;
+            return;
+        }
+        status = mkl_sparse_destroy(eye);
+        if(status != SPARSE_STATUS_SUCCESS)
+        {
+            std::cerr<<"sparse destory error!checkAndScale fail!"<<std::endl;
+            return;
+        }
+        delete [] row_indx2;
+        delete [] col_indx2;
+        delete [] value2;
+        delete [] X;
+        if (s2> 1 + 1e-10){
+            b = b/sqrt(s2);
+            M.getValVector() = M.getValVector()/sqrt(s2);
+        }
         //TODO:以下意义不明,暂定为判断是否为复数存储
         //if ~isreal(M * rand(n,1))
         //  eopts.isreal = false;
         //end
-
+        status = mkl_sparse_destroy(A);
+        if(status != SPARSE_STATUS_SUCCESS)
+        {
+            std::cerr<<"sparse destory error!checkAndScale fail!"<<std::endl;
+            return;
+        }
         bool isreal = M.getValVector().getValueType() == Aurora::Normal;
 
         //TODO：以下操作意义不明
@@ -30,6 +149,7 @@ namespace Recon
         //    M = M ./ sqrt( s2);
         //end
     }
+
     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());
diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
index e870fe8..a2fbf67 100644
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
@@ -1,10 +1,12 @@
 #ifndef __TVAL_H__
 #define __TVAL_H__
-#include "Matrix.h"
 #include "Sparse.h"
 #include "tvalstruct.h"
 
 namespace Recon {
+
+void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n);
+
 Aurora::Matrix callTval3(Aurora::Sparse &M, Aurora::Matrix &b,
                          const Aurora::Matrix &dims, int device,
                          struct TVALOptions &options);
diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp
index a5c9d5c..dbfdcf6 100644
--- a/test/Reconstruction_Test.cpp
+++ b/test/Reconstruction_Test.cpp
@@ -5,11 +5,13 @@
 #include "Function1D.h"
 #include "MatlabReader.h"
 #include "Matrix.h"
+#include "Sparse.h"
 #include "config/config.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/DGradient.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/FMM.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
 #include "transmissionReconstruction/reconstruction/reconstruction.h"
+#include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h"
 
 
 
@@ -37,6 +39,18 @@ protected:
     }
 };
 
+TEST_F(Reconstruction_Test, checkAndScale) {
+    MatlabReader m("/home/sun/testData/buildMatrix.mat");
+    auto i  = m.read("i1");
+    auto j  = m.read("j1");
+    auto s = m.read("s1");
+    MatlabReader m2("/home/sun/testData/checkAndScale.mat");
+    auto b = m2.read("b");
+
+    Aurora::Sparse sparse(i,j,s,709613,1196032);
+    Recon::checkAndScale(sparse,b,1196032);
+}
+
 TEST_F(Reconstruction_Test, calculateMinimalMaximalTransducerPositions) {
     auto senderList  = Aurora::Matrix::fromRawData(new double[6]{1, 2, 3, 1, 2, 4}, 3, 2);
     auto receiverList  = Aurora::Matrix::fromRawData(new double[6]{1, 8, 3, 1, 2, 1}, 3, 2);