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

Discussion in 'C++' started by Hamilton Woods, Nov 30, 2006.

  1. 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
     
    Hamilton Woods, Nov 30, 2006
    #1
    1. Advertising

  2. 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
     
    Oliver Bleckmann, Nov 30, 2006
    #2
    1. Advertising

  3. > 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-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
    };
    // ...
     
    Christoph Flohr, Nov 30, 2006
    #3
  4. 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$-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-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
    > };
    > // ...
    >
     
    Hamilton Woods, Dec 1, 2006
    #4
    1. Advertising

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

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. christopher diggins
    Replies:
    16
    Views:
    785
    Pete Becker
    May 4, 2005
  2. LiDongning
    Replies:
    2
    Views:
    548
    red floyd
    Jul 2, 2009
  3. Phlip
    Replies:
    5
    Views:
    591
    Stefan Behnel
    Jan 13, 2010
  4. A L
    Replies:
    1
    Views:
    530
    Alf P. Steinbach /Usenet
    Aug 25, 2010
  5. abargaddon
    Replies:
    1
    Views:
    218
    clintmazur
    Feb 4, 2008
Loading...

Share This Page