a_mat_sq.cxx

Go to the documentation of this file.

/*
         Copyright 2010-17 Pierre SMARS (smars@yuntech.edu.tw)
         This program is free software: you can redistribute it and/or modify
         it under the terms of the GNU General Public License as published by
         the Free Software Foundation, either version 2 of the License, or
         (at your option) any later version.
 
         This program is distributed in the hope that it will be useful,
         but WITHOUT ANY WARRANTY; without even the implied warranty of
         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
         GNU General Public License for more details.
 
         You should have received a copy of the GNU General Public License
         along with this program.  If not, see <http://www.gnu.org/licenses/>.
         */
#include "a_mat_sq.h"
// .PORTABILITY : ansi C++
 
//***************************************************************************
//member functions
//-a_mat_sq-----------------------------------------------------------------
a_mat_sq::~a_mat_sq()
{
        //a_mat::~a_mat();
}
//-compatible----------------------------------------------------------------
void a_mat_sq::compatible(const a_mat& m) const
{
        if (m.maxi() != m.maxj()) throw compatibility_error();
        if (m.maxi() != maxi_)    throw compatibility_error();
}
//-operator=-----------------------------------------------------------------
a_mat_sq & a_mat_sq::operator=(const a_mat_sq& m)
{
        maxi_ = m.maxi();
        maxj_ = maxi_;
        (*x_) = m.x();
        return (*this);
}
//---------------------------------------------------------------------------
void a_mat_sq::settoI()
{
        a_mat_sq m(this->dim());
        for (UI i = 0; i < maxi_; i++)
                m(i,i) = 1.;
        (*this) = m;
}
//-operator=-----------------------------------------------------------------
a_mat_sq & a_mat_sq::operator=(const a_mat& m)
{
        this->compatible(m);
        (*x_) = m.x();
        return (*this);
}
//-operator*=----------------------------------------------------------------
a_mat_sq & a_mat_sq::operator*=(const a_mat_sq m)
{
        this->compatible(m);
        a_mat r(maxi_,maxj_);
        double v;
        for (UI i = 0; i < maxi_; i++)
        {
                for (UI j = 0; j < maxi_; j++)
                {
                        v = 0;
                        for (UI k = 0; k < maxi_; k++)
                                v += (*this)(i,k)*m(k,j);
                        r(i,j) = v;
                }
        }
        (*x_) = r.x();
        return (*this);
}
//-trace----------------------------------------------------------------------
double a_mat_sq::trace() const
{
        double v;
        for (int i = 0; i<maxi_; i++)
                v += (*this)(i,i);
        return v;
}
//-determinant----------------------------------------------------------------
double a_mat_sq::determinant() const
{
        if (maxi_==1)
                return (*this)(0,0);
        else if (maxi_==2)
                return (*this)(0,0)*(*this)(1,1)-(*this)(1,0)*(*this)(0,1);
        else if (maxi_==3)
        {
                double t1 = (*this)(0,0)*(*this)(1,1)*(*this)(2,2);
                double t2 = (*this)(0,1)*(*this)(1,2)*(*this)(2,0);
                double t3 = (*this)(0,2)*(*this)(1,0)*(*this)(2,1);
                double t4 = (*this)(0,0)*(*this)(1,2)*(*this)(2,1);
                double t5 = (*this)(0,1)*(*this)(1,0)*(*this)(2,2);
                double t6 = (*this)(0,2)*(*this)(1,1)*(*this)(2,0);
                return t1+t2+t3-t4-t5-t6;
        }
        else
        {
                a_mat_sq M(maxi_-1);
                double d = 0.;
                for (int i = 0; i<maxi_; i++)
                {
                        M = this->sub_matrix(0,i);
                        d += (*this)(0,i)*M.determinant()*(1.-(i%2)*2.); 
                }
                return d;
        }
}