ProteoWizard
Functions
MatrixInverse.hpp File Reference
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix_expression.hpp>
#include <iostream>
#include <fstream>
#include <vector>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>

Go to the source code of this file.

Functions

template<typename M >
bool InvertMatrix (const M &input, M &inverse)
 
template<class T >
boost::numeric::ublas::matrix< T > gjinverse (const boost::numeric::ublas::matrix< T > &m, bool &singular)
 Invert a matrix via gauss-jordan algorithm (PARTIAL PIVOT)
 

Function Documentation

◆ InvertMatrix()

template<typename M >
bool InvertMatrix ( const M &  input,
M &  inverse 
)

Definition at line 32 of file MatrixInverse.hpp.

34{
35 using namespace boost::numeric::ublas;
36
37 typedef permutation_matrix<std::size_t> pmatrix;
38
39 // create a working copy of the input
40 M A(input);
41 // create a permutation matrix for the LU-factorization
42 pmatrix pm(A.size1());
43
44 // perform LU-factorization
45 int res = lu_factorize(A,pm);
46 if( res != 0 ) return false;
47
48 // create identity matrix of "inverse"
49 inverse.assign(identity_matrix<typename M::value_type>(A.size1()));
50
51 // backsubstitute to get the inverse
52 lu_substitute(A, pm, inverse);
53
54 return true;
55}
#define A

References A.

◆ gjinverse()

template<class T >
boost::numeric::ublas::matrix< T > gjinverse ( const boost::numeric::ublas::matrix< T > &  m,
bool &  singular 
)

Invert a matrix via gauss-jordan algorithm (PARTIAL PIVOT)

Parameters
mThe matrix to invert. Must be square.
singularIf the matrix was found to be singular, then this is set to true, else set to false.
Returns
If singular is false, then the inverted matrix is returned. Otherwise it contains random values.

Definition at line 70 of file MatrixInverse.hpp.

72{
73 using namespace boost::numeric::ublas;
74
75 const size_t size = m.size1();
76
77 // Cannot invert if non-square matrix or 0x0 matrix.
78 // Report it as singular in these cases, and return
79 // a 0x0 matrix.
80 if (size != m.size2() || size == 0)
81 {
82 singular = true;
83 matrix<T> A(0,0);
84 return A;
85 }
86
87 // Handle 1x1 matrix edge case as general purpose
88 // inverter below requires 2x2 to function properly.
89 if (size == 1)
90 {
91 matrix<T> A(1, 1);
92 if (m(0,0) == 0.0)
93 {
94 singular = true;
95 return A;
96 }
97 singular = false;
98 A(0,0) = 1/m(0,0);
99 return A;
100 }
101
102 // Create an augmented matrix A to invert. Assign the
103 // matrix to be inverted to the left hand side and an
104 // identity matrix to the right hand side.
105 matrix<T> A(size, 2*size);
106 matrix_range<matrix<T> > Aleft(A,
107 range(0, size),
108 range(0, size));
109 Aleft = m;
110 matrix_range<matrix<T> > Aright(A,
111 range(0, size),
112 range(size, 2*size));
113 Aright = identity_matrix<T>(size);
114
115 // Doing partial pivot
116 for (size_t k = 0; k < size; k++)
117 {
118 // Swap rows to eliminate zero diagonal elements.
119 for (size_t kk = 0; kk < size; kk++)
120 {
121 if ( A(kk,kk) == 0 ) // XXX: test for "small" instead
122 {
123 // Find a row(l) to swap with row(k)
124 int l = -1;
125 for (size_t i = kk+1; i < size; i++)
126 {
127 if ( A(i,kk) != 0 )
128 {
129 l = i;
130 break;
131 }
132 }
133
134 // Swap the rows if found
135 if ( l < 0 )
136 {
137 std::cerr << "Error:" << __FUNCTION__ << ":"
138 << "Input matrix is singular, because cannot find"
139 << " a row to swap while eliminating zero-diagonal.";
140 singular = true;
141 return Aleft;
142 }
143 else
144 {
145 matrix_row<matrix<T> > rowk(A, kk);
146 matrix_row<matrix<T> > rowl(A, l);
147 rowk.swap(rowl);
148
149/*#if defined(DEBUG) || !defined(NDEBUG)
150 std::cerr << __FUNCTION__ << ":"
151 << "Swapped row " << kk << " with row " << l
152 << ":" << A << "\n";
153#endif*/
154 }
155 }
156 }
157
158 /////////////////////////////////////////////////////////////////////////////////////////////////////////
159 // normalize the current row
160 for (size_t j = k+1; j < 2*size; j++)
161 A(k,j) /= A(k,k);
162 A(k,k) = 1;
163
164 // normalize other rows
165 for (size_t i = 0; i < size; i++)
166 {
167 if ( i != k ) // other rows // FIX: PROBLEM HERE
168 {
169 if ( A(i,k) != 0 )
170 {
171 for (size_t j = k+1; j < 2*size; j++)
172 A(i,j) -= A(k,j) * A(i,k);
173 A(i,k) = 0;
174 }
175 }
176 }
177
178/*#if defined(DEBUG) || !defined(NDEBUG)
179 std::cerr << __FUNCTION__ << ":"
180 << "GJ row " << k << " : " << A << "\n";
181#endif*/
182 }
183
184 singular = false;
185 return Aright;
186}
Kernel k

References A, and k.