Commit 02e6ac9c authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

Matrix indentation

parent 358aff20

drivers/closuresTest.cpp

+25 −1
Original line number Diff line number Diff line 
#include "Matrix.h"

#include "Vector.h"



int main(int argc, char** argv)

void testVector()

{

	int n = 4;

	std::vector<double> v1(n,1.0);
@@ -28,3 +29,26 @@ int main(int argc, char** argv) 
//	std::cout<<v3;

}



void testMatrix()

{

	SizeType n = 4;

	PsimagLite::Matrix<double> m1(n,n);

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

		m1(i,i) = i;

		SizeType iPlus1 = i + 1;

		if (iPlus1 >= n) iPlus1 = 0;

		m1(i,iPlus1) = m1(iPlus1,i) = 0.1;

	}



	PsimagLite::Matrix<double> m2;

	// multiplication by scalar

	m2 = m1*1.2;

	std::cout<<m2;

}



int main(int, char**)

{

	//testVector();

	testMatrix();

}

src/Matrix.h

+180 −176
Original line number Diff line number Diff line
@@ -339,7 +339,32 @@ void symbolicPrint(std::ostream& os,const Matrix<T>& A) 
	os<<"\n";

}



template <class T>

template<typename T>

void printNonZero(const Matrix<T>& m,std::ostream& os)

{

	os<<"#size="<<m.n_row()<<"x"<<m.n_col()<<"\n";

	for (SizeType i=0;i<m.n_row();i++) {

		SizeType nonzero = 0;

		for (SizeType j=0;j<m.n_col();j++) {

			const T& val = m(i,j);

			if (val!=static_cast<T>(0)) {

				os<<val<<" ";

				nonzero++;

			}

		}

		if (nonzero>0) os<<"\n";

	}

	os<<"#diagonal non-zero\n";

	for (SizeType i=0;i<m.n_row();i++) {

		const T& val = m(i,i);

		if (val!=static_cast<T>(0)) {

			os<<val<<" ";

		}

	}

	os<<"\n";

}



template<typename T>

std::istream& operator>>(std::istream& is, Matrix<T>& A)

{

	SizeType nrow=0,ncol=0;
@@ -360,6 +385,65 @@ std::istream& operator >> (std::istream& is, Matrix<T>& A) 
	return is;

}



template<typename T>

bool isHermitian(Matrix<T> const &A,bool verbose=false)

{

	SizeType n=A.n_row();

	double eps=1e-6;

	if (n!=A.n_col())

		throw RuntimeError("isHermitian called on a non-square matrix.\n");

	for (SizeType i=0;i<n;i++) for (SizeType j=0;j<n;j++) {

		if (std::norm(A(i,j)-std::conj(A(j,i)))>eps) {

			if (verbose) {

				std::cerr<<"A("<<i<<","<<j<<")="<<A(i,j);

				std::cerr<<" A("<<j<<","<<i<<")="<<A(j,i)<<"\n";

			}

			return false;

		}

	}

	return true;

}



template<typename T>

bool isTheIdentity(Matrix<T> const &a)

{



	for (SizeType i=0;i<a.n_row();i++) {

		for (SizeType j=0;j<a.n_col();j++) {

			if (i!=j && std::norm(a(i,j))>1e-6)  {

				std::cerr<<"a("<<i<<","<<j<<")="<<a(i,j)<<"\n";

				return false;

			}

		}

	}



	for (SizeType i=0;i<a.n_row();i++)

		if (std::norm(a(i,i)-static_cast<T>(1.0))>1e-6) return false;



	return true;

}



template<typename T>

bool isZero(Matrix<std::complex<T> > const &a)

{



	for (SizeType i=0;i<a.n_row();i++)

		for (SizeType j=0;j<a.n_col();j++)

			if (norm(a(i,j))>0) return false;

	return true;

}



template<typename T>

bool isZero(const Matrix<T>& m)

{

	for (SizeType i=0;i<m.n_row();i++)

		for (SizeType j=0;j<m.n_col();j++)

			if (fabs(m(i,j))>0) return false;

	return true;

}



// closures start



template<typename T>

Matrix<T> operator+(const Matrix<T>& a,const Matrix<T>& b)

{
@@ -441,6 +525,82 @@ Matrix<T2> operator*(const Matrix<T2>& a,const T1& val) 
	return val*a;

}



template<typename T>

T norm2(const Matrix<T>& m)

{

	T sum = 0;

	for (SizeType i=0;i<m.n_row();i++)

		for (SizeType j=0;j<m.n_col();j++)

			sum += m(i,j) * m(i,j);

	return sum;

}



template<typename T>

T norm2(const Matrix<std::complex<T> >& m)

{

	T sum = 0;

	for (SizeType i=0;i<m.n_row();i++)

		for (SizeType j=0;j<m.n_col();j++)

			sum += std::real(std::conj(m(i,j)) * m(i,j));

	return sum;

}



template<typename T>

Matrix<T> multiplyTransposeConjugate(

        const Matrix<T>& O1,

        const Matrix<T>& O2,char modifier='C')

{

	SizeType n=O1.n_row();

	Matrix<T> ret(n,n);

	if (modifier=='C') {

		for (SizeType s=0;s<n;s++)

			for (SizeType t=0;t<n;t++)

				for (SizeType w=0;w<n;w++)

					ret(s,t) += std::conj(O1(w,s))*O2(w,t);

	} else {

		for (SizeType s=0;s<n;s++)

			for (SizeType t=0;t<n;t++)

				for (SizeType w=0;w<n;w++)

					ret(s,t) += O1(w,s)*O2(w,t);

	}

	return ret;

}



template<class T>

void transposeConjugate(Matrix<T>& m2,const Matrix<T>& m)

{

	m2.resize(m.n_col(),m.n_row());

	for (SizeType i=0;i<m2.n_row();i++)

		for (SizeType j=0;j<m2.n_col();j++)

			m2(i,j)=std::conj(m(j,i));



}



template<typename T>

Matrix<T> transposeConjugate(const Matrix<T>& A)

{

	Matrix<T> ret(A.n_col(),A.n_row());

	for (SizeType i=0;i<A.n_col();i++)

		for (SizeType j=0;j<A.n_row();j++)

			ret(i,j)=std::conj(A(j,i));

	return ret;

}



template<typename T>

void outerProduct(Matrix<T>& A,const Matrix<T>& B,const Matrix<T>& C)

{

	SizeType ni = B.n_row();

	SizeType nj = B.n_col();



	A.resize(B.n_row()*C.n_row(),B.n_col()*C.n_col());



	for (SizeType i1 = 0; i1 < B.n_row(); ++i1)

		for (SizeType j1 = 0; j1 < B.n_col(); ++j1)

			for (SizeType i2 = 0; i2 < C.n_row(); ++i2)

				for (SizeType j2 = 0; j2 < C.n_col(); ++j2)

					A(i1+i2*ni,j1+j2*nj) = B(i1,j1) * C(i2,j2);

}



template<typename T>

void exp(Matrix<T>& m)

{
@@ -464,6 +624,8 @@ void exp(Matrix<T>& m) 


}



// closures end



template<typename VectorLikeType>

void svd(char jobz,Matrix<double> &a,VectorLikeType& s,Matrix<double>& vt)

{
@@ -537,164 +699,6 @@ void svd(char jobz,Matrix<double> &a,VectorLikeType& s,Matrix<double>& vt) 
#endif

}



template<typename T>

bool isHermitian(Matrix<T> const &A,bool verbose=false)

{

	SizeType n=A.n_row();

	double eps=1e-6;

	if (n!=A.n_col())

		throw RuntimeError("isHermitian called on a non-square matrix.\n");

	for (SizeType i=0;i<n;i++) for (SizeType j=0;j<n;j++) {

		if (std::norm(A(i,j)-std::conj(A(j,i)))>eps) {

			if (verbose) {

				std::cerr<<"A("<<i<<","<<j<<")="<<A(i,j);

				std::cerr<<" A("<<j<<","<<i<<")="<<A(j,i)<<"\n";

			}

			return false;

		}

	}

	return true;

}



template<typename T>

void printNonZero(const Matrix<T>& m,std::ostream& os)

{

	os<<"#size="<<m.n_row()<<"x"<<m.n_col()<<"\n";

	for (SizeType i=0;i<m.n_row();i++) {

		SizeType nonzero = 0;

		for (SizeType j=0;j<m.n_col();j++) {

			const T& val = m(i,j);

			if (val!=static_cast<T>(0)) {

				os<<val<<" ";

				nonzero++;

			}

		}

		if (nonzero>0) os<<"\n";

	}

	os<<"#diagonal non-zero\n";

	for (SizeType i=0;i<m.n_row();i++) {

		const T& val = m(i,i);

		if (val!=static_cast<T>(0)) {

			os<<val<<" ";

		}

	}

	os<<"\n";

}



template<typename T>

bool isTheIdentity(Matrix<T> const &a)

{



	for (SizeType i=0;i<a.n_row();i++) {

		for (SizeType j=0;j<a.n_col();j++) {

			if (i!=j && std::norm(a(i,j))>1e-6)  {

				std::cerr<<"a("<<i<<","<<j<<")="<<a(i,j)<<"\n";

				return false;

			}

		}

	}



	for (SizeType i=0;i<a.n_row();i++)

		if (std::norm(a(i,i)-static_cast<T>(1.0))>1e-6) return false;



	return true;

}



template<typename T>

bool isZero(Matrix<std::complex<T> > const &a)

{



	for (SizeType i=0;i<a.n_row();i++)

		for (SizeType j=0;j<a.n_col();j++)

			if (norm(a(i,j))>0) return false;

	return true;

}



template<typename T>

bool isZero(const Matrix<T>& m)

{

	for (SizeType i=0;i<m.n_row();i++)

		for (SizeType j=0;j<m.n_col();j++)

			if (fabs(m(i,j))>0) return false;

	return true;

}



template<typename T>

T norm2(const Matrix<T>& m)

{

	T sum = 0;

	for (SizeType i=0;i<m.n_row();i++)

		for (SizeType j=0;j<m.n_col();j++)

			sum += m(i,j) * m(i,j);

	return sum;

}



template<typename T>

T norm2(const Matrix<std::complex<T> >& m)

{

	T sum = 0;

	for (SizeType i=0;i<m.n_row();i++)

		for (SizeType j=0;j<m.n_col();j++)

			sum += std::real(std::conj(m(i,j)) * m(i,j));

	return sum;

}



template<typename T>

Matrix<T> multiplyTransposeConjugate(

        const Matrix<T>& O1,

        const Matrix<T>& O2,char modifier='C')

{

	SizeType n=O1.n_row();

	Matrix<T> ret(n,n);

	if (modifier=='C') {

		for (SizeType s=0;s<n;s++)

			for (SizeType t=0;t<n;t++)

				for (SizeType w=0;w<n;w++)

					ret(s,t) += std::conj(O1(w,s))*O2(w,t);

	} else {

		for (SizeType s=0;s<n;s++)

			for (SizeType t=0;t<n;t++)

				for (SizeType w=0;w<n;w++)

					ret(s,t) += O1(w,s)*O2(w,t);

	}

	return ret;

}



template<class T>

void transposeConjugate(Matrix<T>& m2,const Matrix<T>& m)

{

	m2.resize(m.n_col(),m.n_row());

	for (SizeType i=0;i<m2.n_row();i++)

		for (SizeType j=0;j<m2.n_col();j++)

			m2(i,j)=std::conj(m(j,i));



}



template<typename T>

Matrix<T> transposeConjugate(const Matrix<T>& A)

{

	Matrix<T> ret(A.n_col(),A.n_row());

	for (SizeType i=0;i<A.n_col();i++)

		for (SizeType j=0;j<A.n_row();j++)

			ret(i,j)=std::conj(A(j,i));

	return ret;

}



template<typename T>

void outerProduct(Matrix<T>& A,const Matrix<T>& B,const Matrix<T>& C)

{

	SizeType ni = B.n_row();

	SizeType nj = B.n_col();



	A.resize(B.n_row()*C.n_row(),B.n_col()*C.n_col());



	for (SizeType i1 = 0; i1 < B.n_row(); ++i1)

		for (SizeType j1 = 0; j1 < B.n_col(); ++j1)

			for (SizeType i2 = 0; i2 < C.n_row(); ++i2)

				for (SizeType j2 = 0; j2 < C.n_col(); ++j2)

					A(i1+i2*ni,j1+j2*nj) = B(i1,j1) * C(i2,j2);

}



#ifdef USE_MPI

namespace MPI {