-----------------------------------------------------------------------
                            Matrix class Library 
                            Instruction Manual
                            Copyright (c) Gerox
                            Author Kazuyuki Kobayashi
                            Version 2.0  01/03/95
                                         03/22/96
                                         01/10/97 
-----------------------------------------------------------------------

0. Introduction

This class library is aimed at easy operation of matrix calculations by use of the C++ operator overload function. By using template function, different matrix type such as float, double, complex and BCD type matrices can be manipulate in C++.
To use this matrix library, no prior knowledge of the conditions for matrix multiplication is necessary. 
One of powerful capability of this library is we can use similar expression as MATLAB like matrix manipulation. For example, we can express least square is 
x = 1/(~A*A)*(~A)*B;  // Note ~ means transpose matrix

1. Matrix class library
Matrix library can work template supported C++ compiler, such as Turbo C++ Ver. 3.0, Borland C++ Ver. 4.0 and Watcom C++ Ver. 9.5 as well as Visual C++ Version 4. BCD type is only available by using Borland C compilers.

1.1.1 memory check option
If check memory allocation check is required.
#define DEBUG
otherwise
//#define DEBUG
//Commented out of Debugging option may speed up the matrix calculation process.

1.2 Definition of Matrix
In order to avoid terms confusion, here is matrix of reference
(matrix A) which is of size m x n.
Matrix<double> A(m,n) means

                          Column
                Row        0       1       2          n-1
                 0    |   a0,0    a0,1    a0,2   ...  a0,n-1   |
     A =         1    |   a1,0    a1,1    a1,2   ...  a1,n-1   |
                 2    |   a2,0    a2,1    a2,2   ...  a2,n-1   |
                      |   ...     ...     ...    ...  ...      |
                 m-1  |   am-1,0  am-1,1  am-1,2 ...  am-1,n-1 |

A constructor is available. The four constructor types are as following.

(a) Definition of matrix by size
	Matrix<double>     Ad(2,2);// Defines size of matrix Ad as 2*2 for double type
	Matrix<complex>     Ac(2,2);// Defines size of matrix Ac as 2*2 for complex type
	Matrix<bcd> Abcd(2,2);// Defines size of matrix Abcd as 2*2 for BCD type
	//at this time, each element of the matrix is cleared by zero.
	
(b)Definition by using initial substitution
	Matrix<double> B = A;
	// Defines size of Matrix B as equal to the size of matrix A.
	// Provided the size of A was previously defined.

(c) Definition of matrix without size definition
	Matrix<double> C;// Defines C as a matrix but without a specified 
           // Actually matrix size is defined as 0x0.

(d) Definition by using pointers

	Matrix<double> *A = new Matrix<double>(2,2);
	Matrix<dobule> *B = new Matrix<double>;
	A->operator()(0,1) = 2;
	A->operator()(1,1) = 4;
	//or you can define A->elem(1,1,4);
	*B = 1 / (~*A * *A);
	B->mprint();
	delete A;// Pointer type definition requires the delete function.
	delete B;


1.3 data input output method

(a)input method from keyboard
	//input form
	b.minput(1,0);// b matrix of element (1,0)

	a.minput();// each element of a matrix
	cin >> a; // and also available 

(b)input from file
		ifstream fin("datfile.dat");
		a.minput(fin);
		fin >> a;

(c)input method by using program
	//input/output form of each element
	a(0,0) = 3.;
	//input method by using elem function
	a.elem(0,0,3.);// a(0,0) = 3.; 
	//where both expression have the same meaning and are interchangeable.
// elem and () are almost same function, most usage () has advantage.
the elem function is useful following usage.
	cout << (A * B).elem(0,0) << endl;
	
Multiply matrix A and B and show the result of
	matrix element (0,0).() is invalid.
	Example
	cout << (A * B)(0,0) << endl; // invalid


	// << operator is also available for input data
	// input data should be following format are accepted.
	space and ';' will skip 
	Matrix<double> a1(2,2);
	a1 << "1 2 2 3"; // OK
	a1 << "1 2; 2 3"; // OK
	a1 << "1 2   ;   2 3"; // OK
	a1 << "1 2   ;2 3"; // OK
    a1 << "1 2;2 3"; // invalid
              ~~ To distinguish each data should make at least 1 space.
// complex matrix is also same procedure.
	Matrix<complex> ac(2,2);
    a1 << "1.566 2.677 2.898 3.765; 2 3 3 4"; // OK

    a1 << "1.566 2.677i 2.898 3.765i; 2 3i 3 4i"; // OK
    a1 << "1.566 i 2.898 3.765i; 2 3i 3 4i"; // Invalid
                 ~ at this time i recognize 0
    a1 << "1.566 1.0i 2.898 3.765i; 2 3i 3 4i"; // OK
                 ~~~~
	a1 << "1.566 2.677j 2.898 3.765j; 2 3i 3 4j"; // OK
    a1 << "1.566+2.677j 2.898 3.765j; 2 3i 3 4j"; // invalid
                ~ To distinguish each data should make at least 1 space.
// Multiple usage is not supported.
	Matrix<double> b1(1,2);
	b1 << "1" << "2"; // Result will be b1 = |2 0|
// Should be write
	b1 << " 1 2";


(c)Display of matrix 
		Matrix<double> a(2,2);
		...
		a.mprint();
		cout << a << endl;// both expression have the same meaning.
		
		The following is also acceptable.
		(a * a).mprint();// square of Matrix A is displayed.
		cout << (a * a);
		
		cerr << a << endl;
		a.mprint(cout);
		a.mprint(clog);
		a.mprint(cerr);
		
(d) Output to file
		Matrix<double> a(2,2);
		a << "1 2; 2 3";
		fstream fout("datfile.dat");
		a.mprint(fout); // same as fout << a;

output file will produce following data.
C:>type datfile.dat
1 2
2 3

(e) change matrix precision format
	default numeric format is defined "%8.4f" like as C printf format
	so a.mprint(cerr,8,4) is same function as a.mprint().
	a.mprint(cerr,10,2); // "%10.2f" format 
		
1.4 Supported Matrix operator

(a)Redefinition of Matrix size

	Matrix<double> A(2,2);
	Matrix<double> B(2,2);
	Matrix<double> C(2,4);
	...
	C = A * B;

The difference in size between the A*B and C matrix is obvious, but this program will not display the error. In such a situation, the size of Matrix C is cleared and redefined, to be the same size as that of matrix A*B. Memory management of the matrix allocation 
is already taking into account.
(In this example, the matrix size is 2 by 2)

(b)addition and subtraction: allowed on same size row and col.
	Matrix<double> A(2,2);
	Matrix<double> B(2,2);
	Matrix<double> C(2,3);
	Matrix<double> D,E,F,G;
	...
	D = A + B;//OK
	E = A + C;// Invalid because of matrix col size inequality

	F = - A;
	G = + A;// Matrix G will be the same size as that of matrix A

	A += B;// same  meaning like as A = A + B;
	A -= B;// same  meaning like as A = A - B;


(c)multiplication of Matrices: requires that the no. of rows in the first matrix
  is equal to the no. of columns in the second matrix.
	Matrix<double> A(2,1);
	Matrix<double> B(1,2);
	Matrix<double> D(1,3);
	Matrix<double> F(1,3);
	
	Matrix<double> C,E,F;
	...
	C = A * B; // Valid
	E = A * D; // Invalid
	F = A * B; // OK. This  expression is valid. 
	due to the redefinition of the matrix size.
			matrix F is now 2 by 2.

	complex number matrices can be used.
	Matrix<complex> A(2,1);
	Matrix<complex> B;
	...
	B = 3. * A;//OK
	B = complex(3.,0) * A; // same expression
	B = A * 3.;// This expression is legal in mathematics. 
				But in this program it will be detected as an error.

(d)matrix division: is similar to matrix multiplication.
	Matrix<double> A(2,2);
	Matrix<double> B(2,2);
	Matrix<double> C(2,3);
	Matrix<double> D,E,F,G,H;
	...
	D = A / B;
	E = A / C;// Invalid.
	F = 1. / A;// Valid.
	G = A / 2.;// Invalid.
	H = 1. / 2. * A ;//OK Instead of division by 2, A is multiplied by 1/2.

(f)transpose of matrix:
	~ is defined as the transpose of a matrix.

	Matrix<double> Ad(2,2);
	Matrix<double> Dd;
	...
	Dd = ~Ad; // transpose matrix

	Note that the transpose of a complex type matrix, becomes the
	conjugate transpose matrix.
	Matrix<complex> Ac(2,2);
	Matrix<complex> Dc;
	...
	Dc = ~Ac; // Matrix Dc becomes conjugate transpose matrix of Ac

(g)to view matrix row and col size:
	Matrix<double> A(2,2);
	cout << A.Getrow() << endl;
	cout << A.Getcol() << endl;

(h) union operator &
	This operator is different meaning from that of C syntax.
	union operator can unite 2 matrices to one matrix, for example 
	Matrix<double> A(2,3);
	Matrix<double> B(2,5);
	Matrix<double> D;
	D = A & B;
	becomes [D(2,3 + 5)] = [A(2,3) & B(2,5)] 
	& operator can apply only 2 matrices have same row size.
	
	Matrix<double> A(2,3);
	Matrix<double> B(2,5);
	Matrix<double> C(3,5);
	Matrix<double> D = A & B;// D = |A & B|
	Matrix<double> E = A & C;// invalid because of different row size.
    Matrix<double> F = ~(~B & ~C);// OK   F = | B |
                                    | & |
                                    | C |
(i)operator &=
	+=, -= operator are like as normal C syntax,
	&= operator is different meaning.
	&= operator is useful for data collection and data storage.
    Example of usage &= and +=
	Matrix<double> D(1,1);
	Matrix<double> A(1,1);
	Matrix<double> SUM;
	cin >> D;
	SUM = D;
	for(int i = 1;i < 10;i++) {
		cin >> A;
		D &= A;
		SUM += A;
	}
	D.mprint();//  D matrix becomes D(1,10)
	SUM.mprint(); // sum of D matrix (1x1);

(j)  operator != , ==
	Matrix<double> A(2,2);
	Matrix<double> B(2,2);
    *
     if(A == B) *
     if(A != B) *
	==,!= operator check the matrix size as well as inside values of matrix.
	
Logical matrix calculation mistakes(size mismatch) are not checked by the compiler.
Therefore errors may appear in program execution. 
DEBUG option may help to find the some mistakes. 

C++ operator can also be redefined as *=,/=,|,++,--. However these redefinition seems to me no meaning in matrix operations. Therefore, these are NOT supported in this matrix calculation library.

3. sample program

The following programs are sample applications of the matrix class library.
They are only given as examples to demonstrate matrix class usage and these
 are programmed in Borland C++.

Mtest.cpp //example of matrix operation (as shown in appendix A)
Eigenvec.cpp // example of eigen value and eigen vector operation (this is not included in here, if you are interested, please send me E-mail)
Trany.cpp //May be used to simulate a state variable equation.
		this simulation is used by employing pade' method.
Runge.cpp //This is the C++ version of Runge Kutta method. 
		A large matrix size of non-linear state variable equations can be easily calculated.

4. Support
This program is simple and tiny, but It is very useful for digital signal
 processing and simulation as well as real-time signal processing.
If you find any bugs in this program, I will be grateful if you sent them to
me via E-mail.

5. Conclusion
This class library was developed for matrices and it is a part of C++'s
useful function.
The matrix class library is intended to simplify simulation and data 
processing, and to be used as an introduction to C++ programming.

				kobayash@vela.acs.oakland.edu(Oakland Univ)
				PXD00173@niftyserve.or.jp (Niftyserve)
				ikko@wtnb.sc.hosei.ac.jp (Hosei Univ)
				Kazuyuki Kobayashi
				1994 Gerox(c)


------------------------------------------------------------------------------
Appendix A1 executed results of mtest.cpp 
------------------------------------------------------------------------------
A:\matrix>mtest


 incremental data storage test
 Input data 
Matrix[0,0] = 1
Matrix[0,0] = 2
Dd Matrix [1]
 |  1.0000   2.0000 |
Matrix[0,0] = 4
Dd Matrix [2]
 |  1.0000   2.0000   4.0000 |
Matrix[0,0] = 5
Dd Matrix [3]
 |  1.0000   2.0000   4.0000   5.0000 |
Matrix[0,0] = 6
Dd Matrix [4]
 |  1.0000   2.0000   4.0000   5.0000   6.0000 |
push any key
Matrix D(1,5)
 |  1.0000   2.0000   4.0000   5.0000   6.0000 |
 Advanced Usage of Union operator
 E(3,2) = |A(1,1)|C(1,2)|
          |B(2,1)|D(1,1)|
push any key
A Matrix
 |  1.0000 |
B Matrix
 |  3.0000 |
 |  5.0000 |
C Matrix
 |  2.0000 |
 |  4.0000 |
D Matrix
 |  6.0000 |
E Matrix ~(~A&~B)&~(~C&~D)
 |  1.0000   2.0000 |
 |  3.0000   4.0000 |
 |  5.0000   6.0000 |
push any key
Complex type matrix demonstration
Input matrix size:2
Matrix[0,0]:
real part = 1
imag part = 6
Matrix[0,1]:
real part = 2
imag part = 8
Matrix[1,0]:
real part = 9
imag part = 5
Matrix[1,1]:
real part = 2
imag part = 5
Original Matrix(complex type): cout << a1
 |  1.0000+j  6.0000   2.0000+j  8.0000 |
 |  9.0000+j  5.0000   2.0000+j  5.0000 |
Inverse Matrix(complex type): cout << (1./a1)
 | -0.0791+j  0.0235   0.1249+j -0.0192 |
 |  0.0889+j -0.1303  -0.0929+j  0.0068 |
Unit Matrix(complex type): cout << (a1 * (1. / a1))
 |  1.0000+j  0.0000   0.0000+j  0.0000 |
 |  0.0000+j  0.0000   1.0000+j  0.0000 |

 | 1.0+j 0.0  0.0+j 0.0 |
 | 0.0+j 0.0  1.0+j 0.0 |

push any key
Double type matrix demonstration
Matrix[0,0] = 2
Matrix[0,1] = 3
Matrix[1,0] = 56
Matrix[1,1] = 6
Original Matrix(double type): cout << a1d
 |  2.0000   3.0000 |
 | 56.0000   6.0000 |

Inverse Matrix(double type): cout << (1./a1d)
 | -0.0385   0.0192 |
 |  0.3590  -0.0128 |
Unit Matrix(double type): cout << (a1d / ad1)
 |  1.0000   0.0000 |
 |  0.0000   1.0000 |
Eigenvalue: cout << eig.eigen(ILL)
 | -9.1149   0.0000 |
 | 17.1149   0.0000 |

------------------------------------------------------------------------------
Appendix A2 mtest.cpp
------------------------------------------------------------------------------
// Programed by Gerox(c) 1994.12.27
// Matrix test program
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <complex.h>

#ifdef __BORLANDC__
	#include <bcd.h>
#endif

#include "matrixt.h"
#include "eigenvec.h"
void main(void) {
	cout << "\n incremental data storage test" << endl; 
	Matrix<double> Dd(1,1);
	Matrix<double> Ad(1,1);
	Matrix<double> SUMd(1,1);
	cout << " Input data " << endl;
	cin >> Dd;
	for(int i = 1;i < 5;i++) {
		cin >> Ad;
		Dd &= Ad;
		SUMd += Ad;
		cout << "Dd Matrix [" << i << "]" << endl;
		Dd.mprint();
	}
	cout << "push any key"<< endl;
	getch();
	cout << "Matrix D(" << Dd.Getrow()<< "," << Dd.Getcol() <<")"<< endl;
	Dd.mprint();
	getch();
	cout << endl;
	cout << " Advanced Usage of Union operator" << endl;
	cout << " E(3,2) = |A(1,1)|C(1,2)|" << endl;
	cout << "          |B(2,1)|D(1,1)|" << endl;
	cout << "push any key"<< endl;
	Matrix<double> A(1,1);A(0,0) = 1;
	cout << "A Matrix" << A;
	Matrix<double> B(2,1);B(0,0) = 3;B(1,0) = 5;
	cout << "B Matrix" << B;
	Matrix<double>  C(2,1);C(0,0) = 2;C(1,0) = 4;
	cout << "C Matrix" << C;
	Matrix<double>  D(1,1);D(0,0) = 6;
	cout << "D Matrix" << D;
	Matrix<double> E = ~(~A&~B)&~(~C&~D);
	cout << "E Matrix ~(~A&~B)&~(~C&~D)";
	E.mprint();
	cout << "push any key"<< endl;
	cout << "Complex type matrix demonstration" << endl;
	getch();
	cout << "Input matrix size:";
	int n;
	cin >> n;
	Matrix<complex> a1(n,n);
	cin >> a1;
	cout << "Original Matrix(complex type): cout << a1\n";
	cout << a1;
	cout << "Inverse Matrix(complex type): cout << (1./a1)\n";
	cout << (1./a1);
	cout << "Unit Matrix(complex type): cout << (a1 * (1. / a1))\n";
	cout << (a1 * (1. / a1));
	(a1/a1).mprint(cerr,4,1);
	cout << "push any key\n";
	cout << "Double type matrix demonstration" << endl;
	getch();
	Matrix<double>  a1d(n,n);
	cin >> a1d;
	cout << "Original Matrix(double type): cout << a1d\n";
	cout << a1d;
	cout << "Inverse Matrix(double type): cout << (1./a1d)\n";
	cout << (1./a1d);
	cout << "Unit Matrix(double type): cout << (a1d / ad1)\n";
	cout << (a1d / a1d);
	Eigen eig(a1d);
	Matrix<double> ILL;
	cout << "Eigenvalue: cout << eig.eigen(ILL)\n";
	cout << eig.eigen(ILL);

#if defined(__BCD_H)
	cout << "Bcd type matrix demonstration" << endl;
	cout << "push any key\n";
	getch();
	Matrix<bcd> a1bcd(n,n);
	cin >> a1bcd;
	cout << "Original Matrix(bcd type): a1bcd.mprint(coutm10,5)\n";
	a1bcd.mprint(cout,10,5);
	cout << "Inverse Matrix(bcd type): (1./a1bcd).mprint(cout,10,5)\n";
	(1./a1bcd).mprint(cout,10,5);
	cout << "Unit Matrix(bcd type): (a1bcd / a1bcd).mprint(cout,10,5)\n";
	(a1bcd / a1bcd).mprint(cout,10,5);
#endif
}
-----------------------------------------------------------------------
end of mtest.cpp
-----------------------------------------------------------------------


-----------------------------------------------------------------------
Appendix B matrixt.h
-----------------------------------------------------------------------
// Programmed by Gerox(c) 1993.12.5 for C++ Version
// template version double and complex 1994.12.22
// template class float, double, complex, bcd
// Borland C++ Ver.4, TurboC++ Ver.3, Watcom C++ Ver9.5
// iostream support
// speed optimize and (union)| operator support 1995.01.02
// Kazuyuki Kobayashi E-mail : ikko@oakland.edu, ikko@wtnb.sc.hosei.ac.jp
#if !defined(GEROX_MATRIXT_H)
#define GEROX_MATRIXT_H
#endif
#if !defined(__IOSTREAM_H) || !defined(_IOSTREAM_H_INCLUDED)
#include <iostream.h>
#endif

// If you want to check memory allocation.
#define DEBUG

#define LU

#if !defined(_CLASSTYPE)
// Compatibility for Watcom C
#define _CLASSTYPE
#endif

template <class T> class _CLASSTYPE Matrix;

#if defined(__MATH_H) || defined(_MATH_H_INCLUDED)
//////////////////////////////////////////////////////
//typedef Matrix<float> Matf; 
typedef Matrix<double> Matd;
//typedef Matrix<int> Mati;
/////////////////////////////////////////////////////
#endif
#if defined(__BCD_H)
/////////////////////////////////////////////////////
typedef Matrix<bcd> Matbcd;
/////////////////////////////////////////////////////
#endif
#if defined(__COMPLEX_H) || defined(_COMPLEX_H_INCLUDED)
/////////////////////////////////////////////////////
typedef Matrix<complex> Matc;
/////////////////////////////////////////////////////
#endif

#define print cerr<<__FILE__<<" "<<__LINE__<<"\n"

const double eps=1e-10;
template <class T> class _CLASSTYPE Matrix
{
private:
	void dimcheck(int i,int j,int k) const;
	void operr(const Matrix<T>&a,char* op) const;
protected:
	T* m;
	int row;
	int col;
public:
	Matrix(void);
	Matrix(const Matrix<T>& a);
	Matrix(int row,int col);
	~Matrix(void) { delete m;}

	void minput(istream& o = cin);
	void minput(int i,int j);
	void mprint(ostream& o = cerr,int width = 8,int point = 4) const;
	double ffabs(const T& val) const;
	Matrix<T> operator+(void) const;
	Matrix<T> operator-(void) const;
	Matrix<T>& operator =(const Matrix<T>& a);
	Matrix<T>& operator +=(const Matrix<T>& a);
	Matrix<T>& operator -=(const Matrix<T>& a);
	Matrix<T>& operator &=(const Matrix<T>& a);
	friend int operator ==(const Matrix<T>& a,const Matrix<T>& b);
	friend int operator !=(const Matrix<T>& a,const Matrix<T>& b);
	friend Matrix<T> operator &(const Matrix<T>& a,const Matrix<T>& b);
	friend Matrix<T> operator *(const T& t,const Matrix<T>& a);
	friend Matrix<T> operator *(const Matrix<T>& a,const Matrix<T>& b);
	friend Matrix<T> operator +(const Matrix<T>& a,const Matrix<T>& b);
	friend Matrix<T> operator -(const Matrix<T>& a,const Matrix<T>& b);
	friend Matrix<T> operator /(const T& t,const Matrix<T>& b);
	friend Matrix<T> operator /(const Matrix<T>& a,const Matrix<T>& b);
	friend Matrix<T> operator ~(const Matrix<T>& a);
	friend ostream& operator<<(ostream& o, Matrix<T>& t);
	friend istream& operator>>(istream& o, Matrix<T>& t);
	friend Matrix<T>& operator<<(Matrix<T>& t, char* line);
	

	int Getrow(void) const {return row;}
	int Getcol(void) const {return col;}
#if !defined(DEBUG)
	T& operator ()(int i,int j) const { return m[col * i + j];}
	T elem(int i,int j) const { return m[col * i + j];}
	T elem(int i,int j,T t) const { return (m[col * i + j] = t);}
#else
	T& operator ()(int i,int j) const 
	{
		dimcheck(i,j,0);
		return m[col * i + j];
	}
	T elem(int i,int j) const
	{
		dimcheck(i,j,1);
		return m[col * i + j];
	}
	T elem(int i,int j,T t) const 
	{
		dimcheck(i,j,2);
		return (m[col * i + j] = t);
	}
#endif
};
template <class T> void Matrix<T>::dimcheck(int i,int j,int k) const 
{
	if((i >= row) || (j >= col) || (i < 0) || (j < 0)) {
		cerr << "////////////////////////////////////////////////" << endl;
		cerr << "dimension over "<< k << "\n" << endl;
		cerr << "0:(),1:elem(i,j),2:elem(i,j,t)\n";
		cerr << "i = " << i << " j = " << j << endl;
		mprint(cerr);
		cerr << "////////////////////////////////////////////////" << endl;
		exit(1);
	}
}
template <class T> void Matrix<T>::operr(const Matrix<T>& a,char* op) const 
{
		cerr << "////////////////////////////////////////////////" << endl;
		cerr << "Mat("<<row<<","<<col<<")";
		mprint(cerr);
		cerr << "Mat("<<a.row<<","<<a.col<<")";
		a.mprint(cerr);
		cerr << "Syntax Error operator " << op << " (Mat("; 
		cerr <<row<<","<<col<<"),Mat("<<a.row<<","<<a.col<<"))\n";
		cerr << "////////////////////////////////////////////////" << endl;
		exit(1);
}
#if defined(__COMPLEX_H) || defined(_COMPLEX_H_INCLUDED)
double Matrix<complex>::ffabs(const complex& val) const
{ 
		complex tmp = val;
		return abs(tmp);
}
#endif
#if defined(__BCD_H)
double Matrix<bcd>::ffabs(const bcd& val) const
{ 
		bcd tmp = val;
		return real(abs(tmp));
}
#endif
template <class T> double Matrix<T>::ffabs(const T& val) const
{
	return fabs((double)val);
}
template <class T> void Matrix<T>::minput(int i,int j)
{
	char dummy[256];
	double t;
	cin >> dummy;
	t = atof(dummy);
	(*this).elem(i,j,(T)t);
}
#if defined(__COMPLEX_H) || defined(_COMPLEX_H_INCLUDED)
void Matrix<complex>::minput(int i,int j)
{
	char dummy[256];
	complex t;
	cerr << endl;
	cerr << "real part = ";
	cin >>dummy;
	double tmp = atof(dummy);
	cerr << "imag part = ";
	cin >> dummy;
	t = complex(tmp,atof(dummy));
	(*this).elem(i,j,(complex)t);
}
#endif
#if defined(__COMPLEX_H) || defined(_COMPLEX_H_INCLUDED)
void Matrix<complex>::minput(istream& o)
{
	char dummy[256];
	for(int i = 0;i < row;i++) {
		for(int j = 0;j < col;j++) {
			cerr << "Matrix["<<i<<","<<j<<"]:" << endl;
			complex t;
			cerr <<"real part = ";
			o >> dummy;
			double tmp = atof(dummy);
			cerr <<"imag part = ";
			o >> dummy;
			t = complex(tmp,atof(dummy));
			m[i * col + j] = t;
		}
	}
}
#endif
template <class T> void Matrix<T>::minput(istream& o)
{
	char dummy[256];
	for(int i = 0;i < row;i++) {
		for(int j = 0;j < col;j++) {
			cerr << "Matrix["<<i<<","<<j<<"] = ";
			double t;
			o >> dummy;
			t = atof(dummy);
			m[i * col + j] = (T)t;
		}
	}
}
template <class T> Matrix<T>::Matrix(int row,int col)
{
	Matrix<T>::row = row;
	Matrix<T>::col = col;
	m = new T[row * col];
	for(register int i = 0;i < row;i++){
		for(register int j = 0;j < col;j++) {
			m[i * col + j] = (T)0.0;
		}
	}
}
template <class T> Matrix<T>::Matrix(const Matrix<T>& a)
{
	Matrix<T>::row = a.row;
	Matrix<T>::col = a.col;
	m = new T[row * col];
	for(register int i = 0;i < row;i++){
		for(register int j = 0;j < col;j++) {
			m[col * i + j] = a.m[a.col * i + j];
		}
	}
}
template <class T> Matrix<T>::Matrix(void)
{
	Matrix<T>::row = 0;
	Matrix<T>::col = 0;
	m = new T[0];
}
template<class T> Matrix<T> operator +(const Matrix<T>& a,const Matrix<T>& b)
{
	Matrix<T> c(a.row,a.col);
	if ((a.row != b.row) || (a.col != b.col)) a.operr(b,"+");
	else {
		for(register int i = 0;i < c.row;i++) {
			for(register int j = 0;j < c.col;j++) {
				c.m[c.col * i + j] = a.m[a.col * i + j] + b.m[b.col * i + j];
			}
		}
	}
	return c;
}
template <class T> Matrix<T> operator -(const Matrix<T>& a,const Matrix<T>& b)
{
	Matrix<T> c(a.row,a.col);
	if ((a.row != b.row) || (a.col != b.col)) a.operr(b,"-");
	else {
		for(register int i = 0;i < c.row;i++) {
			for(register int j = 0;j < c.col;j++) {
				c.m[c.col * i + j] = a.m[a.col * i + j] - b.m[b.col * i + j];
			}
		}
	}
	return c;
}
template <class T> Matrix<T>& Matrix<T>::operator =(const Matrix<T>& a)
{
	row = a.row;
	col = a.col;
	delete m;
	m = new T[row * col];
	for(register int i = 0;i < a.row;i++) {
		for(register int j = 0;j < a.col;j++) {
			m[col * i + j] = a.m[a.col * i + j];
		}
	}
	return *this;
}
template <class T> Matrix<T>& Matrix<T>::operator +=(const Matrix<T>& a)
{
	if ((a.row != row) || (a.col != col)) operr(a,"+=");
	for(register int i = 0;i < row;i++) {
		for(register int j = 0;j < col;j++) {
			m[col * i + j] += a.m[a.col * i + j];
		}
	}
	return *this;
}
template <class T> Matrix<T>& Matrix<T>::operator -=(const Matrix<T>& a)
{
	if ((a.row != row) || (a.col != col)) operr(a,"-=");
	for(register int i = 0;i < row;i++) {
		for(register int j = 0;j < col;j++) {
			m[col * i + j] -= a.m[a.col * i + j];
		}
	}
	return *this;
}
template <class T> Matrix<T>& Matrix<T>::operator &=(const Matrix<T>& a)
{
	if (a.row != row) operr(a,"&=");
	Matrix<T> tmp = (*this);
	delete m;
	row = tmp.row;
	col = tmp.col + a.col;
	m = new T[row * col];
	for(register int i = 0;i < row;i++){
		for(register int j = 0;j < tmp.col;j++) {
			m[col * i + j] = tmp.m[tmp.col * i + j];
		}
		for(j = tmp.col;j < col;j++) {
			m[col * i + j] = a.m[a.col * i + j - tmp.col];
		}
	}
	return *this;
}
template <class T> Matrix<T> Matrix<T>::operator +(void) const
{
	return *this;
}
template <class T> Matrix<T> Matrix<T>::operator -(void) const
{
	Matrix<T> tmp = *this;
	for(register int i = 0;i < row;i++) {
		for(register int j = 0;j < col;j++) {
			tmp(i,j) = -tmp(i,j);
		}
	}
	return tmp;
}
template <class T> Matrix<T> operator *(const T& t,const Matrix<T>& a)
{
	Matrix<T> b(a.row,a.col);
	for(register int i = 0;i < a.row;i++) {
		for(register int j = 0;j < a.col;j++) {
			b.m[b.col * i + j] = (T)t * a.m[a.col * i + j];
		}
	}
	return b;
}
template <class T> Matrix<T> operator *(const Matrix<T>& a,const Matrix<T>& b)
{
	Matrix<T> c(a.row,b.col);
	if (a.col != b.row) a.operr(b,"*");
	for(register int i = 0;i < a.row;i++) {
		for(register int j = 0;j < b.col;j++) {
			T d = (T)0.;
			for (register int k = 0;k < a.col;k++) {
				d += (a.m[a.col * i + k] * b.m[b.col * k + j]);
			}
			c.m[c.col * i + j] = d;
		}
	}
	return c;
}
template <class T> Matrix<T> operator &(const Matrix<T>& a,const Matrix<T>& b)
{
	Matrix<T> c(a.row,a.col + b.col);
	if (a.row != b.row) a.operr(b,"&");
	else {
		for(register int i = 0;i < a.row;i++) {
			for(register int j = 0;j < a.col;j++) {
				c.m[c.col * i + j] = a.m[a.col * i + j];
			}
			for(j = a.col;j < a.col + b.col;j++) {
				c.m[c.col * i + j] = b.m[b.col * i + j - a.col];
			}
		}
	}
	return c;
}
template <class T> Matrix<T> operator ~(const Matrix<T>& a)
{
	Matrix<T> c(a.col,a.row);
	for(register int i = 0;i < c.row;i++) {
		for(register int j = 0;j < c.col;j++) {
			c.m[c.col * i + j] = a.m[a.col * j + i];
		}
	}
	return c;
}
#if defined(__COMPLEX_H) || defined(_COMPLEX_H_INCLUDED)
Matrix<complex> operator ~(const Matrix<complex>& a)
{
	Matrix<complex> c(a.col,a.row);
	for(register int i = 0;i < c.row;i++) {
		for(register int j = 0;j < c.col;j++) {
			c.m[c.col * i + j] = conj(a.m[a.col * j + i]);
		}
	}
	return c;
}
#endif
template <class T> Matrix<T> operator /(const Matrix<T>& a,const Matrix<T>& b)
{
	Matrix<T> c = 1.0 / b;
	return a * c;
}
#ifdef LU
template <class T> Matrix<T> operator /(const T& t,const Matrix<T>& b)
{
	if(b.col != b.row) {
		Matrix<T> tmp(1,1);
		tmp(0,0) = t;
		tmp.operr(b,"/");
	}
	int j;
	int n = b.Getrow();
	int *ip = new int[n];
	double ttt,u;
	double *weight = new double[n];
	Matrix<T> a = b;
	Matrix<T> a_inv = b;
	T tt;
	T det = 0;
	for(int i = 0;i < n;i++) {
		weight[i] = 0.;
		ip[i] = i;
	}
	for(int k = 0;k < n;k++) {
		u = 0.;
		for(j = 0;j < n;j++) {
			ttt = a.ffabs(a(k,j));
			if(ttt > u) u = ttt;
		}
		if(a.ffabs(u) < eps) goto EXIT;
		weight[k] = 1. / u;
	}
	det = 1.;
	for(k = 0;k < n;k++) {
		u = -1.;
		for(i = k;i < n;i++) {
			int ii = ip[i];
			ttt = a.ffabs(a(ii,k)) * weight[ii];
			if(ttt > u) {
				u = ttt;
				j = i;
			}
			int ik = ip[j];
			if(j != k) {
				ip[j] = ip[k];
				ip[k] = ik;
				det = -det;
			}
			T uu = a(ik,k);
			det *= uu;
			if(a.ffabs(uu) < eps) goto EXIT;
			for(i = k + 1;i < n;i++) {
				int ii = ip[i];
				a(ii,k) /= uu;
				tt = a(ii,k);
				for(j = k + 1;j < n;j++) {
					a(ii,j) -= tt * a(ik,j);
				}
			}
		}
	}
EXIT:
	delete weight;
	if(a.ffabs(det) != 0) {
		for(k = 0;k < n;k++) {
			for(i = 0;i < n;i++) {
				int ii = ip[i];
				if(ii == k) tt = 1.;
				else tt = 0;
				for(j = 0;j < i;j++) {
					tt -= a(ii,j) * a_inv(j,k);
				}
				a_inv(i,k) = tt;
			}
			for(i = n - 1;i >= 0;i--) {
				tt = a_inv(i,k);
				int ii = ip[i];
				for(j = i + 1;j < n;j++) {
					tt -= a(ii,j) * a_inv(j,k);
				}
				a_inv(i,k) = tt / a(ii,i);
			}
		}
	} else cerr << "\n Warning Matrix calculation may not accurate.\n";
	delete ip;
	return t * a_inv;
}
#else
template <class T> Matrix<T> operator /(const T& t,const Matrix<T>& b)
{
	Matrix<T> c = b;
	T p,aik;
	if(b.col != b.row) {
		Matrix<T> tmp(1,1);
		tmp(0,0) = t;
		tmp.operr(b,"/");
	}
	int n = c.row;
	for (register int k = 0;k < n;k++) {
		p = c.m[c.col * k + k];
		c.m[c.col * k + k] = 1.0;
		for (register int j = 0;j < n ; j++) {
			if(c.ffabs(p) < eps) p = eps;
			c.m[c.col * k + j] /= p;
		}
		for (register int i = 0;i < n;i++) {
			if (i != k) {
				aik = c.m[c.col * i + k];
				c.m[c.col * i + k] = 0.0;
				for (j = 0;j < n;j++) {
					c.m[c.col * i + j] -= (aik * c.m[c.col * k + j]);
				}
			}
		}
	}
	return t * c;
}
#endif
template <class T> void Matrix<T>::mprint(ostream& o,int width,int point) const
{
	o << flush;
	cerr << flush;
	cerr << endl;
	for(int i = 0;i < row; i++) {
		cerr <<" |";
		for(int j = 0;j < col;j++) {
			o.width(width);
			o.precision(point);
			o.flags(cout.flags() | ios::fixed | ios::showpoint);
			o << m[col * i + j];
			o << " ";
		}
		cerr <<"|";
		o << "\n";
	}
	o << "\n";
}
#if defined(__COMPLEX_H) || defined(_COMPLEX_H_INCLUDED)
void Matrix<complex>::mprint(ostream& o,int width,int point) const
{
	o << flush;
	cerr << flush;
	cerr << endl;
	for(int i = 0;i < row; i++) {
		cerr <<" |";
		for(int j = 0;j < col;j++) {
			o.width(width);
			o.precision(point);
			o.flags(cout.flags() | ios::fixed | ios::showpoint);
			o << real(m[col * i + j]);
			o << "+j";
			o.width(width);
			o.precision(point);
			o.flags(cout.flags() | ios::fixed | ios::showpoint);
			o << imag(m[col * i + j]);
			o << " ";
		}
		cerr << "|";
		o << "\n";
	}
	o << "\n";
}
#endif
template <class T> ostream & operator<<(ostream & o, Matrix<T> & t)
{
	t.mprint(o);
	return o;
}
template <class T> istream & operator>>(istream & o, Matrix<T> & t)
{
	t.minput(o);
	return o;
}
template <class T> Matrix<T>& operator<<(Matrix<T>& t,char* line)
{
	char *p = line;
	char dummy2[256];
	for(int i = 0;i < t.row * t.col;i++,p += strlen(dummy2)) {
		while(*p == ' ' || *p == ';') p++;
		if (sscanf(p,"%s",dummy2) == EOF) {
			cerr << "Warning: data is not enough"<< endl;
			break;
		}
		T tmp = (T) atof(dummy2);
		t.m[i] = tmp;
	}
	return t;
}
#if defined(__COMPLEX_H) || defined(_COMPLEX_H_INCLUDED)
Matrix<complex>& operator<<(Matrix<complex>& t,char* line)
{
	char *p = line;
	char dummy2[256];
	for(int i = 0;i < t.row * t.col;i++,p += strlen(dummy2)) {
		while(*p == ' ') p++;
		if (sscanf(p,"%s",dummy2) == EOF) {
			cerr << "Warning: data is not enough"<< endl;
			break;
		}
		double tmp1 = atof(dummy2);
		print << p;
		p += strlen(dummy2);
		while(*p == ' ' || *p == ';') p++;
		if (sscanf(p,"%s",dummy2) == EOF) {
			cerr << "Warning: data is not enough"<< endl;
			break;
		}
		complex tmp = complex(tmp1,atof(dummy2));
		t.m[i] = tmp;
		print << p;
	}
	return t;
}
#endif
template <class T> int operator ==(const Matrix<T>& a,const Matrix<T>& b)
{
	if((a.col != b.col) || (a.col != b.col)) {
		return 0;
	}
	for(register int i = 0;i < a.row;i++) {
		for(register int j = 0;j < a.col;j++) {
			T tmp = a.m[a.col * i + j] - b.m[b.col * i + j];
			if(a.ffabs(tmp) > eps) return 0;
		}
	}
	return 1;
}
template <class T> int operator !=(const Matrix<T>& a,const Matrix<T>& b)
{
	if(a == b) return 0;
	else return 1;
}
-----------------------------------------------------------------------
end of matrixt.h
-----------------------------------------------------------------------
This is previous version of matrix library which I wrote.
Non-template type support complex number 
-----------------------------------------------------------------------
// Programed by Gerox(c) 1993.12.5 for C++ Version
// Complex Version
class Matrix
{
protected:
	complex *m;
	int row;
	int col;
public:
	Matrix(int row,int col);
	Matrix(Matrix& a);
	Matrix(void);
	~Matrix() { free(m);}
	void minput(void);
	void minput(int i,int j);
	void cminput(void);
	void cminput(int i,int j);
	Matrix& operator =(const Matrix& a);
	friend Matrix operator *(const complex& t,const Matrix& a);
	friend Matrix operator *(const Matrix& a,const Matrix& b);
	friend Matrix operator +(const Matrix& a,const Matrix& b);
	friend Matrix operator -(const Matrix& a,const Matrix& b);
	friend Matrix operator /(const complex& t,const Matrix& b);
	friend Matrix operator /(const Matrix& a,const Matrix& b);
	friend Matrix operator ~(const Matrix& a);
	complex& operator ()(int i,int j) { return m[col * i + j];}
	void mprint(void);
	void cmprint(void);
	double elem(int i,int j) { return real(m[col * i + j]);}
	complex celem(int i,int j) { return m[col * i + j];}
	complex elem(int i,int j,complex t) { return (m[col * i + j] = t);}
	int Getrow(void) {return row;}
	int Getcol(void) {return col;}

};

void Matrix::minput(int i,int j)
{
	if(i >= row || j >= col) {
		fprintf(stderr,"%d %d %d %d\n",i,j,col,row);
		exit(1);
	}
	char dummy[256];
	scanf("%s",dummy);
	complex t = complex(atof(dummy),0);
	m[col * i + j] = t;
}
void Matrix::cminput(int i,int j)
{
	if(i >= row || j >= col) {
		fprintf(stderr,"%d %d %d %d\n",i,j,col,row);
		exit(1);
	}
	char dummy[256];
	fprintf(stderr,"Input real part =");
	scanf("%s",dummy);
	double tr = atof(dummy);
	fprintf(stderr,"Input imag part =");
	scanf("%s",dummy);
	complex t = complex(tr,atof(dummy));
	m[col * i + j] = t;
}
Matrix::Matrix(int row,int col)
{
	if ((m = (complex *)malloc(sizeof(complex) * row * col)) == NULL) {
		fprintf(stderr,"%^%H%j%/%9$N%a%b%j$,3NJ]$G$-$^$;$s !!");
		fprintf(stderr,"matrix(int,int) row = %d col = %d",row,col);
		exit(0);
	}
	Matrix::row = row;
	Matrix::col = col;
	for(int i = 0;i < row;i++){
		for(int j = 0;j < col;j++) m[i * col + j] = complex(0,0);
	}
}
Matrix::Matrix(Matrix& a)
{
	Matrix::row = a.row;
	Matrix::col = a.col;
	if ((m = (complex *)malloc(sizeof(complex) * row * col)) == NULL) {
		fprintf(stderr,"%^%H%j%/%9$N%a%b%j$,3NJ]$G$-$^$;$s !!");
		fprintf(stderr,"matrix(a) row = %d col = %d",row,col);
		exit(0);
	}
	for(int i = 0;i < row;i++){
		for(int j = 0;j < col;j++) m[col * i + j] = a.m[a.col * i + j];
	}
}
Matrix::Matrix(void)
{
	row = 1;
	col = 1;
	if ((m = (complex *)malloc(sizeof(complex) * row * col)) == NULL) {
		fprintf(stderr,"%^%H%j%/%9$N%a%b%j$,3NJ]$G$-$^$;$s !!");
		fprintf(stderr,"matrix(void)row = %d col = %d",row,col);
		exit(0);
	}
	for(int i = 0;i < row;i++){
		for(int j = 0;j < col;j++) m[i * col + j] = complex(0,0);
	}
}

void Matrix::cminput(void)
{
	for(int i = 0;i < row;i++) {
		for(int j = 0;j < col;j++) {
			fprintf(stderr,"Matrix[%d,%d] = ",i,j);
			char dummy[256];
			fprintf(stderr,"Input real part =");
			scanf("%s",dummy);
			double tr = atof(dummy);
			fprintf(stderr,"Input imag part =");
			scanf("%s",dummy);
			complex t = complex(tr,atof(dummy));
			m[i * col + j] = t;
		}
	}
}
void Matrix::minput(void)
{
	for(int i = 0;i < row;i++) {
		for(int j = 0;j < col;j++) {
			fprintf(stderr,"Matrix[%d,%d] = ",i,j);
			double val;
			scanf("%lg",&val);
			m[i * col + j] = complex(val,0);
		}
	}
}

Matrix operator +(const Matrix& a,const Matrix& b)
{
	Matrix c(a.row,a.col);

	if ((a.row != b.row) || (a.col != b.col)) {
		fprintf(stderr,"Not add !!\n");
		exit(1);
	} else {
		c.row = a.row;
		c.col = a.col;
		for(int i = 0;i < c.row;i++) {
			for(int j = 0;j < c.col;j++) {
				c.m[c.col * i + j] = a.m[a.col * i + j] + b.m[b.col * i + j];
			}
		}
	}
	return c;
}
Matrix operator -(const Matrix& a,const Matrix& b)
{
	Matrix c(a.row,a.col);

	if ((a.row != b.row) || (a.col != b.col)) {
		fprintf(stderr,"Not add !!\n");
		exit(1);
	} else {
		c.row = a.row;
		c.col = a.col;
		for(int i = 0;i < c.row;i++) {
			for(int j = 0;j < c.col;j++) {
				c.m[c.col * i + j] = a.m[a.col * i + j] - b.m[b.col * i + j];
			}
		}
	}
	return c;
}
Matrix& Matrix::operator =(const Matrix& a)
{
	int i,j;
	if(a.row != row || a.col != col) {
		free(m);
		row = a.row;
		col = a.col;
		if ((m = (complex *)malloc(sizeof(complex) * row * col)) == NULL) {
			fprintf(stderr,"%^%H%j%/%9$N%a%b%j$,3NJ]$G$-$^$;$s !!");
			fprintf(stderr,"operator = row = %d col = %d",row,col);
			exit(0);
		}
	}
	for(i = 0;i < a.row;i++){
		for(j = 0;j < a.col;j++)
			m[col * i + j] = a.m[a.col * i + j];
	}
	return *this;
}
Matrix operator *(const complex& t,const Matrix& a)
{
	Matrix b(a.row,a.col);
	if(a.row == b.row && a.col == b.col) {
	for(int i = 0;i < a.row;i++){
		for(int j = 0;j < a.col;j++)
			b.m[b.col * i + j] = (complex)t * a.m[a.col * i + j];
		}
	}
	else fprintf(stderr,"2C;;$G$-$^$;$s!*\n");
	return b;
}
Matrix operator ~(const Matrix& a)
{
	Matrix c(a.col,a.row);
	for(int i = 0;i < c.row;i++) {
		for(int j = 0;j < c.col;j++) {
			c.m[c.col * i + j] = conj(a.m[a.col * j + i]);
		}
	}
	return c;
}
Matrix operator *(const Matrix& a,const Matrix& b)
{
	Matrix c(a.row,b.col);
	int i,j,k;

	if (a.col == b.row) {
		for(i = 0;i < a.row;i++) {
			for(j = 0;j < b.col;j++) {
				complex d = 0;
				for (k = 0;k < a.col;k++)
					d += (a.m[a.col * i + k] * b.m[b.col * k + j]);
					c.m[c.col * i + j] = d;
			}
		}
	} else {
		fprintf(stderr,"Not multiple");
		exit(1);
	}
	return c;
}
Matrix operator /(const Matrix& a,const Matrix& b)
{
	Matrix c = 1 / b;
	return a * c;
}
Matrix operator /(const complex& t,const Matrix& b)
{
	Matrix c;
	c = b;
	int i,j,k;
	complex p,aik;
	if(c.col != c.row) {
		fprintf(stderr,"Not divide col = %d row = %d",c.col,c.row);
		exit(1);
	}
	int n = c.row;
	for (k = 0;k < n;k++) {
		p = c.m[c.col * k + k];
		c.m[c.col * k + k] = 1.0;
		for (j = 0;j < n ; j++) {
			if(abs(p) < 0.0001) p = 0.0001;
			c.m[c.col * k + j] /= p;
		}
		for (i = 0;i < n;i++) {
			if (i != k) {
				aik = c.m[c.col * i + k];
				c.m[c.col * i + k] = 0.0;
				for (j = 0;j < n;j++) {
					c.m[c.col * i + j] -= (aik * c.m[c.col * k + j]);
				}
			}
		}
	}
	return t * c;
}

void Matrix::cmprint(void)
{
	int i,j;
	for(i = 0;i < row; i++) {
		fprintf(stderr," | ");
		for(j = 0;j < col;j++) {
			fprintf(stderr,"%8.4f+j%8.4f ",real(m[col * i + j]),imag(m[col * i + j]));
		}
		fprintf(stderr," |\n");
	}
	fprintf(stderr,"\n");
}
void Matrix::mprint(void)
{
	int i,j;
	for(i = 0;i < row; i++) {
		fprintf(stderr," | ");
		for(j = 0;j < col;j++) {
			fprintf(stderr,"%8.4f ",real(m[col * i + j]));
		}
		fprintf(stderr," |\n");
	}
	fprintf(stderr,"\n");
}