svd_quda.h

Go to the documentation of this file.

#ifndef _SVD_QUDA_H_
#define _SVD_QUDA_H_

#define DEVICEHOST __device__ __host__
#define SVDPREC 1e-11
#define LOG2 0.69314718055994530942


//namespace quda{

template<class Cmplx> 
  inline DEVICEHOST
typename RealTypeId<Cmplx>::Type cabs(const Cmplx & z)
{
  typename RealTypeId<Cmplx>::Type max, ratio, square;
  if(fabs(z.x) > fabs(z.y)){ max = z.x; ratio = z.y/max; }else{ max=z.y; ratio = z.x/max; }
  square = max*max*(1.0 + ratio*ratio);
  return sqrt(square);
}


/*
   template<class Cmplx>
   inline DEVICEHOST
   typename RealTypeId<Cmplx>::Type cabs(const Cmplx & z)
   {
   return sqrt(z.x*z.x + z.y*z.y);
   }
   */

template<class T, class U> 
inline DEVICEHOST typename PromoteTypeId<T,U>::Type quadSum(const T & a, const U & b){
  typename PromoteTypeId<T,U>::Type ratio, square, max;
  if(fabs(a) > fabs(b)){ max = a; ratio = b/a; }else{ max=b; ratio = a/b; }
  square = max*max*(1.0 + ratio*ratio);
  return sqrt(square);
}

/*
   template<class T, class U>
   inline DEVICEHOST
   T quadSum(const T & a, const U & b)
   {
   return sqrt(a*a + b*b);
   }
   */



// In future iterations of the code, I would like to use templates to compute the norm
DEVICEHOST
float getNorm(const Array<float2,3>& a){
  float temp1, temp2, temp3;
  temp1 = cabs(a[0]);
  temp2 = cabs(a[1]);
  temp3 = quadSum(temp1,temp2);
  temp1 = cabs(a[2]);
  return quadSum(temp1,temp3);
}


DEVICEHOST
double getNorm(const Array<double2,3>& a){
  double temp1, temp2, temp3;
  temp1 = cabs(a[0]);
  temp2 = cabs(a[1]);
  temp3 = quadSum(temp1,temp2);
  temp1 = cabs(a[2]);
  return quadSum(temp1,temp3);
}


template<class T>
  DEVICEHOST 
void constructHHMat(const T & tau, const Array<T,3> & v, Matrix<T,3> & hh)
{
  Matrix<T,3> temp1, temp2;
  outerProd(v,v,&temp1);

  temp2 = Conj(tau)*temp1;    

  setIdentity(&temp1);

  hh =  temp1 - temp2; 
  return;
}


template<class Real>
DEVICEHOST
void getLambdaMax(const Matrix<Real,3> & b, Real & lambda_max){

  Real m11 = b(1,1)*b(1,1) + b(0,1)*b(0,1);
  Real m22 = b(2,2)*b(2,2) + b(1,2)*b(1,2);
  Real m12 = b(1,1)*b(1,2);
  Real dm  = (m11 - m22)/2.0;

  Real norm1 = quadSum(dm, m12);
  if( dm >= 0.0 ){ 
    lambda_max = m22 - (m12*m12)/(norm1 + dm);
  }else{
    lambda_max = m22 + (m12*m12)/(norm1 - dm);
  }
  return;
}


template<class Real>
  DEVICEHOST
void getGivensRotation(const Real & alpha, const Real & beta, Real & c, Real & s)
{
  Real ratio;
  if( beta == 0.0 ){
    c = 1.0; s = 0.0;
  }else if( fabs(beta) > fabs(alpha) ){
    ratio = -alpha/beta;
    s = rsqrt(1.0 + ratio*ratio);
    c = ratio*s;
  }else{
    ratio = -beta/alpha;
    c = rsqrt(1.0 + ratio*ratio);
    s = ratio*c; 
  }
  return;
}

template<class Real>
inline DEVICEHOST
void accumGivensRotation(int index, const Real & c, const Real & s, Matrix<Real,3> & m){
  int index_p1 = index+1;
  Real alpha, beta;
  for(int i=0; i<3; ++i){
    alpha = m(i,index);
    beta  = m(i,index_p1);
    m(i,index) = alpha*c - beta*s;
    m(i,index_p1) = alpha*s + beta*c;
  }
  return;
}

template<class Real>
inline DEVICEHOST
void assignGivensRotation(const Real & c, const Real & s, Matrix<Real,2> & m){
  m(0,0) = c;
  m(0,1) = s;
  m(1,0) = -s;
  m(1,1) =  c;
  return;
}

template<class Real>
inline DEVICEHOST
void swap(Real & a, Real & b){
  Real temp = a;
  a = b;
  b = temp;
  return;
}

template<class Real>
inline DEVICEHOST
void smallSVD(Matrix<Real,2> & u, Matrix<Real,2> & v, Matrix<Real,2> & m){
  // set u and v to the 2x2 identity matrix
  setIdentity(&u);
  setIdentity(&v);

  Real c, s;

  if( m(0,0) == 0.0 ){

    getGivensRotation(m(0,1), m(1,1), c, s);

    m(0,0) = m(0,1)*c - m(1,1)*s;
    m(0,1) = 0.0;
    m(1,1) = 0.0;

    // exchange columns in v 
    v(0,0) = 0.0;
    v(0,1) = 1.0;
    v(1,0) = 1.0;
    v(1,1) = 0.0;

    // u is a Givens rotation 
    assignGivensRotation(c, s, u);

  }else if( m(1,1) == 0.0 ){

    getGivensRotation(m(0,0), m(0,1), c, s);  

    m(0,0) = m(0,0)*c - m(0,1)*s;
    m(0,1) = 0.0;
    m(1,1) = 0.0;

    assignGivensRotation(c,s,v);

  }else if( m(0,1) != 0.0 ){

    // need to calculate (m(1,1)**2 + m(0,1)**2 - m(0,0)**2)/(2*m(0,0)*m(0,1))
    Real abs01 = fabs(m(0,1));
    Real abs11 = fabs(m(1,1));
    Real min, max;
    if( abs01 > abs11 ){ min = abs11; max = abs01; } 
    else { min = abs01; max = abs11; }


    Real ratio = min/max;
    Real alpha = 2.0*log(max) + log(1.0 + ratio*ratio);

    Real abs00 = fabs(m(0,0));
    Real beta = 2.0*log(abs00);

    int sign;
    Real temp;

    if( alpha > beta ){
      sign = 1;
      temp = alpha + log(1.0 - exp(beta-alpha));
    }else{
      sign = -1;        
      temp = beta + log(1.0 - exp(alpha-beta));
    }
    temp -= LOG2 + log(abs00) + log(abs01);
    temp = sign*exp(temp);

    if( m(0,0) < 0.0 ){ temp *= -1.0; }
    if( m(0,1) < 0.0 ){ temp *= -1.0; }

    beta = quadSum(1.0, temp);

    if( temp >= 0.0 ){
      temp = 1.0/(temp + beta);
    }else{
      temp = 1.0/(temp - beta);
    }

    // Calculate beta = sqrt(1 + temp**2)
    beta = quadSum(1.0, temp);

    c = 1.0/beta;
    s = temp*c;

    Matrix<Real,2> p;

    p(0,0) = c*m(0,0) - s*m(0,1);
    p(1,0) =          - s*m(1,1);
    p(0,1) = s*m(0,0) + c*m(0,1);
    p(1,1) = c*m(1,1);

    assignGivensRotation(c, s, v);

    // Make the column with the largest norm the first column
    alpha = quadSum(p(0,0),p(1,0));
    beta  = quadSum(p(0,1),p(1,1));

    if(alpha < beta){
      swap(p(0,0),p(0,1));
      swap(p(1,0),p(1,1));

      swap(v(0,0),v(0,1));
      swap(v(1,0),v(1,1));
    }    

    getGivensRotation(p(0,0), p(1,0), c, s);

    m(0,0) = p(0,0)*c - s*p(1,0);
    m(1,1) = p(0,1)*s + c*p(1,1);
    m(0,1) = 0.0;

    assignGivensRotation(c,s,u);        
  }

  return;
}

// Change this so that the bidiagonal terms are not returned 
// as a complex matrix
template<class Cmplx>
  DEVICEHOST
void getRealBidiagMatrix(const Matrix<Cmplx,3> & mat, 
    Matrix<Cmplx,3> & u,
    Matrix<Cmplx,3> & v)
{
  Matrix<Cmplx,3> p;
  Array<Cmplx,3> vec;
  Matrix<Cmplx,3> temp;


  const Cmplx COMPLEX_UNITY = makeComplex<Cmplx>(1,0);
  const Cmplx COMPLEX_ZERO  = makeComplex<Cmplx>(0,0);
  // Step 1: build the first left reflector v,
  //          calculate the first left rotation
  //          apply to the original matrix
  typename RealTypeId<Cmplx>::Type x = cabs(mat(1,0));
  typename RealTypeId<Cmplx>::Type y = cabs(mat(2,0));
  typename RealTypeId<Cmplx>::Type norm1 = quadSum(x,y);
  typename RealTypeId<Cmplx>::Type beta;
  Cmplx w, tau, z;

  if(norm1 == 0 && mat(0,0).y == 0){
    p = mat;
  }else{
    Array<Cmplx,3> temp_vec;
    copyColumn(mat, 0, &temp_vec);

    beta = getNorm(temp_vec); 

    if(mat(0,0).x > 0.0){ beta = -beta; }

    w.x = mat(0,0).x - beta; // work around for LLVM
    w.y = mat(0,0).y;
    norm1 = cabs(w);
    w = Conj(w)/norm1; 

    // Assign the vector elements
    vec[0] = COMPLEX_UNITY; 
    vec[1] = mat(1,0)*w/norm1;
    vec[2] = mat(2,0)*w/norm1;

    // Now compute tau
    tau.x =  (beta - mat(0,0).x)/beta;
    tau.y =  mat(0,0).y/beta;
    // construct the Householder matrix
    constructHHMat(tau, vec, temp);
    p = conj(temp)*mat;    
    u = temp;
  }

  // Step 2: build the first right reflector
  typename RealTypeId<Cmplx>::Type norm2 = cabs(p(0,2));
  if(norm2 != 0.0 || p(0,1).y != 0.0){
    norm1 = cabs(p(0,1));
    beta  = quadSum(norm1,norm2);
    vec[0] = COMPLEX_ZERO;
    vec[1] = COMPLEX_UNITY;  

    if( p(0,1).x > 0.0 ){ beta = -beta; }
    w.x = p(0,1).x-beta; // work around for LLVM
    w.y = p(0,1).y;
    norm1 = cabs(w);
    w = Conj(w)/norm1; 
    z = Conj(p(0,2))/norm1;
    vec[2] = z*Conj(w);

    tau.x = (beta - p(0,1).x)/beta; // work around for LLVM
    tau.y = (- p(0,1).y)/beta;
    // construct the Householder matrix
    constructHHMat(tau, vec, temp);
    p = p*temp;
    v = temp;
  }

  // Step 3: build the second left reflector
  norm2 = cabs(p(2,1));
  if(norm2 != 0.0 || p(1,1).y != 0.0){
    norm1 = cabs(p(1,1));
    beta  = quadSum(norm1,norm2);

    // Set the first two elements of the vector
    vec[0] = COMPLEX_ZERO;
    vec[1] = COMPLEX_UNITY;

    if( p(1,1).x > 0 ){ beta = -beta; }
    w.x = p(1,1).x - beta; // work around for LLVM
    w.y = p(1,1).y;
    norm1 = cabs(w);
    w = Conj(w)/norm1;
    z = p(2,1)/norm1;
    vec[2] = z*w;

    tau.x  = (beta - p(1,1).x)/beta;
    tau.y  = p(1,1).y/beta;
    // I could very easily change the definition of tau
    // so that we wouldn't need to call conj(temp) below.
    // I would have to be careful to make sure this change 
    // is consistent with the other parts of the code.
    constructHHMat(tau, vec, temp);
    p = conj(temp)*p;
    u = u*temp;
  }

  // Step 4: build the second right reflector
  setIdentity(&temp);
  if( p(1,2).y != 0.0 ){
    beta = cabs(p(1,2));
    if( p(1,2).x > 0.0 ){ beta = -beta; }
    temp(2,2) = Conj(p(1,2))/beta;
    p(2,2) = p(2,2)*temp(2,2); // This is the only element of p needed below
    v = v*temp;
  }


  // Step 5: build the third left reflector
  if( p(2,2).y != 0.0 ){
    beta = cabs(p(2,2));
    if( p(2,2).x > 0.0 ){ beta = -beta; }
    temp(2,2) = p(2,2)/beta;
    u = u*temp;
  }
  return;
}



template<class Real>
  DEVICEHOST
void bdSVD(Matrix<Real,3>& u, Matrix<Real,3>& v, Matrix<Real,3>& b, int max_it)
{

  Real c,s;

  // set u and v matrices equal to the identity
  setIdentity(&u);
  setIdentity(&v);

  Real alpha, beta;

  int it=0;
  do{  

    if( fabs(b(0,1)) < SVDPREC*( fabs(b(0,0)) + fabs(b(1,1)) ) ){ b(0,1) = 0.0; } 
    if( fabs(b(1,2)) < SVDPREC*( fabs(b(0,0)) + fabs(b(2,2)) ) ){ b(1,2) = 0.0; }

    if( b(0,1) != 0.0 && b(1,2) != 0.0 ){
      if( b(0,0) == 0.0 ){

        getGivensRotation(-b(0,1), b(1,1), s, c);

        for(int i=0; i<3; ++i){
          alpha = u(i,0);
          beta = u(i,1);
          u(i,0) = alpha*c - beta*s;
          u(i,1) = alpha*s + beta*c;
        }

        b(1,1) = b(0,1)*s + b(1,1)*c;
        b(0,1) = 0.0;

        b(0,2) = -b(1,2)*s;
        b(1,2) *= c;

        getGivensRotation(-b(0,2), b(2,2), s, c);

        for(int i=0; i<3; ++i){
          alpha = u(i,0);
          beta = u(i,2);
          u(i,0) = alpha*c - beta*s;
          u(i,2) = alpha*s + beta*c;
        }
        b(2,2) = b(0,2)*s + b(2,2)*c;
        b(0,2) = 0.0;

      }else if( b(1,1) == 0.0 ){
        // same block
        getGivensRotation(-b(1,2), b(2,2), s, c);
        for(int i=0; i<3; ++i){
          alpha = u(i,1);
          beta = u(i,2);
          u(i,1) = alpha*c - beta*s;
          u(i,2) = alpha*s + beta*c;
        }
        b(2,2) = b(1,2)*s + b(2,2)*c;
        b(1,2) = 0.0;
        // end block
      }else if( b(2,2) == 0.0 ){

        getGivensRotation(b(1,1), -b(1,2), c, s);
        for(int i=0; i<3; ++i){
          alpha = v(i,1);
          beta = v(i,2);
          v(i,1) = alpha*c + beta*s;
          v(i,2) = -alpha*s + beta*c;
        }
        b(1,1) = b(1,1)*c + b(1,2)*s;
        b(1,2) = 0.0;

        b(0,2) = -b(0,1)*s;
        b(0,1) *= c;

        // apply second rotation, to remove b02
        getGivensRotation(b(0,0), -b(0,2), c, s);
        for(int i=0; i<3; ++i){
          alpha = v(i,0);
          beta = v(i,2);
          v(i,0) = alpha*c + beta*s;
          v(i,2) = -alpha*s + beta*c;
        }
        b(0,0) = b(0,0)*c + b(0,2)*s;
        b(0,2) = 0.0;

      }else{ // Else entering normal QR iteration

        Real lambda_max; 
        getLambdaMax(b, lambda_max); // defined above

        alpha = b(0,0)*b(0,0) - lambda_max;
        beta  = b(0,0)*b(0,1);

        // c*beta + s*alpha = 0
        getGivensRotation(alpha, beta, c, s);
        // Multiply right v matrix 
        accumGivensRotation(0, c, s, v);
        // apply to bidiagonal matrix, this generates b(1,0)
        alpha = b(0,0);
        beta  = b(0,1);  

        b(0,0) = alpha*c - beta*s;
        b(0,1) = alpha*s + beta*c;
        b(1,0) = -b(1,1)*s;
        b(1,1) = b(1,1)*c;

        // Calculate the second Givens rotation (this time on the left)
        getGivensRotation(b(0,0), b(1,0), c, s);

        // Multiply left u matrix
        accumGivensRotation(0, c, s, u);

        b(0,0) = b(0,0)*c - b(1,0)*s;
        alpha  = b(0,1);
        beta   = b(1,1);
        b(0,1) = alpha*c - beta*s;
        b(1,1) = alpha*s + beta*c;
        b(0,2) = -b(1,2)*s;
        b(1,2) = b(1,2)*c;
        // Following from the definition of the Givens rotation, b(1,0) should be equal to zero.
        b(1,0) = 0.0; 

        // Calculate the third Givens rotation (on the right)
        getGivensRotation(b(0,1), b(0,2), c, s);

        // Multiply the right v matrix
        accumGivensRotation(1, c, s, v);

        b(0,1) = b(0,1)*c - b(0,2)*s;
        alpha = b(1,1);
        beta  = b(1,2);

        b(1,1) = alpha*c - beta*s;
        b(1,2) = alpha*s + beta*c;
        b(2,1) = -b(2,2)*s;
        b(2,2) *= c;
        b(0,2) = 0.0;


        // Calculate the fourth Givens rotation (on the left)
        getGivensRotation(b(1,1), b(2,1), c, s);
        // Multiply left u matrix
        accumGivensRotation(1, c, s, u);
        // Eliminate b(2,1)
        b(1,1) = b(1,1)*c - b(2,1)*s;
        alpha = b(1,2);
        beta  = b(2,2);
        b(1,2) = alpha*c - beta*s;
        b(2,2) = alpha*s + beta*c;
        b(2,1) = 0.0;

      } // end of normal QR iteration


    }else if( b(0,1) == 0.0 ){ 
      Matrix<Real,2> m_small, u_small, v_small;

      m_small(0,0) = b(1,1);
      m_small(0,1) = b(1,2);
      m_small(1,0) = b(2,1);
      m_small(1,1) = b(2,2);

      smallSVD(u_small, v_small, m_small); 

      b(1,1) = m_small(0,0);    
      b(1,2) = m_small(0,1);    
      b(2,1) = m_small(1,0);    
      b(2,2) = m_small(1,1);    


      for(int i=0; i<3; ++i){
        alpha = u(i,1);
        beta  = u(i,2);
        u(i,1) = alpha*u_small(0,0) + beta*u_small(1,0);
        u(i,2) = alpha*u_small(0,1) + beta*u_small(1,1);

        alpha = v(i,1);
        beta  = v(i,2);
        v(i,1) = alpha*v_small(0,0) + beta*v_small(1,0);
        v(i,2) = alpha*v_small(0,1) + beta*v_small(1,1);
      }


    }else if( b(1,2) == 0.0 ){
      Matrix<Real,2> m_small, u_small, v_small;

      m_small(0,0) = b(0,0);
      m_small(0,1) = b(0,1);
      m_small(1,0) = b(1,0);
      m_small(1,1) = b(1,1);    

      smallSVD(u_small, v_small, m_small); 

      b(0,0) = m_small(0,0);
      b(0,1) = m_small(0,1);
      b(1,0) = m_small(1,0);
      b(1,1) = m_small(1,1);

      for(int i=0; i<3; ++i){
        alpha = u(i,0);
        beta  = u(i,1);
        u(i,0) = alpha*u_small(0,0) + beta*u_small(1,0);
        u(i,1) = alpha*u_small(0,1) + beta*u_small(1,1);

        alpha = v(i,0);
        beta  = v(i,1);
        v(i,0) = alpha*v_small(0,0) + beta*v_small(1,0);
        v(i,1) = alpha*v_small(0,1) + beta*v_small(1,1);
      }

    } // end if b(1,2) == 0
    it++;
  } while( (b(0,1) != 0.0 || b(1,2) != 0.0) && it < max_it);


  for(int i=0; i<3; ++i){
    if( b(i,i) < 0.0) {
      b(i,i) *= -1;
      for(int j=0; j<3; ++j){
        v(j,i) *= -1;
      }
    }
  }
  return;
}



template<class Cmplx>
  DEVICEHOST
void computeSVD(const Matrix<Cmplx,3> & m, 
    Matrix<Cmplx,3>&  u,
    Matrix<Cmplx,3>&  v,
    typename RealTypeId<Cmplx>::Type singular_values[3])
{

  getRealBidiagMatrix<Cmplx>(m, u, v);
  Matrix<typename RealTypeId<Cmplx>::Type,3> bd, u_real, v_real;
  // Make real
  for(int i=0; i<3; ++i){
    for(int j=0; j<3; ++j){
      bd(i,j) = (conj(u)*m*(v))(i,j).x;
    }
  }

  bdSVD(u_real, v_real, bd, 500);
  for(int i=0; i<3; ++i){
    singular_values[i] = bd(i,i);
  }

  u = u*u_real;
  v = v*v_real;

  return;
}

//} // end namespace quda


#endif // _SVD_QUDA_H