diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp index 715ace1..9867575 100644 --- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp +++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp @@ -54,141 +54,13 @@ 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!"< 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!"< 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!"< 1 + 1e-10 - // b = b ./ sqrt(s2); - // M = M ./ sqrt( s2); - //end + auto s2 = eigs(M); + if (s2> 1 + 1e-10) + { + b = b / sqrt(s2); + M.getValVector() = M.getValVector() / sqrt( s2); + } } Aurora::Matrix callTval3(Aurora::Sparse& M, Aurora::Matrix& b,const Aurora::Matrix& dims,int device, TVALOptions& opt)