0

C++のLapackパッケージを使用して線形方程式系を解きたい。ここからのルーチン、つまりdgesvを使用して、このように実装する予定です。

これは私のコードです：

unsigned short int n=profentries.size();
double **A, *b, *x;
A= new double*[n];
b=new double[n];
x=new double[2];
// Generate the matrices for my equation
for(int i = 0; i < n; ++i)
{
    A[i] = new double[2];
}

for(int i=0; i<n; i++)
{
    A[i][0]=5;
    A[i][1]=1;
    b[i]=3;
}

std::cout<< "printing matrix A:"<<std::endl;
for(int i=0; i<n; i++)
{
    std::cout<< A[i][0]<<" " <<A[i][1]<<std::endl;
}
// Call the LAPACK solver
    x = new double[n];//probably not necessary
dgesv(A, b, 2, x); //wrong result for x!
std::cout<< "printing vector x:"<<std::endl;

/*prints 
3
3
but that's wrong, the solution is (0.6, 0)!
*/
for(int i=0; i<2; i++) 
{
    std::cout<< x[i]<<std::endl;
}

私は次の問題を抱えています：

dgesvが要素{3、3}を使用してベクトルxを計算するのはどうしてですか？解は{0.6、0}である必要があります（matlabで確認）。

ご挨拶

編集：dgesvは正方行列に対して機能する可能性があります。以下の私の解決策は、dgelsを使用して過剰決定系を解決する方法を示しています。

4

2 に答える 2

1

問題は、行優先/列優先の混同だけでなく、明らかに dgesv が非正方行列には適していないという事実でもありました。

過剰決定された線形方程式系は、dgels (#include "lapacke.h") を使用して、最小二乗法で解くことができます。

注目に値するのは、一見しただけではわからないことです。解ベクトル (通常は x で示されます) が b に格納されます。

私の例のシステムは、最初の列に値 1 ～ 10 を保持する行列 A と、{i|1<=i<=10} の値 i^2 を持つベクトル b で構成されます。

    double a[10][2]={1,1,2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1,10,1};
double b[10][1]={1,4,9,16,25,36,49,64,81,100};
lapack_int info,m,n,lda,ldb,nrhs;
int i,j;

// description here: http://www.netlib.org/lapack/double/dgesv.f
m = 10;
n = 2;
nrhs = 1;
lda = 2;
ldb = 1;

for(int i=0; i<m; i++)
{
    for(int k=0; k<lda; k++)
    {
        std::cout << a[i][k]<<" ";
    }
    std::cout <<"" << std::endl;
}
std::cout<< "printing vector b:"<<std::endl;
for(int i=0; i<m; i++) 
{
    std::cout<< *b[i]<<std::endl;
}
std::cout<< "\nStarting calculation..."<<std::endl;

info = LAPACKE_dgels(LAPACK_ROW_MAJOR,'N',m,n,nrhs,*a,lda,*b,ldb);

std::cout<< "Done."<<std::endl;

std::cout<< " "<<std::endl;
std::cout<< "printing vector b:"<<std::endl;
for(int i=0; i<2; i++) 
{
    std::cout<< *b[i]<<std::endl;
}
std::cout<< "Used values: "<< m << ", " <<  n << ", " << nrhs << ", " << lda << ", " << ldb << std::endl;
于 2012-10-10T11:25:38.477 に答える
1

Lapack は、行列が列優先形式で格納されていることを前提としています。行列 A のポインター ツー ポインター形式では機能しません。行列のすべての要素は、メモリに継続的に格納する必要があります。

A を次のように宣言してみてください。

double *A;
A = new double[2*n];

そして for ループ:

for(int i=0; i<n; i++)
{
    A[i+0*n]=5;
    A[i+1*n]=1;
    b[i]=3;
}
于 2012-10-09T10:28:10.390 に答える