nrutil.h

Schrodinger Equation 数值计算中的方程

#pragma once
#include<string>
#include<cmath>
#include<complex>
//#include<iostream>
#include "assert.h"
//using namespace std;

typedef double DP;

    template<class T>
		inline const T SQR(const T &a){return a*a;}
	
	template<class T>
		inline const T MAX(const T &a,const T &b)
	{return b>a?(b):(a);}
	inline float MAX(const double &a,const float &b)
	{return b>a?float(b):(a);}
	inline float MAX(const float &a,const double &b)
	{return b>a?float(b):(a);}

	template<class T>
		inline const T MIN(const T &a,const T &b)
	{return b<a ? (b):(a);}
	inline float MIN(const double &a,const float &b)
	{return b<a ? (b):float(a);}
	inline float MIN(const float &a,const double &b)
	{return b<a ? float(b):(a);}

	template<class T>
		inline const T SIGN(const T &a, const T &b)
	{
		return b>=0 ? (a>=0 ? a:-a):(a>=0 ? -a:a);
	}
	inline float SIGN(const float &a,const double &b)
	{
		return b>=0 ? (a>=0 ? a:-a):(a>=0 ?-a:a);
	}
    inline float SIGN(const double &a,const float &b)
	{
		return b>=0 ? (a>=0 ?a:-a):(a>=0 ? -a:a);
	}

	template<class T>
		inline void SWAP(T &a,T &b)
	{
		T dum=a;a=b;b=dum;
	}

    namespace NR
	{
		//inline void nrerror(const string error_text)
		//{
			//std::cerr<<"Numerical Recipes run-time error..."<<endl;
			//std::cerr<<error_text<<endl;
			//std::cerr<<"...now exiting to system..."<<endl;
		//	exit(1);
		//}
	}
//////////////////////////////////////////////////////////////
template<class T>
class NRVec  
{
private:
	int nn;//size of array. upper index is nn-1
	T *v;

public:
	NRVec();
	explicit NRVec(int n);
	NRVec(const T &a, int n);
	NRVec(const T *a, int n);
	NRVec(const NRVec &rhs);
	NRVec &operator=(const NRVec &rhs);
	NRVec &operator=(const T &a);
    inline T & operator [](const int i)
	{
	return v[i];
	}
	inline const T & operator[](const int i) const
	{
	return v[i];
	}


	inline int size() const
	{
	return nn;
	}

	~NRVec();
};
//////////////////////////////////////////////////////////////
template<class T>
NRVec<T>::NRVec():nn(0),v(0){}

template<class T>
NRVec<T>::NRVec(int n):nn(n),v(new T[n]){}

template<class T>
NRVec<T>::NRVec(const T& a, int n):nn(n),v(new T[n])
{
	for(int i=0;i<n;i++)
		v[i]=a;
}


template<class T>
NRVec<T>::NRVec(const T *a, int n):nn(n),v(new T[n])
{
	for(int i=0;i<nn;i++)
		v[i]=*a++;
}

template<class T>
NRVec<T>::NRVec(const NRVec<T> &rhs):nn(rhs.nn),v(new T[nn])
{
	for(int i=0;i<nn;i++)
		v[i]=rhs[i];
}

template<class T>
NRVec<T> & NRVec<T>::operator =(const NRVec<T> &rhs)
{
	if(this !=&rhs)
	{
		if(nn!=rhs.nn){
			if(v!=0)delete [](v);
			nn=rhs.nn;
			v=new T[nn];
		}
		for(int i=0;i<nn;i++)
			v[i]=rhs[i];
	}
	return *this;
}

template <class T>
NRVec<T> & NRVec<T>::operator =(const T &a)
{
	for(int i=0;i<nn;i++)
		v[i]=a;
	return *this;
}

template <class T>
NRVec<T>::~NRVec()
{
	if(v!=0)
		delete [] (v);
}

//////////////////////////////////////////////////////
template<class T>
class NRMat  
{
private:
	int nn;
	int mm;
	T **v;
public:
	NRMat();
	NRMat(int n, int m);
	NRMat(const T &a, int n, int m);
	NRMat(const T *a, int n, int m);
	NRMat(const NRMat &rhs);
	NRMat & operator=(const NRMat &rhs);
	NRMat & operator=(const T &a);
    inline T* operator [](const int i)
		{return v[i];}
	inline const T* operator [](const int i) const
		{return v[i];}
	inline int nrows() const
	{return nn;}
	inline int ncols() const
		{return mm;}
	~NRMat();

};
//////////////////////////////////////////////////////
template<class T>
NRMat<T>::NRMat():nn(0),mm(0),v(0){}

template <class T>
NRMat<T>::NRMat(int n,int m):nn(n),mm(m),v(new T*[n])
//n--row;m--bollow;
{
	v[0]=new T[m*n];//the start point of pointer;
	
	for(int i=1;i<n;i++)
		v[i]=v[i-1]+m;
	//v[0],which is the 0-th row,  has m elements.
	//v[1],which is the 1-th row,  also has m elements consisting of
	// a row.
}

template <class T>
NRMat<T>::NRMat(const T &a, int n, int m):nn(n),mm(m),v(new T*[n])
{
	int i,j;
	v[0]=new T[m*n];
	for(i=1;i<n;i++)
		v[i]=v[i-1]+m;
	for(i=0;i<n;i++)
		for(j=0;j<m;j++)
			v[i][j]=a;
}

template <class T>
NRMat<T>::NRMat(const T*a, int n, int m):nn(n),mm(m),v(new T*[n])
{
	int i,j;
	v[0]=new T[m*n];
	for(i=1;i<n;i++)
		v[i]=v[i-1]+m;
	for(i=0;i<n;i++)
		for(j=0;j<m;j++)
			v[i][j]=*a++;
}

template<class T>
NRMat<T>::NRMat(const NRMat &rhs):nn(rhs.nn),mm(rhs.mm),v(new T*[nn])
{
	int i,j;
	v[0]=new T[mm*nn];
	for(i=1;i<nn;i++)
		v[i]=v[i-1]+mm;
	for(i=0;i<nn;i++)
		for(j=0;j<mm;j++)
			v[i][j]=rhs[i][j];
}

template <class T>
NRMat<T> & NRMat<T>::operator =(const NRMat<T> &rhs)
{
	if(this !=&rhs){
		int i,j;
		if(nn!=rhs.nn||mm!=rhs.mm){
			if(v!=0){
				delete[](v[0]);
				delete[] (v);
			}
			nn=rhs.nn;
			mm=rhs.mm;
			v=new T*[nn];
			v[0]=new T[mm*nn];
		}
		for(i=1;i<nn;i++)
			v[i]=v[i-1]+mm;
		for(i=0;i<nn;i++)
			for(j=0;j<mm;j++)
				v[i][j]=rhs[i][j];
	}
	return *this;
}

template <class T>
NRMat<T> & NRMat<T>::operator =(const T &a)
{
	for(int i=0;i<nn;i++)
		for(int j=0;j<mm;j++)
			v[i][j]=a;
		return *this;
}


template <class T>
NRMat<T>::~NRMat()
{
	if(v!=0){
		delete [](v[0]);
		delete [](v);
	}
}