LAPACK dgetrs vs dgesv

В документации LAPACK,

dgesv Решает общую систему линейных уравнений AX = B.

dgetrs Решает общую систему линейных уравнений AX = B, AT X = B или AH X = B, используя факторизацию LU, вычисленную DGETRF.

В чем разница ? Кроме того, я сделал следующую программу на C / C ++, и оба дают разные результаты. Почему это так ?

#include <iostream>
#include <iomanip>
#include <vector>
#include <math.h>
#include <time.h>
#include "stdafx.h"
using namespace std;

extern "C" {
// LU decomoposition of a general matrix
void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
//// generate inverse of a matrix given its LU decomposition
//void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int*    lwork, int* INFO);
void dgetrs_(char* C, int* N, int* NRHS, double* A, int* LDA, int* IPIV, double* B, int* LDB, int* INFO);
void dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
}

void solvelineq(double* A, double* B, int N)
{
int *IPIV = new int[N + 1];
int LWORK = N*N;
double *WORK = new double[LWORK];
int INFO;
char C = 'N';
int NRHS = 1;
dgetrf_(&N, &N, A, &N, IPIV, &INFO);
dgetrs_(&C, &N, &NRHS, A, &N, IPIV, B, &N, &INFO);

//dgetri_(&N, A, &N, IPIV, WORK, &LWORK, &INFO);

delete IPIV;
delete WORK;
}

double comparematrices(double* A, double* B, int N)
{
double sum = 0.;
for (int i = 0; i < N; ++i)
sum += fabs(A[i] - B[i]);
return sum;
}

int main() {
int dim;
std::cout << "Enter Dimensions of Matrix : \n";
std::cin >> dim;
clock_t c1, c2;
c1 = clock();

vector<double> a(dim*dim);
vector<double> b(dim);
vector<double> c(dim);
vector<int> ipiv(dim);
srand(1);              // seed the random # generator with a known value
double maxr = (double)RAND_MAX;
for (int r = 0; r < dim; r++) {  // set a to a random matrix, i to the identity
for (int c = 0; c < dim; c++) {
a[r + c*dim] = rand() / maxr;
}
b[r] = rand() / maxr;
c[r] = b[r];
}
c2 = clock();
cout << endl << dim << " x " << dim << " matrix generated in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
// C is identical matrix to B because we are solving by two methods.

c1 = clock();
solvelineq(&*a.begin(), &*c.begin(), dim);
c2 = clock();
cout << endl << " dgetrf_ and dgetrs_ completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
vector<double> a1(a);
vector<double> b1(b);
int info;
int one = 1;
c1 = clock();
dgesv_(&dim, &one, &*a.begin(), &dim, &*ipiv.begin(), &*b.begin(), &dim, &info);
c2 = clock();
cout << endl << " dgesv_ completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "info is " << info << endl;
double eps = 0.;
c1 = clock();
for (int i = 0; i < dim; ++i)
{
double sum = 0.;
for (int j = 0; j < dim; ++j)
sum += a1[i + j*dim] * b[j];
eps += fabs(b1[i] - sum);
}
c2 = clock();
cout << endl << " dgesv_ check completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "check is " << eps << endl;
cout << "\nFinal Matrix 1 : " << endl;
for (int i = 0; i < dim; ++i)
cout << b[i] << endl;

cout << "\nFinal Matrix 2 : " << endl;
for (int i = 0; i < dim; ++i)
cout << c[i] << endl;
double diff;
c1 = clock();
diff = comparematrices(&*b.begin(), &*c.begin(), dim);
c2 = clock();
cout << endl << " Both functions compared in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "Sum of difference is : " << diff << endl;
return 0;
}

Мой результат:
Введите размеры матрицы:
5

Матрица 5 x 5, сгенерированная за: 0 секунд

dgetrf_ и dgetrs_ выполнено за: 0,001 секунды

dgesv_ выполнено за: 0 секунд
информация 0

dgesv_ проверка завершена: 0 секунд
чек 2.02009e-15

Финальная матрица 1:
-2,97629
4,83948
2,00769
-1,98887
-5,15561

Финальная матрица 2:
-1,40622
2,29029
0.480829
-1,63597
0,71962

Обе функции сравниваются за: 0 секунд

Сумма разницы: 11,8743