fix bug in SVD, update comments in sph_harm.hpp

include/sctl/mat_utils.txx

@@ -107,7 +107,7 @@ template <> inline void cublasgemm<double>(char TransA, char TransB, int M, int
 //#define SCTL_SVD_DEBUG
 
 template <class ValueType> static inline void GivensL(Iterator<ValueType> S_, StaticArray<Long, 2> &dim, Long m, ValueType a, ValueType b) {
-  auto S = [S_,dim](Long i, Long j){ return S_[(i) * dim[1] + (j)]; };
+  auto S = [S_,dim](Long i, Long j) -> ValueType& { return S_[(i) * dim[1] + (j)]; };
 
   ValueType r = sqrt<ValueType>(a * a + b * b);
   if (r == 0) return;
@@ -127,7 +127,7 @@ template <class ValueType> static inline void GivensL(Iterator<ValueType> S_, St
 }
 
 template <class ValueType> static inline void GivensR(Iterator<ValueType> S_, StaticArray<Long, 2> &dim, Long m, ValueType a, ValueType b) {
-  auto S = [S_,dim](Long i, Long j){ return S_[(i) * dim[1] + (j)]; };
+  auto S = [S_,dim](Long i, Long j) -> ValueType& { return S_[(i) * dim[1] + (j)]; };
 
   ValueType r = sqrt<ValueType>(a * a + b * b);
   if (r == 0) return;
@@ -147,9 +147,9 @@ template <class ValueType> static inline void GivensR(Iterator<ValueType> S_, St
 }
 
 template <class ValueType> static inline void SVD(StaticArray<Long, 2> &dim, Iterator<ValueType> U_, Iterator<ValueType> S_, Iterator<ValueType> V_, ValueType eps = -1) {
-  auto U = [U_,dim](Long i, Long j){ return U_[(i) * dim[0] + (j)]; };
-  auto S = [S_,dim](Long i, Long j){ return S_[(i) * dim[1] + (j)]; };
-  auto V = [V_,dim](Long i, Long j){ return V_[(i) * dim[1] + (j)]; };
+  auto U = [U_,dim](Long i, Long j) -> ValueType& { return U_[(i) * dim[0] + (j)]; };
+  auto S = [S_,dim](Long i, Long j) -> ValueType& { return S_[(i) * dim[1] + (j)]; };
+  auto V = [V_,dim](Long i, Long j) -> ValueType& { return V_[(i) * dim[1] + (j)]; };
 
   assert(dim[0] >= dim[1]);
 #ifdef SCTL_SVD_DEBUG

include/sctl/matrix.txx

@@ -17,7 +17,7 @@ template <class ValueType> std::ostream& operator<<(std::ostream& output, const
   for (Long i = 0; i < M.Dim(0); i++) {
     for (Long j = 0; j < M.Dim(1); j++) {
       float f = ((float)M(i, j));
-      if (fabs<ValueType>(f) < 1e-25) f = 0;
+      if (fabs<float>(f) < 1e-25) f = 0;
       output << std::setw(10) << ((double)f) << ' ';
     }
     output << ";\n";

include/sctl/sph_harm.hpp

@@ -36,8 +36,8 @@ template <class Real> class SphericalHarmonics{
     /**
      * \brief Compute spherical harmonic coefficients from grid values.
      * \param[in] X Grid values {X(t0,p0), X(t0,p1), ... , X(t1,p0), X(t1,p1), ... }, where, {cos(t0), cos(t1), ... } are the Gauss-Legendre nodes of order (Nt-1) in the interval [-1,1] and {p0, p1, ... } are equispaced in [0, 2*pi].
-     * \param[in] Nt Number of grid points \theta \in (1,pi).
-     * \param[in] Np Number of grid points \phi \in (1,2*pi).
+     * \param[in] Nt Number of grid points \theta \in (0,pi).
+     * \param[in] Np Number of grid points \phi \in (0,2*pi).
      * \param[in] p Order of spherical harmonic expansion.
      * \param[in] arrange Arrangement of the coefficients.
      * \param[out] S Spherical harmonic coefficients.
@@ -49,8 +49,8 @@ template <class Real> class SphericalHarmonics{
      * \param[in] S Spherical harmonic coefficients.
      * \param[in] arrange Arrangement of the coefficients.
      * \param[in] p Order of spherical harmonic expansion.
-     * \param[in] Nt Number of grid points \theta \in (1,pi).
-     * \param[in] Np Number of grid points \phi \in (1,2*pi).
+     * \param[in] Nt Number of grid points \theta \in (0,pi).
+     * \param[in] Np Number of grid points \phi \in (0,2*pi).
      * \param[out] X Grid values {X(t0,p0), X(t0,p1), ... , X(t1,p0), X(t1,p1), ... }, where, {cos(t0), cos(t1), ... } are the Gauss-Legendre nodes of order (Nt-1) in the interval [-1,1] and {p0, p1, ... } are equispaced in [0, 2*pi].
      * \param[out] X_theta \theta derivative of X evaluated at grid points.
      * \param[out] X_phi \phi derivative of X evaluated at grid points.
@@ -77,8 +77,8 @@ template <class Real> class SphericalHarmonics{
     /**
      * \brief Compute vector spherical harmonic coefficients from grid values.
      * \param[in] X Grid values {X(t0,p0), X(t0,p1), ... , X(t1,p0), ... , Y(t0,p0), ... , Z(t0,p0), ... }, where, {cos(t0), cos(t1), ... } are the Gauss-Legendre nodes of order (Nt-1) in the interval [-1,1] and {p0, p1, ... } are equispaced in [0, 2*pi].
-     * \param[in] Nt Number of grid points \theta \in (1,pi).
-     * \param[in] Np Number of grid points \phi \in (1,2*pi).
+     * \param[in] Nt Number of grid points \theta \in (0,pi).
+     * \param[in] Np Number of grid points \phi \in (0,2*pi).
      * \param[in] p Order of spherical harmonic expansion.
      * \param[in] arrange Arrangement of the coefficients.
      * \param[out] S Vector spherical harmonic coefficients.
@@ -90,8 +90,8 @@ template <class Real> class SphericalHarmonics{
      * \param[in] S Vector spherical harmonic coefficients.
      * \param[in] arrange Arrangement of the coefficients.
      * \param[in] p Order of spherical harmonic expansion.
-     * \param[in] Nt Number of grid points \theta \in (1,pi).
-     * \param[in] Np Number of grid points \phi \in (1,2*pi).
+     * \param[in] Nt Number of grid points \theta \in (0,pi).
+     * \param[in] Np Number of grid points \phi \in (0,2*pi).
      * \param[out] X Grid values {X(t0,p0), X(t0,p1), ... , X(t1,p0), X(t1,p1), ... , Y(t0,p0), ... , Z(t0,p0), ... }, where, {cos(t0), cos(t1), ... } are the Gauss-Legendre nodes of order (Nt-1) in the interval [-1,1] and {p0, p1, ... } are equispaced in [0, 2*pi].
      */
     static void VecSHC2Grid(const Vector<Real>& S, SHCArrange arrange, Long p, Long Nt, Long Np, Vector<Real>& X);
@@ -130,6 +130,12 @@ template <class Real> class SphericalHarmonics{
 
     static void StokesEvalKSelf(const Vector<Real>& S, SHCArrange arrange, Long p, const Vector<Real>& coord, bool interior, Vector<Real>& U);
 
+    /**
+     * \brief Nodes and weights for Gauss-Legendre quadrature rule
+     */
+    static const Vector<Real>& LegendreNodes(Long p1);
+    static const Vector<Real>& LegendreWeights(Long p1);
+
     static void test_stokes() {
       int p = 6;
       int dof = 3;
@@ -423,8 +429,6 @@ template <class Real> class SphericalHarmonics{
     static void LegPolyDeriv(Vector<Real>& poly_val, const Vector<Real>& X, Long degree);
     static void LegPolyDeriv_(Vector<Real>& poly_val, const Vector<Real>& X, Long degree);
 
-    static const Vector<Real>& LegendreNodes(Long p1);
-    static const Vector<Real>& LegendreWeights(Long p1);
     static const Vector<Real>& SingularWeights(Long p1);
 
     static const Matrix<Real>& MatFourier(Long p0, Long p1);

include/sctl/sph_harm.txx

@@ -205,7 +205,7 @@ template <class Real> void SphericalHarmonics<Real>::WriteVTK(const char* fname,
     }
 
     if(v_ptr) {
-      Long data__dof = V.Dim() / (2*p1*(p1+1));
+      Long data__dof = V.Dim() / (N_ves * 2*p1*(p1+1));
       for(Long k=0;k<N_ves;k++){ // Set point_value
         for(Long i=0;i<p1+1;i++){
           for(Long j=0;j<2*p1;j++){