19#include <cusp/array1d.h>
20#include <cusp/array2d.h>
21#include <cusp/linear_operator.h>
30template <
typename IndexType,
typename ValueType,
typename MemorySpace,
typename OrientationA>
32 cusp::array1d<IndexType,MemorySpace>& pivot)
34 const int n = A.num_rows;
37 for (
int k = 0; k < n; k++)
41 ValueType max = std::fabs(A(k,k));
43 for (
int j = k + 1;
j < n;
j++)
45 if (max < std::fabs(A(
j,k)))
47 max = std::fabs(A(
j,k));
55 for (
int j = 0;
j < n;
j++)
56 std::swap(A(k,
j), A(pivot[k],
j));
63 for (
int i = k + 1; i < n; i++)
67 for (
int i = k + 1; i < n; i++)
68 for (
int j = k + 1;
j < n;
j++)
69 A(i,
j) -= A(i,k) * A(k,
j);
76template <
typename IndexType,
typename ValueType,
typename MemorySpace,
77 typename OrientationA,
typename OrientationB>
78int block_lu_solve(
const cusp::array2d<ValueType,MemorySpace,OrientationA>& A,
79 const cusp::array1d<IndexType,MemorySpace>& pivot,
80 const cusp::array2d<ValueType,MemorySpace,OrientationB>& b,
81 cusp::array2d<ValueType,MemorySpace,OrientationB>& x,
84 const int n = A.num_rows;
85 const int numRHS = b.num_cols;
90 for (
int k = 0; k < n; k++)
93 for (
int j = 0;
j < numRHS;
j++)
94 std::swap(x(k,
j),x(pivot[k],
j));
97 for (
int i = 0; i < k; i++){
98 for (
int j = 0;
j< numRHS;
j++)
99 x(k,
j) -= A(k,i) * x(i,
j);
105 for (
int k = n - 1; k >= 0; k--)
107 for (
int j = 0;
j< numRHS;
j++){
108 for (
int i = k + 1; i < n; i++){
109 x(k,
j) -= A(k,i) * x(i,
j);
123template <
typename ValueType,
typename MemorySpace>
126 cusp::array2d<ValueType,cusp::host_memory>
lu;
127 cusp::array1d<int,cusp::host_memory>
pivot;
132 template <
typename MatrixType>
134 linear_operator<ValueType,MemorySpace>(A.num_rows, A.num_cols, A.num_entries)
136 CUSP_PROFILE_SCOPED();
139 pivot.resize(A.num_rows);
143 template <
typename VectorType1,
typename VectorType2>
146 CUSP_PROFILE_SCOPED();
block_lu_solver(const MatrixType &A)
void operator()(const VectorType1 &b, VectorType2 &x) const
cusp::array1d< int, cusp::host_memory > pivot
cusp::array2d< ValueType, cusp::host_memory > lu
int block_lu_factor(cusp::array2d< ValueType, MemorySpace, OrientationA > &A, cusp::array1d< IndexType, MemorySpace > &pivot)
int block_lu_solve(const cusp::array2d< ValueType, MemorySpace, OrientationA > &A, const cusp::array1d< IndexType, MemorySpace > &pivot, const cusp::array2d< ValueType, MemorySpace, OrientationB > &b, cusp::array2d< ValueType, MemorySpace, OrientationB > &x, cusp::array2d_format)