Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C++ > Old template class works in VC++ Not in C++ Builder 5

Reply
Thread Tools

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

 
 
Hamilton Woods
Guest
Posts: n/a
 
      11-30-2006
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(;

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[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>: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[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>: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>:aste(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



 
Reply With Quote
 
 
 
 
Oliver Bleckmann
Guest
Posts: n/a
 
      11-30-2006
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(;
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[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>: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[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>: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>:aste(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


 
Reply With Quote
 
 
 
 
Christoph Flohr
Guest
Posts: n/a
 
      11-30-2006
> 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>: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
};
// ...

 
Reply With Quote
 
Hamilton Woods
Guest
Posts: n/a
 
      12-01-2006
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@invalid> wrote in message
news:eknlug$up6$00$(E-Mail Removed)-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>: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
> };
> // ...
>



 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Declaring a template class with two template params a friend in anon-template class A L C++ 1 08-25-2010 07:25 AM
Template function within class template not not resolving when passedfunction local class? jrwats C++ 1 03-01-2010 08:30 AM
Where's a DOM builder that uses the Builder Pattern to ... buildDOMs? Phlip Python 5 01-13-2010 12:48 PM
class nested in template class: Dev-C works while VC6 not LiDongning C++ 2 07-02-2009 03:44 PM
A parameterized class (i.e. template class / class template) is not a class? christopher diggins C++ 16 05-04-2005 12:26 AM



Advertisments