LAPACKE C ++ комплексная матричная инверсия - PullRequest
Новая
       73

LAPACKE C ++ комплексная матричная инверсия

1 голос
/ 21 января 2020

Я пытаюсь инвертировать сложную матрицу в C ++, используя пакет LAPACKE. Мой код выглядит следующим образом: 

void inv(std::vector<std::vector<std::complex<double>>> &ans, std::vector<std::vector<std::complex<double>>> MAT){
    int N = MAT.size();
    int *IPIV = new int[N];
    std::complex<double> ** arr = new std::complex<double>*[N];
    for (int i = 0; i<N; i++){
        arr[i] = new std::complex<double>[N];
        for (int j = 0; j<N; j++){
            arr[i][j] = MAT[i][j];
        }
    }

    LAPACKE_cgetrf(LAPACK_ROW_MAJOR, N, N, arr[0], N, IPIV); 
    LAPACKE_cgetri(LAPACK_ROW_MAJOR, N, arr[0], N, IPIV);
    for (int i = 0; i<N; i++){
        for (int j = 0; j<N; j++){                                                                                                                                                      
            ans[i][j] = arr[i][j];                                                                                                                                                  
        }                                                                                                                                                                       
    }                                                                                                                                                                                                                                                                                                                                                       
    delete IPIV;                                                                                                                                                                
    delete arr;                                                                                                                                                             
}

, который выдает ошибку: 

utils.cc:253:59: error: cannot convert 'std::complex<double>*' to '__complex__ 
float*' for argument '4' to 'int LAPACKE_cgetrf(int, int, int, __complex__ float*, 
int, int* '                                                                                                                                                                                
LAPACKE_cgetrf(LAPACK_ROW_MAJOR, N, N, arr[0], N, IPIV);                                                                                                                                                                          ^                                                                                                                
utils.cc:254:56: error: cannot convert 'std::complex<double>*' to '__complex__ 
float*' for argument '3' to 'int LAPACKE_cgetri(int, int, __complex__ float*, int, 
const int*)'                                                                                                                                                                               
LAPACKE_cgetri(LAPACK_ROW_MAJOR, N, arr[0], N, IPIV);

Я пытался сделать: 

__complex__ float ** arr = new __complex__ float *[N];                                                                                                                      
for (int i = 0; i<N; i++){                                                                                                                                                      
    arr[i] = new __complex__ float [N];                                                                                                                                         
    for (int j = 0; j<N; j++){                                                                                                                                                      
        arr[i][j] = reinterpret_cast<float(&)[2]>(MAT[i][j]);                                                                                                                   
    }                                                                                                                                                                       
}

для первого l oop в функция, но она дает: 

 error: cannot convert 'float [2]' to '__complex__ float' in assignment                                                                                                  
 arr[i][j] = reinterpret_cast<float(&)[2]>(MAT[i][j]);

что я делаю не так? У меня есть другая функция inv, которая отлично работает с double. Это сложные типы, с которыми у меня проблемы.

Спасибо

1 Ответ

0 голосов
/ 21 января 2020

LAPACK может ожидать непрерывного блока данных. Вдоль линий 

void inv(std::vector<std::vector<std::complex<double>>> &ans, std::vector<std::vector<std::complex<double>>> MAT){
    using cfloat = __complex__ float;
    int N = MAT.size();

    cfloat* arr = new cfloat [N*N];
    for (int i = 0; i != N; ++i) {            
        for (int j = 0; j != N; ++j){
            int idx = i*N + j;
            arr[idx] = (MAT[i][j]).real() + _Complex_I*(MAT[i][j]).imag();
        }
    }

    int* IPIV = new int[N];

    LAPACKE_cgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV); 
    LAPACKE_cgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);
    for (int i = 0; i != N; ++i){
        for (int j = 0; j != N; ++j){
            int idx = i*N + j;
            ans[i][j] = arr[idx];                                                                                                 
        }                                                                                                                                                                       
    }                                                                                                                                                                                                                                                                                                                                                       
    delete [] IPIV;                                                                                                                                                                
    delete [] arr;                                                                                                                                                             
}
...