// 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 >> <<