Matrix.cpp

  }
  // remove old memory
  const size_t tx=nx;
  const size_t ty=ny;
  deleteMem();  // resets nx,ny
  // replace memory
  V=Vt;
  nx=ty;
  ny=tx;

  return *this;
}

template<>
int
Matrix<int>::GaussJordan(Kernel::Matrix<int>&)
  /**
    Not valid for Integer
    @return zero
  */
{
  return 0;
}

template<typename T>
int
Matrix<T>::GaussJordan(Matrix<T>& B)
  /**
    Invert this Matrix and solve the 
    form such that if A.x=B then  solve to generate x.  
    This requires that B is B[A.nx][Any]
    The result is placed back in B
    @return the calculation result
   */
{
  // check for input errors
  if (nx!=ny || B.nx!=nx)
    return -1;
  
  // pivoted rows 
  std::vector<int> pivoted(nx);
  fill(pivoted.begin(),pivoted.end(),0);

  std::vector<int> indxcol(nx);   // Column index
  std::vector<int> indxrow(nx);   // row index

  size_t irow(0),icol(0);
  for(size_t i=0;i<nx;i++)
  {
    // Get Biggest non-pivoted item
    double bigItem=0.0;         // get point to pivot over
    for(size_t j=0;j<nx;j++)
        {
            if (pivoted[j]!= 1)  //check only non-pivots
            {
              for(size_t k=0;k<nx;k++)
        {
                if (!pivoted[k])
                      {
                        if (fabs(V[j][k]) >=bigItem)
                        {
                                bigItem=fabs(V[j][k]);
                                irow=j;
                                icol=k;
                        }
                  }
        }
      }
            else if (pivoted[j]>1)
      {
              throw std::runtime_error("Error doing G-J elem on a singular matrix");
      }
    }
    pivoted[icol]++;
    // Swap in row/col to make a diagonal item 
    if (irow != icol)
        {
            swapRows(irow,icol);
            B.swapRows(irow,icol);
          }
    indxrow[i]=static_cast<int>(irow);
    indxcol[i]=static_cast<int>(icol);
    
    if (V[icol][icol] == 0.0)
          {
            std::cerr<<"Error doing G-J elem on a singular matrix"<<std::endl;
            return 1;
    }
    const T pivDiv= T(1.0)/V[icol][icol];
    V[icol][icol]=1;
    for(size_t l=0;l<nx;l++) 
    {
      V[icol][l] *= pivDiv;
    }
    for(size_t l=0;l<B.ny;l++)
    {
      B.V[icol][l] *=pivDiv;
    }
    
    for(size_t ll=0;ll<nx;ll++)
    {
      if (ll!=icol)
      {
        const T div_num=V[ll][icol];
        V[ll][icol]=0.0;
        for(size_t l=0;l<nx;l++) 
        {
          V[ll][l] -= V[icol][l]*div_num;
        }
        for(size_t l=0;l<B.ny;l++)
        {
          B.V[ll][l] -= B.V[icol][l]*div_num;
        }
      }
    }
  }

  // Un-roll interchanges
  if( nx > 0 )
  {
    for(int l=static_cast<int>(nx)-1;l>=0;l--)
    {
      if (indxrow[l] !=indxcol[l])
      {
        swapCols(indxrow[l],indxcol[l]);
      }
    }
  }
  return 0;  
}

template<typename T>
std::vector<T>
Matrix<T>::Faddeev(Matrix<T>& InvOut)
  /**
    Return the polynominal for the matrix
    and the inverse. 
    The polynomial is such that
    \f[
      det(sI-A)=s^n+a_{n-1}s^{n-1} \dots +a_0
    \f]
    @param InvOut ::: output 
    @return Matrix self Polynomial (low->high coefficient order)
  */
{
  if (nx!=ny)
    throw Kernel::Exception::MisMatch<size_t>(nx,ny,"Matrix::Faddev(Matrix)");

  Matrix<T>& A(*this);
  Matrix<T> B(A);
  Matrix<T> Ident(nx,ny);
  Ident.identityMatrix();

  T tVal=B.Trace();                     // Trace of the matrix
  std::vector<T> Poly;
  Poly.push_back(1);
  Poly.push_back(tVal);

  for(size_t i=0;i<nx-2;i++)   // skip first (just copy) and last (to keep B-1)
    {
      B=A*B - Ident*tVal;
      tVal=B.Trace();
      Poly.push_back(tVal/static_cast<T>(i+1));
    }
  // Last on need to keep B;
  InvOut=B;
  B=A*B - Ident*tVal;
  tVal=B.Trace();
  Poly.push_back(tVal/static_cast<T>(nx));
  
  InvOut-= Ident* (-Poly[nx-1]);
  InvOut/= Poly.back();
  return Poly;
}

template<typename T>
T
Matrix<T>::Invert()
  /**
    If the Matrix is square then invert the matrix
    using LU decomposition
    @return Determinant (0 if the matrix is singular)
  */
{
  if (nx!=ny && nx<1)
    return 0;
  
  int* indx=new int[nx];   // Set in lubcmp

  double* col=new double[nx];
  int d;
  Matrix<T> Lcomp(*this);
  Lcomp.lubcmp(indx,d);
        
  double det=static_cast<double>(d);
  for(size_t j=0;j<nx;j++)
    det *= Lcomp.V[j][j];
  
  for(size_t j=0;j<nx;j++)
    {
      for(size_t i=0;i<nx;i++)
        col[i]=0.0;
      col[j]=1.0;
      Lcomp.lubksb(indx,col);
      for(size_t i=0;i<nx;i++)
        V[i][j]=static_cast<T>(col[i]);
    } 
  delete [] indx;
  delete [] col;
  return static_cast<T>(det);
}

template<typename T>
T
Matrix<T>::determinant() const
  /**
    Calculate the derminant of the matrix
    @return Determinant of matrix.
  */
{
  if (nx!=ny)
    throw Kernel::Exception::MisMatch<size_t>(nx,ny,
        "Determinant error :: Matrix is not NxN");

  Matrix<T> Mt(*this);         //temp copy
  T D=Mt.factor();
  return D;
}

template<typename T>
T
Matrix<T>::factor()
  /** 
     Gauss jordan diagonal factorisation 
     The diagonal is left as the values, 
     the lower part is zero.
     @return the factored matrix
  */
{
  if (nx!=ny || nx<1)
    throw std::runtime_error("Matrix::factor Matrix is not NxN");

  double Pmax;
  double deter=1.0;
  for(int i=0;i<static_cast<int>(nx)-1;i++)       //loop over each row
    {
      int jmax=i;
      Pmax=fabs(V[i][i]);
      for(int j=i+1;j<static_cast<int>(nx);j++)     // find max in Row i
        {
          if (fabs(V[i][j])>Pmax)
            {
              Pmax=fabs(V[i][j]);
              jmax=j;
            }
        }
      if (Pmax<1e-8)         // maxtrix signular 
        {
//          std::cerr<<"Matrix Singular"<<std::endl;
          return 0;
        }
      // Swap Columns
      if (i!=jmax) 
        {
          swapCols(i,jmax);
          deter*=-1;            //change sign.
        }
      // zero all rows below diagonal
      Pmax=V[i][i];
      deter*=Pmax;
      for(int k=i+1;k<static_cast<int>(nx);k++)  // row index
        {
          const double scale=V[k][i]/Pmax;
          V[k][i]=static_cast<T>(0);
          for(int q=i+1;q<static_cast<int>(nx);q++)   //column index
            V[k][q]-=static_cast<T>(scale*V[i][q]);
        }
    }
  deter*=V[nx-1][nx-1];
  return static_cast<T>(deter);
}

template<typename T>
void 
Matrix<T>::normVert()
  /**
    Normalise EigenVectors
    Assumes that they have already been calculated
  */
{
  for(size_t i=0;i<nx;i++)
  {
    T sum=0;
    for(size_t j=0;j<ny;j++)
    {
      sum += V[i][j]*V[i][j];
    }
    sum=static_cast<T>(std::sqrt(static_cast<double>(sum)));
    for(size_t j=0;j<ny;j++)
    {
      V[i][j]/=sum;
    }
  }
  return;
}

template<typename T>
T
Matrix<T>::compSum() const
  /**
    Add up each component sums for the matrix
    @return \f$ \sum_i \sum_j V_{ij}^2 \f$
   */
{
  T sum(0);
  for(size_t i=0;i<nx;i++)
  {
    for(size_t j=0;j<ny;j++)
    {
      sum+=V[i][j]*V[i][j];
    }
  }
  return sum;
}

template<typename T>
void 
Matrix<T>::lubcmp(int* rowperm,int& interchange)
  /**
    Find biggest pivot and move to top row. Then
    divide by pivot.
    @param interchange :: odd/even nterchange (+/-1)
    @param rowperm :: row permuations [nx values]
  */
{
  int imax(0),j,k;
  double sum,dum,big,temp;

  if (nx!=ny || nx<2)
    {
      std::cerr<<"Error with lubcmp"<<std::endl;
      return;
    }
  double *vv=new double[nx];
  interchange=1;
  for(int i=0;i<static_cast<int>(nx);i++)
    {
      big=0.0;
      for(j=0;j<static_cast<int>(nx);j++)
      if ((temp=fabs(V[i][j])) > big) 
          big=temp;

      if (big == 0.0) 
        {
          delete [] vv;
          return;
        }
      vv[i]=1.0/big;
    }

  for (j=0;j<static_cast<int>(nx);j++)
    {
      for(int i=0;i<j;i++)
        {
          sum=V[i][j];
          for(k=0;k<i;k++)
            sum-= V[i][k] * V[k][j];
          V[i][j]=static_cast<T>(sum);
        }
      big=0.0;
      imax=j;
      for (int i=j;i<static_cast<int>(nx);i++)
        {
          sum=V[i][j];
          for (k=0;k<j;k++)
            sum -= V[i][k] * V[k][j];
          V[i][j]=static_cast<T>(sum);
          if ( (dum=vv[i] * fabs(sum)) >=big)
            {
              big=dum;
              imax=i;
            }
        }

      if (j!=imax)
        {
          for(k=0;k<static_cast<int>(nx);k++)
            {                     //Interchange rows
              dum=V[imax][k];
              V[imax][k]=V[j][k];
              V[j][k]=static_cast<T>(dum);
            }
          interchange *= -1;
          vv[imax]=static_cast<T>(vv[j]);
        }
      rowperm[j]=imax;
      
      if (V[j][j] == 0.0) 
        V[j][j]=static_cast<T>(1e-14);
      if (j != static_cast<int>(nx)-1)
        {
          dum=1.0/(V[j][j]);
          for(int i=j+1;i<static_cast<int>(nx);i++)
            V[i][j] *= static_cast<T>(dum);
        }
    }
  delete [] vv;
  return;
}

template<typename T>
void 
Matrix<T>::lubksb(const int* rowperm,double* b)
  /**
    Impliments a separation of the Matrix
    into a triangluar matrix
  */
{
  int ii= -1;
 
  for(int i=0;i<static_cast<int>(nx);i++)
    {
      int ip=rowperm[i];
      double sum=b[ip];
      b[ip]=b[i];
      if (ii != -1) 
        for (int j=ii;j<i;j++) 
          sum -= V[i][j] * b[j];
      else if (sum) 
        ii=i;
      b[i]=sum;
    }

  for (int i=static_cast<int>(nx)-1;i>=0;i--)
    {
      double sum=static_cast<T>(b[i]);
      for (int j=i+1;j<static_cast<int>(nx);j++)
        sum -= V[i][j] * b[j];
      b[i]=sum/V[i][i];
    }
  return;
}

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

template<typename T> 
std::vector<T>
Matrix<T>::Diagonal() const
  /**
    Returns the diagonal form as a vector
    @return Diagonal elements
  */
{
  const size_t Msize=(ny>nx) ? nx : ny;
  std::vector<T> Diag(Msize);
  for(size_t 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
    @return Trace of matrix 
  */
{
  const size_t Msize=(ny>nx) ? nx : ny;
  T Trx = 0;
  for(size_t 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;
          throw(std::invalid_argument(" Matrix is not in an eigenvalue format"));
    }
  std::vector<int> index;
  std::vector<T> X=DiagMatrix.Diagonal();
  indexSort(X,index);
  Matrix<T> EigenVec(*this);
  for(size_t Icol=0;Icol<nx;Icol++)
  {
    for(size_t 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 
    @return :: 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(size_t i=0;i<nx;i++)
    for(size_t 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(size_t 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(size_t ip=0;ip<nx-1;ip++)
        for(size_t 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;
          for(size_t 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/static_cast<int>(nx*nx) : 0.0;

      for(int ip=0;ip<static_cast<int>(nx)-1;ip++)
        {
          for(int iq=ip+1;iq<static_cast<int>(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<static_cast<int>(nx);j++)
                    A.rotate(tau,sinAngle,ip,j,iq,j);
                  for(int j=0;j<static_cast<int>(nx);j++)
                    EigenVec.rotate(tau,sinAngle,j,ip,j,iq);
                  iteration++;
                }
            }
        }   
      for(size_t 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>
bool Matrix<T>::isRotation() const
/** Check if a matrix represents a proper rotation
@ return :: true/false
*/
{
  if (this->nx != this->ny) throw(std::invalid_argument("matrix is not square"));
//  std::cout << "Matrix determinant-1 is " << (this->determinant()-1) << std::endl;
  if (fabs(this->determinant()-1) >1e-5) 
  {
    return false;
  }
  else
  {
    Matrix<T> prod(nx,ny),ident(nx,ny,true);
    prod= this->operator*(this->Tprime());
//    std::cout << "Matrix * Matrix' = " << std::endl << prod << std::endl;
    return prod.equals(ident, 1e-5);
  }
}

template<typename T>
bool Matrix<T>::isOrthogonal() const
/** Check if a matrix is orthogonal. Same as isRotation, but allows determinant to be -1
@ return :: true/false
*/
{
  if (this->nx != this->ny) throw(std::invalid_argument("matrix is not square"));
  if (fabs(fabs(this->determinant())-1.) >1e-5) 
  {
    return false;
  }
  else
  {
    Matrix<T> prod(nx,ny),ident(nx,ny,true);
    prod= this->operator*(this->Tprime());
    return prod.equals(ident, 1e-7);
  }
}



template<typename T>
std::vector<T> Matrix<T>::toRotation() 
/**
  Transform the matrix to a rotation matrix, by normalizing each column to 1
  @return :: a vector of scaling factors
  @throw :: std::invalid_argument if the absolute value of the determinant is less then 1e-10 or not square matrix
*/
{
  if (this->nx != this->ny) throw(std::invalid_argument("matrix is not square"));
  if (fabs(this->determinant()) <1e-10) throw(std::invalid_argument("Determinant is too small"));
  // step 1: orthogonalize the matrix
  double spself,spother;
  for (size_t i=0; i<this->ny;++i)
  {
    spself=0.;
    for (size_t j=0; j<this->nx;++j) spself+=(V[j][i]*V[j][i]);
    for (size_t k=i+1; k<this->ny;++k) 
    {
      spother=0;
      for (size_t j=0; j<this->nx;++j) spother+=(V[j][i]*V[j][k]);
      for (size_t j=0; j<this->nx;++j) V[j][k]-=static_cast<T>(V[j][i]*spother/spself);
    }
    
  }
  // step 2: get scales and rescsale the matrix
  std::vector<T> scale(this->nx);
  T currentScale;
  for (size_t i=0; i<this->ny;++i)
  {
    currentScale=T(0.);
    for (size_t j=0; j<this->nx;++j)
      currentScale+=(V[j][i]*V[j][i]);
    currentScale=static_cast<T>(sqrt(static_cast<double>(currentScale)));
    if (currentScale <1e-10) throw(std::invalid_argument("Scale is too small"));
    scale[i]=currentScale;
  }
  Matrix<T> scalingMatrix(nx,ny),change(nx,ny,true);
  for (size_t i=0; i<this->ny;++i) scalingMatrix[i][i]=static_cast<T>(1.0/scale[i]);
  *this = this->operator*(scalingMatrix);
  if (this->determinant() <0.)
  {
    scale[0]=-scale[0];
    change[0][0]=(T)(-1);
    *this = this->operator*(change);
  }
  return scale;
}

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 size_t blockNumber((blockCnt>0) ? blockCnt : ny);  
  size_t BCnt(0);
  do
    {
      const size_t ACnt=BCnt;
      BCnt += blockNumber;
      if (BCnt>ny)
      {
          BCnt=ny;
      }

      if (ACnt)
      {
        Fh<<" ----- "<<ACnt<<" "<<BCnt<<" ------ "<<std::endl;
      }
      for(size_t i=0;i<nx;i++)
      {
        for(size_t 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 
    @return String value of output
  */
{
  std::ostringstream cx;
  for(size_t i=0;i<nx;i++)
  {
    for(size_t j=0;j<ny;j++)
    {
      cx<<std::setprecision(6)<<V[i][j]<<" ";
    }
  }
  return cx.str();
}

///\cond TEMPLATE

// Symbol definitions for common types
template class MANTID_KERNEL_DLL Matrix<double>;
template class MANTID_KERNEL_DLL Matrix<int>;
template class MANTID_KERNEL_DLL Matrix<float>;

template MANTID_KERNEL_DLL std::ostream& operator<<(std::ostream&, const DblMatrix&);
template MANTID_KERNEL_DLL std::istream& operator>>(std::istream&, DblMatrix&);
template MANTID_KERNEL_DLL std::ostream& operator<<(std::ostream&, const Matrix<float>&);
template MANTID_KERNEL_DLL std::istream& operator>>(std::istream&, Matrix<float>&);
template MANTID_KERNEL_DLL std::ostream& operator<<(std::ostream&, const IntMatrix&);
template MANTID_KERNEL_DLL std::istream& operator>>(std::istream&, IntMatrix&);
///\endcond TEMPLATE

} // namespace Kernel
} // namespace