28 #ifndef __TERRALIB_COMMON_INTERNAL_MATRIXUTILS_H
29 #define __TERRALIB_COMMON_INTERNAL_MATRIXUTILS_H
38 #include <boost/numeric/ublas/io.hpp>
39 #include <boost/numeric/ublas/lu.hpp>
40 #include <boost/numeric/ublas/matrix.hpp>
41 #include <boost/numeric/ublas/triangular.hpp>
42 #include <boost/numeric/ublas/vector.hpp>
43 #include <boost/numeric/ublas/vector_proxy.hpp>
60 if( ( inputMatrix.size1() == 0 ) || ( inputMatrix.size2() == 0 ) )
67 boost::numeric::ublas::matrix<T>
A( inputMatrix );
69 const unsigned int size1 = (
unsigned int)A.size1();
71 boost::numeric::ublas::permutation_matrix<std::size_t> pm( size1 );
73 if( boost::numeric::ublas::lu_factorize( A, pm ) != 0.0 )
80 for (
unsigned int pmi = 0; pmi < size1; ++pmi)
81 if ( pmi != pm( pmi ) )
86 for(
unsigned int i = 0 ; i < size1 ; i++ )
87 determinant *=
A(i,i);
89 determinant = determinant * pmSign;
105 boost::numeric::ublas::matrix<T>& outputMatrix)
107 assert( inputMatrix.size1() == inputMatrix.size2() );
110 boost::numeric::ublas::matrix<T>
A( inputMatrix );
113 boost::numeric::ublas::permutation_matrix<std::size_t> pm( A.size1() );
116 if( boost::numeric::ublas::lu_factorize( A, pm ) != 0 )
123 outputMatrix = boost::numeric::ublas::identity_matrix<T>( A.size1() );
128 boost::numeric::ublas::lu_substitute( A, pm, outputMatrix );
153 boost::numeric::ublas::matrix<T>& outputMatrix)
155 if( inputMatrix.size1() > inputMatrix.size2() )
157 boost::numeric::ublas::matrix<T> trans( boost::numeric::ublas::trans(
160 boost::numeric::ublas::matrix<T> aux1( boost::numeric::ublas::prod( trans,
163 boost::numeric::ublas::matrix<T> aux1Inv;
167 outputMatrix = boost::numeric::ublas::prod( aux1Inv, trans );
175 else if( inputMatrix.size1() < inputMatrix.size2() )
177 boost::numeric::ublas::matrix<T> trans( boost::numeric::ublas::trans(
180 boost::numeric::ublas::matrix<T> aux1( boost::numeric::ublas::prod(
181 inputMatrix, trans ) );
183 boost::numeric::ublas::matrix<T> aux1Inv;
187 outputMatrix = boost::numeric::ublas::prod( trans, aux1Inv );
211 bool EigenVectors(
const boost::numeric::ublas::matrix<T>& inputMatrix, boost::numeric::ublas::matrix<T> &eigenVectorsMatrix, boost::numeric::ublas::matrix<T> &eigenValuesMatrix)
218 int dim = inputMatrix.size1();
227 double range, fdim, anorm, anrmx,
229 sinx, cosx, sinx2, cosx2,
235 eigenVectorsMatrix.resize(dim, dim);
236 eigenVectorsMatrix.clear();
237 for (i = 0; i < dim; i++)
238 eigenVectorsMatrix(i, i) = 1.0;
240 int fat = (dim * dim + dim) / 2;
244 cov =
new double[fat];
245 e_vec =
new double[dim * dim];
246 e_val =
new double[fat];
248 if( cov == NULL || e_vec == NULL || e_val == NULL )
252 for (i = 0; i < dim; i++)
253 for (j = 0; j <= i; j++)
254 cov[k++] = inputMatrix(i, j);
256 for (i = 0; i < ((dim * dim + dim) / 2); i++)
260 for (i = 0; i < dim; i++)
263 for (j = 0; j < dim; j++)
275 for (j = 0; j < dim; j++)
277 for (i = 0; i <= j; i++)
280 ia = i + (j * j + j) / 2;
281 anorm = anorm + e_val[ia] * e_val[ia];
287 anorm = 1.414 * sqrt((
double)anorm);
288 anrmx = anorm * range / fdim;
311 mq = (m * m + m) / 2;
312 lq = (l * l + l) / 2;
315 if (fabs((
double)(e_val[lm])) >= thr)
321 x = 0.5 * (e_val[ll] - e_val[mm]);
322 z = e_val[lm] * e_val[lm] + x * x;
323 y = - e_val[lm] / sqrt((
double)(z));
328 z = sqrt( (
double)(1.0 - y * y) );
329 sinx = y / sqrt( (
double)(2.0 * (1.0 + z)) );
332 cosx = sqrt( (
double)(1.0 - sinx2) );
340 for (i = 0; i < dim; i++)
342 iq = (i * i + i) / 2;
357 x = e_val[il] * cosx - e_val[im] * sinx;
358 e_val[im] = e_val[il] * sinx + e_val[im] * cosx;
366 x = e_vec[ilr] * cosx - e_vec[imr] * sinx;
367 e_vec[imr] = e_vec[ilr] * sinx + e_vec[imr] * cosx;
371 x = 2.0 * e_val[lm] * sincs;
372 y = e_val[ll] * cosx2 + e_val[mm] * sinx2 - x;
373 x = e_val[ll] * sinx2 + e_val[mm] * cosx2 + x;
375 e_val[lm] = (e_val[ll] - e_val[mm]) * sincs + e_val[lm] * (cosx2 - sinx2);
405 for (i = 0; i < dim; i++)
408 ll = i + (i * i + i) / 2;
411 for (j = i; j < dim; j++)
414 mm = j + (j * j + j) / 2;
416 if (e_val[ll] < e_val[mm])
419 e_val[ll] = e_val[mm];
422 for (k = 0; k < dim; k++)
427 e_vec[ilr] = e_vec[imr];
435 eigenValuesMatrix.resize(fat, 1);
436 eigenValuesMatrix.clear();
438 for (i = 0; i < dim; i++)
439 for (j = 0; j < dim; j++)
440 eigenVectorsMatrix(j, i) = e_vec[k++];
441 for (i = 0; i < fat; i++)
442 eigenValuesMatrix(i, 0) = e_val[i];
454 #endif //__TERRALIB_COMMON_INTERNAL_MATRIXUTILS_H
bool GetDeterminant(const boost::numeric::ublas::matrix< T > &inputMatrix, T &determinant)
Get the Matrix determinant value.
Configuration flags for the TerraLib Common Runtime module.
bool GetPseudoInverseMatrix(const boost::numeric::ublas::matrix< T > &inputMatrix, boost::numeric::ublas::matrix< T > &outputMatrix)
Pseudo matrix inversion.
bool EigenVectors(const boost::numeric::ublas::matrix< T > &inputMatrix, boost::numeric::ublas::matrix< T > &eigenVectorsMatrix, boost::numeric::ublas::matrix< T > &eigenValuesMatrix)
Computes the eigenvectors of a given matrix.
bool GetInverseMatrix(const boost::numeric::ublas::matrix< T > &inputMatrix, boost::numeric::ublas::matrix< T > &outputMatrix)
Matrix inversion.