I am trying to use LAPACK & armadillo libraries in my c++ code My problem is that the solution I get from LAPACK does not match the one I get using MATLAB.
I thought I must have messed up some argument in launching LAPACK, but I have tried all possible combinations and none of them match the MATLAB result.
C++ code (with the different tries of LAPACK)
#include <armadillo>
#include <lapacke.h>
int main()
{
const arma::mat AA = {
{-2.2554e+01, 4.0744e-03, 5.5231e-12, 5.0063e-01, 5.5051e-03, 1.4626e-08, 1.3794e-10, 7.0855e-01, 3.2284e-03, 6.7785e-01, 2.0654e+01},
{0, -5.1240e+02, 1.0125e-08, 6.6856e-08, 2.0260e-03, 5.7915e-03, 4.1114e-01, 1.6642e+02, 3.3348e-11, 3.4556e+02, 1.8098e-10},
{0, 0, -5.2337e+03, 1.9156e-12, 1.4217e-06, 3.4750e-06, 5.1531e+03, 5.1900e-08, 3.1478e+01, 1.8905e+01, 3.0223e+01},
{0, 0, 0, -1.4089e+03, 1.8106e-09, 1.6573e-03, 7.0136e-11, 1.4084e+03, 3.7770e-04, 2.9686e-11, 8.0124e-09},
{0, 0, 0, 0, -4.2774e+02, 8.1569e-02, 1.9177e-02, 1.5263e-11, 3.6179e+01, 1.0643e+02, 2.8503e+02},
{0, 0, 0, 0, 0, -1.7360e+01, 4.3968e-02, 5.4283e-07, 7.2528e+00, 9.9706e+00, 3.6747e-03},
{0, 0, 0, 0, 0, 0, -5.1551e+03, 1.8164e-09, 4.8217e-09, 1.6071e-12, 1.5191e+00},
{0, 0, 0, 0, 0, 0, 0, -1.5836e+03, 7.1831e+00, 3.4879e-09, 8.2123e-01},
{0, 0, 0, 0, 0, 0, 0, 0, -8.2139e+01, 9.9345e-09, 4.2855e-02},
{0, 0, 0, 0, 0, 0, 0, 0, 0, -4.8154e+02, 2.3236e-08},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3.3829e+02}};
const arma::vec bb = {-9.8760e-14, 5.0000e-02, 2.3451e-04, 3.9929e-04, 6.9937e-05, 4.4099e-04, 4.5265e-05, 2.2968e-03, 1.4573e-03, 2.1959e-03, 2.7900e-03};
arma::mat A1 = AA;
arma::vec b1 = bb;
int lda = 11;
int ldb = 11;
int ipiv[11];
int info = LAPACKE_dsysv(LAPACK_COL_MAJOR, 'L', 11, 1, A1.memptr(), lda, ipiv, b1.memptr(), ldb);
std::cout << "Result1 " << info << std::endl;
b1.print(std::cout);
arma::mat A2 = AA;
arma::vec b2 = bb;
lda = 11;
ldb = 11;
ipiv[11];
info = LAPACKE_dsysv(LAPACK_ROW_MAJOR, 'L', 11, 1, A2.memptr(), lda, ipiv, b2.memptr(), ldb);
std::cout << "Result2 " << info << std::endl;
b2.print(std::cout);
arma::mat A3 = AA;
arma::vec b3 = bb;
lda = 11;
ldb = 11;
ipiv[11];
info = LAPACKE_dsysv(LAPACK_COL_MAJOR, 'U', 11, 1, A3.memptr(), lda, ipiv, b3.memptr(), ldb);
std::cout << "Result3 " << info << std::endl;
b3.print(std::cout);
arma::mat A4 = AA;
arma::vec b4 = bb;
lda = 11;
ldb = 11;
ipiv[11];
info = LAPACKE_dsysv(LAPACK_ROW_MAJOR, 'U', 11, 1, A4.memptr(), lda, ipiv, b4.memptr(), ldb);
std::cout << "Result4 " << info << std::endl;
b4.print(std::cout);
return 0;
}
Matlab code
A= [ -2.2554e+01, 4.0744e-03, 5.5231e-12, 5.0063e-01, 5.5051e-03, 1.4626e-08, 1.3794e-10, 7.0855e-01, 3.2284e-03, 6.7785e-01, 2.0654e+01;
0, -5.1240e+02, 1.0125e-08, 6.6856e-08, 2.0260e-03, 5.7915e-03, 4.1114e-01, 1.6642e+02, 3.3348e-11, 3.4556e+02, 1.8098e-10;
0, 0, -5.2337e+03, 1.9156e-12, 1.4217e-06, 3.4750e-06, 5.1531e+03, 5.1900e-08, 3.1478e+01, 1.8905e+01, 3.0223e+01;
0, 0, 0, -1.4089e+03, 1.8106e-09, 1.6573e-03, 7.0136e-11, 1.4084e+03, 3.7770e-04, 2.9686e-11, 8.0124e-09;
0, 0, 0, 0, -4.2774e+02, 8.1569e-02, 1.9177e-02, 1.5263e-11, 3.6179e+01, 1.0643e+02, 2.8503e+02;
0, 0, 0, 0, 0, -1.7360e+01, 4.3968e-02, 5.4283e-07, 7.2528e+00, 9.9706e+00, 3.6747e-03;
0, 0, 0, 0, 0, 0, -5.1551e+03, 1.8164e-09, 4.8217e-09, 1.6071e-12, 1.5191e+00;
0, 0, 0, 0, 0, 0, 0, -1.5836e+03, 7.1831e+00, 3.4879e-09, 8.2123e-01;
0, 0, 0, 0, 0, 0, 0, 0, -8.2139e+01, 9.9345e-09, 4.2855e-02;
0, 0, 0, 0, 0, 0, 0, 0, 0, -4.8154e+02, 2.3236e-08;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3.3829e+02;]
A=A+A'
b=[ -9.8760e-14,
5.0000e-02,
2.3451e-04,
3.9929e-04,
6.9937e-05,
4.4099e-04,
4.5265e-05,
2.2968e-03,
1.4573e-03,
2.1959e-03,
2.7900e-03]
A\b
And the outputs of matlab (which I assume is right because it's harder messed the computation up) is:
ans =
-3.5245e-06
-5.7627e-05
-1.6269e-07
-2.6132e-06
-5.7499e-06
-2.1939e-05
-8.9090e-08
-4.9435e-06
-1.1355e-05
-2.3825e-05
-6.6679e-06
While the result from c++ is:
Result1 0 4.3788e-15
-9.7580e-05
-4.4808e-08
-2.8341e-07
-1.6350e-07
-2.5403e-05
-8.7806e-09
-1.4504e-06
-1.7742e-05
-4.5602e-06
-8.2474e-06
Result2 0
-2.8856e-02
5.0000e-02
2.3451e-04
3.9929e-04
6.9937e-05
4.4099e-04
4.5265e-05
2.2968e-03
1.4573e-03
2.1959e-03
2.7900e-03
Result3 0
-0.6278
-0.6280
-0.6277
-0.6277
-0.6278
-0.6277
-0.6277
-0.6277
-0.6274
-0.6279
-0.6278
Result4 0
4.3788e-15
5.0000e-02
2.3451e-04
3.9929e-04
6.9937e-05
4.4099e-04
4.5265e-05
2.2968e-03
-2.8858e-02
2.1959e-03
2.7900e-03
What am I getting wrong?
I checked again, and can get correct results making AA explicitly symmetrical
AA = AA + AA.t();
But in from the documentation, it seemed that the second argument is used to access just the upper or lower triangle of the matrix. So this should not be needed.
Documentation:
UPLO (input) = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.
N (input) The number of linear equations, i.e., the order of the matrix A. N >= 0.
NRHS (input) The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.
A (input/output) On entry, the symmetric matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A con- tains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower tri- angular part of the matrix A, and the strictly upper triangular part of A is not referenced.
