Matrix.cpp

      b[i]=sum/V[i][i];
    }
  return;
}

template<typename T>
void
Matrix<T>::averSymmetric()
  /*!
    Simple function to take an average symmetric matrix
    out of the Matrix
  */
{
  const int minSize=(nx>ny) ? ny : nx;
  for(int i=0;i<minSize;i++)
    for(int j=i+1;j<minSize;j++)
      {
	V[i][j]=(V[i][j]+V[j][i])/2;
	V[j][i]=V[i][j];
      }
  return;
}

template<typename T> 
std::vector<T>
Matrix<T>::Diagonal() const
  /*!
    Returns the diagonal form as a vector
    \returns Diagonal elements
  */
{
  const int Msize=(ny>nx) ? nx : ny;
  std::vector<T> Diag(Msize);
  for(int i=0;i<Msize;i++)
    Diag[i]=V[i][i];
  return Diag;
}

template<typename T> 
T
Matrix<T>::Trace() const
  /*!
    Calculates the trace of the matrix
    \returns Trace of matrix 
  */
{
  const int Msize=(ny>nx) ? nx : ny;
  T Trx=0;
  for(int i=0;i<Msize;i++)
    Trx+=V[i][i];
  return Trx;
}

template<typename T> 
void
Matrix<T>::sortEigen(Matrix<T>& DiagMatrix) 
  /*!
    Sorts the eigenvalues into increasing
    size. Moves the EigenVectors correspondingly
    \param DiagMatrix :: matrix of the EigenValues
  */
{
  if (ny!=nx || nx!=DiagMatrix.nx || nx!=DiagMatrix.ny)
    {
      std::cerr<<"Matrix not Eigen Form"<<std::endl;
    }
  std::vector<int> index;
  std::vector<T> X=DiagMatrix.Diagonal();
  indexSort(X,index);
  Matrix<T> EigenVec(*this);
  for(int Icol=0;Icol<nx;Icol++)
    {
      for(int j=0;j<nx;j++)
	V[j][Icol]=EigenVec[j][index[Icol]];
      DiagMatrix[Icol][Icol]=X[index[Icol]];
    }

  return;
}

template<typename T>
int 
Matrix<T>::Diagonalise(Matrix<T>& EigenVec,Matrix<T>& DiagMatrix) const
  /*!
    Attempt to diagonalise the matrix IF symmetric
    \param EigenVec :: (output) the Eigenvectors matrix 
    \param DiagMatrix  :: the diagonal matrix of eigenvalues 
    \returns :: 1  on success 0 on failure
  */
{
  double theta,tresh,tanAngle,cosAngle,sinAngle;
  if(nx!=ny || nx<1)
    {
      std::cerr<<"Matrix not square"<<std::endl;
      return 0;
    }
  for(int i=0;i<nx;i++)
    for(int j=i+1;j<nx;j++)
      if (fabs(V[i][j]-V[j][i])>1e-6)
        {
	  std::cerr<<"Matrix not symmetric"<<std::endl;
	  std::cerr<< (*this);
	  return 0;
	}
    
  Matrix<T> A(*this);
  //Make V an identity matrix
  EigenVec.setMem(nx,nx);
  EigenVec.identityMatrix();
  DiagMatrix.setMem(nx,nx);
  DiagMatrix.zeroMatrix();

  std::vector<double> Diag(nx);
  std::vector<double> B(nx);
  std::vector<double> ZeroComp(nx);
  //set b and d to the diagonal elements o A
  for(int i=0;i<nx;i++)
    {
      Diag[i]=B[i]= A.V[i][i];
      ZeroComp[i]=0;
    }

  int iteration=0;
  double sm; 
  for(int i=0;i<100;i++)        //max 50 iterations
    {
      sm=0.0;           // sum of off-diagonal terms
      for(int ip=0;ip<nx-1;ip++)
	for(int iq=ip+1;iq<nx;iq++)
	  sm+=fabs(A.V[ip][iq]);

      if (sm==0.0)           // Nothing to do return...
	{
	  // Make OUTPUT -- D + A
	  // sort Output::
	  std::vector<int> index;
//	  indexSort(Diag,index);
//	  for(int ii=0;ii<nx;ii++)
//	    {
//	      DiagMatrix.V[ii][ii]=static_cast<T>(Diag[index[ii]]);
//	      EigenVector

	  for(int ix=0;ix<nx;ix++)
	    DiagMatrix.V[ix][ix]=static_cast<T>(Diag[ix]);
	  return 1;
	}

      // Threshold large for first 5 sweeps
      tresh= (i<6) ? 0.2*sm/(nx*nx) : 0.0;

      for(int ip=0;ip<nx-1;ip++)
	{
	  for(int iq=ip+1;iq<nx;iq++)
	    {
	      double g=100.0*fabs(A.V[ip][iq]);
	      // After 4 sweeps skip if off diagonal small
	      if (i>6 && 
		  static_cast<float>(fabs(Diag[ip]+g))==static_cast<float>(fabs(Diag[ip])) &&
		  static_cast<float>(fabs(Diag[iq]+g))==static_cast<float>(fabs(Diag[iq])))
		A.V[ip][iq]=0;

	      else if (fabs(A.V[ip][iq])>tresh)
		{
		  double h=Diag[iq]-Diag[ip];
		  if (static_cast<float>((fabs(h)+g)) == static_cast<float>(fabs(h)))
		    tanAngle=A.V[ip][iq]/h;      // tanAngle=1/(2theta)
		  else
		    {
		      theta=0.5*h/A.V[ip][iq];
		      tanAngle=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
		      if (theta<0.0)
			tanAngle = -tanAngle;
		    }
		  cosAngle=1.0/sqrt(1+tanAngle*tanAngle);
		  sinAngle=tanAngle*cosAngle;
		  double tau=sinAngle/(1.0+cosAngle);
		  h=tanAngle*A.V[ip][iq];
		  ZeroComp[ip] -= h;
		  ZeroComp[iq] += h;
		  Diag[ip] -= h;
		  Diag[iq] += h;
		  A.V[ip][iq]=0;
		  // Rotations 0<j<p
		  for(int j=0;j<ip;j++)
		    A.rotate(tau,sinAngle,j,ip,j,iq);
		  for(int j=ip+1;j<iq;j++)
		    A.rotate(tau,sinAngle,ip,j,j,iq);
		  for(int j=iq+1;j<nx;j++)
		    A.rotate(tau,sinAngle,ip,j,iq,j);
		  for(int j=0;j<nx;j++)
		    EigenVec.rotate(tau,sinAngle,j,ip,j,iq);
		  iteration++;
		}
	    }
	}   
      for(int j=0;j<nx;j++)
	{
	  B[j]+=ZeroComp[j];
	  Diag[j]=B[j];
	  ZeroComp[j]=0.0;
	}
	    
    }
  std::cerr<<"Error :: Iterations are a problem"<<std::endl;
  return 0;
}

template<typename T>
void
Matrix<T>::print() const
  /*! 
    Simple print out routine 
   */
{
  write(std::cout,10);
  return;
}

template<typename T>
void
Matrix<T>::write(std::ostream& Fh,const int blockCnt) const
  /*!
    Write out function for blocks of 10 Columns 
    \param Fh :: file stream for output
    \param blockCnt :: number of columns per line (0 == full)
  */
{

  std::ios::fmtflags oldFlags=Fh.flags();
  Fh.setf(std::ios::floatfield,std::ios::scientific);
  const int blockNumber((blockCnt>0) ? blockCnt : ny);  
  int BCnt(0);

  do
    {
      const int ACnt=BCnt;
      BCnt+=blockNumber;
      if (BCnt>ny)
	BCnt=ny;

      if (ACnt)
	Fh<<" ----- "<<ACnt<<" "<<BCnt<<" ------ "<<std::endl;
      for(int i=0;i<nx;i++)
        {
	  for(int j=ACnt;j<BCnt;j++)
	    Fh<<std::setw(10)<<V[i][j]<<"  ";
	  Fh<<std::endl;
	}
    } while(BCnt<ny);

  Fh.flags(oldFlags);
  return;
}

template<typename T>
std::string
Matrix<T>::str() const
  /*!
    Convert the matrix into a simple linear string expression 
    \returns String value of output
  */
{
  std::ostringstream cx;
  for(int i=0;i<nx;i++)
    for(int j=0;j<ny;j++)
      {
	cx<<std::setprecision(6)<<V[i][j]<<" ";
      }
  return cx.str();
}

template class Matrix<double>;
template class Matrix<int>;

}   // NAMESPACE Geometry

}   // NAMESPACE MAntid