Commit 3e49e15d authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

CrsMAtrix: tests and improvements

parent 74781f6a

drivers/testCRSMatrix.cpp

+44 −7
Original line number Diff line number Diff line 
#include "SampleCRSMatrix.h"

//#include "CrsMatrix.h"

#include "CrsMatrix.h"

#include <cstdlib>

#include <fstream>

using namespace PsimagLite;
@@ -14,6 +15,14 @@ std::ostream& operator<<(std::ostream& os,const typename Vector<T>::Type& v) 
	return os;

}



void fillRandomVector(Vector<SizeType>::Type& x, SizeType maxValue)

{

	unsigned int long seed = 7334211;

	srand48(seed);

	for (SizeType i=0;i<x.size();i++)

		x[i] = drand48()*maxValue;

}



void fillRandomVector(Vector<RealType>::Type& x,RealType maxValue)

{

	unsigned int long seed = 7334211;
@@ -23,9 +32,9 @@ void fillRandomVector(Vector<RealType>::Type& x,RealType maxValue) 
}



template<typename T>

void testMultiply(const SampleCRSMatrix<T>& m,RealType maxValue)

void testMultiply(const CrsMatrix<T>& m,RealType maxValue)

{

	typename Vector<RealType>::Type x(m.rank(),0.0),y(m.rank());

	typename Vector<RealType>::Type x(m.rows(),0.0),y(m.rows());

	fillRandomVector(y,maxValue);

	std::cout<<"initial vector:\n";

	std::cout<<y;
@@ -34,6 +43,30 @@ void testMultiply(const SampleCRSMatrix<T>& m,RealType maxValue) 
	std::cout<<x;

}



template<typename T>

CrsMatrix<T> createRandomCrs(SizeType rank, SizeType seed, SizeType nonZeros, T maxValue)

{

	typename Vector<SizeType>::Type rows;

	typename Vector<SizeType>::Type cols;

	typename Vector<T>::Type vals;



	srand48(seed);

	for (SizeType i = 0; i < nonZeros; ++i) {

		// pick a row

		const SizeType row = SizeType(drand48()*rank);

		// and a column

		const SizeType col = SizeType(drand48()*rank);

		// and a value

		const T val = drand48()*maxValue;

		rows.push_back(row);

		cols.push_back(col);

		vals.push_back(val);

	}



	// fill the matrix with this data:

	return CrsMatrix<T>(rank, rows, cols, vals);

}



int main(int argc,char *argv[])

{
@@ -43,18 +76,22 @@ int main(int argc,char *argv[]) 
		RealType ratio = std::atof(argv[2]);

		SizeType nonZeros = SizeType(ratio * rank *rank);

		RealType maxValue = 10.0;

		SampleCRSMatrix<RealType> m(rank,seed,nonZeros,maxValue);

		CrsMatrix<RealType> m = createRandomCrs(rank, seed, nonZeros, maxValue);

		std::cout<<m;



		testMultiply(m,maxValue);

	} else if (argc==2) {

		std::ifstream fin(argv[1]);

		SampleCRSMatrix<RealType> m(fin);

		Matrix<RealType> mdense(fin);

		fin.close();

		std::cout<<m;

		std::cout<<mdense;



		CrsMatrix<RealType> m(mdense);

		RealType maxValue = 10.0;

		testMultiply(m, maxValue);

		std::cout<<m;

		std::cout<<"----------\n";

		std::cout<<m.toDense();

	} else {

		throw RuntimeError("Wrong number of arguments\n");

	}

src/CrsMatrix.h

+124 −93
Original line number Diff line number Diff line
@@ -87,6 +87,8 @@ DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS. 
#include "loki/TypeTraits.h"

#include "Mpi.h"

#include "Io/IoSerializerStub.h"

#include <fstream>

#include "Sort.h"



namespace PsimagLite {
@@ -189,6 +191,35 @@ public: 
		setRow(a.rows(),counter);

	}



	CrsMatrix(SizeType rank, // square matrix ONLY for now

	          const Vector<SizeType>::Type& rows2,

	          const Vector<SizeType>::Type& cols,

	          const typename Vector<T>::Type& vals)

	    : rowptr_(rank + 1), nrow_(rank), ncol_(rank)

	{

		Sort<Vector<SizeType>::Type > s;

		Vector<SizeType>::Type iperm(rows2.size());

		Vector<SizeType>::Type rows = rows2;

		s.sort(rows, iperm);

		SizeType counter = 0;

		SizeType prevRow = rows[0]+1;

		for (SizeType i=0;i<rows.size();i++) {

			SizeType row = rows[i];

			if (prevRow!=row) {

				// add new row

				rowptr_[row] = counter++;

				prevRow = row;

			}



			colind_.push_back(cols[iperm[i]]);

			values_.push_back(vals[iperm[i]]);

		}



		SizeType lastNonZeroRow = rows[rows.size()-1];

		for (SizeType i = lastNonZeroRow + 1; i <= rank; ++i)

			rowptr_[i] = counter;

	}



	// start closure ctors



	CrsMatrix(const std::ClosureOperator<CrsMatrix,
@@ -321,7 +352,6 @@ public: 
		values_ *= x;

	}





	bool operator==(const CrsMatrix<T>& op) const

	{

		return (nrow_ == op.nrow_ &&
@@ -437,6 +467,7 @@ public: 
			colind_[i]=i;

			rowptr_[i]=i;

		}



		rowptr_[row]=row;

	}

src/Matrix.h

+14 −0
Original line number Diff line number Diff line
@@ -29,6 +29,7 @@ Please see full open source license included in file LICENSE. 
#include "TypeToString.h"

#include "Mpi.h"

#include "Io/IoSerializerStub.h"

#include <fstream>



namespace PsimagLite {
@@ -91,6 +92,19 @@ public: 
				data_[i+j*nrow_] = m(i,j);

	}



	Matrix(std::ifstream& io)

	{

		io>>nrow_;

		io>>ncol_;

		if (nrow_ <= 0 || ncol_ <= 0)

			throw RuntimeError("Matrix::nrows or ncols negative: cannot construct\n");



		data_.resize(nrow_*ncol_);

		for (SizeType i = 0; i < nrow_; ++i)

			for (SizeType j = 0; j < ncol_; ++j)

				io >> data_[i + j*nrow_];

	}



	// ctor closures

	Matrix(const std::ClosureOperator<Matrix,Matrix,std::ClosureOperations::OP_MULT>& c)

	{