39 lines
No EOL
918 B
C++
39 lines
No EOL
918 B
C++
#include "factorization_solve.h"
|
|
#include <vector>
|
|
#include <math.h>
|
|
#include <assert.h>
|
|
#include <lapack.h>
|
|
#include <iostream>
|
|
|
|
using namespace std;
|
|
|
|
|
|
void factorization_solve(vector<double> &A, vector<double> &b, const int32_t &n_rhs)
|
|
{
|
|
int32_t nelem = A.size();
|
|
int32_t N = sqrt(nelem);
|
|
assert(N == static_cast<int32_t>(b.size())/n_rhs);
|
|
|
|
vector<int> IPIV(N); // pivot indices
|
|
int info_factorization;
|
|
|
|
dgetrf_(&N, &N, A.data(), &N, IPIV.data(), &info_factorization);
|
|
|
|
const char transp = 'N';
|
|
int info_solve;
|
|
dgetrs_(&transp, &N, &n_rhs, A.data(), &N, IPIV.data(), b.data(), &N, &info_solve, 1); // 1 is length of parameter 'N'
|
|
|
|
}
|
|
|
|
void print_matrix(vector<double> &A, size_t M, size_t N)
|
|
{
|
|
for (size_t i = 0; i < M; ++i)
|
|
{
|
|
for (size_t j = 0; j < N; ++j)
|
|
{
|
|
cout << A[N*i + j] << "\t";
|
|
}
|
|
cout << endl;
|
|
}
|
|
cout << endl;
|
|
} |