solveur/src/matpleine.cpp

#include "gestionversion.h"

//---------------------------------------------------------------------------

#pragma hdrstop
#include <iomanip>
#include "matpleine.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)

#include <math>


//------------------------------------
//   Implementation of Matrice Class
//------------------------------------


#include <fstream>
#include <cmath>


Matrice::Matrice()
{
        Data = NULL;
        NbCols = NbRows = 0;
}


Matrice::Matrice( int Rows, int Cols,double InitVal)
{
        NbRows = Rows;
        NbCols = Cols;

        Data = new double*[NbRows];
        for(int i=0; i<NbRows; ++i) {
                Data[i] = new double[NbCols];
                for(int j=0; j<NbCols; ++j) {
                        Data[i][j] = InitVal;
                }
        }
}


 Matrice::Matrice(double* m_pData, int Rows, int Cols)
{
        NbRows = Rows;
        NbCols = Cols;

        Data = new double*[NbRows];
        for(int i=0; i<NbRows; ++i) {
                Data[i] = new double[NbCols];
                for(int j=0; j<NbCols; ++j) {
                        Data[i][j] = m_pData[i*NbCols+j];
                }
        }
}

Matrice::Matrice(double** m_nData, int Rows, int Cols)
{
        NbRows = Rows;
        NbCols = Cols;

        Data = new double*[NbRows];
        for(int i=0; i<NbRows; ++i) {
                Data[i] = new double[NbCols];
                for(int j=0; j<NbCols; ++j) {
                        Data[i][j] = m_nData[i][j];
                }
        }
}


Matrice::Matrice(const Matrice& obj)
{
        NbRows = obj.NbRows;
        NbCols = obj.NbCols;

        Data = new double*[NbRows];
        for(int i=0; i<NbRows; ++i) {
                Data[i] = new double[NbCols];
                for(int j=0; j<NbCols; ++j) {
                        Data[i][j] = obj.Data[i][j];
                }
        }
}


Matrice::~Matrice()
{
if(Data){
        for(int i=0; i<NbRows; ++i) {
                delete[] Data[i];
                Data[i] = NULL;
        }
        delete[] Data;
        }
        Data = NULL;
}

//assignment operator
Matrice& Matrice::operator =(const Matrice& obj)
{
        if(&obj == this) return *this;

        if(NbRows || NbCols) {
                for(int i=0; i<NbRows; ++i) {
                        delete[] Data[i];
                }
                delete[] Data;
        }

        NbRows = obj.NbRows;
        NbCols = obj.NbCols;

        Data = new double*[NbRows];
        for(int i=0; i<NbRows; ++i) {
                Data[i] = new double[NbCols];
                for(int j=0; j<NbCols; ++j) {
                        Data[i][j] = obj.Data[i][j];
                }
        }

        return *this;
}

Matrice& Matrice::operator +(const Matrice& obj) const
{
        assert((NbRows == obj.NbRows)&& (NbCols == obj.NbCols)) ;


        static Matrice temp;
        temp = Matrice( NbRows, NbCols,0.0);

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[i][j] = Data[i][j] + obj.Data[i][j];
                }
        }

        return temp;
}


 Matrice& Matrice::operator -(const Matrice& obj) const
{
        assert((NbRows == obj.NbRows)&& (NbCols == obj.NbCols)) ;

        static Matrice temp;
        temp = Matrice(NbRows, NbCols,0.);

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[i][j] = Data[i][j] - obj.Data[i][j];
                }
        }

        return temp;
}


 Matrice& Matrice::operator *(const Matrice& obj) const
{
        assert(NbCols == obj.NbRows);

        double sum = double(0);
        double prod = double(1);

        static Matrice temp;
        temp = Matrice( NbRows, obj.NbCols,0.0);

        for(int i=0; i<temp.NbRows; ++i) {
                for(int j=0; j<temp.NbCols; ++j) {
                        sum = double(0);
                        for(int q=0; q<NbCols; ++q) {
                                prod = Data[i][q] * obj.Data[q][j];
                                sum += prod;
                        }
                        temp.Data[i][j] = sum;
                }
        }

        return temp;
}

//multiply a Matrice by a double scalar
Matrice& Matrice::operator *(const double _d) const
{
        static Matrice temp;
        temp = Matrice(Data, NbRows,NbCols);

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[i][j] *= _d;
                }
        }

return temp;
}


//divide a Matrice by a double scalar
 Matrice& Matrice::operator /(const double _d) const
{
        assert(_d != 0) ;

        static Matrice temp;
        temp = Matrice(Data, NbCols, NbRows);

        for(int i=0; i<temp.NbRows; ++i) {
                for(int j=0; j<temp.NbCols; ++j) {
                        temp.Data[i][j] /= _d;
                }
        }

        return temp;
}


//add two matrices and assign
Matrice& Matrice::operator +=(const Matrice& obj)
{
        return(*this = *this + obj);
}

//subtract two matrices and assign
 Matrice& Matrice::operator -=(const Matrice& obj)
{
        return(*this = *this - obj);
}

//multiply two matrices and assign
 Matrice& Matrice::operator *=(const Matrice& obj)
{
        return(*this = *this * obj);
}

//multiply a Matrice by a double scalar and assign
 Matrice& Matrice::operator *=(const double _d)
{
        return(*this = *this * _d);
}


//divide a Matrice by a double scalar and assign
 Matrice& Matrice::operator /=(const double _d)
{
        return(*this = *this / _d);
}


//equality operator
 bool Matrice::operator ==(const Matrice& obj) const
{
        if(NbRows != obj.NbRows) return false;
        if(NbCols != obj.NbCols) return false;

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        if(Data[i][j] != obj.Data[i][j]) return false;
                }
        }

        return true;
}

//inequality operator
 bool Matrice ::operator !=(const Matrice& obj) const
{
        if(NbRows != obj.NbRows) return true;
        if(NbCols != obj.NbCols) return true;

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        if(Data[i][j] != obj.Data[i][j]) return true;
                }
        }

        return false;
}

//renvoie la ieme ligne
double* Matrice ::operator [](const int _i) const
{
        assert((_i < NbRows)&& (_i >= 0));

        return Data[_i];
}


//another indexing operator, but takes both rows and columns
//remember, a zero (0) element DOES exist in this Matrice, although not in a real Matrice
double& Matrice ::operator ()(const int _i, const int _j) const
{
        assert((_i < NbRows) && (_j < NbCols));
        assert((_i >= 0) && (_j >= 0));

        return Data[_i][_j];
}

//returns true if the calling Matrice is an identity Matrice
 bool Matrice::IsIdentity() const
{
        if(NbCols != NbRows) return false;

        for(int i=0; i<NbCols; ++i) {
                for(int j=0; j<NbRows; ++j) {
                        if(i == j) {
                                if(Data[i][j] != 1) return false;
                        }
                        else if(Data[i][j] != 0) return false;
                }
        }

        return true;
}

//returns true if every element in the calling Matrice is zero (0)
 bool Matrice::IsZeroMatrice() const
{
        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        if(Data[i][j] != 0) return false;
                }
        }
        
        return true;
}


//renvoie la somme des tous les elemnets de la matrice
//sur les lignes et sur les colones

double Matrice::Get_SumAll() const
{
        double sum = 0;

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        sum += Data[i][j];
                }
        }

        return sum;
}


//Renvoie la somme de tous les termes d une ligne dans la matrice

double Matrice::Get_SumRow(const int Row) const
{
        double sum = double(0.);

        for(int j=0; j<NbCols; ++j) {
                sum += Data[Row][j];
        }

        return sum;
}

//Renvoie la somme de tous les ternes d<une colonnne dans la matrice

double Matrice::Get_SumColumn(const int Col) const
{
        double sum = double(0);

        for(int i=0; i<NbRows; ++i) {
                sum += Data[i][Col];
        }

        return sum;
}


//Renvoie le terme le plus grand dans la matrice

double Matrice::Get_Max() const
{
        double max = Data[0][0];

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        if(Data[i][j] > max) max = Data[i][j];
                }
        }

        return max;
}
// Renvoie le termes le plus grand en valeur absolue dans la matrice

double Matrice::Get_AbMax() const
{
        double max = fabs(Data[0][0]);

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        if(fabs(Data[i][j]) > max) max = fabs(Data[i][j]);
                }
        }

        return max;
}

//Renvoie le terme le plus petit dans la matrice

double Matrice ::Get_Min() const
{
        double min = Data[0][0];

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        if(Data[i][j] < min) min = Data[i][j];
                }
        }

        return min;
}

//returns the largest value in row 'Row' in the Matrice
double Matrice::Get_RowMax(const int Row) const
{
        double max = Data[Row][0];

        for(int j=0; j<NbCols; ++j) {
                if(Data[Row][j] > max) max = Data[Row][j];
        }

        return max;
}

//returns the smallest value in row 'Row' in the Matrice
double Matrice::Get_RowMin(const int Row) const
{
        double min = Data[Row][0];

        for(int j=0; j<NbCols; ++j) {
                if(Data[Row][j] < min) min = Data[Row][j];
        }

        return min;
}

//returns the largest value in column 'Col' in the Matrice
double Matrice::Get_ColumnMax(const int Col) const
{
        double max = Data[0][Col];

        for(int i=0; i<NbRows; ++i) {
                if(Data[i][Col] > max) max = Data[i][Col];
        }

        return max;
}

//returns the smallest value in column 'Col' in the Matrice
double Matrice::Get_ColumnMin(const int Col) const
{
        double min = Data[0][Col];

        for(int i=0; i<NbRows; ++i) {
                if(Data[i][Col] < min) min = Data[i][Col];
        }

        return min;
}


// Renvoie la matrice dans un tableau a une dimension;

double* Matrice ::Get_DataOneDimen() const
{
        double* newData;

        newData = new double[NbRows * NbCols];

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        newData[(i*NbRows)+j] = Data[i][j];
                }
        }

        return newData;
}

// Renvoie la matrice dans un tableu a deux dimension

double** Matrice ::Get_DataTwoDimen() const
{
        double** newData;

        newData = new double*[NbRows];
        for(int i=0; i<NbRows; ++i) {
                newData[i] = new double[NbCols];
                for(int j=0; j<NbCols; ++j) {
                        newData[i][j] = Data[i][j];
                }
        }

        return newData;
}

//returns number of rows of the Matrice
 int Matrice ::Get_NbRows() const
{
        return NbRows;
}

//returns number of columns of the Matrice
 int Matrice::Get_NbCols() const
{
        return NbCols;
}

//clears every entry in the calling Matrice
 Matrice& Matrice::Clear()
{
        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        Data[i][j] = 0;
                }
        }

        return *this;
}

//Mettre tous les termes d<une ligne ala valeur zero

Matrice& Matrice::ClearRow(const int Row)
{
        for(int j=0; j<NbCols; ++j) {
                Data[Row][j] = 0;
        }

        return *this;
}

//Metrre tous les termes d' une aolonne en zero

 Matrice& Matrice::ClearColumn(const int Col)
{
        for(int i=0; i<NbRows; ++i) {
                Data[i][Col] = 0;
        }

        return *this;
}

//Metrre tous les termes D<une matrice a la veleur '_d'
Matrice& Matrice::Fill(const double _d)
{
        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        Data[i][j] = _d;
                }
        }

        return *this;
}

//fills every entry in the row 'Row' of the calling Matrice to '_d'
 Matrice& Matrice::FillRow(const int Row, const double _d)
{
        for(int j=0; j<NbCols; ++j) {
                Data[Row][j] = _d;
        }

        return *this;
}

//fills every entry in the column 'Col' of the calling Matrice to '_d'
 Matrice& Matrice::FillColumn(const int Col, const double _d)
{
        for(int i=0; i<NbRows; ++i) {
                Data[i][Col] = _d;
        }

        return *this;
}


 Matrice& Matrice::Get_SubMatrice(const int RowSpot, const int ColSpot, const int RowLen, const int ColLen) const
{
        static Matrice temp;
        temp = Matrice( RowLen, ColLen,0.0);

        for(int i=RowSpot, k=0; i<(RowLen+RowSpot); i++, k++) {
                for(int j=ColSpot, l=0; j<(ColLen+ColSpot); j++, l++) {
                        temp(k,l) = Data[i][j];
                }
        }

        return temp;
}

//changes the calling Matrice into a subMatrice starting at spot (RowSpot, ColSpot)
//and with lengths of 'RowLen' rows and 'ColLen' columns
Matrice& Matrice::Set_SubMatrice(Matrice& sub,const int RowSpot, const int ColSpot, const int RowLen, const int ColLen)
{

          for(int i=RowSpot, k=0;i<RowSpot+RowLen;i++,++k)
            {
              for(int j=ColSpot, l=0;j< ColSpot+ColLen;j++,++l)
               Data[i][j]=sub(k,l);
            }

        return *this;
}

Matrice Matrice::Get_SubMatrice(vecteur& lgn,vecteur& col)
{
int t_tmp=lgn.getDim();
Matrice tmp(t_tmp,t_tmp,0.);
 for(int i=0;i<t_tmp;i++)
 {
 int I= lgn[i];
  for(int j=0;j<t_tmp;j++)
   {
    int J= col[j];
    tmp.Data[i][j] = Data[I][J];
   }
 }
return tmp;
}

//swaps two rows, Row1 and Row2, from the calling Matrice
 Matrice& Matrice::Get_SwapRows(const int Row1, const int Row2)
{
        double* temp;

        temp = new double[NbCols];

        for(int j=0; j<NbCols; ++j) {
                temp[j] = Data[Row1][j];
                Data[Row1][j] = Data[Row2][j];
                Data[Row2][j] = temp[j];
        }

        delete[] temp;

        return *this;
}

//swaps two columns, Col1 and Col2, from the calling Matrice
 Matrice& Matrice::Get_SwapCols(const int Col1, const int Col2)
{
        double* temp;

        temp = new double[NbRows];

        for(int i=0; i<NbRows; ++i) {
                temp[i] = Data[i][Col1];
                Data[i][Col1] = Data[i][Col2];
                Data[i][Col2] = temp[i];
        }

        delete[] temp;

        return *this;
}

//returns the transposition of the calling Matrice
 Matrice& Matrice::Get_Transpose() const
{
        static Matrice temp;
        temp = Matrice(NbCols, NbRows,0.0);

        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[j][i] = Data[i][j];
                }
        }

        return temp;
}

//transposes the calling Matrice
 Matrice& Matrice::Transpose()
{
        *this = this->Get_Transpose();

        return *this;
}


//CMAR = Concatenate Matrice As Rows
//concatenates Matrice 'obj' on to the right of the calling Matrice
 Matrice& Matrice::Get_ConMatAsRow(const Matrice& obj)
{
        assert(NbCols == obj.NbCols);// ErrorMsg("mismatched matrices in row concatenation", true);

        Matrice temp((NbRows + obj.NbRows), NbCols,0.0);
        int i;
        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[i][j] = Data[i][j];
                }
        }
        for(int k=0; i<temp.NbRows; ++i, ++k) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[i][j] = obj.Data[k][j];
                }
        }

        *this = temp;

        return *this;
}

//CMAC = Concatenate Matrice As Columns
//concatenates Matrice 'obj' on to the bottom of the calling Matrice
 Matrice& Matrice::Get_ConMatAsCol(const Matrice& obj)
{
        assert(NbRows == obj.NbRows);// ErrorMsg("mismatched matrices in column concatenation", true);

        Matrice temp(NbRows, (NbCols + obj.NbCols),0.0);
        int i,j;
        for(i=0; i<NbRows; ++i) {
                for( j=0; j<NbCols; ++j) {
                        temp.Data[i][j] = Data[i][j];
                }
                for(int l=0; l<obj.NbCols; ++l, ++j) {
                        temp.Data[i][j] = obj.Data[i][l];
                }
        }

        *this = temp;

        return *this;
}

//CMAR = Concatenate Matrice As Rows
//returns a new Matrice that is the calling object + 'obj' on the right
 Matrice& Matrice::GetCMAR(const Matrice& obj) const
{
        static Matrice temp;

        temp = *this;
        temp.Get_ConMatAsRow(obj);

        return temp;
}

//CMAC = Concatenate Matrice As Columns
//returns a new Matrice that is the valling object + 'obj' on the bottom
 Matrice& Matrice::GetCMAC(const Matrice& obj) const
{
        static Matrice temp;

        temp = *this;
        temp.Get_ConMatAsCol(obj);

        return temp;
}

//adds a row onto the right of the calling Matrice
 Matrice& Matrice::ConcatenateRow(const double* RowData)
{
        Matrice temp( NbRows+1, NbCols,0.0);
        int i;
        for(i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[i][j] = Data[i][j];
                }
        }
        int j;
        for(i=NbRows, j=0; j<NbCols; ++j) {
                temp.Data[i][j] = RowData[j];
        }

        *this = temp;

        return *this;
}

//adds a column onto the bottom of the calling Matrice
 Matrice& Matrice::ConcatenateColumn(const double* ColumnData)
{
        
        Matrice temp( NbRows, NbCols+1,0.0);
        int i;
        for(i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[i][j] = Data[i][j];
                }
        }
        int j;
        for(j=NbCols, i=0; i<NbRows; ++i) {
                temp.Data[i][j] = ColumnData[i];
        }

        *this = temp;

        return *this;
}

//adds a row into the Matrice in the spot 'RowSpot'
 Matrice& Matrice::Add_Row(const vecteur& RowData, const int RowSpot)
{
        int New_NbRows=NbRows;
        int New_NbCols=NbCols;

        Matrice temp( NbRows+1, NbCols,0.0);

        //if(NbRows==0)
        New_NbRows=NbRows+1;


        for(int i=0, k=0; i<New_NbRows; ++i, ++k) {
                if(i == RowSpot) {
                        for(int j=0; j<New_NbCols; ++j) {
                                temp.Data[i][j] += RowData[j];
                        }
                        ++k;
                }

                if(k<New_NbRows){
                   for(int j=0; j<New_NbCols; ++j)
                        temp.Data[k][j] += Data[i][j];

                }
        }

        *this = temp;

        return *this;
}

//Ajouter une Colonne dans la Matrice in the spot 'ColumnSpot'

Matrice& Matrice::Add_Column(const vecteur& ColumnData, const int ColumnSpot)
{
        int New_NbRows=NbRows;
        int New_NbCols=NbCols;

        Matrice temp(NbRows, NbCols+1,0.0);

       // if(NbCols==0)
         New_NbCols=NbCols+1;

        for(int i=0; i<NbRows; ++i) {
                for(int j=0, k=0; j<New_NbCols; ++j, ++k) {
                        if(j == ColumnSpot) {
                                temp.Data[i][k] += ColumnData[i];
                                ++k;
                        }
                        if(k<New_NbCols)
                        temp.Data[i][k] += Data[i][j];
                }
        }

        *this = temp;

        return *this;
}

//extraire une ligne

vecteur Matrice::Get_Row(int RowSpot)
{
vecteur RowData(NbCols);

        for(int j=0; j<NbCols; ++j) {
              //        if(i == RowSpot) {
              //                for(int j=0; j<NbCols; ++j) {
                               RowData[j]=Data[RowSpot][j];
             //         }
            //  break;
            //  }

        }
return RowData;
}

// Extraire une colone
vecteur Matrice::Get_Col(int ColSpot)
{
vecteur ColData(NbRows);

        for(int i=0; i<NbRows; ++i) {
                //for(int j=0; j<NbCols; ++j) {
                //if(j == ColSpot) {
                               ColData[i]=Data[i][ColSpot];
                //      }
              //        break;
           //   }

        }
return ColData;

}


//removes the specified row from the calling Matrice
 Matrice& Matrice::RemoveRow(const int Row)
{
        Matrice temp(NbRows-1, NbCols,0.0);

        for(int i=0, k=0; i<NbRows; ++i, ++k) {
                if(i == Row) ++i;
                for(int j=0; j<NbCols; ++j) {
                        temp.Data[k][j] = Data[i][j];
                }
        }

        *this = temp;

        return *this;
}

//removes the specified column from the calling Matrice
 Matrice& Matrice::RemoveColumn(const int Column)
{
        Matrice temp(0.0, NbRows, NbCols-1);

        for(int i=0; i<NbRows; ++i) {
                for(int j=0, l=0; j<NbCols; ++j, ++l) {
                        if(j == Column) ++j;
                        temp.Data[i][l] = Data[i][j];
                }
        }

        *this = temp;

        return *this;
}


//returns the covariant of the calling Matrice (transposed(obj) * obj)
 Matrice& Matrice::Get_Covariant() const
{
        Matrice temp;

        temp = *this;
        temp.Transpose();

        return(*this * temp);
}

//turns this Matrice into its covariant(obj = transposed(obj) * obj)
 Matrice& Matrice::Set_Covariant()
{
        *this = this->Get_Covariant();

        return *this;
}


//returns an identity Matrice of size Diagonal
 Matrice& Matrice::IdentityMatrice(int Diagonal)
{
        static Matrice temp;
        temp = Matrice( Diagonal, Diagonal,0.0);

        for(int q=0; q<Diagonal; ++q) {
                temp.Data[q][q] = 1;
        }

        return temp;
}


void Matrice:: affiche()
{
        cout << "[";
        for(int i=0; i<NbRows; ++i) {
                for(int j=0; j<NbCols; ++j) {
                        cout << "[" << Data[i][j] << "]";
                }
                if(NbRows-i-1) cout << "\n ";
        }
        cout << "]\n";


}

//easily display the Matrice, usually to the console
//formatted output


Matrice& Matrice::Get_LU(vecteur& perm,double&d)
{
        static const double petit=double(1.0e-20);
        int taille=NbRows;
        int imax,i,j,k;
        vector<double> v1(taille);
        double maxi,temp,som,dum;
        d=1;
        if(taille&&taille==NbCols)
         {
             for(i=0;i<taille;i++)
                {
                        maxi=double(0);
                        for(j=0;j<taille;j++)
                                if((temp=(double)fabs(Data[i][j]))>maxi)
                                        maxi=temp;
                        if(maxi==0)
                        {
                         this->Clear();
                         return *this;
                        }
                        v1[i]=1./maxi;
                }

              for(j=0;j<taille;j++)
                {
                //cherche les colonnnes  de lamatrice L (triangulaire inferieur
                        for(i=0;i<j;i++)
                                {
                                        som=Data[i][j];
                                        for(k=0;k<i;k++)
                                                som-=Data[i][k]*Data[k][j];
                                        Data[i][j]=som;
                                }
                        maxi=0;

               //cherches les lignes de la matrice triangulaire superieur
                        for(i=j;i<taille;i++)
                                {
                                som=Data[i][j];
                                for(k=0;k<j;k++)
                                        som-=Data[i][k]*Data[k][j];
                                Data[i][j]=som;
                        //v�rificstion de pivot
                                if((dum=v1[i]*fabs(som))>=maxi)
                                        {
                                                maxi=dum;
                                                imax=i;
                                        }
                                }
                        //permuation des lignes imax et j
                        if(j!=imax)
                                {
                                           for(k=0;k<taille;k++)
                                                {
                                                        dum=Data[imax][k];
                                                        Data[imax][k]=Data[j][k];
                                                        Data[j][k]=dum;
                                                }
                                        d*=-1;
                                        v1[imax]=v1[j];
                                }
                       // stoker dans perm la permuattion
                        perm[j]=imax;
                       // si le pivot est null 
                        if(Data[j][j]==0.0)
                                Data[j][j]=petit;
                        if(j!=taille-1)
                                {
                                        dum=1./Data[j][j];
                                        for(i=j+1;i<taille;i++)
                                                Data[i][j]*=dum;
                                }
                }
 }
 return *this;
}


vecteur Matrice::Get_Solution(vecteur& perm,const vecteur&b)
{
        vecteur r1;
        int taille=NbRows,i,j,ip,ii=-1;
        double som;

        if(taille&&taille==NbCols&&taille==(int)b.getDim())//   Get_Size())
                {
                        r1=b;
                        for(i=0;i<taille;i++)
                                {
                                        ip=perm[i];
                                        som=r1[ip];
                                        r1[ip]=r1[i];
                                        if(ii>=0)
                                                for(j=ii;j<=i-1;j++)
                                                        som-=Data[i][j]*r1[j];
                                        else if(som) ii=i;
                                        r1[i]=som;
                                }
                        for(i=taille-1;i>=0;i--)
                                {
                                        som=r1[i];
                                        for(j=i+1;j<taille;j++)
                                                som-=Data[i][j]*r1[j];
                                        r1[i]=som/Data[i][i];
                                }
                }
        return r1;
}


Matrice Matrice ::Get_Inverse(Matrice& K)

{

 double d;

 Matrice inv(NbRows,0,0.);

 vecteur perm(NbRows);

 K= K. Get_LU(perm,d);


 for(int i=0;i<NbRows;i++)

 {

  vecteur ei(NbRows);

  ei[i]=1.;

  if (i>0) ei[i-1]=0.;

  vecteur& xi=K.Get_Solution(perm,ei);

  inv.Add_Column(xi,i);


 }

return inv;


}


void Matrice::Get_Solution(vecteur& perm, vecteur& b, vecteur&  x)

{

        int taille=NbRows,i,j,ip,ii=-1;
        double som;


                       for(i=0;i<taille;i++) x[i]=b[i];

                        for(i=0;i<taille;i++)
                                {
                                        ip=perm[i];
                                        som=x[ip];
                                        x[ip]=x[i];
                                        if(ii>=0)
                                                for(j=ii;j<=i-1;j++)
                                                        som-=Data[i][j]*x[j];
                                        else if(som) ii=i;
                                        x[i]=som;
                                }
                        for(i=taille-1;i>=0;i--)
                                {
                                        som=x[i];
                                        for(j=i+1;j<taille;j++)
                                                som-=Data[i][j]*x[j];
                                        x[i]=som/Data[i][i];
                                }


}


vecteur Matrice:: Multiply(const vecteur& v)

{

assert (NbCols==v.getDim());

vecteur z(NbRows);

        for(int i=0; i<NbRows; ++i) {

                for(int j=0; j<NbCols; ++j) {

                        z[i] += Data[i][j]*v[j];
                }
        }


return z;


}


//extraire une matrice de la ligne 1 jusqua iigne 2 et de la colone 1 jusqua la colonne2

 /*

Matrice Matrice ::Get_SubMatrice(int ln1,int ln2, int cln1,int cln2)

{

int nb_lgne=ln2-ln1+1;

int nb_colne=cln2-cln1+1;

Matrice m(nb_lgne, nb_colne,0.);


   for(int i=ln1;i<=ln2;i++)

    {

    int I=i-ln1;

    for(int j=cln1;j<=cln2;j++)

      {

       int J=j-cln1;

       m(I,J)=Data[i][j] ;

      }


    }


 return m;

}

 */


std::ostream& operator <<(std::ostream& os, const Matrice& obj)

{
//ostr<<"MATRIX :"<<"#"<<obj.Get_NbRows()<<setw(2)<< "#"<<obj.Get_NbCols()<< endl;
os<<"------------------------------------------------------------------------------------"<<endl;
os<<"#  MATRICE:   DIM1=  "<< obj.Get_NbRows()<<"    , DIM2:  "<<obj.Get_NbCols()<<"  #"<<endl;
os<<"------------------------------------------------------------------------------------"<<endl;
  for(int i=0;i<obj.Get_NbRows();i++) {
     for(int j=0;j<obj.Get_NbCols();j++) {
         os<<setw(18)<<obj(i,j);
                       }
            os<<endl;
            }
os<<"------------------------------------------------------------------------------------"<<endl;            
            return os;
        /*
        for(int i=0; i<obj.Get_NbRows(); ++i) {
               for(int j=0; j<obj.Get_NbCols(); ++j) {
               if (i==0)ostr<<"["<<setw(9);
                        ostr <<setw(10)<<obj[i][j];
                }
                if(obj.Get_NbRows()-i-1) ostr << "\n ";
        }
        ostr<<"]"<<endl;

        return ostr; */
}

vecteur operator * (const Matrice& m,const vecteur& v)
{

assert (m.Get_NbCols()==v.getDim()) ;

vecteur z(m.Get_NbRows());

        for(int i=0; i<m.Get_NbRows(); ++i) {

                for(int j=0; j<m.Get_NbCols(); ++j) {

                        z[i] += m(i,j)*v[j];
                }
        }


return z;


}

/*

void Matrice::convert_Mat(const Matrice &mat,double valmin,MatCr &mat1)

        {
                int i,j,k,taille=mat.Get_NbRows();

                if(mat.Get_NbCols()==taille)
                        {
                                mat1._index.resize((unsigned int)(taille+1));
                                mat1._values.resize((unsigned int)(taille+1));
                                for(i=0;i<taille;mat1._values[(unsigned int)i]=mat(i,i),i++);
                                mat1._index[0]=taille+1;
                                k=taille;
                                for(i=0;i<taille;i++)
                                        {
                                                for(j=0;j<taille;j++)
                                                        if(fabs(mat(i,j))>valmin&&i!=j)
                                                                {
                                                                        k++;
                                                                        mat1._index.push_back(j);
                                                                        mat1._values.push_back(mat(i,j));
                                                                }
                                                mat1._index[i+1]=k+1;
                                        }
                        }
                else
                        {
                                mat1._index.resize(0);
                                mat1._values.resize(0);
                        }
        }


 */


Revision:	166
Committed:	Thu Oct 23 21:23:28 2008 UTC (16 years, 6 months ago) by souaissa
Original Path:	*magic/lib/solveur/solveur/src/matpleine.cpp*
File size:	25055 byte(s)
Error occurred while calculating annotation data.
Log Message:	mise a jour du solveur : ajout de la resolution partie plastique