// Copyright (C) 1996 DIMACS Center, Rutgers, The State University of New Jersey
// Author(s): ??, Jonathan Berry

// This software is copyrighted by the DIMACS Center at Rutgers, The State
// University of New Jersey.  IT IS PROVIDED AS IS, AND THE AUTHORS, DIMACS, AND
// RUTGERS, THE STATE UNIVERSITY OF NEW JERSEY  DISCLAIM
// ALL LIABILITY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
// DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY
// DERIVATIVES THEREOF, EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

// THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,
// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE
// IS PROVIDED ON AN "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE
// NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
// MODIFICATIONS.

// The authors hereby grant permission to use, copy, modify, distribute,
// and license this software and its documentation for any purpose, provided
// that existing copyright notices are retained in all copies and that this
// notice is included verbatim in any distributions. No written agreement,
// license, or royalty fee is required for any of the authorized uses.
// Modifications to this software may be copyrighted by their authors
// and need not follow the licensing terms described here, provided that
// the new terms are clearly indicated on the first page of each file where
// they apply.

// Last File Update: 31-Jul-1996
// 

#ifndef MATRIX_CC
#define MATRIX_CC

#include <fstream.h>
#include <stdlib.h>
#include <stdio.h> 

/////////////////////////////////////////////////////////////////////////////
//
//  asymmetric matrix constructors and destructor 
//
/////////////////////////////////////////////////////////////////////////////

template <class Item>
AsymMatrix<Item>::AsymMatrix(int numRows, int numCols) : 
				Matrix<Item>(numRows, numCols) 
{
    if ( (numRows <= 0) || (numCols <= 0) ) 
        cerr << "Matrices must have at least 1 row and column.\n";
    Elements = (Item**) calloc(numRows, sizeof(Item*));
    for (int i=0; i < numRows; i++) 
        Elements[i] = (Item*) calloc(numCols, sizeof(Item));
}


template <class Item>
AsymMatrix<Item>::AsymMatrix(const AsymMatrix<Item>& A) : 
				Matrix<Item>(A.numrows, A.numcols)
{
    Elements = (Item**) calloc(A.numrows, sizeof(Item*));
    int i, j;
    for (i=0; i < A.numrows; i++) 
        Elements[i] = (Item*) calloc(A.numcols, sizeof(Item));
    Item*ptr, *temp;
    for (i=0; i < numrows; i++) {
        ptr = Elements[i];
        temp = A.Elements[i];
        for (j=0; j < numcols; j++, ptr++, temp++) 
            *ptr = *temp;
    }
}


template <class Item>
AsymMatrix<Item>::~AsymMatrix()
{
    for (int i=0; i < numrows; i++) 
        free((char *) Elements[i]);
    free((char *) Elements);
} 


//
// Saves an AsymMatrix to a file specified by the user 
//
template <class Item>
void
AsymMatrix<Item>::save(char * filename)
{
    ofstream curFile;
    curFile.open(filename, ios::out);
    if (!curFile) {
        cerr << "Error saving to file: " << filename << ".\n";
        exit(-1);
    }
    curFile << "A " << numrows << ", " << numcols << "\n";
    Item*ptr;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        for (int j=0; j < numcols; j++, ptr++) {
            curFile << *ptr;
            if (j != numcols-1)   
                curFile << ",";
        }
        curFile << ";\n";
    }
}    

template <class Item>
ostream&
AsymMatrix<Item>::display(ostream& os) const
{
    Item*ptr;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        for (int j=0; j < numcols; j++, ptr++) {
            os << *ptr;
            if (j != numcols-1)   
                os << " ";
        }
        os << "\n";
    }
    return os;
}


//
// Loads an AsymMatrix from a file specified by the user    
//
template <class Item>
void
AsymMatrix<Item>::load(char * filename)
{
    ifstream curFile;
    curFile.open(filename, ios::in);
    if (!curFile) {
        cerr << "Error loading from file: " << filename << ".\n";
        exit(-1);
    }
    char ch;
    int numRows, numCols;
    curFile.get(ch);
    if (ch != 'A') { 
        cerr << "Matrices are of different types.  We cannot complete load.\n";
        exit(-1);
    }
    curFile >> numRows;
    curFile.get(ch);
    curFile >> numCols;
    curFile.get(ch);
    if ((numRows != numrows) || (numCols != numcols)) { 
        cerr << "Matrices are of different types.  We cannot complete load.\n";
        exit(-1);
    }
    Item*ptr;
    for (int i=0; i < numRows; i++) {
        ptr = Elements[i];
        for (int j=0; j < numCols; j++, ptr++) {
            curFile >> *ptr;
            if (j != numCols-1) 
                curFile.get(ch);
        }
        curFile.get(ch).get(ch);
    }
}


//
// Returns element i,j.  Can also be used to set element i,j 
//
template <class Item>
Item&
AsymMatrix<Item>::operator() (int row, int col) const
{
    if (( row >= numrows) || (col >= numcols)) 
        cerr << "Out of bounds error.\n";
    return Elements[row][col];
}


//
// returns element i,j
//
//template <class Item>
//Item
//AsymMatrix<Item>::operator() (int row, int col) const
//{
//    return Elements[row][col];
//}


//
//  multiplies two asymmetric matrices
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::operator*(const AsymMatrix& A)
{
    if (numcols != A.numrows) 
        cerr << "Matrices not compatable.\n";
    AsymMatrix<Item> result(numrows, A.numcols);
    Item*rptr;
    for (int i=0; i < numrows; i++) {
        rptr = result.Elements[i]; 
        for (int j=0; j < A.numcols; j++, rptr++) 
            for (int k=0; k < numcols; k++) 
                *rptr += (*this)(i,k) * A(k,j);
    }
    return result;
}


//
// multiplies an asymmetric matrix by a symmetric matrix     
// returns an asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::operator*(const SymMatrix<Item>& S)
{
    if (numcols != S.numrows) 
        cerr << "Matrices not compatable.\n";
    AsymMatrix<Item> result(numrows, S.numcols);
    Item*rptr; 
    for (int i=0; i < numrows; i++) {
        rptr = result.Elements[i];
        for (int j=0; j < S.numcols; j++, rptr++) 
            for (int k=0; k < numcols; k++) 
                *rptr += (*this)(i,k) * S(k,j);
    }
    return result;
}


//
// multiplies an asymmetric matrix by a scalar
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::operator*(const Item& scalar)
{
    AsymMatrix<Item> result(numrows, numcols);
    Item*ptr, *rptr;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        rptr = result.Elements[i];
        for (int j=0; j < numcols; j++, ptr++, rptr++) 
            *rptr = (*ptr) * scalar;
    }
    return result;
}


//
// adds two asymmetric matrices, returns and asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::operator+(const AsymMatrix<Item>& A)
{     
    if ((numrows != A.numrows) || (numcols != A.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot add them.\n";
    AsymMatrix<Item> result(numrows, numcols);
    Item*ptr, *rptr, *temp;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        rptr = result.Elements[i];
        temp = A.Elements[i];
        for (int j=0; j < numcols; j++, ptr++, rptr++, temp++)  
            *rptr = *ptr + *temp;
    }
    return result;
}


//
// adds an asymmetric matrix to a symmetric matrix
// return an asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::operator+(const SymMatrix<Item>& S)
{ 
    if ((numrows != S.numrows) || (numcols != S.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot add.\n";
    AsymMatrix<Item> result(numrows, numcols);
    Item*ptr, *rptr;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        rptr = result.Elements[i]; 
        for (int j=0; j < numcols; j++, ptr++, rptr++) 
            *rptr = *ptr + S(i,j);
    }
    return result;
}


//
// substracts two asymetric matrices, returns an asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::operator-(const AsymMatrix<Item>& A)
{
    if ((numrows != A.numrows) || (numcols != A.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot sub them.\n";
    AsymMatrix<Item> result(numrows, numcols);
    Item*ptr, *rptr, *temp;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        rptr = result.Elements[i];
        temp = A.Elements[i];
        for (int j=0; j < numcols; j++, ptr++, rptr++, temp++) 
            *rptr = *ptr - *temp;
    }
    return result;
}


//
// subtracts a symmetric matrix from an asymmetric one
// returns an asymmetric matrix 
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::operator-(const SymMatrix<Item>& S)
{
    if ((numrows != S.numrows) || (numcols != S.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot sub them.\n";
    AsymMatrix<Item> result(numrows, numcols);
    Item*ptr, *rptr;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        rptr = result.Elements[i];
        for (int j=0; j < numcols; j++, ptr++, rptr++) 
            *rptr = *ptr - S(i,j);
    }
    return result;
}


//
// one asymmetric matrix is assigned the values of another asymmetric matrix
// 
template <class Item>
AsymMatrix<Item>&
AsymMatrix<Item>::operator=(const AsymMatrix<Item>& A)
{
    if ((numrows != A.numrows) || (numcols != A.numcols)) 
        cerr << "Matrices are of different sizes.  Assignment not allowed.\n";
    if (this == &A) return *this;
    Item*ptr, *temp;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        temp = A.Elements[i];
        for (int j=0; j < numcols; j++, ptr++, temp++) 
            *ptr = *temp;
    }
    return *this;
}


//
// returns TRUE if the two asymmetric matrices are equal,
// else returns FALSE
//
template <class Item>
Bool  
AsymMatrix<Item>::operator==(const AsymMatrix<Item>& A)
{
    if ((numrows != A.numrows) || (numcols != A.numcols)) 
        return FALSE; 
    Item*ptr, *temp;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        temp = A.Elements[i]; 
        for (int j=0; j < numcols; j++, ptr++, temp++) 
            if ((*ptr) != (*temp))
                return FALSE;
    }
    return TRUE;                     // if everything is equal, return TRUE 
}


//
// returns TRUE if the asymmetric and symmetric matrices are equal, 
// else returns FALSE
//
template <class Item>
Bool
AsymMatrix<Item>::operator==(const SymMatrix<Item>& S)
{
    if ((numrows != S.numrows) || (numcols != S.numcols)) 
        return FALSE;
    Item*ptr;
    for (int i=0; i < numrows; i++) {
        ptr = Elements[i];
        for (int j=0; j < numcols; j++, ptr++)   
            if ( *ptr != S(i,j) )
                return FALSE;
    }
    return TRUE;
}

 
//
// returns TRUE if the two asymmetric matrices are not equal, 
// else returns FALSE
//
template <class Item>
Bool
AsymMatrix<Item>::operator!=(const AsymMatrix<Item>& A) { return !(*this == A);}


//
// returns TRUE if the asymmetric and symmetric matrices are not equal,
// else returns FALSE
//
template <class Item>
Bool
AsymMatrix<Item>::operator!=(const SymMatrix<Item>& S)  { return !(*this == S);}


//
// returns the transpose of the asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
AsymMatrix<Item>::transpose()
{
    AsymMatrix<Item> result(numcols,numrows);
    for (int i=0; i < numrows; i++) 
        for (int j=0; j < numcols; j++) 
            result.Elements[j][i] = Elements[i][j];
    return result;
}

 
/////////////////////////////////////////////////////////////////////////////
//
// symmetric matrix constructors and destructor
//
template <class Item>
SymMatrix<Item>::SymMatrix(int numRows, int numCols) //numCols not used
		: Matrix<Item>(numRows, numRows)
{
    if (numRows <= 0) {
        cerr << "Error in size of Symmetric Matrix<Item>.\n";
        cerr << "dimension must be > 0.\n";
    }  
    size = ((numRows * numRows + numRows)/2);
    LowTriangMat = new Item[size];
#if 1
    for (int i=0; i < size; i++) 
        LowTriangMat[i] = 0;
#endif
}


template <class Item>
SymMatrix<Item>::SymMatrix(const SymMatrix<Item>& S) : 
					Matrix<Item>(numrows, numcols)
{
    size = S.size;
    Item *LowTriangMat = new Item[size];
    Item*ptr = LowTriangMat;
    Item*temp = S.LowTriangMat;
    for (int i=0; i < size; i++, ptr++, temp++) 
        *ptr = *temp;
}


template <class Item>
SymMatrix<Item>::~SymMatrix()   {delete [] LowTriangMat;}

//
// saves a SymMatrix<Item> to a file specified by the user
//
template <class Item>
void
SymMatrix<Item>::save(char * filename)
{
    ofstream curFile;
    curFile.open(filename, ios::out);
    if (!curFile) {
        cerr << "Error saving to file: " << filename << ".\n";
        exit(-1);
    }
    curFile << "S " << numrows << ", " << numcols << "\n";
    Item*ptr = LowTriangMat;
    for (int i=0; i < size; i++, ptr++) {
        curFile << *ptr;
        if (i != size-1) 
            curFile << ",";
    }
    curFile << ";\n";
}

template <class Item>
ostream&
SymMatrix<Item>::display(ostream& os) const
{
    Item*ptr = LowTriangMat;
    for (int i=0; i < size; i++) {
    	for (int j=0; j <=i; j++)
        	os << *ptr++;
	os << endl;
    }
    return os;
}

template <class Item>
ostream& operator<<(ostream& os, const Matrix<Item>& m)
{
	return m.display(os); 
}

//
// loads a SymMatrix<Item> from a file specified by the user
//
template <class Item>
void
SymMatrix<Item>::load(char * filename)
{
    ifstream curFile;
    curFile.open(filename, ios::in);
    if (!curFile) {
        cerr << "Error loading from file: " << filename << ".\n";
        exit(-1);
    }
    char ch;
    int numRows, numCols;
    curFile.get(ch);
    if (ch != 'S') {
        cerr << "Matrices are of different types.  We cannot complete load.\n";
        exit(-1);
    } 
    curFile >> numRows;
    curFile.get(ch);
    curFile >> numCols;
    curFile.get(ch);
    if ((numRows != numrows) || (numCols != numcols)) { 
        cerr << "Matrices are of different sizes. We cannot complete load.\n";
        exit(-1);
    }
    Item*ptr = LowTriangMat;
    for (int i=0; i < size; i++, ptr++) {
        curFile >> *ptr;
        if (i != size-1) 
            curFile.get(ch);
    }
}


//
// returns the element i,j.  Can also set the element i,j
//
template <class Item>
Item&
SymMatrix<Item>::operator() (int row, int col) const
{
    if ((row >= numrows) || (col >= numcols)) 
        cerr << "Out of bounds error.\n";
    if (row >= col) 
        return LowTriangMat[((row * row) + row)/2 + col];
    else 
        return LowTriangMat[((col * col) + col)/2 + row];
}


//
// returns element i,j   
//
//template <class Item>
//Item
//SymMatrix<Item>::operator() (int row, int col) const
//{
//    if (row >= col) 
//        return LowTriangMat[((row*row) + row)/2 + col];
//    else
//        return LowTriangMat[((col*col) + col)/2 + row];
//}


//
// multiplies two symmetric matrices, returns an asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
SymMatrix<Item>::operator*(const SymMatrix<Item>& S)
{
    if (numcols != S.numrows) 
        cerr << "Matrices are not compatable.\n";
    AsymMatrix<Item> result(numrows, S.numcols);
    Item*rptr;
    for (int i=0; i < numrows; i++) { 
        rptr = result.Elements[i];
        for (int j=0; j < S.numcols; j++, rptr++) 
            for (int k=0; k < numcols; k++) 
                *rptr += (*this)(i,k) * S(k,j);
    }
    return result;
}


//
// multiplies a symmetric matrix by an asymmetric matrix
// returns an asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
SymMatrix<Item>::operator*(const AsymMatrix<Item>& A)
{
    if (numcols != A.numrows)  
        cerr << "Matrices are not compatable.\n";
    AsymMatrix<Item> result(numrows, A.numcols);
    Item*rptr;
    for (int i=0; i < numrows; i++) {
        rptr = result.Elements[i];
        for (int j=0; j < A.numcols; j++, rptr++) 
            for (int k=0; k < numcols; k++) 
                *rptr += (*this)(i,k) * A(k,j);
    }
    return result;
}


//
// mutiplies a symmetric matrix by a scalar, returns a symmetric matrix
//
template <class Item>
SymMatrix<Item>
SymMatrix<Item>::operator*(const Item& scalar)
{
    SymMatrix<Item> result(numrows, numcols);
    Item*ptr = LowTriangMat;
    Item*rptr= result.LowTriangMat;
    for (int i=0; i < size; i++, ptr++, rptr++) 
        *rptr = (*ptr) * scalar;
    return result;
}


//
// Adds two symmetric matrices, returns a symmetric matrix
//
template <class Item>
SymMatrix<Item>
SymMatrix<Item>::operator+(const SymMatrix<Item>& S)
{
    if ((numrows != S.numrows) || (numcols != S.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot add them.\n";
    SymMatrix<Item> result(numrows, numcols);
    Item*ptr = LowTriangMat;
    Item*rptr = result.LowTriangMat;
    Item*temp = S.LowTriangMat;
    for (int i=0; i < size; i++, ptr++, rptr++, temp++)
        *rptr = *ptr + *temp;
    return result;
}


//
// adds a symmetric matrix and an asymmetric matrix
// returns an asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
SymMatrix<Item>::operator+(const AsymMatrix<Item>& A)
{
    if ((numrows != A.numrows) || (numcols != A.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot add them.\n";
    AsymMatrix<Item> result(numrows, numcols);
    Item*rptr, *temp;
    for (int i=0; i < numrows; i++) {
        rptr = result.Elements[i];
        temp = A.Elements[i];
        for (int j=0; j < numcols; j++, rptr++, temp++) 
            *rptr = (*this)(i,j) + *temp;
    }
    return result;
}


//
// subtracts two symmetric matrices, returns a symmetric matrix
//
template <class Item>
SymMatrix<Item>
SymMatrix<Item>::operator-(const SymMatrix<Item>& S)
{
    if ((numrows != S.numrows) || (numcols != S.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot sub them.\n";
    SymMatrix<Item> result(numrows, numcols);      
    Item*ptr = LowTriangMat;
    Item*rptr = result.LowTriangMat;
    Item*temp = S.LowTriangMat;
    for (int i=0; i < size; i++, ptr++, rptr++, temp++) 
        *rptr = *ptr - *temp;
    return result;
}


//
// subtracts an asymmetric matrix from a symmetric matrix
// returns an asymmetric matrix
//
template <class Item>
AsymMatrix<Item>
SymMatrix<Item>::operator-(const AsymMatrix<Item>& A)
{
    if ((numrows != A.numrows) || (numcols != A.numcols)) 
        cerr << "Matrices are of different sizes.  We cannot sub them.\n";
    AsymMatrix<Item> result(numrows, numcols);
    Item*rptr, *temp;
    for (int i=0; i < numrows; i++) { 
        rptr = result.Elements[i];
        temp = A.Elements[i];
        for (int j=0; j < numcols; j++, rptr++, temp++) 
            *rptr = (*this)(i,j) - *temp;
    }
    return result;
}


//
// one symmetric matrix is assigned the values of another symmetric matrix
//
template <class Item>
SymMatrix<Item>&
SymMatrix<Item>::operator=(const SymMatrix<Item>& S)
{
    if ((numrows != S.numrows) || (numcols != S.numcols)) 
        cerr << "Matrices are of different sizes.  Assignment not allowed.\n";
    if (this == &S) return *this;
    Item*ptr = LowTriangMat;
    Item*temp = S.LowTriangMat;
    for (int i=0; i < size; i++, ptr++, temp++) 
        *ptr = *temp;
    return *this;
}


//
// returns TRUE if the two symmetric matrices are equal, 
// else returns FALSE
//
template <class Item>
Bool
SymMatrix<Item>::operator==(const SymMatrix<Item>& S)
{
    if ((numrows != S.numrows) || (numcols != S.numcols)) 
        return FALSE;
    Item*ptr = LowTriangMat;
    Item*temp = S.LowTriangMat;
    for (int i=0; i < size; i++, ptr++, temp++) 
        if ((*ptr) != (*temp))
            return FALSE;
    return TRUE;
}


// 
// returns TRUE if the symmetric matrix equals the asymmetric matrix, 
// else returns FALSE
//
template <class Item>
Bool 
SymMatrix<Item>::operator==(const AsymMatrix<Item>& A)
{
    if ((numrows != A.numrows) || (numcols != A.numcols))
        return FALSE;
    for (int i=0; i < numrows; i++)
        for (int j=0; j < numcols; j++)
            if ((*this)(i,j) != A(i,j)) 
                return FALSE;
    return TRUE;
}


//
// returns TRUE if the two symmetric matrices are not equal, 
// else returns FALSE
//
template <class Item>
Bool 
SymMatrix<Item>::operator!=(const SymMatrix<Item>& S)  { return !(*this == S);}


//
// returns TRUE if the symmetric matrix is not equal to the asymmetric matrix,
// else returns FALSE
//
template <class Item>
Bool
SymMatrix<Item>::operator!=(const AsymMatrix<Item>& A) { return !(*this == A); }


//
// returns the transpose of the symmetric matrix
// the transpose of a symmetric matrix is the matrix itself: S == S.transpose()
//
template <class Item>
SymMatrix<Item>
SymMatrix<Item>::transpose()  { return *this;}

#endif