73 using namespace boost::numeric::ublas;
75 const size_t size = m.size1();
80 if (size != m.size2() || size == 0)
105 matrix<T>
A(size, 2*size);
106 matrix_range<matrix<T> > Aleft(
A,
110 matrix_range<matrix<T> > Aright(
A,
112 range(size, 2*size));
113 Aright = identity_matrix<T>(size);
116 for (
size_t k = 0;
k < size;
k++)
119 for (
size_t kk = 0; kk < size; kk++)
125 for (
size_t i = kk+1; i < size; i++)
137 std::cerr <<
"Error:" << __FUNCTION__ <<
":"
138 <<
"Input matrix is singular, because cannot find"
139 <<
" a row to swap while eliminating zero-diagonal.";
145 matrix_row<matrix<T> > rowk(
A, kk);
146 matrix_row<matrix<T> > rowl(
A, l);
160 for (
size_t j =
k+1; j < 2*size; j++)
165 for (
size_t i = 0; i < size; i++)
171 for (
size_t j =
k+1; j < 2*size; j++)
172 A(i,j) -=
A(
k,j) *
A(i,
k);