// testovaci program pro template matrix
#include <stdlib.h> // 'cause of "random"
#include "tmatrix.h"
void definice(void)
{
MATRIX <float> a(3,3);
MATRIX <float> b=a,c;
MATRIX <float> d(2,2);
a.Print(" nulova 3x3");
b.Print(" totez ");
c=a;
c.Print(" jeste jednou totez");
d=c;
d.Print(" jeste jednou totez");
d=2;
d.Print(" sama dvojka v 3x3");
a.Unitary(5);
a.Print(" jednotkova 5x5 ");
b.Unitary(a.Rows());
b.Print(" jdnotkova ");
}
void secitani(void)
{
MATRIX <float> a(5,5);
MATRIX <float> b(5,5);
MATRIX <float> c(4,4);
MATRIX <float> d;
+a;
d = a + c;
d.Print(" chyba ");
a = 4;
b = 3;
d = a + b;
d.Print(" sama 7 ");
d += b;
d.Print(" sama 10 ");
}
void odecitani(void)
{
MATRIX <float> a(5,5);
MATRIX <float> b(5,5);
MATRIX <float> c(4,4);
MATRIX <float> d;
d = a - c;
d.Print(" chyba ");
a = -7;
b = 2;
d = -a - b;
d.Print(" sama 5 ");
d -= -a;
d.Print(" sama -2 ");
}
void nasobeni(void)
{
int i,j;
MATRIX <int> a(3,4);
MATRIX <int> b(4,5);
MATRIX <int> c;
// a.Fill(1,3,4);
a.Print("");
for (i=0;i<4;i++)
for (j=0;j<5;j++)
b(i,j) = i + j;
c = b * a;
c.Print(" chyba ");
c = a * b;
b.Print(" sikma ");
c.Print(" vysledek a (jednicky) * b ");
c = a;
c *= b;
/* c.Print(" totez ");
c *= 2;
c.Print(" dvojnasobek ");
c *= 2;
c.Print( " jeste dvojnasobek ");
c = b;
c *= a;
c.Print(" chyba ");
a = 2;
b = a * 7;
b.Print(" sama 14 ");
a = 2;
b = 7 * a ;
b.Print(" sama 14 ");
*/
}
void inverze(void)
{
int i,j;
MATRIX <double> b(10,10);
MATRIX <double> a(5, 7);
MATRIX <double> c(5,5);
MATRIX <double> d(4,4);
d.Inv();
d.Print(" chyba ");
c = 7;
c.Inv();
c.Print(" chyba ");
a.Inv();
a.Print(" chyba ");
for (i=0;i<10;i++)
for (j=0;j<10;j++)
// b(i,j) = random(10);// borland
b(i,j) = rand() / RAND_MAX; // microsoft
b.Print(" b sikma ");
a = b;
a.Print(" totez ");
a.Inv();
a.Print (" inverzni ");
c = a * b;
c.Print (" jednotkova ");
for (i=0;i<10;i++)
for (j=0;j<10;j++)
// b(i,j) = random(10);// borland
b(i,j) = rand() / RAND_MAX; // microsoft
a = b;
a.Inv();
c = a * b;
c.Print (" jednotkova ");
for (i=0;i<10;i++)
for (j=0;j<10;j++)
// b(i,j) = random(10); // borland
b(i,j) = rand() / RAND_MAX; // microsoft
a = b;
a.Inv();
c = a * b;
c.Print (" jednotkova ");
for (i=0;i<10;i++)
for (j=0;j<10;j++)
// b(i,j) = random(10);// borland
b(i,j) = rand() / RAND_MAX; // microsoft
a = b;
a.Inv();
c = a * b;
c.Print (" jednotkova ");
for (i=0;i<10;i++)
for (j=0;j<10;j++)
// b(i,j) = random(10);// borland
b(i,j) = rand() / RAND_MAX; // microsoft
a = b;
c = Inv(b) * b;
c.Print (" jednotkova ");
for (i=0;i<10;i++)
for (j=0;j<10;j++)
// b(i,j) = random(10);// borland
b(i,j) = rand() / RAND_MAX; // microsoft
}
void transponovani(void)
{
MATRIX <float> b(6,6),c;
int i,j;
for (i=0;i<6;i++)
for (j=0;j<6;j++)
b(i,j) = i+10*j;
b.Print("original");
c = b;
c.Transp();
c.Print("transponovana");
c = Transp(b);
c.Print("transponovana");
}
void deleni()
{
MATRIX<double > a(5,5),b;
a = 3;
a /= 2;
a.Print(" 5x5 plna 1.5 ");
a = 3;
a /= 0;
a.Print(" nulou se nedeli ");
a.Fill(3,5,5);
b = a / 2.;
b.Print(" 5x5 plna 1.5 ");
b = a / 0.;
b.Print(" nulou se nedeli ");
}
void reseni_rovnic(void)
{
MATRIX <double> y,x(5,5),a(5,1),v,x1(10,5);
int i,j;
for (i=0;i<5;i++)
{
a(i,0) = 2.1 * i - i * i ;
for (j=0;j<5;j++)
// x(i,j) = random(30); // borland
x(i,j) = rand() / RAND_MAX; // microsoft
}
x.Print(" Matice x");
a.Print(" Matice koeficientu a");
y = x * a;
y.Print(" Matice vysledku y");
v = SolveLinEq(x,y);
v.Print(" vysledna matice koeficientu ");
for (i=0;i<10;i++)
{
for (j=0;j<5;j++)
// x1(i,j) = random(50);// borland
x(i,j) = rand() / RAND_MAX; // microsoft
}
x1.Print(" Matice x1 pro metodu nejmensich ctvercu");
// vice rovnic nez koeficientu
a.Print(" Matice koeficientu a");
y = x1 * a;
y.Print(" Matice vysledku y");
v = SolveLinEq(x1,y);
v.Print(" vysledna matice koeficientu ");
}
void zmena_rozmeru(void)
{
MATRIX <double> m;
m.Fill( 2,3,4);
m.Print(" vytvoreno Fill (3,4,5) ");
m.AddRow(5);
m.Print(" add row (5) ");
m.AddColumn(6);
m.Print(" add column(6) ");
m.AddColumn(8);
m.Print(" add column(8) ");
m.AddRow();
m.Print(" add row () ");
m.AddRowCol(12);
m.Print(" add row column(12) ");
}
void vymena(void)
{
MATRIX <float> b(6,6);
int i,j;
for (i=0;i<6;i++)
for (j=0;j<6;j++)
b(i,j) = i+10*j;
b.Print(" original ");
b.ChangeRows(1,3);
b.Print(" Change Rows (1,3)");
b.Print(" original ");
b.ChangeColumns(4,2);
b.SetName(" b je tato matice");
b.Print(" Change Columns (4,2)");
}
// priklad na ziskani matice jako navratove hodnoty z funkce
MATRIX <int> VratMatici(MATRIX <int> &a,MATRIX<int> &b)
{
MATRIX <int> pom;
pom = a * b;
return pom;
}
void Return(void)
{
MATRIX <int> a(3,4);
MATRIX <int> b(4,5);
MATRIX <int> c;
c = VratMatici(a,b);
c.Print(" vysledek ");
}
void main(void)
{
// definice();
// secitani();
// odecitani();
// nasobeni();
// inverze();
// deleni();
// transponovani();
// reseni_rovnic();
// zmena_rozmeru();
// vymena();
Return();
printf("====================== THE END =======================");
}
==================== hlavička ================
// Template for matrix operations
// first argument row, second column
// definition: MATRIX <type> name....; (MATRIX <float> a,b(5,5),c(6,32,2);
// work is the same like for another class a.Inv();b.Unitary();
// e += d = a * b + c; Inv(a);
// changes from previous version:
// name of funcion Transposition changed to Transp
// both InsertMatrix & CutMatrix become InserMatrix
// Print function has a new parameter, it is the output stream (FILE *...)
// 1.01
// Unitary matrix has implicit parameter now
// new SolveLinEq (previouse version was very wrong)
// 1.02
// deleted duplicity in implicit parameters
#ifndef MATRIX_H
#define MATRIX_H
#include <stdio.h>
//#include <mem.h>
#include <conio.h>
#include <math.H>
#include <string.h>
#define ERR_NOTHING 0
#define ERR_NO_MEMORY -1
#define ERR_BAD_SIZE -2
#define ERR_COUNTINGS -3
#define ERR_CANNOT_CHANGE -4
#define ERR_CANNOT_INSERT -5
#define ERR_SMALL_PIVOT -6
#define ERR_DIVIDE_BY_ZERO -7
#define ERR_INDEX -8
#define MATRIX_ALOCATION 1
#define MATRIX_INSIDE_RET 2
#define MATRIX_EQUAL 3
#define MATRIX_EQUAL_SIZE 4
#define MATRIX_DEALOCATE 5
#define MATRIX_MULTIPLICATION 6
#define MATRIX_INSIDE 7
#define MATRIX_SQUARE 8
#define LenghtofErrTxt 25
#define NofErr 10
#define textchyba0 {" bez chyby "}
#define textchyba1 {" malo pameti "}
#define textchyba2 {" spatne rozmery "}
#define textchyba3 {" chyba pri vypoctu "}
#define textchyba4 {" nelze zmenit "}
#define textchyba5 {" nelze vlozit "}
#define textchyba6 {" maly pivot "}
#define textchyba7 {" deleni nulou "}
#define textchyba8 {" spatny index pole "}
#define textchyba9 {" nezavedena "}
//template <class T> class MATRIX;
//MATRIX<T> Inv(MATRIX<T> &orig,double ZeroRange=0.);
template <class T>
class MATRIX {
// mat is the pointer on column of pointers on rows
// if used first operand is row, second column
int rows; // number of rows - first argument
int columns; // number of columns - second
int error; // error code
T **mat; // matrix data
T null; // null value, this value is returned if something
// must be returned and apropriate value isn't available
char *Name; // name of matrix
public:
//================== Constructor & destructor ====================
// default constructor
MATRIX() {rows=0;columns=0;mat=NULL;error=0;null=0;Name=NULL;}; //1.01
// constructor for matrix with size - rows2, columns2
// if data are present matrix is filled with this ones
// else matrix is filled by zero value
// if (dealoc==1) memory with data is dealocated
MATRIX(int rows2,int columns2,T *data=NULL,int dealoc=0); //1.01
// copy constructor
MATRIX(MATRIX &m); //1.01
// destructor
~MATRIX(void){DealocateMatrix();if (Name!=NULL) delete[] Name;
Name = NULL;} //1.01
//================== functoins ====================
void AddRow(double initvalue = 0); //1.01
void AddColumn(double initvalue = 0); //1.01
void AddName(char *name); //1.01
void AddRowCol(double initvalue = 0)
{ AddColumn(initvalue); AddRow(initvalue);} //1.01
// matrix memory alocation
int AlocateMatrix(T **&mat,int rows,int columns); //1.0
int AlocateMatrix(void); //1.0
void ChangeColumns(int col1,int col2);
void ChangeRows(int row1,int row2);
// check indexes for different types of situation
// type of checking index, olderr - error state of second matrix
// rows2, columns2 compared indexes (second matrix)
//1.0
int CheckIndex(int type,int olderr=ERR_NOTHING,int rows2=0,int columns2=0);
//1.01
void ClearName(void) { if (Name!=NULL) delete[] Name;Name = NULL;}
// return number of columns
int Columns(void) {return columns;} //1.0
// dealocate memory for matrix
void DealocateMatrix(T **&mat, int rows); //1.0
void DealocateMatrix(void); //1.0
// fill matrix vith given value, if needed new size is generated
// better (if matrix don't change size) is operator = ("= value")
void Fill(double value,int rows2=-1,int columns2=-1); //1.0
//1.01
char *GetName() { return Name;}
// matrix m is inserted to 'this' matrix
MATRIX& InsertMatrix(MATRIX &m,int fromrow,int fromcolumn,
int torow,int tocolumn,int crows,int ccolumns);//1.0
// inversion of matrix
// range for zero pivot - matrix is treated as singular matrix
// i.e. err is generated if pivot element is less then ZeroRange
void Inv(double ZeroRange = 0); //1.0
// friend MATRIX<T> Inv(MATRIX<T> &orig,double ZerRange=0.); //1.02
// return inverse matrix
// function is called every time when any error occured
void MatrixError(int e);
// print matrix on screen
void Print(char text[] = NULL,FILE *f = stdout); //1.0
int Rows(void) {return rows;} //1.0
void SetName(char *NewName); //1.01
// return matrix dimensions //1.0
void Size(int &rows2,int &columns2) { rows2 = rows; columns2 = columns;}
// transposition of matrix
void Transp(void); //1.0
friend MATRIX<T> Transp(MATRIX<T> &orig); //1.0
// generate Unitary matrix - square with dimension dim
void Unitary(int dim = 0); //1.01
//================== Operators ===========================
// return reference on element row,column in matrix
T& operator () (int row,int column); //1.0
// unary plus
MATRIX& operator + (void); //1.0
// unary minus, return negative matrix - don't change original matrix
MATRIX operator - (void); //1.0
// operator equal
MATRIX& operator = (const MATRIX &m); //1.02
MATRIX& operator = (double v); //1.0
// add two matrices
MATRIX operator + (MATRIX &m); //1.0
MATRIX& operator += (const MATRIX &m); //1.02
// substraction of two matrix
MATRIX operator - (MATRIX &m); //1.0
MATRIX& operator -= (const MATRIX &m); //1.02
// multiplications & dividions
MATRIX operator * (MATRIX &m); //1.0
MATRIX& operator *= (const MATRIX &m); //1.02
MATRIX& operator *= (T hodn); //1.0
friend MATRIX<T> operator * (T hodn,MATRIX<T> &m); //1.0
friend MATRIX<T> operator * (MATRIX<T> &m,T hodn); //1.0
MATRIX& operator /= (T hodn); //1.0
friend MATRIX<T> operator / (MATRIX<T> &m,T hodn); //1.0
};
//================== Constructor & destructor ====================
// constructor for matrix with size - rows2, columns2
// if input data are present, matrix is filled with this ones,
// else is filled by zero value
// if (dealoc==1) memory with data is dealocated
//1.0
template <class T> MATRIX<T>::
MATRIX(int rows2,int columns2,T *data,int dealoc)
{
int i;
long int lenght;
null=0;
rows=rows2;
columns=columns2;
if (CheckIndex(MATRIX_ALOCATION)!=ERR_NOTHING)
return; // err occur
if (data!=NULL) // alocation of matrix is done
{ // copy data
lenght=(long int)sizeof(T)*columns;
for(i=0;i<rows;i++)
memcpy(mat[i],&data[i*columns],lenght);
}
if (dealoc) // dealocate data
delete[] data;
Name = NULL;
}
// copy konstruktor
//1.0
template <class T> MATRIX<T>::
MATRIX(MATRIX<T> &m)
{
int i;
long int lenght;
rows=m.rows;
columns=m.columns;
if (CheckIndex(MATRIX_ALOCATION)!=ERR_NOTHING)
return; // if errorr
lenght=(long int)sizeof(T)*columns;
for(i=0;i<rows;i++)
memcpy(mat[i],m.mat[i],lenght); // copy data
Name = NULL;
SetName(m.Name);
}
//================== functoins ====================
// add column to the matrix
//1.01
template <class T> void MATRIX<T>::
AddColumn(double initvalue)
{
T *row;
int i,j;
for (i=0;i<rows;i++)
{
row = (T*) new T[columns+1];
if (row == NULL)
{
MatrixError(ERR_NO_MEMORY);
return;
}
for(j=0;j<columns;j++)
row[j] = mat[i][j];
delete[] mat[i];
mat[i] = row;
mat[i][columns] = initvalue;
}
columns ++;
}
// add row to the matrix
//1.01
template <class T> void MATRIX<T>::
AddRow(double initvalue)
{
T **col;
int i;
col = (T**) new T*[rows+1];
if (col == NULL)
{
MatrixError(ERR_NO_MEMORY);
return;
}
col[rows] = (T*) new T[columns];
if (col[rows] == NULL)
{
MatrixError(ERR_NO_MEMORY);
delete[] col;
return;
}
for(i=0;i<rows;i++)
col[i] = mat[i];
delete[] mat;
mat = col;
for (i=0;i<columns;i++)
mat[rows][i] = initvalue;
rows++;
}
//1.0
template <class T> void MATRIX<T>::
AddName(char *NewName)
{
int delka;
if (NewName==NULL)
return;
delka = strlen(NewName);
if (Name)
delka += strlen(Name);
if (delka == 0)
return;
char *n = (char *) new char[delka+1];
if (n == NULL)
return;
n[0] = '\0';
if (Name)
{
strcpy(n,Name);
delete[] Name;
Name = NULL;
}
strcat(n,NewName);
Name = n;
}
// alocate memory for matrix
//1.0
template <class T> int MATRIX<T>::
AlocateMatrix(T **&mat,int rows,int columns)
{
int i,j;
mat = (T**) new T*[rows]; // column
if (mat == NULL)
return (ERR_NO_MEMORY);
for (i=0;i<rows;i++)
mat[i] = NULL;
for (i=0;i<rows;i++)
{
mat[i] = (T*) new T[columns]; // rows
if (mat[i] == NULL)
{ // errorr
DealocateMatrix();
return(ERR_NO_MEMORY);
}
for (j=0;j<columns;j++)
mat[i][j] = 0; // all elements are zero
}
return(ERR_NOTHING);
}
// alocate memory for matrix
//1.0
template <class T> int MATRIX<T>::
AlocateMatrix(void)
{
return(AlocateMatrix(mat,rows,columns));
}
// exchange two columns
//1.01
template <class T> void MATRIX<T>::
ChangeColumns(int col1,int col2)
{
T p;
int i;
if ((col1>=0) && (col1<columns) && (col2>=0) && (col2 < columns))
{
for (i=0;i<rows;i++)
{
p = mat[i][col1];
mat[i][col1] = mat[i][col2];
mat[i][col2] = p;
}
}
else
MatrixError(ERR_CANNOT_CHANGE);
}
// exchange two rows
//1.01
template <class T> void MATRIX<T>::
ChangeRows(int row1,int row2)
{
T p;
int i;
if ((row1>=0) && (row1<rows) && (row2>=0) && (row2 < rows))
{
for (i=0;i<columns;i++)
{
p = mat[row1][i];
mat[row1][i] = mat[row2][i];
mat[row2][i] = p;
}
}
else
MatrixError(ERR_CANNOT_CHANGE);
}
// check indexes for different types of situation
// type of checking index, olderr - error state of second matrix
// rows2, columns2 compared indexes (second matrix)
//1.01
template <class T> int MATRIX<T>::
CheckIndex(int type,int olderr,int rows2,int columns2)
{
null=0; // zero element of class MATRIX
if (olderr!=ERR_NOTHING)
error=olderr; // old data are bad - everything is bad
// but think of this this way is probably also bad
else
{
switch(type){
case MATRIX_DEALOCATE : // only dealocation of matrix
error=ERR_COUNTINGS;
break; // 0
case MATRIX_ALOCATION: // alocation of matrix
mat=NULL; // first definition
if ((rows>0)&&(columns>0))
error=AlocateMatrix(); // size of matrix is good
else
error=ERR_BAD_SIZE; // is bad
break;
case MATRIX_SQUARE: // square matrix test
if (columns!=rows)
error=ERR_BAD_SIZE;
break;
case MATRIX_EQUAL: // if different sizes, new matrix alocation
if ((columns!=columns2)||(rows!=rows2))
{ // matrix have different size
DealocateMatrix(); // delete old matrix
columns = columns2;rows = rows2;
error = AlocateMatrix(); // alocate proper size of matrix
}
break;
case MATRIX_EQUAL_SIZE: // if two matrices are the same size
if (error!=ERR_NOTHING)
break;
if ((columns!=columns2)||(rows!=rows2)||(columns<=0)||(rows<=0))
return(ERR_BAD_SIZE);
else
return(ERR_NOTHING);
case MATRIX_MULTIPLICATION: // multiplication
if (error!=ERR_NOTHING)
break;
if ((columns!=rows2)||(columns<=0)||(rows<=0)||(columns2<=0))
return ERR_BAD_SIZE;
else
return ERR_NOTHING;
case MATRIX_INSIDE: // index columns2,rows2 must lay inside given matrix
if (error!=ERR_NOTHING)
break;
if ((columns2>=columns)||(rows2>=rows)||(columns2<0)||(rows2<0))
error=ERR_BAD_SIZE;
break;
case MATRIX_INSIDE_RET: // if rows2,columns2 are inside given matrix
if (error!=ERR_NOTHING)
break;
if ((columns2>=columns)||(rows2>=rows)||(columns<0)||(rows2<0))
return (ERR_BAD_SIZE);
else
return (ERR_NOTHING);
default: MatrixError(ERR_INDEX);
getch();
}
}
if (error!=ERR_NOTHING)
{
MatrixError(error);
DealocateMatrix();
}
return(error);
}
// dealocate memory for matrix
//1.0
template <class T> void MATRIX<T>::
DealocateMatrix(T **&mat, int rows)
{
int i;
if (mat!=NULL)
{
for (i=0;i<rows;i++)
{
if (mat[i]!=NULL)
delete[] mat[i];
}
delete[] mat;
mat=NULL;
}
}
// dealocate memory for matrix
//1.0
template <class T> void MATRIX<T>::
DealocateMatrix(void)
{
DealocateMatrix(mat,rows);
rows=0;
columns=0;
}
// fill matrix vith given value, if needed new size is generated
//1.0
template <class T> void MATRIX<T>::
Fill(double value,int rows2,int columns2)
{
int i,j,lenght;
if (columns2>0) // is new size needed
{
if (CheckIndex(MATRIX_EQUAL,ERR_NOTHING,rows2,columns2)!=ERR_NOTHING)
return; // if sizes are different - redefine size
}
else
if (CheckIndex(MATRIX_INSIDE,ERR_NOTHING,rows2,columns2)!=ERR_NOTHING)
return;
lenght=sizeof(T)*columns;
for (j=0;j<columns;j++)
mat [0][j]=value;
for (i=1;i<rows;i++)
memcpy(mat[i],mat[0],lenght);
}
// insert matrix to another one
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
InsertMatrix(MATRIX<T> &m,int fromrow,int fromcolumn,
int torow,int tocolumn,int crows,int ccolumns)
{
int i,j;
if ((fromcolumn+ccolumns)>m.columns)
ccolumns = m.columns-fromcolumn;
if ((tocolumn+ccolumns)>columns)
ccolumns = columns - tocolumn;
if ((fromrow+crows)>m.rows)
crows = m.rows-fromrow;
if ((torow+crows)>rows)
crows = rows - torow;
if ((torow>=rows)||(torow<0)||(tocolumn>=columns)||(tocolumn<0)
||(crows==0)||(ccolumns==0))
{
MatrixError(ERR_CANNOT_INSERT);
return *this;
}
for (i=0;i<crows;i++)
for (j=0;j<ccolumns;j++)
mat[torow+i][tocolumn+j] = m.mat[fromrow+i][fromcolumn+j];
return *this;
}
// inversion of matrix
//1.0
template <class T> void MATRIX<T>::
Inv(double ZeroRange)
{
int i,j,k,ip,jp;
double pivot,pvt,*r;
MATRIX a=*this,b;
int n=columns;
if (CheckIndex(MATRIX_SQUARE,a.error))
{
error=ERR_BAD_SIZE;
return;
}
b.Unitary(columns);
if (columns <= 0)
return;
r=new double [columns];
if ((r==NULL)||(b.error!=ERR_NOTHING))
{
CheckIndex(MATRIX_DEALOCATE);
MatrixError(ERR_COUNTINGS);
return ;
}
for(i=0;i<n;i++)
r[i]=1;
for(k=0;k<n;k++)
{
pivot=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
pvt=a.mat[i][j]*r[i];
if(pvt<0.0)
pvt=-pvt;
if(pvt>pivot)
{
pivot=pvt;
ip=i;
jp=j;
}
}
if (fabs(pivot)<=ZeroRange)
{
delete[] r;
CheckIndex(MATRIX_DEALOCATE);
error=ERR_COUNTINGS;
MatrixError(ERR_SMALL_PIVOT);
return ;
}
r[ip]=0;
pivot=a.mat[ip][jp];
for(j=0;j<n;j++)
{
a.mat[ip][j]/=pivot;
b.mat[ip][j]/=pivot;
}
for (i=0;i<n;i++)
if(i!=ip)
{
for(j=0;j<n;j++)
b.mat[i][j]-=a.mat[i][jp]*b.mat[ip][j];
for(j=0;j<n;j++)
if(j!=jp)
a.mat[i][j]-=a.mat[i][jp]*a.mat[ip][j];
}
for(i=0;i<n;i++)
if (i!=ip)
a.mat[i][jp]=0;
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
if (a.mat[i][j]!=0.0)
for(k=0;k<n;k++)
mat[j][k]=b.mat[i][k];
delete[] r;
return ;
}
//1.0
template <class T> MATRIX<T>
Inv(MATRIX<T> &mat,double ZeroRange=0.)
{
MATRIX<T> pom=mat;
pom.Inv(ZeroRange);
return pom;
}
// function is called when any error occured
//1.01
template <class T> void MATRIX<T>::
MatrixError(int ec)
{
// you can add your code here
// you have "Name" of matrix & error code "ec" so you can print error message
//--- All errorr messages are not directed to this procedure yet.
}
// print matrix on screen
//1.0
template <class T>void MATRIX<T>::
Print(char text[],FILE * f)
{
int i,j;
int err;
char errtexts[NofErr][LenghtofErrTxt]= {
textchyba0,
textchyba1,
textchyba2,
textchyba3,
textchyba4,
textchyba5,
textchyba6,
textchyba7,
textchyba8,
textchyba9
};
if (f == NULL)
return;
if (text!=NULL)
{
if (Name!=NULL)
fprintf(f,"\n============= %s ; %s ==================\n",Name,text);
else
fprintf(f,"\n============= %s ==================\n",text);
}
else
{
if (Name!=NULL)
fprintf(f,"\n============= %s ==================\n",Name);
else
fprintf(f,"\n====================================\n");
}
if (error<0)
err = -error;
else
err = error;
fprintf(f," sloupcu = %d radku = %d chyba = %s \n",
columns,rows,&errtexts[err][0]);
for (i=0;i<rows;i++)
{
for (j=0;j<columns;j++)
fprintf(f," %f ",(double)mat[i][j]);
fprintf (f,"\n");
}
if (f==stdout)
getch();
}
// sets matrix name
//1.01
template <class T> void MATRIX<T>::
SetName(char *NewName )
{
if (NewName == NULL)
return;
if (Name != NULL)
delete[] Name;
Name = (char *) new char[strlen(NewName)+1];
if (Name == NULL)
return;
strcpy(Name,NewName);
// memcpy(Name,NewName,(strlen(NewName)+1)*sizeof(char));
}
// transposition of matrix
//1.0
template <class T> void MATRIX<T>::
Transp(void)
{
T **newmat;
T pom;
int i,j;
if (error!=ERR_NOTHING)
return;
if (columns==rows)
{
for(i=0;i<columns;i++)
for (j=i+1;j<columns;j++)
{
pom=mat[i][j];
mat[i][j]=mat[j][i];
mat[j][i]=pom;
}
}
else
{
if (AlocateMatrix(newmat,columns,rows)!=ERR_NOTHING)
{
CheckIndex(MATRIX_DEALOCATE);
return;
}
for (i=0;i<columns;i++)
for (j=0;j<rows;j++)
newmat[i][j]=mat[j][i];
DealocateMatrix(mat,rows);
mat=newmat;
i=rows;
rows=columns;
columns=i;
}
}
//1.0
template <class T> MATRIX<T>
Transp(MATRIX<T> &mat)
{
MATRIX<T> pom=mat;
pom.Transp();
return pom;
}
// Unitary matrix
//1.01
template <class T> void MATRIX<T>::
Unitary(int rad)
{
int i,j;
if (rad==0)
rad=columns;
if ((columns!=rad)||(rows!=rad))
{
if (mat!=NULL)
DealocateMatrix();
}
else
{
for (i=0;i<rad;i++)
{
for (j=0;j<rad;j++)
mat[i][j] = 0;
mat [i][i] = 1;
}
return;
}
columns=rad;rows=rad;
if (CheckIndex(MATRIX_ALOCATION)==ERR_NOTHING)
{
for (i=0;i<rad;i++)
mat[i][i]=1;
}
}
//================== Operators ===========================
// return reference on element (row,column) of matrix
//1.0
template <class T> T& MATRIX<T>::
operator ()(int row,int column)
{
if (CheckIndex(MATRIX_INSIDE_RET,ERR_NOTHING,row,column)==ERR_NOTHING)
return mat[row][column];
else // wrong size
return null; // return reference on zero element (sets up in CheckIndex)
}
// operator =
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
operator = (const MATRIX<T> &m)
{
int i;
int lenght;
if (this == &m)
return *this;
if (CheckIndex(MATRIX_EQUAL,m.error,m.rows,m.columns)!=ERR_NOTHING)
return *this; // errorr
lenght = sizeof(T)*columns;
for(i=0;i<rows;i++)
memcpy(mat[i],m.mat[i],lenght); // copy data to new matrix
return *this;
}
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
operator = (double v)
{
int i,j;
for(i=0;i<rows;i++)
for (j=0;j<columns;j++)
mat[i][j] = (T)v;
return *this;
}
// unary plus
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
operator + (void)
{
return *this;
}
// unary minus - return negative matrix - don't change original matrix
//1.0
template <class T> MATRIX<T> MATRIX<T>::
operator - (void)
{
int i,j;
MATRIX<T> pom(rows,columns);
if (pom.error==ERR_NOTHING)
for (i=0;i<rows;i++)
for (j=0;j<columns;j++)
pom.mat[i][j]=-mat[i][j];
pom.ClearName();
pom.AddName(" - ");
pom.AddName(Name);
return MATRIX<T>(pom);
}
// add two matrices
//1.0
template <class T> MATRIX<T> MATRIX<T>::
operator + (MATRIX<T> &m)
{
int i,j;
MATRIX<T> pom;
if (CheckIndex(MATRIX_EQUAL_SIZE,m.error,m.rows,m.columns)
==ERR_NOTHING)
{
pom=m;
for (i=0;i<rows;i++)
for (j=0;j<columns;j++)
pom.mat[i][j]+=mat[i][j];
pom.SetName(Name);
pom.AddName(" + ");
pom.AddName(m.Name);
}
else
pom.error = ERR_BAD_SIZE;
return MATRIX<T>(pom);
}
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
operator += (const MATRIX<T> &m)
{
int i,j;
if (CheckIndex(MATRIX_EQUAL_SIZE,m.error,m.rows,m.columns)==ERR_NOTHING)
{
for (i=0;i<rows;i++)
for (j=0;j<columns;j++)
mat[i][j]+=m.mat[i][j];
}
else
CheckIndex(MATRIX_DEALOCATE);
return *this;
}
// substraction of two matrix
//1.0
template <class T> MATRIX<T> MATRIX<T>::
operator - (MATRIX<T> &m)
{
int i,j;
MATRIX<T> pom ;
if (CheckIndex(MATRIX_EQUAL_SIZE,m.error,m.rows,m.columns)==ERR_NOTHING)
{
pom=*this;
for (i=0;i<rows;i++)
for (j=0;j<columns;j++)
pom.mat[i][j]-=m.mat[i][j];
pom.SetName(Name);
pom.AddName(" - ");
pom.AddName(m.Name);
}
else
pom.error = ERR_BAD_SIZE;
return MATRIX<T>(pom);
}
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
operator -= (const MATRIX<T> &m)
{
int i,j;
if (CheckIndex(MATRIX_EQUAL_SIZE,m.error,m.rows,m.columns)==ERR_NOTHING)
{
for (i=0;i<rows;i++)
for (j=0;j<columns;j++)
mat[i][j]-=m.mat[i][j];
}
else
CheckIndex(MATRIX_DEALOCATE);
return *this;
}
// multiplaing two matrixes
//1.0
template <class T> MATRIX<T> MATRIX<T>::
operator * (MATRIX<T> &m)
{
int i,j,k;
if (CheckIndex(MATRIX_MULTIPLICATION,m.error,m.rows,m.columns)==ERR_NOTHING)
{
MATRIX pom(rows,m.columns);
pom.SetName(Name);
pom.AddName(" * ");
pom.AddName(m.Name);
pom.columns=m.columns;
pom.rows=rows;
for (i=0;i<rows;i++)
for (j=0;j<m.columns;j++)
for (k=0;k<columns;k++)
pom.mat[i][j]+=mat[i][k]*m.mat[k][j];
return MATRIX<T>(pom);
}
else
{
MATRIX pom;
pom.error = ERR_BAD_SIZE;
return MATRIX<T>(pom);
}
}
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
operator *= (const MATRIX<T> &m)
{
int i,j,k;
int px,py;
T **pmat=NULL;
px=AlocateMatrix(pmat,rows,m.columns);
if (px==ERR_NOTHING)
{
if (CheckIndex(MATRIX_MULTIPLICATION,m.error,m.rows,m.columns)==ERR_NOTHING)
{
px=m.columns;
py=rows;
for (i=0;i<rows;i++)
for (j=0;j<m.columns;j++)
for (k=0;k<columns;k++)
pmat[i][j]+=mat[i][k]*m.mat[k][j];
DealocateMatrix(mat,rows);
mat=pmat;
columns=px;
rows=py;
return *this;
}
}
if (pmat!=NULL)
DealocateMatrix(pmat,rows);
CheckIndex(MATRIX_DEALOCATE);
return *this;
}
//1.0
template <class T> MATRIX<T>& MATRIX<T>::
operator *=(T hodn)
{
int i,j;
for (i=0;i<rows;i++)
for (j=0;j<columns;j++)
mat[i][j]*=hodn;
return *this;
}
//1.01
template <class T> MATRIX<T>& MATRIX<T>::
operator /=(T hodn)
{
int i,j;
if (hodn != 0)
for (i=0;i<rows;i++)
for (j=0;j<columns;j++)
mat[i][j]/=hodn;
else
{
CheckIndex(MATRIX_DEALOCATE);
MatrixError(ERR_DIVIDE_BY_ZERO );
error = ERR_DIVIDE_BY_ZERO;
}
return *this;
}
//1.0
template <class T> MATRIX<T>
operator * (MATRIX<T> &m,T hodn)
{
MATRIX<T> pom=m;
pom.AddName("*value, ");
pom*=hodn;
return pom;
}
//1.0
template <class T> MATRIX<T>
operator * (T hodn,MATRIX<T> &m)
{
MATRIX<T> pom=m;
pom.AddName("*value, ");
pom*=hodn;
return pom;
}
//1.01
template <class T> MATRIX<T>
operator / (MATRIX<T> &m,T hodn)
{
MATRIX<T> pom=m;
pom.AddName("/value, ");
if (hodn == 0)
{
m.DealocateMatrix();
m.MatrixError(ERR_DIVIDE_BY_ZERO );
return pom;
}
pom/=hodn;
return pom;
}
// linear equation solution
// it is solved by pseudoinversion (least square solution if there is
// more equations than coeficients)
// y matrix of results
// x values
// a coeficients we are looking for
// y = x . a
template <class T> MATRIX<T>
SolveLinEq(MATRIX<T> &x, MATRIX<T> &y)
{
MATRIX<T> a = x;
a.Transp();
a *= x;
a.Inv();
a *= Transp(x);
a *= y;
return a;
}
#endif
// methods which will be added
// operators >> <<