У меня есть код C ++, который должен вычислять инверсию различных матриц (ковариационных матриц).
На данный момент существует первая версия, которая явно выполняет разложение Холецкого с помощью вычисление факторизации матриц tri angular и после вычисления обратной матрицы.
Эта явная подпрограмма использует директивы OpenMP #pragma для оптимизации, просто улучшая ее, например, выполняя двойное l oop , но это не topi c здесь моей проблемы:
#pragma omp parallel for schedule(dynamic, num_threads)
for(int i=0; i<F_matrix_A.size(); i++){
for(int j=i; j<F_matrix_A.size(); j++){
Fisher_new[i][j] = Fisher_new[j][i];
}
Ниже приведен код этой подпрограммы с именем matrix_inverse, который принимает 2 "матрицы" в качестве аргументов (F_Matrix_A и F_previous) типа, подобного "матрице", то есть двумерный массив, объявленный как (Мне нужно глубже изучить работу этой подпрограммы, чтобы улучшить gr asp: я подозреваю, что процедура Cramer инвертирует матрицу в этой явной подпрограмме, но не уверен на 100%):
vector<vector<double>> F_matrix
Часть I) Текущий код процедуры матрицы инверсии по Разложение Холецкого (что хорошо работает):
vector<vector<double>> matrix_inverse(vector<vector<double>> F_matrix_A, vector<vector<double>> F_previous){
cout<<"Begin the LT matrix inversion"<<endl;
vector<double> F_temp(F_matrix_A.size()); vector<double> temp_prev(F_matrix_A.size());
for(int i=0; i < F_matrix_A.size(); i++){
for(int j=0; j < F_matrix_A.size(); j++){
if(i == j){
for(int k=0; k < F_matrix_A.size(); k++) {
if(F_previous[j][k] != 0){
F_previous[j][k] = F_previous[j][k]/F_matrix_A[j][i];
F_matrix_A[j][k] = F_matrix_A[j][k]/F_matrix_A[j][i];
}
}
}
else{
if(F_matrix_A[j][i] != 0){
for(int k=0; k < F_matrix_A.size(); k++) {
if(F_previous[i][k] != 0){
F_previous[j][k] = F_previous[j][k] - F_matrix_A[j][i]*F_previous[i][k];
F_matrix_A[j][k] = F_matrix_A[j][k] - F_matrix_A[j][i]*F_matrix_A[i][k];
}
}
}
}
}
}
vector<vector<double>> Fisher_new(F_matrix_A.size(), vector<double>(F_matrix_A.size(), 0));
vector<vector<double>> F_previous_T(F_matrix_A.size(), vector<double>(F_matrix_A.size(), 0));
#pragma omp parallel for schedule(dynamic, num_threads)
for(int i=0; i<F_matrix_A.size(); i++){
for(int j=0; j<F_matrix_A.size(); j++){
F_previous_T[i][j] = F_previous[j][i];
}
}
#pragma omp parallel for schedule(dynamic, num_threads)
for(int i=0; i<F_matrix_A.size(); i++){
for(int k=0; k<F_matrix_A.size(); k++){
for(int j=0; j<=i; j++){
if(F_previous_T[i][k] != 0 && F_previous[k][j] !=0){
Fisher_new[i][j] += F_previous_T[i][k]*F_previous[k][j];
}
}
}
}
cout<<"Building the inverse covariance matrix"<<endl;
#pragma omp parallel for schedule(dynamic, num_threads)
for(int i=0; i<F_matrix_A.size(); i++){
for(int j=i; j<F_matrix_A.size(); j++){
Fisher_new[i][j] = Fisher_new[j][i];
}
}
return Fisher_new;
}
Часть II) Телефонный код:
И, например, эта подпрограмма называется так:
vector<vector<double>> CO_CL(lsize*(2*Dim_x+Dim_y), vector<double>(lsize*(2*Dim_x+Dim_y), 0));
vector<vector<double>> CO_CL_AB(lsize*(2*Dim_x+Dim_y), vector<double>(lsize*(2*Dim_x+Dim_y), 0));
vector<vector<double>> CO_I(lsize*(2*Dim_x+Dim_y), vector<double>(lsize*(2*Dim_x+Dim_y), 0));
vector<vector<double>> CO_CL_D(lsize*(2*Dim_x+Dim_y), vector<double>(lsize*(2*Dim_x+Dim_y), 0));
for(int i=0; i<CO_I.size(); i++){
CO_I[i][i] = 1.;
}
for(int i=0; i<CO_WL_I.size(); i++){
CO_WL_I[i][i] = 1.;
}
for(int FX=0; FX<Fisher_M.size(); FX++){
for(int FY=FX; FY<Fisher_M.size(); FY++){
for(int lll=0; lll<lsize; lll++){
for(int i=0;i<zrange.size();i++){
...
}
... // Computing of CO_CL matrixes
if(probe == "GCp"){
if(FX == 0 && FY == 0){
vector<vector<double>> CO_CL_temp(lsize*Dim_x, vector<double>(lsize*Dim_x, 0));
CO_CL_AB = CO_CL_temp; CO_I = CO_CL_temp;
for(int i=0; i<lsize*Dim_x; i++){
CO_I[i][i] = 1.;
for(int j=0; j<lsize*Dim_x; j++){
CO_CL_temp[i][j] = CO_CL[i][j];
}
}
CO_CL = CO_CL_temp;
}
vector<vector<double>> CO_CL_temp_D(lsize*Dim_x, vector<double>(lsize*Dim_x, 0));
for(int i=0; i<lsize*Dim_x; i++){
for(int j=0; j<lsize*Dim_x; j++){
CO_CL_temp_D[i][j] = CO_CL_D[i][j];
}
}
CO_CL_D = CO_CL_temp_D;
}
if(probe == "WL"){
if(FX == 0 && FY == 0){
vector<vector<double>> CO_CL_temp(lsize*Dim_x, vector<double>(lsize*Dim_x, 0));
CO_CL_AB = CO_CL_temp; CO_I = CO_CL_temp;
for(int i=lsize*Dim_x; i<2*lsize*Dim_x; i++){
CO_I[i-lsize*Dim_x][i-lsize*Dim_x] = 1.;
for(int j=lsize*Dim_x; j<2*lsize*Dim_x; j++){
CO_CL_temp[i-lsize*Dim_x][j-lsize*Dim_x] = CO_CL[i][j];
}
}
CO_CL = CO_CL_temp;
}
vector<vector<double>> CO_CL_temp_D(lsize*Dim_x, vector<double>(lsize*Dim_x, 0));
for(int i=lsize*Dim_x; i<2*lsize*Dim_x; i++){
for(int j=lsize*Dim_x; j<2*lsize*Dim_x; j++){
CO_CL_temp_D[i-lsize*Dim_x][j-lsize*Dim_x] = CO_CL_D[i][j];
}
}
CO_CL_D = CO_CL_temp_D;
}
if(FX==0 && FY == 0){
cout<<"Begin Cholesky decomposition"<<endl;
for(int i=0; i<CO_CL.size(); i++){
for(int j=0; j<=i; j++){
double sum=0;
if(i==j){
for(int k=0; k<j; k++){
sum += pow(CO_CL_AB[j][k], 2);
}
CO_CL_AB[i][j] = pow(CO_CL[j][j]-sum,0.5);
}
else{
for(int k=0; k<j; k++){
sum += CO_CL_AB[i][k]*CO_CL_AB[j][k];
}
CO_CL_AB[i][j] = (CO_CL[i][j] - sum)/CO_CL_AB[j][j];
}
}
}
// Original version which give good results
CO_CL = matrix_inverse(CO_CL_AB, CO_I);
// Lapacke version
cout << "here lapack" << endl;
// BAD RESULTS if I replace below by matrix_inverse above
//CO_CL = matrix_inverse_lapack(CO_CL_AB);
}
for(int i=0; i<CO_CL.size(); i++){
for(int j=0; j<CO_CL.size(); j++){
Fisher_M[FX][FY] += CO_CL[i][j]*CO_CL_D[i][j];
}
}
Fisher_M[FY][FX] = Fisher_M[FX][FY];
CO_CL_D = CO_CL_ref;
if(probe != "GCp"){
for(int lll=lsize; lll<l.size(); lll++){
if(FX == 0 && FY == 0){
if(lll==lsize){
cout<<"Computing high l's covariance matrix"<<endl;
}
ifstream ifile_ABCD_LL(C_folder[1][1]+"/C_fid"+"/COVAR_fid_"+to_string(lll));
for(int i=0;i<zrange.size();i++){
for(int j=0;j<zrange.size();j++){
ifile_ABCD_LL>>C_ij_ABCD_LL[i][j];
}
}
ifile_ABCD_LL.close();
}
ifstream ifile_ABCD_LL_up_PX(C_folder[1][1]+"/C_"+param_chain[FX]+"_up"+"/COVAR_up_"+to_string(lll));
ifstream ifile_ABCD_LL_dw_PX(C_folder[1][1]+"/C_"+param_chain[FX]+"_dw"+"/COVAR_dw_"+to_string(lll));
ifstream ifile_ABCD_LL_up_PY(C_folder[1][1]+"/C_"+param_chain[FY]+"_up"+"/COVAR_up_"+to_string(lll));
ifstream ifile_ABCD_LL_dw_PY(C_folder[1][1]+"/C_"+param_chain[FY]+"_dw"+"/COVAR_dw_"+to_string(lll));
for(int i=0;i<zrange.size();i++){
for(int j=0;j<zrange.size();j++){
ifile_ABCD_LL_up_PX>>C_ij_ABCD_LL_up_PX[i][j]; ifile_ABCD_LL_dw_PX>>C_ij_ABCD_LL_dw_PX[i][j];
ifile_ABCD_LL_up_PY>>C_ij_ABCD_LL_up_PY[i][j]; ifile_ABCD_LL_dw_PY>>C_ij_ABCD_LL_dw_PY[i][j];
}
}
ifile_ABCD_LL_up_PX.close(); ifile_ABCD_LL_dw_PX.close(); ifile_ABCD_LL_up_PY.close(); ifile_ABCD_LL_dw_PY.close();
I_55=0;
for(int I1=0; I1 < zrange.size(); I1++){
for(int I2=0; I2 < zrange.size(); I2++){
for(int I3=0; I3 < zrange.size(); I3++){
for(int I4=0; I4 < zrange.size(); I4++){
if(I2 <= I1 && I4 <= I3){
if(FX == 0 && FY == 0){
CC_LLLL[I_55] = (C_ij_ABCD_LL[I1][I3]*C_ij_ABCD_LL[I2][I4] + C_ij_ABCD_LL[I1][I4]*C_ij_ABCD_LL[I2][I3])/((2*l[lll]+1)*fsky*delta_l); //LLLL
}
if(FX == relat_index && FY == relat_index){
CC_LLLL_D[I_55] = (C_ij_ABCD_LL_up_PX[I1][I2]-C_ij_ABCD_LL_dw_PX[I1][I2])/(2*steps_all[FX]) * (C_ij_ABCD_LL_up_PY[I3][I4]-C_ij_ABCD_LL_dw_PY[I3][I4])/(2*steps_all[FY]);
}
else if(FX == relat_index && FY != relat_index){
CC_LLLL_D[I_55] = (C_ij_ABCD_LL_up_PX[I1][I2]-C_ij_ABCD_LL_dw_PX[I1][I2])/(2*steps_all[FX]) * (C_ij_ABCD_LL_up_PY[I3][I4]-C_ij_ABCD_LL_dw_PY[I3][I4])/(2*fid_all[FY]*steps_all[FY]);
}
else if(FX != relat_index && FY == relat_index){
CC_LLLL_D[I_55] = (C_ij_ABCD_LL_up_PX[I1][I2]-C_ij_ABCD_LL_dw_PX[I1][I2])/(2*fid_all[FX]*steps_all[FX]) * (C_ij_ABCD_LL_up_PY[I3][I4]-C_ij_ABCD_LL_dw_PY[I3][I4])/(2*steps_all[FY]);
}
else{
CC_LLLL_D[I_55] = (C_ij_ABCD_LL_up_PX[I1][I2]-C_ij_ABCD_LL_dw_PX[I1][I2])/(2*fid_all[FX]*steps_all[FX]) * (C_ij_ABCD_LL_up_PY[I3][I4]-C_ij_ABCD_LL_dw_PY[I3][I4])/(2*fid_all[FY]*steps_all[FY]);
}
I_55++;
}
}
}
}
}
k_vec=0;
for(int i=0; i<Dim_x; i++){
for(int j=0; j<Dim_x; j++){
if(FX == 0 && FY == 0){
CC_LLLL_R[i][j] = CC_LLLL[k_vec];
}
CC_LLLL_DR[i][j] = CC_LLLL_D[k_vec];
k_vec++;
}
}
for(int z1=0; z1<Dim_x; z1++){
for(int z2=0; z2<Dim_x; z2++){
if(FX == 0 && FY == 0){
CO_CL_WL[z1*lsize2+lll-lsize][z2*lsize2+lll-lsize] = CC_LLLL_R[z1][z2];
}
CO_CL_WL_D[z1*lsize2+lll-lsize][z2*lsize2+lll-lsize] = CC_LLLL_DR[z1][z2];
}
}
}
if(probe == "WL"){
if(FX == 0 && FY == 0){
vector<vector<double>> CO_CL_temp(lsize2*Dim_x, vector<double>(lsize2*Dim_x, 0));
CO_CL_WL_AB = CO_CL_temp; CO_WL_I = CO_CL_temp;
for(int i=0; i<lsize2*Dim_x; i++){
CO_WL_I[i][i] = 1.;
for(int j=0; j<lsize2*Dim_x; j++){
CO_CL_temp[i][j] = CO_CL_WL[i][j];
}
}
CO_CL_WL = CO_CL_temp;
}
vector<vector<double>> CO_CL_temp_D(lsize2*Dim_x, vector<double>(lsize2*Dim_x, 0));
for(int i=0; i<lsize2*Dim_x; i++){
for(int j=0; j<lsize2*Dim_x; j++){
CO_CL_temp_D[i][j] = CO_CL_WL_D[i][j];
}
}
CO_CL_WL_D = CO_CL_temp_D;
}
if (FX==0 && FY == 0){
cout<<"Begin Cholesky decomposition"<<endl;
for(int i=0; i<CO_CL_WL.size(); i++){
for(int j=0; j<=i; j++){
double sum=0;
if(i==j){
for(int k=0; k<j; k++){
sum += pow(CO_CL_WL_AB[j][k], 2);
}
CO_CL_WL_AB[i][j] = pow(CO_CL_WL[j][j]-sum,0.5);
}
else{
for(int k=0; k<j; k++){
sum += CO_CL_WL_AB[i][k]*CO_CL_WL_AB[j][k];
}
CO_CL_WL_AB[i][j] = (CO_CL_WL[i][j] - sum)/CO_CL_WL_AB[j][j];
}
}
}
// Original version : VALID RESULTS
CO_CL_WL = matrix_inverse(CO_CL_WL_AB, CO_WL_I);
// Lapacke version : BAD RESULTS
cout << "here lapack << endl;
//CO_CL_WL = matrix_inverse_lapack(CO_CL_WL_AB);
cout<<"Computing the Fisher matrix elements"<<endl;
}
for(int i=0; i<CO_CL_WL.size(); i++){
for(int j=0; j<CO_CL_WL.size(); j++){
Fisher_M[FX][FY] += CO_CL_WL[i][j]*CO_CL_WL_D[i][j];
}
}
Fisher_M[FY][FX] = Fisher_M[FX][FY];
}
CO_CL_WL_D = CO_CL_ref;
}
}
ofstream outF; outF.open("output/"+F_mat_name);
for(int FX=0; FX<Fisher_M.size(); FX++){
for(int FY=0; FY<Fisher_M.size(); FY++){
outF<<setprecision(12)<<scientific<<Fisher_M[FX][FY]<<" ";
}
outF<<endl;
}
outF.close();
cout<<probe+" Fisher matrix saved"<<endl;
int DE_X = Fisher_M.size()+1; int DE_Y = Fisher_M.size()+1;
for(int i=0; i<param_chain.size(); i++){
if(param_chain[i] == "w0"){
DE_X = i;
}
if(param_chain[i] == "wa"){
DE_Y = i;
}
}
if(DE_X != Fisher_M.size()+1 && DE_Y != Fisher_M.size()+1){
vector<vector<double>> F_I(Fisher_M.size(), vector<double>(Fisher_M.size(), 0));
vector<vector<double>> Fisher_AB(Fisher_M.size(), vector<double>(Fisher_M.size(), 0));
for(int i=0; i<Fisher_M.size(); i++){
F_I[i][i] = 1;
}
for(int i=0; i<Fisher_M.size(); i++){
for(int j=0; j<=i; j++){
double sum=0;
if(i==j){
for(int k=0; k<j; k++){
sum += pow(Fisher_AB[j][k], 2);
}
Fisher_AB[i][j] = pow(Fisher_M[j][j]-sum,0.5);
}
else{
for(int k=0; k<j; k++){
sum += Fisher_AB[i][k]*Fisher_AB[j][k];
}
Fisher_AB[i][j] = (Fisher_M[i][j] - sum)/Fisher_AB[j][j];
}
}
}
Fisher_M = matrix_inverse(Fisher_AB, F_I);
//Fisher_M = matrix_inverse_lapack(Fisher_AB);
ofstream invFile; invFile.open("invMatrix_Fisher_M.dat");
for(int aa=0; aa < Fisher_M.size(); aa++){
for(int bb=0; bb < Fisher_M.size(); bb++){
invFile << setprecision(12) << scientific << Fisher_M[aa][bb] << " ";
}
invFile << endl;
}
invFile.close();
double FoM = 1./pow(Fisher_M[DE_X][DE_X]*Fisher_M[DE_Y][DE_Y] - Fisher_M[DE_X][DE_Y]*Fisher_M[DE_X][DE_Y],0.5);
cout<<"value = "<<FoM;
}
Вы можете видеть выше всех вызовов функции matrix_inverse. Проблема обычно в следующих параметрах:
// Original version : VALID RESULTS
CO_CL = matrix_inverse(CO_CL_AB, CO_I);
// Lapacke version : hesitate for inversion between CO_CL or CO_CL_AB matrixes : BOTH give BAD RESULTS at the end of code on produced final matrixes
// as argument of matrix_inverse_lapack ??
//CO_CL = matrix_inverse_lapack(CO_CL_AB);
// OR
//CO_CL = matrix_inverse_lapack(CO_CL);
// SO, I WONDER WHAT I HAVE GOT TO TRY ?
}
Как вы можете видеть в конце этого фрагмента кода, я хотел бы использовать LAPACKE функции dgetrf и dgetri для выполнения этого Разложение Холецкого. Я сомневаюсь, вызывая основанную LAPACKE функцию matrix_inverse_lapack между матрицей CO_CL или CO_CL_AB. Я запрограммировал эту процедуру под названием «matrix_inverse_lapack» следующим образом:
Часть III) Моя собственная обратная процедура с LAPACKE
vector<vector<double>> matrix_inverse_lapack(vector<vector<double>> F_matrix) {
// Index for loop and arrays
int i, j, ip, idx;
// Size of F_matrix
int N = F_matrix.size();
int *IPIV = new int[N];
// Statement of main array to inverse
double *arr = new double[N*N];
// Statement of returned matrix of size N
vector<vector<double> > F_final(N, vector<double>(N));
for (i = 0; i<N; i++){
for (j = 0; j<N; j++){
idx = i*N + j;
arr[idx] = F_matrix[i][j];
}
}
// LAPACKE
int info1 = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);
int info2 = LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);
cout << "info1 =" << info1 << endl;
cout << "info2 =" << info2 << endl;
for (i = 0; i<N; i++){
for (j = 0; j<N; j++){
idx = i*N + j;
F_final[i][j] = arr[idx];
}
}
delete[] IPIV;
delete[] arr;
// Return inverse matrix
return F_final;
}
Но, к сожалению, я не может получить ту же самую окончательную инвертированную матрицу между использованием первой явной подпрограммы matrix_inverse и моей собственной функцией matrix_inverse_lapack (если я выберу вызов CO_CL = matrix_inverse_lapack(CO_CL), большинство значений будут закрыты (ниже 0,1% от относительная разница) но есть расхождения по некоторым значениям (до 40% относительной разности).
ОБНОВЛЕНИЕ 1:
Я не уверен, если , используя dgetrf и dgetri, которые должны обрабатывать общий случай матрицы симметрии c, я могу использовать разложение Холецкого ??
Обычно я могу использовать разложение Холецкого для матрицы симметрии c к обратному.
ОБНОВЛЕНИЕ 2: Я связался с человеком, который закодировал эту процедуру matrix_inverse (не версия LAPACKE, другая).
Он сказал мне, что Второй аргумент F_Previous представляет на самом деле в начале матр ix идентичность, и он используется как «опора Гаусса» в исключении Гаусса. В начале кода F_Previous соответствует матрице идентичности.
А в конце процедуры matrix_inverse я умножаю нижнюю матрицу tri angular на верхнюю матрицу tri angular, как this:
Fisher_new[i][j] += F_previous_T[i][k]*F_previous[k][j];
Чтобы обосновать это, я проиллюстрирую эту операцию выше из статьи Разложение Холецкого
, в которой говорится:
Не знаю, если этот метод, использованный выше в статье, соответствует этому продукту между R и R * (Является ли R * обязательным транспонированием R в пространство Реала?).
Теперь проблема в том, чтобы узнать, где я могу ввести функции LAPACKE_dgetrf и LAPACKE_dgetrf в коде.
Наивно, на первом этапе я подумал, что замена всех вызовов matrix_inverse routine на matrix_inverse_lapack решит проблему, но я понял, что это не сработает или иным образом, используя неправильные аргументы (ковариационная матрица).
Но у меня все еще есть проблема с кодом вызова, который вызывает функцию inverse_matrix. Он возвращает инверсию своего первого аргумента (Matrix), но мой собственный inverse_matrix_lapacke (который действительно использует dgetrf и dgetri) не возвращает то же самое.
ОБНОВЛЕНИЕ 3 : ВАЖНО
Я помечаю это обновление как ВАЖНОЕ, поскольку я осознал, что перед выполнением всех этих операций обращения матриц у меня нет одинаковых входных матриц между g++-9 gnu compiler и icpc Intel compiler. Таким образом, очевидно, что результаты будут совершенно другими, когда я применю инверсии (с помощью явного метода Холецкого или процедур LAPACKE).
Действительно, существует высокая чувствительность относительно значений во входных матрицах. И, к сожалению, при построении этих входных матриц, учитывая тот факт, что есть некоторые значения, равные нулю с g++-9 gnu compiler, а соответствующие значения с icpc Intel compiler, например, равны 1e-24 или 1e-26 (я имею в виду не нулевое значение в строгих терминах), расхождения возникают после инверсии этих матриц.
Поэтому моя проблема сейчас заключается в работе с одинаковыми входными матрицами для компиляторов g++-9 gnu и icpc Intel. Эта часть строительных матриц является первой частью кода.
Вот параметры, которые я использую для этих 2-х компиляторов:
g++-9 GNU (установлен с brew на MacOS Catalina):
CXX = g++-9 -O3
CXXFLAGS = -Wall -c
icpc Intel (с веб-сайта Intel и находится в /opt/compilers_and_libraries_2020.1.216 /)
CXX = icpc -std=c++11 -O3
CXXFLAGS = -Wall -c -I${MKLROOT}/include
LDFLAGS = -L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core
На данный момент я отключил библиотеку OpenMP в обеих компиляциях: один шаг за шагом для моего проблема, поэтому сначала устранение OpenMP. Мы увидим, что после.
Должны ли мои опции / флаги 2-х компиляторов генерировать одинаковые входные матрицы?
Действительно, я проверяю это и, как я сказал, иногда, это делает округление до нуля или нет.
Вы можете увидеть этап построения моего кода по следующей ссылке:
код для построения входных матриц для инверсии
Как видите, для построения входных матриц Я просто выполняю операции c для их генерации (сложение, деление, умножение): поэтому я не делаю Я понимаю, почему я не получаю одинаковые значения между двумя компиляторами.
Я подозреваю, что icp c Intel использует неявные флаги, которые я не использую с g ++ - 9 GNU , но Я не знаю достаточно подробно параметры по умолчанию компилятора Intel, особенно о поведении округления для очень маленьких значений.
ОБНОВЛЕНИЕ 4: ПРОХОД «Матрицы» двойного BY REFERENCE в подпрограмму LAPACKE
Следуя некоторым советам в комментариях, я хотел бы узнать результаты, если вместо того, чтобы передать значение параметра, я выполню ссылку по ссылке.
Вот моя новая (немного измененная) подпрограмма matrix_inverse_lapack (выходная матрица - * 1173). * F_output"Матрица" (на самом деле vector<vector<double>> тип):
// Passing Matrixes by Reference
void matrix_inverse_lapack(vector<vector<double>> &F_matrix, vector<vector<double>> &F_output) {
// Index for loop and arrays
int i, j, ip, idx;
// Size of F_matrix
int N = F_matrix.size();
int *IPIV = new int[N];
// Statement of main array to inverse
double *arr = new double[N*N];
for (i = 0; i<N; i++){
for (j = 0; j<N; j++){
idx = i*N + j;
arr[idx] = F_matrix[i][j];
}
}
// LAPACKE routines to get inverse of F_Matrix
int info1 = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);
int info2 = LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);
cout << "info1 =" << info1 << endl;
cout << "info2 =" << info2 << endl;
for (i = 0; i<N; i++){
for (j = 0; j<N; j++){
idx = i*N + j;
F_output[i][j] = arr[idx];
}
}
delete[] IPIV;
delete[] arr;
}
Матрица, сгенерированная F_output, не содержит nan значений, но все еще имеются расхождения в окончательных результатах между этими LAPACK (icpc+MKL) version и GNU (g++-9) version (с процедурой classicla, описанной выше, то есть: vector<vector<double>> matrix_inverse(vector<vector<double>> F_matrix_A, vector<vector<double>> F_previous)).
Итак, на данный момент я могу сказать, что LAPACK Intel из MKL не дает таких же результатов, как Явная версия (я имею в виду с нуля).
Обходной путь должен был бы использовать другие альтернативные процедуры LAPACK для инвертирования Матрицы (я напоминаю, фактически тип vector<vector<double>>). Есть ли альтернативы? и если да, то как их использовать для получения инверсии?
В начале этого поста я использовал combination of dgetrf and dgetri, чтобы найти обратную матрицу.
Я видел сейчас по этой ссылке LAPACK ROUTINES LINEAR ALGEBRA
Я провел тест с подпрограммами, реализованными в matrix_inverse_lapack:
int info1 = LAPACKE_dsytrf_rk(LAPACK_ROW_MAJOR, 'U', N, arr, N, diag, IPIV);
int info2 = LAPACKE_dsytri_3(LAPACK_ROW_MAJOR, 'U', N, arr, N, diag, IPIV);
в конце программы (это последнее производит матрицу 20x20, так что эта последняя инверсия мала по размеру). Для всех других инверсий (с большими матрицами, такими как 7000x7000), я сохранил классическую функцию matrix_inverse.
Таким образом, в результате сгенерированная матрица 20x20 такая же, как в GNU g ++ - 9.
Это уже хороший результат, но если я применю matrix_inverse_lapack в первой большой инверсии, конечные матрицы 20x20 будут плохими (по сравнению с правильными матрицами 20x20, полученными с помощью GNU gcc ++ - 9 версия).
Я продолжаю свои исследования, чтобы сделать равными большие матрицы, инвертированные и расположенные в середине моего кода.
pardiso Lapack рутина кажется интересной, но я не знаю, как использовать ее для инверсии Любой пример приветствуется.