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

H

Hamilton Woods

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<Type> transpose(const Matrix<Type> &a);
template <class Type> // element by element square root (only if all
elements are positive)
Matrix<Type> sqrt(const Matrix<Type> &source);
template <class Type> // normalize matrix
Matrix<Type> colnormalize(const Matrix<Type> &source);
template <class Type> // normalize matrix
Matrix<Type> colunitize(const Matrix<Type> &source, int urow = 1);
template <class Type> // solve aX = b using lu decomposition w/
backsubstitution
Matrix<Type> solve(const Matrix<Type> &a, const Matrix<Type> &b);
template <class Type> // Matrix Inverse. aX = I
Matrix<Type> inv(const Matrix<Type> &a);
template <class Type> // lower upper backsubstitution
Matrix<Type> lubksb (const Matrix<Type> &a, Matrix<int> &indx,
const Matrix<Type> &b);
template <class Type> // lower upper decomposition
Matrix<Type> ludcmp(const Matrix<Type> &asource, bool &ierr,
Matrix<int> &indx, int &d);
template <class Type> // indexsort. Works on column vectors only, for
now.
Matrix<int> indexsort(const Matrix<Type> &source);
template <class Type> // heapsort. Works on column vectors only, for now.
Matrix<Type> heapsort(const Matrix<Type> &source, Matrix<int> &index);
template <class Type> // heapsort. Does not return sort index.
Matrix<Type> heapsort(const Matrix<Type> &source);
template <class Type> // Row Swap
Matrix<Type> rowswap(const Matrix<Type> &source, const Matrix<int>
&swapvec);
template <class Type> // Column Swap
Matrix<Type> colswap(const Matrix<Type> &source, const Matrix<int>
&swapvec);
template <class Type> // Row/Col Swap. Calls rowswap and then colswap
Matrix<Type> rowcolswap(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<Type> eig_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<Type> chol(const Matrix<Type> &source);
template <class Type> // Cholesky decomposition, with p vector
Matrix<Type> chol(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 = 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 = 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;
Array = 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>::eek:perator() (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[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>::eek:perator=(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<Type> operator+(const Matrix<Type> &left, const Matrix<Type> &right)
{
int i, j;
Matrix<Type> sum(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<Type> operator-(const Matrix<Type> &left, const Matrix<Type> &right)
{
int i, j;
Matrix<Type> diff(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<Type> operator-(const Matrix<Type> &source)
{
int i, j;
Matrix<Type> target(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<Type> operator*(const Matrix<Type> &left, const Matrix<Type> &right)
{
int i, j, k;
Matrix<Type> product(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<Type> operator*(double left, const Matrix<Type> &right)
{
int i, j;
Matrix<Type> product(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<Type> operator*(const Matrix<Type> &left, double right)
{
Matrix<Type> product(nrows(left),ncols(left));

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

// A / c
template <class Type>
Matrix<Type> operator/(const Matrix<Type> &left, double right)
{
Matrix<Type> target(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<Type> transpose(const Matrix<Type> &source)
{
int i, j;
Matrix<Type> target(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<Type> sqrt(const Matrix<Type> &source)
{
Matrix<Type> target(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<Type> colnormalize(const Matrix<Type> &source)
{
Matrix<Type> target(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<Type> colunitize(const Matrix<Type> &source, int urow)
{
Matrix<Type> target(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<Type> solve(const Matrix<Type> &a, const Matrix<Type> &b)
{
Matrix<Type> target(nrows(b), ncols(b));
Matrix<Type> lu(nrows(a), ncols(a));
Matrix<int> index(nrows(b));
Matrix<Type> y(nrows(b));
Matrix<Type> X(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<Type> inv(const Matrix<Type> &source)
{
Matrix<Type> target(nrows(source), nrows(source));
Matrix<Type> ident(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<Type> lubksb(const Matrix<Type> &a, Matrix<int> &indx,
const Matrix<Type> &b)
{
int i, ii, j, ll;
double sum;
Matrix<Type> btemp(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<Type> ludcmp(const Matrix<Type> &asource, bool &ierr,
Matrix<int> &indx, int &d)
{
Matrix<Type> a(nrows(asource), nrows(asource));
int imax, i, j, k;
double tiny, aamax, dum, sum;
Matrix<Type> vv(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<int> indexsort(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<int> index(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<Type> heapsort(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<Type> heapsort(const Matrix<Type> &source)
{
Matrix<int> no_index(nrows(source));

return heapsort(source, no_index);
}

// Row swap
template <class Type>
Matrix<Type> rowswap(const Matrix<Type> &source, const Matrix<int> &swapvec)
{
int i, j, k;
bool hit;
Matrix<int> swapndx(nrows(source));
Matrix<Type> target(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<Type> colswap(const Matrix<Type> &source, const Matrix<int> &swapvec)
{
int i, j, k;
bool hit;
Matrix<int> swapndx(nrows(source));
Matrix<Type> target(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<Type> rowcolswap(const Matrix<Type> &source, const Matrix<int>
&swapvec)
{
Matrix<Type> target(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<Type> e(nrows(A));
Matrix<Type> LinvT(nrows(A), ncols(A));
Matrix<Type> H(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<Type> no_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<Type> e(nrows(A));
Matrix<int> index(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<Type> no_eigvec(nrows(A), ncols(A));

eig(A, eigval, no_eigvec);
}

// Mass normalize eigenvectors
template <class Type>
Matrix<Type> eig_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<Type> eigmn(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<Type> chol(const Matrix<Type> &source)
{
int i,j, n;
n = nrows(source);
Matrix<Type> cholesky(n, n);
Matrix<Type> p(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<Type> chol(const Matrix<Type> &source, Matrix<Type> &p)
{
int i,j,k, n;
double sum;
Matrix<Type> target(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
 
O

Oliver Bleckmann

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<Type> transpose(const Matrix<Type> &a);
template <class Type> // element by element square root (only if all
elements are positive)
Matrix<Type> sqrt(const Matrix<Type> &source);
template <class Type> // normalize matrix
Matrix<Type> colnormalize(const Matrix<Type> &source);
template <class Type> // normalize matrix
Matrix<Type> colunitize(const Matrix<Type> &source, int urow = 1);
template <class Type> // solve aX = b using lu decomposition w/
backsubstitution
Matrix<Type> solve(const Matrix<Type> &a, const Matrix<Type> &b);
template <class Type> // Matrix Inverse. aX = I
Matrix<Type> inv(const Matrix<Type> &a);
template <class Type> // lower upper backsubstitution
Matrix<Type> lubksb (const Matrix<Type> &a, Matrix<int> &indx,
const Matrix<Type> &b);
template <class Type> // lower upper decomposition
Matrix<Type> ludcmp(const Matrix<Type> &asource, bool &ierr,
Matrix<int> &indx, int &d);
template <class Type> // indexsort. Works on column vectors only, for
now.
Matrix<int> indexsort(const Matrix<Type> &source);
template <class Type> // heapsort. Works on column vectors only, for now.
Matrix<Type> heapsort(const Matrix<Type> &source, Matrix<int> &index);
template <class Type> // heapsort. Does not return sort index.
Matrix<Type> heapsort(const Matrix<Type> &source);
template <class Type> // Row Swap
Matrix<Type> rowswap(const Matrix<Type> &source, const Matrix<int>
&swapvec);
template <class Type> // Column Swap
Matrix<Type> colswap(const Matrix<Type> &source, const Matrix<int>
&swapvec);
template <class Type> // Row/Col Swap. Calls rowswap and then colswap
Matrix<Type> rowcolswap(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<Type> eig_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<Type> chol(const Matrix<Type> &source);
template <class Type> // Cholesky decomposition, with p vector
Matrix<Type> chol(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 = 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 = 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;
Array = 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>::eek:perator() (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[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>::eek:perator=(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<Type> operator+(const Matrix<Type> &left, const Matrix<Type> &right)
{
int i, j;
Matrix<Type> sum(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<Type> operator-(const Matrix<Type> &left, const Matrix<Type> &right)
{
int i, j;
Matrix<Type> diff(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<Type> operator-(const Matrix<Type> &source)
{
int i, j;
Matrix<Type> target(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<Type> operator*(const Matrix<Type> &left, const Matrix<Type> &right)
{
int i, j, k;
Matrix<Type> product(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<Type> operator*(double left, const Matrix<Type> &right)
{
int i, j;
Matrix<Type> product(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<Type> operator*(const Matrix<Type> &left, double right)
{
Matrix<Type> product(nrows(left),ncols(left));

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

// A / c
template <class Type>
Matrix<Type> operator/(const Matrix<Type> &left, double right)
{
Matrix<Type> target(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<Type> transpose(const Matrix<Type> &source)
{
int i, j;
Matrix<Type> target(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<Type> sqrt(const Matrix<Type> &source)
{
Matrix<Type> target(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<Type> colnormalize(const Matrix<Type> &source)
{
Matrix<Type> target(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<Type> colunitize(const Matrix<Type> &source, int urow)
{
Matrix<Type> target(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<Type> solve(const Matrix<Type> &a, const Matrix<Type> &b)
{
Matrix<Type> target(nrows(b), ncols(b));
Matrix<Type> lu(nrows(a), ncols(a));
Matrix<int> index(nrows(b));
Matrix<Type> y(nrows(b));
Matrix<Type> X(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<Type> inv(const Matrix<Type> &source)
{
Matrix<Type> target(nrows(source), nrows(source));
Matrix<Type> ident(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<Type> lubksb(const Matrix<Type> &a, Matrix<int> &indx,
const Matrix<Type> &b)
{
int i, ii, j, ll;
double sum;
Matrix<Type> btemp(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<Type> ludcmp(const Matrix<Type> &asource, bool &ierr,
Matrix<int> &indx, int &d)
{
Matrix<Type> a(nrows(asource), nrows(asource));
int imax, i, j, k;
double tiny, aamax, dum, sum;
Matrix<Type> vv(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<int> indexsort(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<int> index(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<Type> heapsort(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<Type> heapsort(const Matrix<Type> &source)
{
Matrix<int> no_index(nrows(source));

return heapsort(source, no_index);
}

// Row swap
template <class Type>
Matrix<Type> rowswap(const Matrix<Type> &source, const Matrix<int> &swapvec)
{
int i, j, k;
bool hit;
Matrix<int> swapndx(nrows(source));
Matrix<Type> target(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<Type> colswap(const Matrix<Type> &source, const Matrix<int> &swapvec)
{
int i, j, k;
bool hit;
Matrix<int> swapndx(nrows(source));
Matrix<Type> target(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<Type> rowcolswap(const Matrix<Type> &source, const Matrix<int>
&swapvec)
{
Matrix<Type> target(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<Type> e(nrows(A));
Matrix<Type> LinvT(nrows(A), ncols(A));
Matrix<Type> H(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<Type> no_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<Type> e(nrows(A));
Matrix<int> index(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<Type> no_eigvec(nrows(A), ncols(A));

eig(A, eigval, no_eigvec);
}

// Mass normalize eigenvectors
template <class Type>
Matrix<Type> eig_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<Type> eigmn(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<Type> chol(const Matrix<Type> &source)
{
int i,j, n;
n = nrows(source);
Matrix<Type> cholesky(n, n);
Matrix<Type> p(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<Type> chol(const Matrix<Type> &source, Matrix<Type> &p)
{
int i,j,k, n;
double sum;
Matrix<Type> target(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
 
C

Christoph Flohr

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
Click to expand...

Hello,

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

check also:
http://www.parashift.com/c++-faq-lite/templates.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>::eek:perator=(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
};
// ...
 
H

Hamilton Woods

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 said:
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
Click to expand...

Hello,

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

check also:
http://www.parashift.com/c++-faq-lite/templates.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>::eek:perator=(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
};
// ...
Click to expand...
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Similar Threads

Difference between using "let" in a "for" loop 0
How to multiply two matrices of size in using inline assembly in C++ 2
Need Help: Program to Accept 2 Matrices and Show their Sum 0
Java matrix problem 3
Custom matrix multiplication produces different results to glm 0
Filter sober in c++ don't pass test 0
matrix mult with const in c++ 0
Could you fix my code please? (i get no output after inputing the number) 1

Members online

No members online now.
Total: 34 (members: 1, guests: 33)
Robots: 599

Forum statistics

Threads
473,769
Messages
2,569,578
Members
45,052
Latest member
LucyCarper

Latest Threads

Top