|| || || || || || || ||

---

MatInverse.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatInverse_h
#define Impala_Core_Matrix_MatInverse_h

#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/MatSet.h"
#include <iostream>
#include <math.h>

namespace Impala
{
namespace Core
{
namespace Matrix
{


static const double Matrix_EPS = 1e-12;

template<class C> 
inline int
Ludcmp(C* a, int n, short *indx, double *d)
{
    double *vv = new double [n];

    *d=1.0;
    int i, j, k;
    for ( i=0;i<n;i++)
    {
        double big=0.0;
        for (j=0;j<n;j++)
        {
            double temp = fabs(a[i*n+j]);
            if (temp > big)
                big=temp;
        }
        if (big == 0.0)
            return(0); // Singular matrix 
        vv[i]=1.0/big;
    }
    for (j=0;j<n;j++)
    {
        for (i=0;i<j;i++)
        {
            double sum=a[i*n+j];
            for (k=0;k<i;k++)
                sum -= a[i*n+k]*a[k*n+j];
            a[i*n+j]=sum;
        }
        double big=0.0;
        int imax = 0;
        for (i=j;i<n;i++)
        {
            double sum=a[i*n+j];
            for (k=0;k<j;k++)
                sum -= a[i*n+k]*a[k*n+j];
            a[i*n+j]=sum;

            double dum = vv[i]*fabs(sum);
            if ( dum >= big)
            {
                big=dum;
                imax=i;
            }
        }
        if (j != imax)
        {
            for (k=0;k<n;k++)
            {
                double dum=a[imax*n+k];
                a[imax*n+k]=a[j*n+k];
                a[j*n+k]=dum;
            }
            *d = -(*d);
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (a[j*n+j] == 0.0)
            a[j*n+j]=Matrix_EPS;
        if (j != n-1)
        {
            double dum=1.0/a[j*n+j];
            for (i=j+1;i<n;i++)
                a[i*n+j] *= dum;
        }
    }
    delete [] vv;

    return(1);
}

template<class C> 
inline void
Lubksb(C* a, int n, short *indx, double *b)
{
    int ii = -1;

    int i;
    for (i=0;i<n;i++)
    {
        int ip=indx[i];
        double sum=b[ip];
        b[ip]=b[i];
        if (ii>=0)
            for (int j=ii;j<=i-1;j++)
                sum -= a[i*n+j]*b[j];
        else if (sum) ii=i;
        b[i]=sum;
    }
    for (i=n-1;i>=0;i--)
    {
        double sum=b[i];
        for (int j=i+1;j<n;j++)
            sum -= a[i*n+j]*b[j];
        b[i]=sum/a[i*n+i];
    }
}

template<class ArrayT>
inline ArrayT*
MatInverse(ArrayT* a)
{
    if (MatNrCol(a) != MatNrRow(a))
    {
        std::cerr << "MatInverse: matrix is not square!" << std::endl;
        return 0;
    }

    int size = MatNrRow(a);
    ArrayT* m = MatCreate<ArrayT>(size, size);
    ArrayT* tmp = MatCreate<ArrayT>(MatNrRow(a), MatNrCol(a));
    MatSet(tmp, a);
    short* idx = new short[size];
    double d;
    if (!Ludcmp(MatE(tmp, 0, 0), size, idx, &d))
    {
        std::cerr << "MatInverse: singular matrix can't be inverted." << std::endl;
        delete [] idx;
        return 0;
    }
    
    double* res = new double[size];
    for (int j=0 ; j<size ; j++)
    {
        for (int i=0 ; i<size; i++)
            res[i] = 0.0;
        res[j] = 1.0;
        Lubksb(MatE(tmp, 0, 0), size, idx, res);
        for (int k=0 ; k<size ; k++)
            *MatE(m, k, j) = res[k];
    }

    delete [] res;
    delete [] idx;
    delete tmp;

    return m;
}

} // namespace Matrix
} // namespace Core
} // namespace Impala

#endif

---