The program output is also avaliable.

//
//	Test file - TFMatrix class and some TMatrix member functions 
//
#include <iostream.h>
#include "tfullmat.h"
#include "tsfulmat.h"

int main() 
{
	REAL coef,ST1[4],ST2[4];
	int i;
	TSFMatrix D(4);		// D is a 4x4 symmetric zero matrix  
	TFMatrix B(4,5),	// B: 4x5 zero matrix
		 A(4,4,10.),	// A,X,Y: 4x4 matrices so that
		 X(4,4,1.),	// A(i,j) = 10, X(i,j) = 1 and 
		 Y(4,4,50),	// Y(i,j) = 50 for all i and j. 
	         Z(D),		// Z is initialized as a copy of D
 		 M(2,2,ST1,4),	// M and P have a special memory
		 P(2,2,ST2,4),	// allocation, associated with ST  
	         b(2,1);	// b: bidimensional column
	         		// vector with zero components
	
	// handling with components of a matrix
	
	for(i=0;i<4;i++) 	// Modify the first line of A:
	   A.PutVal(0,i,i+1);	// PutVal performs  A(0,i) = i+1	
	A.Print("A");		// Shows the matrix A
	coef = A.GetVal(2,2);	// GetVal performs coef = A(2,2)
	B(1,3) = coef;		// It's equivalent to PutVal(1,3)
	B.Print("B, with B(1,3) = A(2,2)");
	
	// Operations, in details
	
	
	A.Multiply(X,Z);		// Product of Matrices
	Z.Print("Z = A*X");
	A.Add(X,Z);			// Sum of Matrices
	Z.Print("Z = A + X");
	A.MultAdd(X,Y,Z,2.,2.2);	// Combine Sum and Product
	Z.Print("Z = A*(2.0*x) + 2.2*Y");
	A.ZAXPY(3,X);
	A.Print("A = A + 3*X");
	A += X*Y - 100.;		// Overloaded operators
	A.Print("A = A + X*Y - C, C(i,j) = 100	 for all i,j");
	Y.TimesBetaPlusZ(.02,X);
	Y.Print("Y = 0.02*Y + X");	
	
	// Changing the matrix structure
	 
	A.Resize(3,6);
	A.Print("A redimensioned");
	Z.Transpose(&Y);	// The argument will be modified
	Y.Print("Y = Transpose of Z");	

//	Solving a Linear System	Mx = b, using LU decomposition
	
	M(0,0) = 2. ;		// We will enter M(i,j) values
	M(1,0) = 1. ;
	M(1,1) = 2. ;
	b(0,0) = 1. ;		
	cout << M << b	;	// Operator << can be used
	out << M << b	;
	P = M;
	P.Decompose_LU();	// LU decomposition of M
	P.Substitution(&b);	// Back Substitution
	b.Print("b: solution of Mx = b");
				
	return(0);	
}

Output:

Writing matrix 'A' (4 x 4):
	1  2  3  4  
	10  10  10  10  
	10  10  10  10  
	10  10  10  10  

Writing matrix 'B, with B(1,3) = A(2,2)' (4 x 5):
	0  0  0  0  0  
	0  0  0  10  0  
	0  0  0  0  0  
	0  0  0  0  0  

Writing matrix 'Z = A*X' (4 x 4):
	10  10  10  10  
	40  40  40  40  
	40  40  40  40  
	40  40  40  40  

Writing matrix 'Z = A + X' (4 x 4):
	2  3  4  5  
	11  11  11  11  
	11  11  11  11  
	11  11  11  11  

Writing matrix 'Z = A*(2.0*x) + 2.2*Y' (4 x 4):
	132  133  134  135  
	201  201  201  201  
	201  201  201  201  
	201  201  201  201  

Writing matrix 'A = A + 3*X' (4 x 4):
	4  5  6  7  
	13  13  13  13  
	13  13  13  13  
	13  13  13  13  

Writing matrix 'A = A + X*Y - C, C(i,j) = 100	 for all i,j' (4 x 4):
	104  105  106  107  
	113  113  113  113  
	113  113  113  113  
	113  113  113  113  

Writing matrix 'Y = 0.02*Y + X' (4 x 4):
	2  2  2  2  
	2  2  2  2  
	2  2  2  2  
	2  2  2  2  

Writing matrix 'A redimensioned' (3 x 6):
	104  105  106  107  0  0  
	113  113  113  113  0  0  
	113  113  113  113  0  0  

Writing matrix 'Y = Transpose of Z' (4 x 4):
	132  201  201  201  
	133  201  201  201  
	134  201  201  201  
	135  201  201  201  

Writing matrix 'operator << ' (2 x 2):
	2  0  
	1  2  

Writing matrix 'operator << ' (2 x 1):
	1  
	0  

Writing matrix 'b: solution of Mx = b' (2 x 1):
	0.5  
	-0.25  

Hints:

Use with other classes

The result of an operation involving a fullmat matrix is also a fullmat matrix, doesn't matter the class of the other matrices present at the operation.

PutVal and GetVal use

Row and column index start on zero. To (i,j) position, use PutVal(i-1,j-1,value) and GetVal(i-1,j-1)

When things go wrong

If two or more matrices don't match at some operation, as sum or multiplication, the function which performs that operation will return a error message, and the program will be stopped. The same occurs when we try to solve a linear system whose matrix is singular.

Designing fast programs

Allocate dynamic variables takes a considerable amount of time, so as floating operations. The least you create intermediate objects, the least time you'll lose. It can be done using Multiply instead of the operator "*" and associating a new object to a fixed position of memory (see objects M and P), for example.

About REAL

REAL is a type of data that may be double or float, depending on the choice of precision. This choice counts for all the process of compilation, not only at matrix level.