473,224 Members | 1,646 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,224 software developers and data experts.

Old template class works in VC++ Not in C++ Builder 5

Diehards,

I developed a template matrix class back around 1992 using Borland C++ 4.5
(ancestor of C++ Builder) and haven't touched it until a few days ago. I
pulled it from the freezer and thawed it out. I built a console app using
Microsoft Visual C++ 6 (VC++) and it worked great. Only one line in the
header file had to be commented out.

I built a console app using Borland C++ Builder 5. The linker complained of
references to external functions (my friend functions that were referenced)
that were not in the object file. I then used the conversion tool supplied
with the IDE to convert the VC++ project to a Borland BPR. Then the linker
complained that ODBCCP32.LIB could not be opened. I included the folder for
ODBCCP32.LIB in the linker library path and tried again, unsuccessfully.

I include the source code with output from the VC++ executable. I hope
someone can help me get this fine piece of code :) working again.

Thanks,
Hamilton Woods

---------- Beginning of Xcholsl.cpp ----------
#include <math.h>

#include <iostream.h>

#include <fstream.h>

#include "tmatrix.h"

void main()

{ // xcholsl

int N = 3;

dmat a(N,N), b(N,N),bchol(N,N);

dmat c(N,N);

dmat eigval(N), eigvec(N,N);

cout.setf(ios::scientific, ios::floatfield);

cout.setf(ios::right, ios::adjustfield);

cout.precision(8);

cout.width(17);

a(1,1) = 7.;

a(1,2) = 4.;

a(1,3) = 3.;

a(2,1) = 4.;

a(2,2) = 8.;

a(2,3) = 2.;

a(3,1) = 3.;

a(3,2) = 2.;

a(3,3) = 6.;

b(1,1) = 8.;

b(1,2) = 1.;

b(1,3) = 3.;

b(2,1) = 1.;

b(2,2) = 6.;

b(2,3) = 4.;

b(3,1) = 3.;

b(3,2) = 4.;

b(3,3) = 4.;

cout << endl << "Original 'a' Matrix:" << endl << a << endl;

cout << endl << "Original 'b' Matrix:" << endl << b << endl;

eigGen(a, b, eigval, eigvec);

cout << endl << "Eigenvalues:" << endl << eigval << endl;

cout << endl << "Eigenvectors:" << endl << eigvec << endl;

bchol = chol(b);

cout << endl << "Cholesky factors:" << endl << bchol << endl;

c = a * b;

cout << endl << "c = a * b" << endl << c << endl;

}

---------- Beginning of TMatrix.h ----------
// tmatrix.h
//
// definition of Matrix class
//
#if !defined _TMatrix
#define _TMatrix
#include <math.h>
#include <string.h>
#include <iostream.h>
#include <stdlib.h>
#ifdef _MSC_VER
//**** ghw2006-11-27 Microsoft VC++ now has a bool type
// typedef int bool;
//****
#define true 1
#define false 0
#endif
#define cmat Matrix<char>
#define dmat Matrix<double>
#define fmat Matrix<float>
#define imat Matrix<int>
#define lmat Matrix<long int>
#define smat Matrix<short int>

template <class Type>
class Matrix
{
public:
// Constructors & Destructor
Matrix(int num_rows, int num_cols = 1);
Matrix(const Matrix<Type&a);
~Matrix();
// Accessors
void SetTab(const char *tabstr);
char *GetErr();
void SetErr(const char *error);
friend int nrows(const Matrix<Type&source);
friend int ncols(const Matrix<Type&source);
Type &operator() (int row, int col = 1) const;
// Type &operator() (int row) const;
// Operator Overloads
friend ostream &operator<<(ostream &target, const Matrix<Type&source);
friend istream &operator>>(istream &source, Matrix<Type&target);
Matrix<Type&operator=(const Matrix<Type&a); // Assignment
friend Matrix<Type // Matrix Addition
operator+(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&source);
friend Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*(double left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, double right);
friend Matrix<Type // Matrix Division
operator/(const Matrix<Type&left, double right);
// Member Functions
void cut(Matrix<Type&target, int r1, int c1 = 1) const;
void paste(Matrix<Type&target, int r1, int c1 = 1) const;
protected:
int NumRows;
int NumCols;
Type **Array;
char Tab[5]; // delimiter used for << operator
char Err[21]; // string to hold error message
};

// Functions
template <class Type // Build identity matrix
void eye(Matrix<Type&a);
template <class Type // build zero matrix
void zeros(Matrix<Type&a);
template <class Type // Matrix Transpose
Matrix<Typetranspose(const Matrix<Type&a);
template <class Type // element by element square root (only if all
elements are positive)
Matrix<Typesqrt(const Matrix<Type&source);
template <class Type // normalize matrix
Matrix<Typecolnormalize(const Matrix<Type&source);
template <class Type // normalize matrix
Matrix<Typecolunitize(const Matrix<Type&source, int urow = 1);
template <class Type // solve aX = b using lu decomposition w/
backsubstitution
Matrix<Typesolve(const Matrix<Type&a, const Matrix<Type&b);
template <class Type // Matrix Inverse. aX = I
Matrix<Typeinv(const Matrix<Type&a);
template <class Type // lower upper backsubstitution
Matrix<Typelubksb (const Matrix<Type&a, Matrix<int&indx,
const Matrix<Type&b);
template <class Type // lower upper decomposition
Matrix<Typeludcmp(const Matrix<Type&asource, bool &ierr,
Matrix<int&indx, int &d);
template <class Type // indexsort. Works on column vectors only, for
now.
Matrix<intindexsort(const Matrix<Type&source);
template <class Type // heapsort. Works on column vectors only, for now.
Matrix<Typeheapsort(const Matrix<Type&source, Matrix<int&index);
template <class Type // heapsort. Does not return sort index.
Matrix<Typeheapsort(const Matrix<Type&source);
template <class Type // Row Swap
Matrix<Typerowswap(const Matrix<Type&source, const Matrix<int>
&swapvec);
template <class Type // Column Swap
Matrix<Typecolswap(const Matrix<Type&source, const Matrix<int>
&swapvec);
template <class Type // Row/Col Swap. Calls rowswap and then colswap
Matrix<Typerowcolswap(const Matrix<Type&source, const Matrix<int>
&swapvec);
template <class Type // Householder eigensolver for generalized eigen
problem
void eigGen(const Matrix<Type&A, const Matrix<Type&B,
Matrix<Type&eigval, Matrix<Type&eigvec);
template <class Type // Householder eigensolver for generalized eigen
problem
void eigGen(const Matrix<Type&A, const Matrix<Type&B,
Matrix<Type&eigval);
template <class Type // Householder eigensolver
void eig(const Matrix<Type&A, Matrix<Type&eigval, Matrix<Type>
&eigvec);
template <class Type // Householder eigensolver, without eigenvectors
void eig(const Matrix<Type&A, Matrix<Type&eigval);
template <class Type // mass normalize eigenvectors
Matrix<Typeeig_mass_normalize(const Matrix<Type&eigen,
const Matrix<Type&mass);
template <class Type // Householder reduction of a real, symmetric matrix
void tred2(const Matrix<Type&a, Matrix<Type&d, Matrix<Type&e);
template <class Type // eigenvalues, vectors of a real, symmetric,
tridiagonal matrix
void tqli(Matrix<Type&eigval, Matrix<Type&e, Matrix<Type&eigvec);
template <class Type // Cholesky decomposition
Matrix<Typechol(const Matrix<Type&source);
template <class Type // Cholesky decomposition, with p vector
Matrix<Typechol(const Matrix<Type&source, Matrix<Type&p);
template <class Type // Craig-Bampton reduction
void CraigBampton(const Matrix<Type&k, const Matrix<Type&m,
const Matrix<int&swapvec, int NgeneralizedDOF,
Matrix<Type&phiCB, Matrix<Type&kCB, Matrix<Type>
&mCB);

//Constructors and Destructor:

// Constructor for doubly subscripted array.
template <class Type>
Matrix<Type>::Matrix(int num_rows, int num_cols)
{
int i;
strcpy(Err, "");
strcpy(Tab, "\t");
NumRows = num_rows;
NumCols = num_cols;
Array = new Type *[NumRows+1];
for (i=0; i<=NumRows; ++i)
Array[i] = new Type [NumCols+1];
}

// Copy Constructor. Needed because constructor performs memory allocation
template <class Type>
Matrix<Type>::Matrix(const Matrix<Type&source)
{
int i, j;
strcpy(Err, source.Err);
strcpy(Tab, source.Tab);
NumRows = nrows(source);
NumCols = ncols(source);
Array = new Type *[NumRows+1];
for (i=0; i<=NumRows; ++i)
Array[i] = new Type [NumCols+1];
for (i=0; i<=NumRows; ++i)
for (j=0; j<=NumCols; ++j)
(*this)(i,j) = source(i,j);
}

// Destructor
template <class Type>
Matrix<Type>::~Matrix()
{
int i;
for (i=0; i<=NumRows; ++i)
{
delete [] Array[i];
Array[i] = 0;
}
delete [] Array;
Array = 0;
}

//Accessors:

// Set separator character for use with <<
template <class Type>
void Matrix<Type>::SetTab(const char *tabstr)
{
strcpy(Tab, tabstr);
}

// Get Error Code of Matrix object
template <class Type>
char *Matrix<Type>::GetErr()
{
return Err;
}

// Set Error Code of Matrix object

template <class Type>
void Matrix<Type>::SetErr(const char *error)
{
strcpy(Err, error);
}

// Overload () operator to access elements by reference
template <class Type>
Type &Matrix<Type>::operator() (int row, int col) const
{
int i, j;
if (row NumRows || col NumCols)
{
cerr << "(): out of range" << endl;
i = 0;
j = 0;
}
else
{
i = row;
j = col;
}
return Array[i][j];
}

// Access number of rows
template <class Type>
int nrows(const Matrix<Type&source)
{
return source.NumRows;
}

// Access number of columns
template <class Type>
int ncols(const Matrix<Type&source)
{
return source.NumCols;
}

//Operator Overloads:

// overload operator << for use with Matrix class
template <class Type>
ostream &operator<<(ostream &target, const Matrix<Type&source)
{
int i, j;
for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
target << source.Tab << source(i,j);
}
target << endl;
}
return target;
}

// overload operator >for use with Matrix class
template <class Type>
istream &operator>>(istream &source, Matrix<Type&target)
{
int i, j;
for (i=1; i<=nrows(target); ++i)
{
for (j=1; j<=ncols(target); ++j)
{
source >target(i,j);
}
}
return source;
}

// Overload assignment operator =
template <class Type>
Matrix<Type&Matrix<Type>::operator=(const Matrix<Type&source)
{
int i, j;
if (NumRows != nrows(source) && NumCols != ncols(source))
strcpy(Err, "=: size mismatch");
else
{
NumRows = nrows(source);
NumCols = ncols(source);
strcpy(Err, source.Err);
strcpy(Tab, source.Tab);
for (i=0; i<=NumRows; ++i)
for (j=0; j<=NumCols; ++j)
(*this)(i,j) = source(i,j);
}
return *this;
}

// Overload addition operator + for Matrix class
template <class Type>
Matrix<Typeoperator+(const Matrix<Type&left, const Matrix<Type&right)
{
int i, j;
Matrix<Typesum(nrows(left), ncols(left));

if (nrows(left) != nrows(right) || ncols(left) != ncols(right))
{
sum.SetErr("+: size mismatch");
}
else
{
for (i=0; i<=nrows(left); ++i)
for (j=0; j<=ncols(left); ++j)
sum(i,j) = left(i,j) + right(i,j);
}
return sum;
}

// Overload addition operator - for Matrix class
template <class Type>
Matrix<Typeoperator-(const Matrix<Type&left, const Matrix<Type&right)
{
int i, j;
Matrix<Typediff(nrows(left), ncols(left));

if (nrows(left) != nrows(right) || ncols(left) != ncols(right))
{
diff.SetErr("-: size mismatch");
}
else
{
for (i=1; i<=nrows(left); ++i)
for (j=1; j<=ncols(left); ++j)
diff(i,j) = left(i,j) - right(i,j);
}
return diff;
}

// Overload unary operator - for Matrix class
template <class Type>
Matrix<Typeoperator-(const Matrix<Type&source)
{
int i, j;
Matrix<Typetarget(nrows(source), ncols(source));

for (i=1; i<=nrows(source); ++i)
for (j=1; j<=ncols(source); ++j)
target(i,j) = -source(i,j);
return target;
}

// Overload multiply operator * for Matrix class
// implemented as a friend function in order to overload for pre- and post-
// multiplies: A * B; c * A; A * c; A * V; V * A.
template <class Type>
Matrix<Typeoperator*(const Matrix<Type&left, const Matrix<Type&right)
{
int i, j, k;
Matrix<Typeproduct(nrows(left),ncols(right));

if (ncols(left) != nrows(right))
{
product.SetErr("*: size mismatch");
}
else
{
for (i=0; i<=nrows(product); ++i)
for (k=0; k<=ncols(product); ++k)
product(i,k) = 0.;
for (i=1; i<=nrows(left); ++i)
for (j=1; j<=ncols(left); ++j)
for (k=1; k<=ncols(right); ++k)
product(i,k) += left(i,j) * right(j,k);
}
return product;
}

// c * A
template <class Type>
Matrix<Typeoperator*(double left, const Matrix<Type&right)
{
int i, j;
Matrix<Typeproduct(nrows(right),ncols(right));

for (i=1; i<=nrows(product); ++i)
for (j=1; j<=ncols(product); ++j)
product(i,j) = left * right(i,j);
return product;
}

// A * c
template <class Type>
Matrix<Typeoperator*(const Matrix<Type&left, double right)
{
Matrix<Typeproduct(nrows(left),ncols(left));

product = right * left; // use c * A
return product;
}

// A / c
template <class Type>
Matrix<Typeoperator/(const Matrix<Type&left, double right)
{
Matrix<Typetarget(nrows(left), ncols(left));
int i, j;

if (right == 0.)
{
target.SetErr("/: 0 divisor");
return target;
}
for (i=1; i<=nrows(left); ++i)
{
for (j=1; j<=ncols(left); ++j)
{
target(i,j) = left(i,j) / right;
}
}
return target;
}

// Member Functions

// cut a piece of this matrix starting at (r1, c1) that fills the target
Matrix
template <class Type>
void Matrix<Type>::cut(Matrix<Type&target, int r1, int c1) const
{
int i, j;

for (i=1; i<=nrows(target); ++i)
{
for (j=1; j<=ncols(target); ++j)
{
target(i,j) = (*this)(r1+i-1,c1+j-1);
}
}
}

// paste this matrix into the target starting at (r1, c1). Assumes target
is
// at least as large as this matrix
template <class Type>
void Matrix<Type>::paste(Matrix<Type&target, int r1, int c1) const
{
int i, j;

if (nrows(target) < (NumRows - (r1 - 1)) || ncols(target) < (NumCols -
(c1 - 1)))
{
target.SetErr("paste: Insufficient space in target");
}
else
{
for (i=1; i<=NumRows; ++i)
{
for (j=1; j<=NumCols; ++j)
{
target(r1+i-1,c1+j-1) = (*this)(i,j);
}
}
}
}

//Friend Functions:

// Identity Matrix
template <class Type>
void eye(Matrix<Type&a)
{
int i, j;
for (i=1; i<=nrows(a); ++i)
{
for (j=1; j<=nrows(a); ++j)
a(i,j) = Type(0);
a(i,i) = Type(1);
}
}

// Zero Matrix
template <class Type>
void zeros(Matrix<Type&a)
{
int i, j;
for (i=1; i<=nrows(a); ++i)
{
for (j=1; j<=ncols(a); ++j)
a(i,j) = Type(0);
}
}

// friend function to perform matrix transpose
template <class Type>
Matrix<Typetranspose(const Matrix<Type&source)
{
int i, j;
Matrix<Typetarget(ncols(source), nrows(source));
for (i=1; i<=nrows(source); ++i)
for (j=1; j<=ncols(source); ++j)
target(j,i) = source(i,j);
return target;
}

// Square root of matrix elements
template <class Type>
Matrix<Typesqrt(const Matrix<Type&source)
{
Matrix<Typetarget(nrows(source), ncols(source));
int i, j;

for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
if (source(i,j) < -0.01)
{
target.SetErr("sqrt: Negative value");
target(i,j) = -1.;
}
else if (source(i,j) >= -.01 && source(i,j) <= 1.e-20)
{
target(i,j) = 0.;
}
else
{
target(i,j) = sqrt(source(i,j));
}
}
}
return target;
}

// normalize matrix columns so that root sum square is unity
template <class Type>
Matrix<Typecolnormalize(const Matrix<Type&source)
{
Matrix<Typetarget(nrows(source), ncols(source));
int i, j;
double sumsquare;

for (i=1; i<=nrows(source); ++i)
{
sumsquare = 0;
for (j=1; j<=ncols(source); ++j)
{
sumsquare += (source(i,j) * source(i,j));
}
sumsquare = sqrt(sumsquare);
for (j=1; j<=ncols(source); ++j)
{
target(i,j) = source(i,j) / sumsquare;
}
}
return target;
}

// scale matrix such that row urow becomes unity
template <class Type>
Matrix<Typecolunitize(const Matrix<Type&source, int urow)
{
Matrix<Typetarget(nrows(source), ncols(source));
int i, j;

for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
if (source(urow,j) == 0)
{
}
else
{
target(i,j) = source(i,j) / source(urow,j);
}
}
}
return target;
}

// solve aX = b. Calls ludcmp to find LU(a) and then calls lubksb for each
// column in b.
template <class Type>
Matrix<Typesolve(const Matrix<Type&a, const Matrix<Type&b)
{
Matrix<Typetarget(nrows(b), ncols(b));
Matrix<Typelu(nrows(a), ncols(a));
Matrix<intindex(nrows(b));
Matrix<Typey(nrows(b));
Matrix<TypeX(nrows(b));
int i, j, d;
bool ierr;

if (nrows(a) != nrows(b))
{
target.SetErr("solve: size mismatch");
}
else
{
lu = ludcmp(a, ierr, index, d);
for (j=1; j<=ncols(b); ++j)
{
for (i=1; i<=nrows(b); ++i)
y(i) = b(i,j);
X = lubksb(lu, index, y);
for (i=1; i<=nrows(b); ++i)
target(i,j) = X(i);
}
}
return target;
}

// invert a matrix using lower upper decompostion with backsubstitution
template <class Type>
Matrix<Typeinv(const Matrix<Type&source)
{
Matrix<Typetarget(nrows(source), nrows(source));
Matrix<Typeident(nrows(source), nrows(source));

if (nrows(source) != ncols(source))
{
target.SetErr("inv: size mismatch");
}
else
{
eye(ident);
target = solve(source, ident);
}
return target;
}

// solves the set of linear equations aX=b, where a is LU(a).
template <class Type>
Matrix<Typelubksb(const Matrix<Type&a, Matrix<int&indx,
const Matrix<Type&b)
{
int i, ii, j, ll;
double sum;
Matrix<Typebtemp(nrows(b));

btemp = b;
ii = 0;
for (i=1; i<=nrows(a); ++i)
{
ll = indx(i);
sum = btemp(ll);
btemp(ll) = btemp(i);
if (ii != 0)
for (j=ii; j<i; ++j)
sum -= a(i,j)*btemp(j);
else if (sum != 0.)
ii = i;
btemp(i) = sum;
}
for (i=nrows(a); i>=1; --i)
{
sum = btemp(i);
for (j=i+1; j<=nrows(a); ++j)
sum -= a(i,j)*btemp(j);
btemp(i) = sum / a(i,i);
}
return btemp;
}

// lower upper decomposition of square matrix.
template <class Type>
Matrix<Typeludcmp(const Matrix<Type&asource, bool &ierr,
Matrix<int&indx, int &d)
{
Matrix<Typea(nrows(asource), nrows(asource));
int imax, i, j, k;
double tiny, aamax, dum, sum;
Matrix<Typevv(nrows(asource));
tiny = 1.e-20;

a = asource;
ierr = false;
d=1;
for (i=1; i<=nrows(a); ++i)
{
aamax = 0.;
for (j=1; j<= nrows(a); ++j)
if (fabsl(a(i,j)) aamax) aamax = fabsl(a(i,j));
if (aamax == 0.)
{
ierr = true;
return a;
}
vv(i) = 1./aamax;
}
for (j=1; j<=nrows(a); ++j)
{
for (i=1; i<j; ++i)
{
sum = a(i,j);
for (k=1; k<i; ++k)
sum -= a(i,k)*a(k,j);
a(i,j) = sum;
}
aamax = 0.;
imax = j;
for (i=j; i<=nrows(a); ++i)
{
sum = a(i,j);
for (k=1; k<j; ++k)
sum -= a(i,k) * a(k,j);
a(i,j) = sum;
dum = vv(i) * fabsl(sum);
if (dum >= aamax)
{
imax = i;
aamax = dum;
}
}
if (j != imax)
{
for (k=1; k<=nrows(a); ++k)
{
dum = a(imax,k);
a(imax,k) = a(j,k);
a(j,k) = dum;
}
d = -d;
vv(imax) = vv(j);
}
indx(j) = imax;
if (a(j,j) == 0.) a(j,j) = tiny;
if (j != nrows(a))
{
dum = 1./a(j,j);
for (i=j+1; i<=nrows(a); ++i)
a(i,j) = a(i,j) * dum;
}
}
return a;
}

// indexed version of heapsort. Returns index vector that contains indices
of
// the source matrix such that elements source(index(i)) are sorted in
ascending
// order
template <class Type>
Matrix<intindexsort(const Matrix<Type&source)
/*c
c Heapsort algorithm as presented in "Numerical Recipes" and described in
c "Sorting and Searching", vol. 3 of "The Art of Computer Programming" by
c Donald Knuth.
c
c Works only for column vectors, so far.
c*/
{
int i, j, L, indext, ir;
double q;
int n;
n = nrows(source);
Matrix<intindex(n);

for (j=1; j<=n; ++j)
{
index(j) = j;
}
if (n == 1)
{
return index;
}
L = n/2 + 1;
ir = n;
while (true)
{
if (L 1)
{
L--;
indext = index(L);
q = source(indext);
}
else
{
indext = index(ir);
q = source(indext);
index(ir) = index(1);
ir--;
if (ir == 1)
{
index(1) = indext;
return index;
}
}
i = L;
j = L + L;
while (j <= ir)
{
if (j < ir)
{
if (source(index(j)) < source(index(j+1))) j++;
}
if (q < source(index(j)))
{
index(i) = index(j);
i = j;
j += j;
}
else
{
j = ir + 1;
}
}
index(i) = indext;
}
}

// heapsort with index returns sorted matrix and also index by which source
is sorted
template <class Type>
Matrix<Typeheapsort(const Matrix<Type&source, Matrix<int&index)
{
index = indexsort(source);
return rowswap(source, index);
}

// heapsort w/o index returns only the sorted matrix
template <class Type>
Matrix<Typeheapsort(const Matrix<Type&source)
{
Matrix<intno_index(nrows(source));

return heapsort(source, no_index);
}

// Row swap
template <class Type>
Matrix<Typerowswap(const Matrix<Type&source, const Matrix<int&swapvec)
{
int i, j, k;
bool hit;
Matrix<intswapndx(nrows(source));
Matrix<Typetarget(nrows(source), ncols(source));

for (i=1; i<=nrows(swapvec); ++i)
{
for (j=1; j<=i-1; ++j)
{
if (swapvec(j) == swapvec(i))
{
target.SetErr("rowswap: double entry in swap vector");
return target;
}
}
swapndx(i) = swapvec(i);
}
k = 0;
for (i=1; i<=nrows(source); ++i)
{
hit = false;
for (j=1; j<=nrows(swapvec); ++j)
{
if (i==swapvec(j))
{
hit = true;
break;
}
}
if (!hit)
{
++k;
swapndx(nrows(swapvec)+ k) = i;
}
}
for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
target(i,j) = source(swapndx(i),j);
}
}
return target;
}

// Column Swap
template <class Type>
Matrix<Typecolswap(const Matrix<Type&source, const Matrix<int&swapvec)
{
int i, j, k;
bool hit;
Matrix<intswapndx(nrows(source));
Matrix<Typetarget(nrows(source), ncols(source));

for (i=1; i<=nrows(swapvec); ++i)
{
for (j=1; j<=i-1; ++j)
{
if (swapvec(j) == swapvec(i))
{
target.SetErr("colswap: double entry in swap vector");
return target;
}
}
swapndx(i) = swapvec(i);
}
k = 0;
for (j=1; j<=ncols(source); ++j)
{
hit = false;
for (i=1; i<=nrows(swapvec); ++i)
{
if (j==swapvec(i))
{
hit = true;
break;
}
}
if (!hit)
{
++k;
swapndx(nrows(swapvec)+ k) = j;
}
}
for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
target(i,j) = source(i,swapndx(j));
}
}
return target;
}

// Row/Column Swap
template <class Type>
Matrix<Typerowcolswap(const Matrix<Type&source, const Matrix<int>
&swapvec)
{
Matrix<Typetarget(nrows(source), ncols(source));

target = rowswap(source, swapvec);
target = colswap(target, swapvec);
return target;
}

// Householder Eigensolver. Generalized eigen problem (Ax = lBx) with
eigenvectors.
template <class Type>
void eigGen(const Matrix<Type&A, const Matrix<Type&B, Matrix<Type>
&eigval,
Matrix<Type&eigvec)
{
Matrix<Typee(nrows(A));
Matrix<TypeLinvT(nrows(A), ncols(A));
Matrix<TypeH(nrows(A), ncols(A));

H = inv(chol(B));
LinvT = transpose(H);
H = H * A * LinvT;
eig(H, eigval, eigvec);
eigvec = LinvT * eigvec;
eigvec = eig_mass_normalize(eigvec, B);
}

// Householder Eigensolver. Generalized eigen problem (Ax = lBx) w/o
eigenvectors.
template <class Type>
void eigGen(const Matrix<Type&A, const Matrix<Type&B, Matrix<Type>
&eigval)
{
Matrix<Typeno_eigvec(nrows(A), ncols(A));

eigGen(A, B, eigval, no_eigvec);
}

// Householder Eigensolver. Eigen problem (Ax = lx) with eigenvectors.
// from Numerical Recipes
template <class Type>
void eig(const Matrix<Type&A, Matrix<Type&eigval, Matrix<Type&eigvec)
{
Matrix<Typee(nrows(A));
Matrix<intindex(nrows(A));

eigvec = A;
tred2(eigvec,eigval,e);
tqli(eigval,e,eigvec);
eigval = heapsort(eigval, index);
eigvec = colswap(eigvec, index);
}

// Householder Eigensolver. Eigen problem (Ax = lx) w/o eigenvectors.
// from Numerical Recipes
template <class Type>
void eig(const Matrix<Type&A, Matrix<Type&eigval)
{
Matrix<Typeno_eigvec(nrows(A), ncols(A));

eig(A, eigval, no_eigvec);
}

// Mass normalize eigenvectors
template <class Type>
Matrix<Typeeig_mass_normalize(const Matrix<Type&eigen,
const Matrix<Type&mass)
/*
This subroutine uses the methodology found in Paz, Structural Dynamics,
to scale the general eigenvectors (eigen) such that the resulting
eigenvectors (eigmn) are mass normalized.
*/
{
int n, i, j, k;
n = nrows(eigen);
Matrix<Typeeigmn(n,n);
double denom, temp;

for (j=1; j<=n; ++j)
{
denom = 0.;
for (k=1; k<=n; ++k)
{
temp = 0.;
for (i=1; i<=n; ++i)
{
temp = temp + eigen(i,j) * mass(i,k);
}
denom = denom + temp * eigen(k,j);
}
for (i=1; i<=n; ++i)
{
eigmn(i,j) = eigen(i,j) / denom;
}
}
return eigmn;
}

// Householder reduction
template <class Type>
void tred2(const Matrix<Type&a, Matrix<Type&d, Matrix<Type&e)
{
int i,j,k,l, n;
double f,g,h,hh,scale;

n = nrows(e);
for (i=n; i>=2;--i)
{
l=i-1;
h=0.;
scale=0.;
if(l>1)
{
for (k=1; k<=l; ++k)
{
scale=scale+fabs(a(i,k));
}
if(scale==0.)
{
e(i)=a(i,l);
}
else
{
for (k=1; k<=l; ++k)
{
a(i,k)=a(i,k)/scale;
h=h+a(i,k)*a(i,k);
}
f=a(i,l);
g=-(fabs(f)/f)*sqrt(h);
e(i)=scale*g;
h=h-f*g;
a(i,l)=f-g ;
f=0.;
for (j=1; j<=l; ++j)
{
//C Omit following line if finding only eigenvalues
a(j,i)=a(i,j)/h;
g=0.;
for (k=1; k<=j; ++k)
{
g=g+a(j,k)*a(i,k);
}
for (k=j+1; k<=l; ++k)
{
g=g+a(k,j)*a(i,k);
}
e(j)=g/h;
f=f+e(j)*a(i,j);
}
hh=f/(h+h);
for (j=1; j<=l; ++j)
{
f=a(i,j);
g=e(j)-hh*f;
e(j)=g;
for (k=1; k<=j; ++k)
{
a(j,k)=a(j,k)-f*e(k)-g*a(i,k);
}
}
}
}
else
{
e(i)=a(i,l);
}
d(i)=h;
}
//C Omit following line if finding only eigenvalues.
d(1)=0.;
e(1)=0.;
for (i=1; i<=n; ++i)
{
//C Delete lines from here ...
l=i-1;
if (d(i)!=0.)
{
for (j=1; j<=l; ++j)
{
g=0.;
for (k=1; k<=l; ++k)
{
g=g+a(i,k)*a(k,j);
}
for (k=1; k<=l; ++k)
{
a(k,j)=a(k,j)-g*a(k,i);
}
}
}
//C ... to here when finding only eigenvalues.
d(i)=a(i,i);
//C Also delete lines from here ...
a(i,i)=1.;
for (j=1; j<=l; ++j)
{
a(i,j)=0.;
a(j,i)=0.;
}
//C ... to here when finding only eigenvalues.
}
}

// eigenvalues, vectors of a real, symmetric, tridiagonal matrix
template <class Type>
void tqli(Matrix<Type&eigval, Matrix<Type&e, Matrix<Type&eigvec)
{
//CU USES pythag
int i,iter,k,l,m, n;
double b,c,dd,f,g,p,r,s;

n = nrows(e);
for (i=2; i<=n; ++i)
{
e(i-1)=e(i);
}
e(n)=0.;
for (l=1; l<=n; ++l)
{
iter=0;
label1:
for (m=l; m<=n-1; ++m)
{
dd=fabs(eigval(m))+fabs(eigval(m+1));
if (fabs(e(m))+dd==dd) goto label2;
}
m=n;
label2:
if(m!=l)
{
if(iter==30)
{
cerr << endl << "too many iterations in tqli" << endl;
exit(1);
}
iter=iter+1;
g=(eigval(l+1)-eigval(l))/(2.*e(l));
r=sqrt(g*g + 1.);
g=eigval(m)-eigval(l)+e(l)/(g+(fabs(g)/g)*r);
s=1.;
c=1.;
p=0.;
for (i=m-1; i>=l; --i)
{
f=s*e(i);
b=c*e(i);
r=sqrt(f*f + g*g);
e(i+1)=r;
if (r==0.)
{
eigval(i+1)=eigval(i+1)-p;
e(m)=0.;
goto label1;
}
s=f/r;
c=g/r;
g=eigval(i+1)-p;
r=(eigval(i)-g)*s+2.*c*b;
p=s*r;
eigval(i+1)=g+p;
g=c*r-b;
//C Omit lines from here ...
if (nrows(eigvec) 0)
{
for (k=1; k<=n; ++k)
{
f=eigvec(k,i+1);
eigvec(k,i+1)=s*eigvec(k,i)+c*f;
eigvec(k,i)=c*eigvec(k,i)-s*f;
}
}
//C ... to here when finding only eigenvalues.
}
eigval(l)=eigval(l)-p;
e(l)=g;
e(m)=0.;
goto label1;
}
}
}

// Cholesky decomposition.
template <class Type>
Matrix<Typechol(const Matrix<Type&source)
{
int i,j, n;
n = nrows(source);
Matrix<Typecholesky(n, n);
Matrix<Typep(n);

cholesky = chol(source, p);
for (i=1; i<=n; ++i)
{
for (j=1; j<=n; ++j)
{
if (i>j)
{
// cholesky(i,j)=source(i,j);
}
else if (i==j)
{
cholesky(i,j)=p(i);
}
else
{
cholesky(i,j)=0.;
}
}
}
return cholesky;
}

// Cholesky decomposition, returning p.
template <class Type>
Matrix<Typechol(const Matrix<Type&source, Matrix<Type&p)
{
int i,j,k, n;
double sum;
Matrix<Typetarget(nrows(source), ncols(source));

n = nrows(source);
target = source;
for (i=1; i<=n; ++i)
{
for (j=i; j<=n; ++j)
{
sum = source(i,j);
for (k=i-1; k>=1; --k)
{
sum = sum - target(i,k) * target(j,k);
}
if (i == j)
{
if(sum <= 0.)
{
cout << endl << "choldc failed: i = " << i;
cout << ", source(i,i) = " << source(i,i) << endl;
exit(1);
}
p(i) = sqrt(sum);
}
else
{
target(j,i) = sum / p(i);
}
}
}
return target;
}

// Craig-Bampton reduction
template <class Type>
void CraigBampton(const Matrix<Type&k, const Matrix<Type&m,
const Matrix<int&swapvec, int NgeneralizedDOF,
Matrix<Type&phiCB, Matrix<Type&kCB, Matrix<Type&mCB)
{
int N = nrows(k);
int Nswap = nrows(swapvec);
int Nreduced = Nswap + NgeneralizedDOF;
dmat kswap(N, N);
dmat mswap(N, N);
dmat kpg(N-Nswap, Nswap);
dmat kgg(N-Nswap, N-Nswap);
dmat mgg(N-Nswap, N-Nswap);
dmat I(Nswap, Nswap);
dmat Zero(Nswap, NgeneralizedDOF);
dmat phic(N-Nswap, Nswap);
dmat phin(N-Nswap, N-Nswap);
dmat phin_red(N-Nswap, NgeneralizedDOF);
dmat phival(N-Nswap);

kswap = rowcolswap(k, swapvec);
mswap = rowcolswap(m, swapvec);
kswap.cut(kgg, Nswap+1, Nswap+1);
kswap.cut(kpg, Nswap+1, 1);
mswap.cut(mgg, Nswap+1, Nswap+1);
eye(I);
zeros(Zero);
phic = -inv(kgg) * kpg;
eigGen(kgg, mgg, phival, phin);
phin.cut(phin_red, 1, 1);
phin_red = colunitize(phin_red, Nreduced);
I.paste(phiCB, 1, 1);
Zero.paste(phiCB, 1, Nswap+1);
phic.paste(phiCB, Nswap+1, 1);
phin_red.paste(phiCB, Nswap+1, Nswap+1);
kCB = transpose(phiCB) * kswap * phiCB;
mCB = transpose(phiCB) * mswap * phiCB;
}
#endif

---------- Beginning of Tester.out ----------

Original 'a' Matrix:
7.00000000e+000 4.00000000e+000 3.00000000e+000
4.00000000e+000 8.00000000e+000 2.00000000e+000
3.00000000e+000 2.00000000e+000 6.00000000e+000
Original 'b' Matrix:
8.00000000e+000 1.00000000e+000 3.00000000e+000
1.00000000e+000 6.00000000e+000 4.00000000e+000
3.00000000e+000 4.00000000e+000 4.00000000e+000
Eigenvalues:
4.98743268e-001
1.11676647e+000
1.12511569e+001
Eigenvectors:
3.17174802e-001 -1.47618076e-001 -3.79836432e-001
-2.21648368e-001 -1.28188516e-001 -8.37320949e-001
-1.18811637e-001 -2.40036434e-001 1.22267452e+000
Cholesky factors:
2.82842712e+000 0.00000000e+000 0.00000000e+000
3.53553391e-001 2.42383993e+000 0.00000000e+000
1.06066017e+000 1.49556081e+000 7.98935462e-001
c = a * b
6.90000000e+001 4.30000000e+001 4.90000000e+001
4.60000000e+001 6.00000000e+001 5.20000000e+001
4.40000000e+001 3.90000000e+001 4.10000000e+001

Nov 30 '06 #1
3 3718
tried to fix some not yet finished, but no more time - still linker error

#include <fstream>
#include <stdexcept>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cctype>
#include <stdio.h>
#include <iostream>
#include <string>
#include <fstream// Dateioperationen
#include <math.h>
#include <iostream>
#include "TMatrix.h"

int main()
{ // xcholsl
int N = 3;
dmat a(N,N), b(N,N),bchol(N,N);
dmat c(N,N);
dmat eigval(N), eigvec(N,N);
cout.setf(ios::scientific, ios::floatfield);
cout.setf(ios::right, ios::adjustfield);
cout.precision(8);
cout.width(17);
a(1,1) = 7.;
a(1,2) = 4.;
a(1,3) = 3.;
a(2,1) = 4.;
a(2,2) = 8.;
a(2,3) = 2.;
a(3,1) = 3.;
a(3,2) = 2.;
a(3,3) = 6.;
b(1,1) = 8.;
b(1,2) = 1.;
b(1,3) = 3.;
b(2,1) = 1.;
b(2,2) = 6.;
b(2,3) = 4.;
b(3,1) = 3.;
b(3,2) = 4.;
b(3,3) = 4.;
cout << endl << "Original 'a' Matrix:" << endl << a << endl;
cout << endl << "Original 'b' Matrix:" << endl << b << endl;
eigGen(a, b, eigval, eigvec);
cout << endl << "Eigenvalues:" << endl << eigval << endl;
cout << endl << "Eigenvectors:" << endl << eigvec << endl;
bchol = chol(b);
cout << endl << "Cholesky factors:" << endl << bchol << endl;
c = a * b;
cout << endl << "c = a * b" << endl << c << endl;
}

///////////////////////////////////

#if !defined _TMatrix
#define _TMatrix
#include <math.h>
#include <string.h>
#include <iostream>
//#include <gstream.h>
#include <stdlib.h>
#ifdef _MSC_VER
#define true 1
#define false 0
#endif
#define cmat Matrix<char>
#define dmat Matrix<double>
#define fmat Matrix<float>
#define imat Matrix<int>
#define lmat Matrix<long int>
#define smat Matrix<short int>
using namespace std;
template <class Type>
class Matrix
{
public:
// Constructors & Destructor
Matrix(int num_rows, int num_cols = 1);
Matrix(const Matrix<Type&a);
~Matrix();
// Accessors
void SetTab(const char *tabstr);
char *GetErr();
void SetErr(const char *error);
friend int nrows(const Matrix<Type&source);
friend int ncols(const Matrix<Type&source);
Type &operator() (int row, int col = 1) const;
// Type &operator() (int row) const;
// Operator Overloads
friend ostream& operator << (ostream &target, const Matrix<Type&source);
friend istream& operator >(istream &source, Matrix<Type&target);
Matrix<Type&operator=(const Matrix<Type&a); // Assignment
friend Matrix<Type // Matrix Addition
operator+(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&source);
friend Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*(double left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, double right);
friend Matrix<Type // Matrix Division
operator/(const Matrix<Type&left, double right);
// Member Functions
void cut(Matrix<Type&target, int r1, int c1 = 1) const;
void paste(Matrix<Type&target, int r1, int c1 = 1) const;
protected:
int NumRows;
int NumCols;
Type **Array;
char Tab[5]; // delimiter used for << operator
char Err[21]; // string to hold error message
};

// Functions
template <class Type // Build identity matrix
void eye(Matrix<Type&a);
template <class Type // build zero matrix
void zeros(Matrix<Type&a);
template <class Type // Matrix Transpose
Matrix<Typetranspose(const Matrix<Type&a);
template <class Type // element by element square root (only if all
elements are positive)
Matrix<Typesqrt(const Matrix<Type&source);
template <class Type // normalize matrix
Matrix<Typecolnormalize(const Matrix<Type&source);
template <class Type // normalize matrix
Matrix<Typecolunitize(const Matrix<Type&source, int urow = 1);
template <class Type // solve aX = b using lu decomposition w/
backsubstitution
Matrix<Typesolve(const Matrix<Type&a, const Matrix<Type&b);
template <class Type // Matrix Inverse. aX = I
Matrix<Typeinv(const Matrix<Type&a);
template <class Type // lower upper backsubstitution
Matrix<Typelubksb (const Matrix<Type&a, Matrix<int&indx,
const Matrix<Type&b);
template <class Type // lower upper decomposition
Matrix<Typeludcmp(const Matrix<Type&asource, bool &ierr,
Matrix<int&indx, int &d);
template <class Type // indexsort. Works on column vectors only, for
now.
Matrix<intindexsort(const Matrix<Type&source);
template <class Type // heapsort. Works on column vectors only, for now.
Matrix<Typeheapsort(const Matrix<Type&source, Matrix<int&index);
template <class Type // heapsort. Does not return sort index.
Matrix<Typeheapsort(const Matrix<Type&source);
template <class Type // Row Swap
Matrix<Typerowswap(const Matrix<Type&source, const Matrix<int>
&swapvec);
template <class Type // Column Swap
Matrix<Typecolswap(const Matrix<Type&source, const Matrix<int>
&swapvec);
template <class Type // Row/Col Swap. Calls rowswap and then colswap
Matrix<Typerowcolswap(const Matrix<Type&source, const Matrix<int>
&swapvec);
template <class Type // Householder eigensolver for generalized eigen
problem
void eigGen(const Matrix<Type&A, const Matrix<Type&B,
Matrix<Type&eigval, Matrix<Type&eigvec);
template <class Type // Householder eigensolver for generalized eigen
problem
void eigGen(const Matrix<Type&A, const Matrix<Type&B,
Matrix<Type&eigval);
template <class Type // Householder eigensolver
void eig(const Matrix<Type&A, Matrix<Type&eigval, Matrix<Type>
&eigvec);
template <class Type // Householder eigensolver, without eigenvectors
void eig(const Matrix<Type&A, Matrix<Type&eigval);
template <class Type // mass normalize eigenvectors
Matrix<Typeeig_mass_normalize(const Matrix<Type&eigen,
const Matrix<Type&mass);
template <class Type // Householder reduction of a real, symmetric matrix
void tred2(const Matrix<Type&a, Matrix<Type&d, Matrix<Type&e);
template <class Type // eigenvalues, vectors of a real, symmetric,
tridiagonal matrix
void tqli(Matrix<Type&eigval, Matrix<Type&e, Matrix<Type&eigvec);
template <class Type // Cholesky decomposition
Matrix<Typechol(const Matrix<Type&source);
template <class Type // Cholesky decomposition, with p vector
Matrix<Typechol(const Matrix<Type&source, Matrix<Type&p);
template <class Type // Craig-Bampton reduction
void CraigBampton(const Matrix<Type&k, const Matrix<Type&m,
const Matrix<int&swapvec, int NgeneralizedDOF,
Matrix<Type&phiCB, Matrix<Type&kCB, Matrix<Type>
&mCB);

//Constructors and Destructor:

// Constructor for doubly subscripted array.
template <class Type>
Matrix<Type>::Matrix(int num_rows, int num_cols)
{
int i;
strcpy(Err, "");
strcpy(Tab, "\t");
NumRows = num_rows;
NumCols = num_cols;
Array = new Type *[NumRows+1];
for (i=0; i<=NumRows; ++i)
Array[i] = new Type [NumCols+1];
}

// Copy Constructor. Needed because constructor performs memory allocation
template <class Type>
Matrix<Type>::Matrix(const Matrix<Type&source)
{
int i, j;
strcpy(Err, source.Err);
strcpy(Tab, source.Tab);
NumRows = nrows(source);
NumCols = ncols(source);
Array = new Type *[NumRows+1];
for (i=0; i<=NumRows; ++i)
Array[i] = new Type [NumCols+1];
for (i=0; i<=NumRows; ++i)
for (j=0; j<=NumCols; ++j)
(*this)(i,j) = source(i,j);
}

// Destructor
template <class Type>
Matrix<Type>::~Matrix()
{
int i;
for (i=0; i<=NumRows; ++i)
{
delete [] Array[i];
Array[i] = 0;
}
delete [] Array;
Array = 0;
}

//Accessors:

// Set separator character for use with <<
template <class Type>
void Matrix<Type>::SetTab(const char *tabstr)
{
strcpy(Tab, tabstr);
}

// Get Error Code of Matrix object
template <class Type>
char *Matrix<Type>::GetErr()
{
return Err;
}

// Set Error Code of Matrix object

template <class Type>
void Matrix<Type>::SetErr(const char *error)
{
strcpy(Err, error);
}

// Overload () operator to access elements by reference
template <class Type>
Type &Matrix<Type>::operator() (int row, int col) const
{
int i, j;
if (row NumRows || col NumCols)
{
cerr << "(): out of range" << endl;
i = 0;
j = 0;
}
else
{
i = row;
j = col;
}
return Array[i][j];
}

// Access number of rows
template <class Type>
int nrows(const Matrix<Type&source)
{
return source.NumRows;
}

// Access number of columns
template <class Type>
int ncols(const Matrix<Type&source)
{
return source.NumCols;
}

//Operator Overloads:

// overload operator << for use with Matrix class
template <class Type>
ostream &operator<<(ostream &target, const Matrix<Type&source)
{
int i, j;
for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
target << source.Tab << source(i,j);
}
target << endl;
}
return target;
}

// overload operator >for use with Matrix class
template <class Type>
istream &operator>>(istream &source, Matrix<Type&target)
{
int i, j;
for (i=1; i<=nrows(target); ++i)
{
for (j=1; j<=ncols(target); ++j)
{
source >target(i,j);
}
}
return source;
}

// Overload assignment operator =
template <class Type>
Matrix<Type&Matrix<Type>::operator=(const Matrix<Type&source)
{
int i, j;
if (NumRows != nrows(source) && NumCols != ncols(source))
strcpy(Err, "=: size mismatch");
else
{
NumRows = nrows(source);
NumCols = ncols(source);
strcpy(Err, source.Err);
strcpy(Tab, source.Tab);
for (i=0; i<=NumRows; ++i)
for (j=0; j<=NumCols; ++j)
(*this)(i,j) = source(i,j);
}
return *this;
}

// Overload addition operator + for Matrix class
template <class Type>
Matrix<Typeoperator+(const Matrix<Type&left, const Matrix<Type&right)
{
int i, j;
Matrix<Typesum(nrows(left), ncols(left));

if (nrows(left) != nrows(right) || ncols(left) != ncols(right))
{
sum.SetErr("+: size mismatch");
}
else
{
for (i=0; i<=nrows(left); ++i)
for (j=0; j<=ncols(left); ++j)
sum(i,j) = left(i,j) + right(i,j);
}
return sum;
}

// Overload addition operator - for Matrix class
template <class Type>
Matrix<Typeoperator-(const Matrix<Type&left, const Matrix<Type&right)
{
int i, j;
Matrix<Typediff(nrows(left), ncols(left));

if (nrows(left) != nrows(right) || ncols(left) != ncols(right))
{
diff.SetErr("-: size mismatch");
}
else
{
for (i=1; i<=nrows(left); ++i)
for (j=1; j<=ncols(left); ++j)
diff(i,j) = left(i,j) - right(i,j);
}
return diff;
}

// Overload unary operator - for Matrix class
template <class Type>
Matrix<Typeoperator-(const Matrix<Type&source)
{
int i, j;
Matrix<Typetarget(nrows(source), ncols(source));

for (i=1; i<=nrows(source); ++i)
for (j=1; j<=ncols(source); ++j)
target(i,j) = -source(i,j);
return target;
}

// Overload multiply operator * for Matrix class
// implemented as a friend function in order to overload for pre- and post-
// multiplies: A * B; c * A; A * c; A * V; V * A.
template <class Type>
Matrix<Typeoperator*(const Matrix<Type&left, const Matrix<Type&right)
{
int i, j, k;
Matrix<Typeproduct(nrows(left),ncols(right));

if (ncols(left) != nrows(right))
{
product.SetErr("*: size mismatch");
}
else
{
for (i=0; i<=nrows(product); ++i)
for (k=0; k<=ncols(product); ++k)
product(i,k) = 0.;
for (i=1; i<=nrows(left); ++i)
for (j=1; j<=ncols(left); ++j)
for (k=1; k<=ncols(right); ++k)
product(i,k) += left(i,j) * right(j,k);
}
return product;
}

// c * A
template <class Type>
Matrix<Typeoperator*(double left, const Matrix<Type&right)
{
int i, j;
Matrix<Typeproduct(nrows(right),ncols(right));

for (i=1; i<=nrows(product); ++i)
for (j=1; j<=ncols(product); ++j)
product(i,j) = left * right(i,j);
return product;
}

// A * c
template <class Type>
Matrix<Typeoperator*(const Matrix<Type&left, double right)
{
Matrix<Typeproduct(nrows(left),ncols(left));

product = right * left; // use c * A
return product;
}

// A / c
template <class Type>
Matrix<Typeoperator/(const Matrix<Type&left, double right)
{
Matrix<Typetarget(nrows(left), ncols(left));
int i, j;

if (right == 0.)
{
target.SetErr("/: 0 divisor");
return target;
}
for (i=1; i<=nrows(left); ++i)
{
for (j=1; j<=ncols(left); ++j)
{
target(i,j) = left(i,j) / right;
}
}
return target;
}

// Member Functions

// cut a piece of this matrix starting at (r1, c1) that fills the target
Matrix
template <class Type>
void Matrix<Type>::cut(Matrix<Type&target, int r1, int c1) const
{
int i, j;

for (i=1; i<=nrows(target); ++i)
{
for (j=1; j<=ncols(target); ++j)
{
target(i,j) = (*this)(r1+i-1,c1+j-1);
}
}
}

// paste this matrix into the target starting at (r1, c1). Assumes target
is
// at least as large as this matrix
template <class Type>
void Matrix<Type>::paste(Matrix<Type&target, int r1, int c1) const
{
int i, j;

if (nrows(target) < (NumRows - (r1 - 1)) || ncols(target) < (NumCols -
(c1 - 1)))
{
target.SetErr("paste: Insufficient space in target");
}
else
{
for (i=1; i<=NumRows; ++i)
{
for (j=1; j<=NumCols; ++j)
{
target(r1+i-1,c1+j-1) = (*this)(i,j);
}
}
}
}

//Friend Functions:

// Identity Matrix
template <class Type>
void eye(Matrix<Type&a)
{
int i, j;
for (i=1; i<=nrows(a); ++i)
{
for (j=1; j<=nrows(a); ++j)
a(i,j) = Type(0);
a(i,i) = Type(1);
}
}

// Zero Matrix
template <class Type>
void zeros(Matrix<Type&a)
{
int i, j;
for (i=1; i<=nrows(a); ++i)
{
for (j=1; j<=ncols(a); ++j)
a(i,j) = Type(0);
}
}

// friend function to perform matrix transpose
template <class Type>
Matrix<Typetranspose(const Matrix<Type&source)
{
int i, j;
Matrix<Typetarget(ncols(source), nrows(source));
for (i=1; i<=nrows(source); ++i)
for (j=1; j<=ncols(source); ++j)
target(j,i) = source(i,j);
return target;
}

// Square root of matrix elements
template <class Type>
Matrix<Typesqrt(const Matrix<Type&source)
{
Matrix<Typetarget(nrows(source), ncols(source));
int i, j;

for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
if (source(i,j) < -0.01)
{
target.SetErr("sqrt: Negative value");
target(i,j) = -1.;
}
else if (source(i,j) >= -.01 && source(i,j) <= 1.e-20)
{
target(i,j) = 0.;
}
else
{
target(i,j) = sqrt(source(i,j));
}
}
}
return target;
}

// normalize matrix columns so that root sum square is unity
template <class Type>
Matrix<Typecolnormalize(const Matrix<Type&source)
{
Matrix<Typetarget(nrows(source), ncols(source));
int i, j;
double sumsquare;

for (i=1; i<=nrows(source); ++i)
{
sumsquare = 0;
for (j=1; j<=ncols(source); ++j)
{
sumsquare += (source(i,j) * source(i,j));
}
sumsquare = sqrt(sumsquare);
for (j=1; j<=ncols(source); ++j)
{
target(i,j) = source(i,j) / sumsquare;
}
}
return target;
}

// scale matrix such that row urow becomes unity
template <class Type>
Matrix<Typecolunitize(const Matrix<Type&source, int urow)
{
Matrix<Typetarget(nrows(source), ncols(source));
int i, j;

for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
if (source(urow,j) == 0)
{
}
else
{
target(i,j) = source(i,j) / source(urow,j);
}
}
}
return target;
}

// solve aX = b. Calls ludcmp to find LU(a) and then calls lubksb for each
// column in b.
template <class Type>
Matrix<Typesolve(const Matrix<Type&a, const Matrix<Type&b)
{
Matrix<Typetarget(nrows(b), ncols(b));
Matrix<Typelu(nrows(a), ncols(a));
Matrix<intindex(nrows(b));
Matrix<Typey(nrows(b));
Matrix<TypeX(nrows(b));
int i, j, d;
bool ierr;

if (nrows(a) != nrows(b))
{
target.SetErr("solve: size mismatch");
}
else
{
lu = ludcmp(a, ierr, index, d);
for (j=1; j<=ncols(b); ++j)
{
for (i=1; i<=nrows(b); ++i)
y(i) = b(i,j);
X = lubksb(lu, index, y);
for (i=1; i<=nrows(b); ++i)
target(i,j) = X(i);
}
}
return target;
}

// invert a matrix using lower upper decompostion with backsubstitution
template <class Type>
Matrix<Typeinv(const Matrix<Type&source)
{
Matrix<Typetarget(nrows(source), nrows(source));
Matrix<Typeident(nrows(source), nrows(source));

if (nrows(source) != ncols(source))
{
target.SetErr("inv: size mismatch");
}
else
{
eye(ident);
target = solve(source, ident);
}
return target;
}

// solves the set of linear equations aX=b, where a is LU(a).
template <class Type>
Matrix<Typelubksb(const Matrix<Type&a, Matrix<int&indx,
const Matrix<Type&b)
{
int i, ii, j, ll;
double sum;
Matrix<Typebtemp(nrows(b));

btemp = b;
ii = 0;
for (i=1; i<=nrows(a); ++i)
{
ll = indx(i);
sum = btemp(ll);
btemp(ll) = btemp(i);
if (ii != 0)
for (j=ii; j<i; ++j)
sum -= a(i,j)*btemp(j);
else if (sum != 0.)
ii = i;
btemp(i) = sum;
}
for (i=nrows(a); i>=1; --i)
{
sum = btemp(i);
for (j=i+1; j<=nrows(a); ++j)
sum -= a(i,j)*btemp(j);
btemp(i) = sum / a(i,i);
}
return btemp;
}

// lower upper decomposition of square matrix.
template <class Type>
Matrix<Typeludcmp(const Matrix<Type&asource, bool &ierr,
Matrix<int&indx, int &d)
{
Matrix<Typea(nrows(asource), nrows(asource));
int imax, i, j, k;
double tiny, aamax, dum, sum;
Matrix<Typevv(nrows(asource));
tiny = 1.e-20;

a = asource;
ierr = false;
d=1;
for (i=1; i<=nrows(a); ++i)
{
aamax = 0.;
for (j=1; j<= nrows(a); ++j)
if (fabsl(a(i,j)) aamax) aamax = fabsl(a(i,j));
if (aamax == 0.)
{
ierr = true;
return a;
}
vv(i) = 1./aamax;
}
for (j=1; j<=nrows(a); ++j)
{
for (i=1; i<j; ++i)
{
sum = a(i,j);
for (k=1; k<i; ++k)
sum -= a(i,k)*a(k,j);
a(i,j) = sum;
}
aamax = 0.;
imax = j;
for (i=j; i<=nrows(a); ++i)
{
sum = a(i,j);
for (k=1; k<j; ++k)
sum -= a(i,k) * a(k,j);
a(i,j) = sum;
dum = vv(i) * fabsl(sum);
if (dum >= aamax)
{
imax = i;
aamax = dum;
}
}
if (j != imax)
{
for (k=1; k<=nrows(a); ++k)
{
dum = a(imax,k);
a(imax,k) = a(j,k);
a(j,k) = dum;
}
d = -d;
vv(imax) = vv(j);
}
indx(j) = imax;
if (a(j,j) == 0.) a(j,j) = tiny;
if (j != nrows(a))
{
dum = 1./a(j,j);
for (i=j+1; i<=nrows(a); ++i)
a(i,j) = a(i,j) * dum;
}
}
return a;
}

// indexed version of heapsort. Returns index vector that contains indices
of
// the source matrix such that elements source(index(i)) are sorted in
ascending
// order
template <class Type>
Matrix<intindexsort(const Matrix<Type&source)
/*c
c Heapsort algorithm as presented in "Numerical Recipes" and described in
c "Sorting and Searching", vol. 3 of "The Art of Computer Programming" by
c Donald Knuth.
c
c Works only for column vectors, so far.
c*/
{
int i, j, L, indext, ir;
double q;
int n;
n = nrows(source);
Matrix<intindex(n);

for (j=1; j<=n; ++j)
{
index(j) = j;
}
if (n == 1)
{
return index;
}
L = n/2 + 1;
ir = n;
while (true)
{
if (L 1)
{
L--;
indext = index(L);
q = source(indext);
}
else
{
indext = index(ir);
q = source(indext);
index(ir) = index(1);
ir--;
if (ir == 1)
{
index(1) = indext;
return index;
}
}
i = L;
j = L + L;
while (j <= ir)
{
if (j < ir)
{
if (source(index(j)) < source(index(j+1))) j++;
}
if (q < source(index(j)))
{
index(i) = index(j);
i = j;
j += j;
}
else
{
j = ir + 1;
}
}
index(i) = indext;
}
}

// heapsort with index returns sorted matrix and also index by which source
is sorted
template <class Type>
Matrix<Typeheapsort(const Matrix<Type&source, Matrix<int&index)
{
index = indexsort(source);
return rowswap(source, index);
}

// heapsort w/o index returns only the sorted matrix
template <class Type>
Matrix<Typeheapsort(const Matrix<Type&source)
{
Matrix<intno_index(nrows(source));

return heapsort(source, no_index);
}

// Row swap
template <class Type>
Matrix<Typerowswap(const Matrix<Type&source, const Matrix<int&swapvec)
{
int i, j, k;
bool hit;
Matrix<intswapndx(nrows(source));
Matrix<Typetarget(nrows(source), ncols(source));

for (i=1; i<=nrows(swapvec); ++i)
{
for (j=1; j<=i-1; ++j)
{
if (swapvec(j) == swapvec(i))
{
target.SetErr("rowswap: double entry in swap vector");
return target;
}
}
swapndx(i) = swapvec(i);
}
k = 0;
for (i=1; i<=nrows(source); ++i)
{
hit = false;
for (j=1; j<=nrows(swapvec); ++j)
{
if (i==swapvec(j))
{
hit = true;
break;
}
}
if (!hit)
{
++k;
swapndx(nrows(swapvec)+ k) = i;
}
}
for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
target(i,j) = source(swapndx(i),j);
}
}
return target;
}

// Column Swap
template <class Type>
Matrix<Typecolswap(const Matrix<Type&source, const Matrix<int&swapvec)
{
int i, j, k;
bool hit;
Matrix<intswapndx(nrows(source));
Matrix<Typetarget(nrows(source), ncols(source));

for (i=1; i<=nrows(swapvec); ++i)
{
for (j=1; j<=i-1; ++j)
{
if (swapvec(j) == swapvec(i))
{
target.SetErr("colswap: double entry in swap vector");
return target;
}
}
swapndx(i) = swapvec(i);
}
k = 0;
for (j=1; j<=ncols(source); ++j)
{
hit = false;
for (i=1; i<=nrows(swapvec); ++i)
{
if (j==swapvec(i))
{
hit = true;
break;
}
}
if (!hit)
{
++k;
swapndx(nrows(swapvec)+ k) = j;
}
}
for (i=1; i<=nrows(source); ++i)
{
for (j=1; j<=ncols(source); ++j)
{
target(i,j) = source(i,swapndx(j));
}
}
return target;
}

// Row/Column Swap
template <class Type>
Matrix<Typerowcolswap(const Matrix<Type&source, const Matrix<int>
&swapvec)
{
Matrix<Typetarget(nrows(source), ncols(source));

target = rowswap(source, swapvec);
target = colswap(target, swapvec);
return target;
}

// Householder Eigensolver. Generalized eigen problem (Ax = lBx) with
eigenvectors.
template <class Type>
void eigGen(const Matrix<Type&A, const Matrix<Type&B, Matrix<Type>
&eigval,
Matrix<Type&eigvec)
{
Matrix<Typee(nrows(A));
Matrix<TypeLinvT(nrows(A), ncols(A));
Matrix<TypeH(nrows(A), ncols(A));

H = inv(chol(B));
LinvT = transpose(H);
H = H * A * LinvT;
eig(H, eigval, eigvec);
eigvec = LinvT * eigvec;
eigvec = eig_mass_normalize(eigvec, B);
}

// Householder Eigensolver. Generalized eigen problem (Ax = lBx) w/o
eigenvectors.
template <class Type>
void eigGen(const Matrix<Type&A, const Matrix<Type&B, Matrix<Type>
&eigval)
{
Matrix<Typeno_eigvec(nrows(A), ncols(A));

eigGen(A, B, eigval, no_eigvec);
}

// Householder Eigensolver. Eigen problem (Ax = lx) with eigenvectors.
// from Numerical Recipes
template <class Type>
void eig(const Matrix<Type&A, Matrix<Type&eigval, Matrix<Type&eigvec)
{
Matrix<Typee(nrows(A));
Matrix<intindex(nrows(A));

eigvec = A;
tred2(eigvec,eigval,e);
tqli(eigval,e,eigvec);
eigval = heapsort(eigval, index);
eigvec = colswap(eigvec, index);
}

// Householder Eigensolver. Eigen problem (Ax = lx) w/o eigenvectors.
// from Numerical Recipes
template <class Type>
void eig(const Matrix<Type&A, Matrix<Type&eigval)
{
Matrix<Typeno_eigvec(nrows(A), ncols(A));

eig(A, eigval, no_eigvec);
}

// Mass normalize eigenvectors
template <class Type>
Matrix<Typeeig_mass_normalize(const Matrix<Type&eigen,
const Matrix<Type&mass)
/*
This subroutine uses the methodology found in Paz, Structural Dynamics,
to scale the general eigenvectors (eigen) such that the resulting
eigenvectors (eigmn) are mass normalized.
*/
{
int n, i, j, k;
n = nrows(eigen);
Matrix<Typeeigmn(n,n);
double denom, temp;

for (j=1; j<=n; ++j)
{
denom = 0.;
for (k=1; k<=n; ++k)
{
temp = 0.;
for (i=1; i<=n; ++i)
{
temp = temp + eigen(i,j) * mass(i,k);
}
denom = denom + temp * eigen(k,j);
}
for (i=1; i<=n; ++i)
{
eigmn(i,j) = eigen(i,j) / denom;
}
}
return eigmn;
}

// Householder reduction
template <class Type>
void tred2(const Matrix<Type&a, Matrix<Type&d, Matrix<Type&e)
{
int i,j,k,l, n;
double f,g,h,hh,scale;

n = nrows(e);
for (i=n; i>=2;--i)
{
l=i-1;
h=0.;
scale=0.;
if(l>1)
{
for (k=1; k<=l; ++k)
{
scale=scale+fabs(a(i,k));
}
if(scale==0.)
{
e(i)=a(i,l);
}
else
{
for (k=1; k<=l; ++k)
{
a(i,k)=a(i,k)/scale;
h=h+a(i,k)*a(i,k);
}
f=a(i,l);
g=-(fabs(f)/f)*sqrt(h);
e(i)=scale*g;
h=h-f*g;
a(i,l)=f-g ;
f=0.;
for (j=1; j<=l; ++j)
{
//C Omit following line if finding only eigenvalues
a(j,i)=a(i,j)/h;
g=0.;
for (k=1; k<=j; ++k)
{
g=g+a(j,k)*a(i,k);
}
for (k=j+1; k<=l; ++k)
{
g=g+a(k,j)*a(i,k);
}
e(j)=g/h;
f=f+e(j)*a(i,j);
}
hh=f/(h+h);
for (j=1; j<=l; ++j)
{
f=a(i,j);
g=e(j)-hh*f;
e(j)=g;
for (k=1; k<=j; ++k)
{
a(j,k)=a(j,k)-f*e(k)-g*a(i,k);
}
}
}
}
else
{
e(i)=a(i,l);
}
d(i)=h;
}
//C Omit following line if finding only eigenvalues.
d(1)=0.;
e(1)=0.;
for (i=1; i<=n; ++i)
{
//C Delete lines from here ...
l=i-1;
if (d(i)!=0.)
{
for (j=1; j<=l; ++j)
{
g=0.;
for (k=1; k<=l; ++k)
{
g=g+a(i,k)*a(k,j);
}
for (k=1; k<=l; ++k)
{
a(k,j)=a(k,j)-g*a(k,i);
}
}
}
//C ... to here when finding only eigenvalues.
d(i)=a(i,i);
//C Also delete lines from here ...
a(i,i)=1.;
for (j=1; j<=l; ++j)
{
a(i,j)=0.;
a(j,i)=0.;
}
//C ... to here when finding only eigenvalues.
}
}

// eigenvalues, vectors of a real, symmetric, tridiagonal matrix
template <class Type>
void tqli(Matrix<Type&eigval, Matrix<Type&e, Matrix<Type&eigvec)
{
//CU USES pythag
int i,iter,k,l,m, n;
double b,c,dd,f,g,p,r,s;

n = nrows(e);
for (i=2; i<=n; ++i)
{
e(i-1)=e(i);
}
e(n)=0.;
for (l=1; l<=n; ++l)
{
iter=0;
label1:
for (m=l; m<=n-1; ++m)
{
dd=fabs(eigval(m))+fabs(eigval(m+1));
if (fabs(e(m))+dd==dd) goto label2;
}
m=n;
label2:
if(m!=l)
{
if(iter==30)
{
cerr << endl << "too many iterations in tqli" << endl;
exit(1);
}
iter=iter+1;
g=(eigval(l+1)-eigval(l))/(2.*e(l));
r=sqrt(g*g + 1.);
g=eigval(m)-eigval(l)+e(l)/(g+(fabs(g)/g)*r);
s=1.;
c=1.;
p=0.;
for (i=m-1; i>=l; --i)
{
f=s*e(i);
b=c*e(i);
r=sqrt(f*f + g*g);
e(i+1)=r;
if (r==0.)
{
eigval(i+1)=eigval(i+1)-p;
e(m)=0.;
goto label1;
}
s=f/r;
c=g/r;
g=eigval(i+1)-p;
r=(eigval(i)-g)*s+2.*c*b;
p=s*r;
eigval(i+1)=g+p;
g=c*r-b;
//C Omit lines from here ...
if (nrows(eigvec) 0)
{
for (k=1; k<=n; ++k)
{
f=eigvec(k,i+1);
eigvec(k,i+1)=s*eigvec(k,i)+c*f;
eigvec(k,i)=c*eigvec(k,i)-s*f;
}
}
//C ... to here when finding only eigenvalues.
}
eigval(l)=eigval(l)-p;
e(l)=g;
e(m)=0.;
goto label1;
}
}
}

// Cholesky decomposition.
template <class Type>
Matrix<Typechol(const Matrix<Type&source)
{
int i,j, n;
n = nrows(source);
Matrix<Typecholesky(n, n);
Matrix<Typep(n);

cholesky = chol(source, p);
for (i=1; i<=n; ++i)
{
for (j=1; j<=n; ++j)
{
if (i>j)
{
// cholesky(i,j)=source(i,j);
}
else if (i==j)
{
cholesky(i,j)=p(i);
}
else
{
cholesky(i,j)=0.;
}
}
}
return cholesky;
}

// Cholesky decomposition, returning p.
template <class Type>
Matrix<Typechol(const Matrix<Type&source, Matrix<Type&p)
{
int i,j,k, n;
double sum;
Matrix<Typetarget(nrows(source), ncols(source));

n = nrows(source);
target = source;
for (i=1; i<=n; ++i)
{
for (j=i; j<=n; ++j)
{
sum = source(i,j);
for (k=i-1; k>=1; --k)
{
sum = sum - target(i,k) * target(j,k);
}
if (i == j)
{
if(sum <= 0.)
{
cout << endl << "choldc failed: i = " << i;
cout << ", source(i,i) = " << source(i,i) << endl;
exit(1);
}
p(i) = sqrt(sum);
}
else
{
target(j,i) = sum / p(i);
}
}
}
return target;
}

// Craig-Bampton reduction
template <class Type>
void CraigBampton(const Matrix<Type&k, const Matrix<Type&m,
const Matrix<int&swapvec, int NgeneralizedDOF,
Matrix<Type&phiCB, Matrix<Type&kCB, Matrix<Type&mCB)
{
int N = nrows(k);
int Nswap = nrows(swapvec);
int Nreduced = Nswap + NgeneralizedDOF;
dmat kswap(N, N);
dmat mswap(N, N);
dmat kpg(N-Nswap, Nswap);
dmat kgg(N-Nswap, N-Nswap);
dmat mgg(N-Nswap, N-Nswap);
dmat I(Nswap, Nswap);
dmat Zero(Nswap, NgeneralizedDOF);
dmat phic(N-Nswap, Nswap);
dmat phin(N-Nswap, N-Nswap);
dmat phin_red(N-Nswap, NgeneralizedDOF);
dmat phival(N-Nswap);

kswap = rowcolswap(k, swapvec);
mswap = rowcolswap(m, swapvec);
kswap.cut(kgg, Nswap+1, Nswap+1);
kswap.cut(kpg, Nswap+1, 1);
mswap.cut(mgg, Nswap+1, Nswap+1);
eye(I);
zeros(Zero);
phic = -inv(kgg) * kpg;
eigGen(kgg, mgg, phival, phin);
phin.cut(phin_red, 1, 1);
phin_red = colunitize(phin_red, Nreduced);
I.paste(phiCB, 1, 1);
Zero.paste(phiCB, 1, Nswap+1);
phic.paste(phiCB, Nswap+1, 1);
phin_red.paste(phiCB, Nswap+1, Nswap+1);
kCB = transpose(phiCB) * kswap * phiCB;
mCB = transpose(phiCB) * mswap * phiCB;
}
#endif
Nov 30 '06 #2
I developed a template matrix class back around 1992 using Borland C++ 4.5
(ancestor of C++ Builder) and haven't touched it until a few days ago. I
pulled it from the freezer and thawed it out. I built a console app using
Microsoft Visual C++ 6 (VC++) and it worked great. Only one line in the
header file had to be commented out.

I built a console app using Borland C++ Builder 5. The linker complained of
references to external functions (my friend functions that were referenced)
that were not in the object file. I then used the conversion tool supplied
with the IDE to convert the VC++ project to a Borland BPR. Then the linker
complained that ODBCCP32.LIB could not be opened. I included the folder for
ODBCCP32.LIB in the linker library path and tried again, unsuccessfully.

I include the source code with output from the VC++ executable. I hope
someone can help me get this fine piece of code :) working again.

Thanks,
Hamilton Woods
Hello,

after making some changes to "TMatrix.h", it compiled using
C++ Builder 4.0.

check also:
http://www.parashift.com/c++-faq-lit...html#faq-35.16

Here are the first lines of the matrix class, the rest is
unchanged.

//---------- Beginning of TMatrix.h ----------
// tmatrix.h
// definition of Matrix class

#if !defined _TMatrix
#define _TMatrix
#include <math.h>
#include <string.h>
#include <iostream.h>
#include <stdlib.h>

#ifdef _MSC_VER
//**** ghw2006-11-27 Microsoft VC++ now has a bool type
// typedef int bool;
//****
#define true 1
#define false 0
#endif

#define cmat Matrix<char>
#define dmat Matrix<double>
#define fmat Matrix<float>
#define imat Matrix<int>
#define lmat Matrix<long int>
#define smat Matrix<short int>

//forward declaration
template <class Type>
class Matrix;

//forward declaration of function template
template <class Type>
int nrows(const Matrix<Type&source);

template <class Type>
int ncols(const Matrix<Type&source);

//forward declaration of function template operator<<
template <class Type>
ostream &operator<< (ostream &target, const Matrix<Type&source);
template <class Type>
istream &operator>(istream &source, Matrix<Type&target);

template <class Type>
Matrix<Type // Matrix Addition
operator+(const Matrix<Type&left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&source);
template <class Type>
Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Multiplication
operator*(double left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, double right);
template <class Type>
Matrix<Type // Matrix Division
operator/(const Matrix<Type&left, double right);

template <class Type>
class Matrix
{
public:
// Constructors & Destructor
Matrix(int num_rows, int num_cols = 1);
Matrix(const Matrix<Type&a);
~Matrix();
// Accessors
void SetTab(const char *tabstr);
char *GetErr();
void SetErr(const char *error);
friend int nrows<Type>(const Matrix<Type&source); // added <Type>
friend int ncols<Type>(const Matrix<Type&source);
Type &operator() (int row, int col = 1) const;
// Type &operator() (int row) const;
// Operator Overloads
friend ostream &operator<< <Type>(ostream &target, const Matrix<Type&source);
friend istream &operator><Type>(istream &source, Matrix<Type&target);

//Matrix<Type& Matrix<Type>::operator=(const Matrix<Type&source);
// Assignment // did not compile

Matrix & operator=(const Matrix<Type&source) ; // changed

friend Matrix<Type // Matrix Addition
operator+<Type>(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-<Type>(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-<Type>(const Matrix<Type&source);
friend Matrix<Type // Matrix Multiplication
operator*<Type>(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*<Type>(double left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*<Type>(const Matrix<Type&left, double right);
friend Matrix<Type // Matrix Division
operator/<Type>(const Matrix<Type&left, double right);
// Member Functions
void cut(Matrix<Type&target, int r1, int c1 = 1) const;
void paste(Matrix<Type&target, int r1, int c1 = 1) const;
protected:
int NumRows;
int NumCols;
Type **Array;
char Tab[5]; // delimiter used for << operator
char Err[21]; // string to hold error message
};
// ...

Nov 30 '06 #3
Wow! Awesome!

Thanks for your help. I looked up template friend functions in my 1991
version of Stroustrup's book. It was there, although not exceptionally
clear. I was beginning to think the language had changed along the way. I
guess Microsoft lets a developer be sloppy, Borland doesn't. Anyway, the
old code now works with your corrections.

Thanks,
Hamilton Woods

"Christoph Flohr" <nomail@invalidwrote in message
news:ek*************@news.t-online.com...
I developed a template matrix class back around 1992 using Borland C++
4.5
(ancestor of C++ Builder) and haven't touched it until a few days ago.
I
pulled it from the freezer and thawed it out. I built a console app
using
Microsoft Visual C++ 6 (VC++) and it worked great. Only one line in the
header file had to be commented out.

I built a console app using Borland C++ Builder 5. The linker
complained of
references to external functions (my friend functions that were
referenced)
that were not in the object file. I then used the conversion tool
supplied
with the IDE to convert the VC++ project to a Borland BPR. Then the
linker
complained that ODBCCP32.LIB could not be opened. I included the folder
for
ODBCCP32.LIB in the linker library path and tried again, unsuccessfully.

I include the source code with output from the VC++ executable. I hope
someone can help me get this fine piece of code :) working again.

Thanks,
Hamilton Woods

Hello,

after making some changes to "TMatrix.h", it compiled using
C++ Builder 4.0.

check also:
http://www.parashift.com/c++-faq-lit...html#faq-35.16

Here are the first lines of the matrix class, the rest is
unchanged.

//---------- Beginning of TMatrix.h ----------
// tmatrix.h
// definition of Matrix class

#if !defined _TMatrix
#define _TMatrix
#include <math.h>
#include <string.h>
#include <iostream.h>
#include <stdlib.h>

#ifdef _MSC_VER
//**** ghw2006-11-27 Microsoft VC++ now has a bool type
// typedef int bool;
//****
#define true 1
#define false 0
#endif

#define cmat Matrix<char>
#define dmat Matrix<double>
#define fmat Matrix<float>
#define imat Matrix<int>
#define lmat Matrix<long int>
#define smat Matrix<short int>

//forward declaration
template <class Type>
class Matrix;

//forward declaration of function template
template <class Type>
int nrows(const Matrix<Type&source);

template <class Type>
int ncols(const Matrix<Type&source);

//forward declaration of function template operator<<
template <class Type>
ostream &operator<< (ostream &target, const Matrix<Type&source);
template <class Type>
istream &operator>(istream &source, Matrix<Type&target);

template <class Type>
Matrix<Type // Matrix Addition
operator+(const Matrix<Type&left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Subtraction
operator-(const Matrix<Type&source);
template <class Type>
Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Multiplication
operator*(double left, const Matrix<Type&right);
template <class Type>
Matrix<Type // Matrix Multiplication
operator*(const Matrix<Type&left, double right);
template <class Type>
Matrix<Type // Matrix Division
operator/(const Matrix<Type&left, double right);

template <class Type>
class Matrix
{
public:
// Constructors & Destructor
Matrix(int num_rows, int num_cols = 1);
Matrix(const Matrix<Type&a);
~Matrix();
// Accessors
void SetTab(const char *tabstr);
char *GetErr();
void SetErr(const char *error);
friend int nrows<Type>(const Matrix<Type&source); // added <Type>
friend int ncols<Type>(const Matrix<Type&source);
Type &operator() (int row, int col = 1) const;
// Type &operator() (int row) const;
// Operator Overloads
friend ostream &operator<< <Type>(ostream &target, const Matrix<Type>
&source);
friend istream &operator><Type>(istream &source, Matrix<Type&target);

//Matrix<Type& Matrix<Type>::operator=(const Matrix<Type&source);
// Assignment // did not compile

Matrix & operator=(const Matrix<Type&source) ; // changed

friend Matrix<Type // Matrix Addition
operator+<Type>(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-<Type>(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Subtraction
operator-<Type>(const Matrix<Type&source);
friend Matrix<Type // Matrix Multiplication
operator*<Type>(const Matrix<Type&left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*<Type>(double left, const Matrix<Type&right);
friend Matrix<Type // Matrix Multiplication
operator*<Type>(const Matrix<Type&left, double right);
friend Matrix<Type // Matrix Division
operator/<Type>(const Matrix<Type&left, double right);
// Member Functions
void cut(Matrix<Type&target, int r1, int c1 = 1) const;
void paste(Matrix<Type&target, int r1, int c1 = 1) const;
protected:
int NumRows;
int NumCols;
Type **Array;
char Tab[5]; // delimiter used for << operator
char Err[21]; // string to hold error message
};
// ...

Dec 1 '06 #4

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

17
by: Alexander Stippler | last post by:
Hi, what do I have to do to get this (incorrect) piece of code to work. The specialization is wrong, but how can I do it? template <typename T, typename V> class Mask { public: Mask(int i)...
7
by: Thomas Matthews | last post by:
Hi, I am converting my table and record classes into templates. My issue is the syntax of declaring a friend class within the template. I have searched the C++ FAQ Lite (web), the C++...
6
by: Alex | last post by:
I've a problem with following simple code (not really usefull): //--------------------------------------------------------------------------- #pragma hdrstop #include "stdio.h" #include...
5
by: Levent | last post by:
Hi, Why doesn't this work? (tried with gcc 3.3.3 and VC++ 7.1): #include <iostream> template<class T, unsigned N> struct Foo { void func(); }; template<class T, unsigned N>
5
by: Hari | last post by:
Guys please help me to solve this strange problem what Iam getting as follows.. Trying to instantiate a global instance of a template class as follows :- when i build this code with debug and...
11
by: Niels Dekker - no reply address | last post by:
The following attempt to pass my template "Base" as a template template argument was rejected by Microsoft VC++ 8.0 (2005), while it still works on VC++ 7.1 (2003). Is it correct C++? And is...
1
by: Leslaw Bieniasz | last post by:
Hello, I have the following problem: file a.h --------------- template <class T> class A { // some stuff
17
by: Fabry | last post by:
Hi All, I'm new of this group and I do not know if this is the correct group for my question. I have a DLL with its export library (.lib) wrote in Borland C++ 6. In borland everything is OK and...
4
by: StephQ | last post by:
According to: http://www.parashift.com/c++-faq-lite/templates.html#faq-35.4 , if my understanding is correct, in template<typename T> class Foo { friend void func (const Foo<T>& foo); }; ...
3
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 3 Jan 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). For other local times, please check World Time Buddy In...
0
by: jianzs | last post by:
Introduction Cloud-native applications are conventionally identified as those designed and nurtured on cloud infrastructure. Such applications, rooted in cloud technologies, skillfully benefit from...
0
by: mar23 | last post by:
Here's the situation. I have a form called frmDiceInventory with subform called subfrmDice. The subform's control source is linked to a query called qryDiceInventory. I've been trying to pick up the...
2
by: jimatqsi | last post by:
The boss wants the word "CONFIDENTIAL" overlaying certain reports. He wants it large, slanted across the page, on every page, very light gray, outlined letters, not block letters. I thought Word Art...
2
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 7 Feb 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:30 (7.30PM). In this month's session, the creator of the excellent VBE...
0
by: fareedcanada | last post by:
Hello I am trying to split number on their count. suppose i have 121314151617 (12cnt) then number should be split like 12,13,14,15,16,17 and if 11314151617 (11cnt) then should be split like...
0
by: stefan129 | last post by:
Hey forum members, I'm exploring options for SSL certificates for multiple domains. Has anyone had experience with multi-domain SSL certificates? Any recommendations on reliable providers or specific...
1
by: davi5007 | last post by:
Hi, Basically, I am trying to automate a field named TraceabilityNo into a web page from an access form. I've got the serial held in the variable strSearchString. How can I get this into the...
0
by: MeoLessi9 | last post by:
I have VirtualBox installed on Windows 11 and now I would like to install Kali on a virtual machine. However, on the official website, I see two options: "Installer images" and "Virtual machines"....

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.