sunwen/UR
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add check and scale

Browse Source
This commit is contained in:
kradchen
2023-06-02 14:37:05 +08:00
parent d5cbcc32cc
commit 10fbe2a763
1 changed files with 6 additions and 134 deletions
Download Patch File Download Diff File Expand all files Collapse all files

140
src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
View File

@@ -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!"<<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：以下操作意义不明
//fh = @(x) ((M * x)' * M)';
//s2 = eigs(fh,n,1,'lm',eopts);
//if real(s2) > 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)