text: /*
     Copyright (c) 2009-2016, Jack Poulson
     All rights reserved.
  
     This file is part of Elemental and is under the BSD 2-Clause License,
     which can be found in the LICENSE file in the root directory, or at
     http://opensource.org/licenses/BSD-2-Clause
  */
  #include <El.hpp>
  
  namespace El {
  
  template<typename Real>
  void FoxLi( Matrix<Complex<Real>>& A, Int n, Real omega )
  {
      EL_DEBUG_CSE
      typedef Complex<Real> C;
      const Real pi = 4*Atan( Real(1) );
      const C phi = Sqrt( C(0,omega/pi) );
  
      // Compute Gauss quadrature points and weights
      Matrix<Real> d, e;
      Zeros( d, n, 1 );
      e.Resize( n-1, 1 );
      for( Int j=0; j<n-1; ++j )
      {
          const Real betaInv = 2*Sqrt(1-Pow(j+Real(1),-2)/4);
          e(j) = 1/betaInv;
      }
      Matrix<Real> x, Z;
      HermitianTridiagEigCtrl<Real> ctrl;
      ctrl.sort = UNSORTED;
      HermitianTridiagEig( d, e, x, Z, ctrl );
      auto z = Z( IR(0), ALL );
      Matrix<Real> sqrtWeights( z ), sqrtWeightsTrans;
      for( Int j=0; j<n; ++j )
          sqrtWeights(0,j) = Sqrt(Real(2))*Abs(sqrtWeights(0,j));
      auto sortPairs = TaggedSort( x, ASCENDING );
      for( Int j=0; j<n; ++j )
          x(j) = sortPairs[j].value;
      ApplyTaggedSortToEachRow( sortPairs, sqrtWeights );
      Transpose( sqrtWeights, sqrtWeightsTrans );
  
      // Form the integral operator
      A.Resize( n, n );
      for( Int j=0; j<n; ++j )
      {
          for( Int i=0; i<n; ++i )
          {
              const Real theta = -omega*Pow(x(i)-x(j),2);
              const Real realPart = Cos(theta);
              const Real imagPart = Sin(theta);
              A(i,j) = phi*C(realPart,imagPart);
          }
      }
  
      // Apply the weighting
      DiagonalScale( LEFT, NORMAL, sqrtWeightsTrans, A );
      DiagonalScale( RIGHT, NORMAL, sqrtWeightsTrans, A );
  }
  
  template<typename Real>
  void FoxLi( AbstractDistMatrix<Complex<Real>>& APre, Int n, Real omega )
  {
      EL_DEBUG_CSE
      typedef Complex<Real> C;
      const Real pi = 4*Atan( Real(1) );
      const C phi = Sqrt( C(0,omega/pi) );
  
      DistMatrixWriteProxy<C,C,MC,MR> AProx( APre );
      auto& A = AProx.Get();
  
      // Compute Gauss quadrature points and weights
      const Grid& g = A.Grid();
      DistMatrix<Real,VR,STAR> d(g), e(g);
      Zeros( d, n, 1 );
      e.Resize( n-1, 1 );
      auto& eLoc = e.Matrix();
      for( Int iLoc=0; iLoc<e.LocalHeight(); ++iLoc )
      {
          const Int i = e.GlobalRow(iLoc);
          const Real betaInv = 2*Sqrt(1-Pow(i+Real(1),-2)/4);
          eLoc(iLoc) = 1/betaInv;
      }
      DistMatrix<Real,VR,STAR> x(g);
      DistMatrix<Real,STAR,VR> Z(g);
      HermitianTridiagEigCtrl<Real> ctrl;
      ctrl.sort = UNSORTED;
      HermitianTridiagEig( d, e, x, Z, ctrl );
      auto z = Z( IR(0), ALL );
      DistMatrix<Real,STAR,VR> sqrtWeights( z );
      auto& sqrtWeightsLoc = sqrtWeights.Matrix();
      for( Int jLoc=0; jLoc<sqrtWeights.LocalWidth(); ++jLoc )
          sqrtWeightsLoc(0,jLoc) = Sqrt(Real(2))*Abs(sqrtWeightsLoc(0,jLoc));
      auto sortPairs = TaggedSort( x, ASCENDING );
      for( Int j=0; j<n; ++j )
          x.Set( j, 0, sortPairs[j].value );
      ApplyTaggedSortToEachRow( sortPairs, sqrtWeights );
  
      // Form the integral operator
      A.Resize( n, n );
      DistMatrix<Real,MC,STAR> x_MC( A.Grid() );
      DistMatrix<Real,MR,STAR> x_MR( A.Grid() );
      x_MC.AlignWith( A );
      x_MR.AlignWith( A );
      x_MC = x;
      x_MR = x;
      auto& ALoc = A.Matrix();
      auto& x_MCLoc = x_MC.Matrix();
      auto& x_MRLoc = x_MR.Matrix();
      for( Int jLoc=0; jLoc<A.LocalWidth(); ++jLoc )
      {
          for( Int iLoc=0; iLoc<A.LocalHeight(); ++iLoc )
          {
              const Real diff = x_MCLoc(iLoc)-x_MRLoc(jLoc);
              const Real theta = -omega*Pow(diff,2);
              const Real realPart = Cos(theta);
              const Real imagPart = Sin(theta);
              ALoc(iLoc,jLoc) = phi*C(realPart,imagPart);
          }
      }
  
      // Apply the weighting
      DistMatrix<Real,VR,STAR> sqrtWeightsTrans(g);
      Transpose( sqrtWeights, sqrtWeightsTrans );
      DiagonalScale( LEFT, NORMAL, sqrtWeightsTrans, A );
      DiagonalScale( RIGHT, NORMAL, sqrtWeightsTrans, A );
  }
  
  #define PROTO(Real) \
    template void FoxLi( Matrix<Complex<Real>>& A, Int n, Real omega ); \
    template void FoxLi \
    ( AbstractDistMatrix<Complex<Real>>& A, Int n, Real omega );
  
  #define EL_NO_INT_PROTO
  #define EL_NO_COMPLEX_PROTO
  #define EL_ENABLE_QUAD
  #define EL_ENABLE_DOUBLEDOUBLE
  #define EL_ENABLE_QUADDOUBLE
  #define EL_ENABLE_BIGFLOAT
  #include <El/macros/Instantiate.h>
  
  } // namespace El
  