...EJB

Author: ebylaska

Nwpw/nwpwlib/lattice/PGrid.cpp

@@ -856,12 +856,13 @@ void PGrid::cc_pack_ndot(const int nb, const int nn, double *a, double *b, doubl
 {
    int one = 1;
    // int ng  = 2*(nida[nb]+nidb[nb]);
-   int ng = 2 * (nida[nb] + nidb[nb]);
-   int ng0 = 2 * nida[nb];
+   int ng  = 2*(nida[nb] + nidb[nb]);
+   int ng0 = 2*nida[nb];
 
-   for (int i = 0; i < nn; ++i) {
-     sum[i] = 2.0 * DDOT_PWDFT(ng, &(a[i * ng]), one, b, one);
-     sum[i] -= DDOT_PWDFT(ng0, &(a[i * ng]), one, b, one);
+   for (size_t i=0; i<nn; ++i) 
+   {
+      sum[i] = 2.0*DDOT_PWDFT(ng, &(a[i * ng]), one, b, one);
+      sum[i] -= DDOT_PWDFT(ng0, &(a[i * ng]), one, b, one);
    }
 
    parall->Vector_SumAll(1, nn, sum);
@@ -3879,17 +3880,36 @@ void PGrid::tc_pack_Mul(const int nb, const double *a, double *c) {
  *    PGrid:tcc_pack_aMulAdd    *
  *                              *
  ********************************/
-void PGrid::tcc_pack_aMulAdd(const int nb, const double alpha, const double *a, const double *b, double *c)
+/**
+ * @brief Performs an interleaved multiply-accumulate: c = c + (alpha * a[i] * b[2i:2i+1])
+ *
+ * This function scales the components of the interleaved array 'b' by the scalar 
+ * value 'a[i]' and the global constant 'alpha', adding the result to 'c'.
+ *
+ * @param nb    The index of the band/neighboring block.
+ * @param alpha A global scaling factor.
+ * @param a     Pointer to input real scalars (Size: N).
+ * @param b     Pointer to interleaved input [Re, Im, Re, Im...] (Size: 2N).
+ * @param c     Pointer to output array for accumulation (Size: 2N).
+ *
+ * @note Use of __restrict__ allows the compiler to use SIMD vectorization.
+ */
+void PGrid::tcc_pack_aMulAdd(const int nb, const double alpha, 
+                            const double *__restrict__ a, 
+                            const double *__restrict__ b, 
+                            double *__restrict__ c)
 {
-   int i, ii;
-   int ng = nida[nb] + nidb[nb];
+   const size_t ng = nida[nb] + nidb[nb];
 
-   ii=0;
-   for (i=0; i<ng; ++i)
+   for (size_t i = 0; i < ng; ++i)
    {
-      c[ii]   += alpha*b[ii] * a[i];
-      c[ii+1] += alpha*b[ii+1]*a[i];
-      ii += 2;
+      // Calculate the effective weight once per scalar element
+      const double weight = alpha * a[i];
+      const size_t idx = 2 * i;
+
+      // Perform the accumulation
+      c[idx]     += weight * b[idx];
+      c[idx + 1] += weight * b[idx + 1];
    }
 }
 
@@ -3898,19 +3918,39 @@ void PGrid::tcc_pack_aMulAdd(const int nb, const double alpha, const double *a,
  *      PGrid:tcc_pack_iMul     *
  *                              *
  ********************************/
-void PGrid::tcc_pack_iMul(const int nb, const double *a, const double *b,
-                          double *c) {
-  int i, ii;
-  int ng = nida[nb] + nidb[nb];
+/**
+ * @brief Performs an interleaved real-to-complex expansion and multiplication.
+ *
+ *  Calculates: i * a[i] * (b_re[i] + i * b_im[i]) 
+ *            = (-a[i] * b_im[i]) + i * (a[i] * b_re[i])
+ *
+ * @param nb The index of the band/neighboring block.
+ * @param a  Pointer to input real scalars (Size: N).
+ * @param b  Pointer to interleaved input [Re, Im, Re, Im...] (Size: 2N).
+ * @param c  Pointer to output expanded components (Size: 2N).
+ *
+ * @note The arrays 'b' and 'c' must be of size at least 2*(nida[nb]+nidb[nb]).
+ */
+void PGrid::tcc_pack_iMul(const int nb, const double *__restrict__ a, const double *__restrict__ b, double *__restrict__ c) 
+{ 
+   const size_t ng = nida[nb] + nidb[nb];
 
-  ii = 0;
-  for (i = 0; i < ng; ++i) {
-    c[ii] = -b[ii + 1] * a[i];
-    c[ii + 1] = b[ii] * a[i];
-    ii += 2;
-  }
+   // Using __restrict__ (if using GCC/Clang/MSVC) tells the compiler 
+   // that the pointers do not overlap, enabling much better SIMD optimization.
+
+   for (size_t i = 0; i < ng; ++i)
+   {
+      const double scalar = a[i];
+      const size_t base_idx = 2 * i;
+      
+      // Real part of result: -a * Im(b)
+      c[base_idx]     = -scalar * b[base_idx + 1];
+      // Imaginary part of result: a * Re(b)
+      c[base_idx + 1] =  scalar * b[base_idx];
+   }
 }
 
+
 /*******************************************
  *                                         *
  *     PGrid:tcr_pack_iMul_unpack_fft      *
@@ -3956,19 +3996,37 @@ void PGrid::tcr_pack_iMul_unpack_fft(const int nb, const double *a, const double
  *     PGrid:tc_pack_iMul       *
  *                              *
  ********************************/
-void PGrid::tc_pack_iMul(const int nb, const double *a, double *c) 
+/**
+ * @brief Performs an in-place interleaved complex rotation and scaling.
+ * 
+ * Computes: c[2i : 2i+1] = c[2i : 2i+1] * (i * a[i])
+ * Which expands to:
+ *    c[2i]   = -a[i] * Im(c_old)
+ *    c[2i+1] =  a[i] * Re(c_old)
+ *
+ * @param nb The index of the band/neighboring block.
+ * @param a  Pointer to input real scalars (Size: N).
+ * @param c  Pointer to interleaved complex array to be modified in-place (Size: 2N).
+ * 
+ * @note Use of __restrict__ is critical here to allow the compiler to 
+ *       vectorize despite 'a' and 'c' being different pointers.
+ */
+void PGrid::tc_pack_iMul(const int nb, const double *__restrict__ a, double *__restrict__ c)
 {
-   double x, y;
-   int ng = nida[nb] + nidb[nb];
-   int ii = 0;
-   for (auto i=0; i<ng; ++i) 
+   const size_t ng = nida[nb] + nidb[nb];
+
+   for (size_t i = 0; i < ng; ++i)
    {
-      x = c[ii];
-      y = c[ii+1];
-     
-      c[ii]   = -y*a[i];
-      c[ii+1] =  x*a[i];
-      ii += 2;
+      const size_t idx = 2 * i;
+      const double scale = a[i];
+
+      // We must capture the original values before overwriting c[idx]
+      // This is essentially what your x and y did, but more explicit.
+      const double re = c[idx];
+      const double im = c[idx + 1];
+
+      c[idx]     = -scale * im;
+      c[idx + 1] =  scale * re;
    }
 }
 
@@ -3977,6 +4035,19 @@ void PGrid::tc_pack_iMul(const int nb, const double *a, double *c)
  *    PGrid:ttc_pack_MulSum2    *
  *                              *
  ********************************/
+/**
+ * @brief Performs weighted accumulation into interleaved array: c = c + (b * a[i])
+ * 
+ * For each index i from 0 to ng-1:
+ *    c[2i]   += b[2i]   * a[i]
+ *    c[2i+1] += b[2i+1] * a[i]
+ *
+ * @param nb The index of the band/neighboring block.
+ * @param a  Pointer to input real scalars (Size: N).
+ * @param b  Pointer to interleaved input array (Size: 2N).
+ * @param c  Pointer to output array for accumulation (Size: Must be at least 2N).
+ * 
+ */
 void PGrid::tcc_pack_MulSum2(const int nb, const double *a, const double *b, double *c) 
 {
    int ng = nida[nb] + nidb[nb];
@@ -4035,27 +4106,65 @@ double PGrid::ttt_pack_MulDot(const int nb, const double *a, const double *b, co
  *      PGrid:cc_pack_Sum2      *
  *                              *
  ********************************/
-void PGrid::cc_pack_Sum2(const int nb, const double *a, double *b) 
+/**
+ * @brief Performs vector addition in-place: b = b + a.
+ * 
+ * This function adds elements of array 'a' to array 'b'.
+ * 
+ * @param nb The index of the band/neighboring block.
+ * @param a  Pointer to input array (Size: ng).
+ * @param b  Pointer to output array (Size: ng) - modified in-place.
+ * 
+ * @note Using __restrict__ is critical here to enable SIMD vectorization.
+ */
+void PGrid::cc_pack_Sum2(const int nb, const double *__restrict__ a, double *__restrict__ b)
 {
-   int ng = 2*(nida[nb] + nidb[nb]);
+   // Use size_t to prevent overflow on large datasets
+   const size_t ng = 2 * (static_cast<size_t>(nida[nb]) + static_cast<size_t>(nidb[nb]));
 
-   for (auto i=0; i<ng; ++i)
+   // The compiler can now use AVX/SSE instructions to add 4-8 doubles at a time
+   for (size_t i = 0; i < ng; ++i)
+   {
       b[i] += a[i];
+   }
 }
 
+
 /********************************
  *                              *
  *     PGrid:cccc_pack_Sum      *
  *                              *
  ********************************/
-void PGrid::cccc_pack_Sum(const int nb, const double *a, const double *b, const double *c, double *d) 
+/**
+ * @brief Performs element-wise sum of three arrays into a fourth: d = a + b + c.
+ * 
+ * This function is a high-throughput streaming operation.
+ *
+ * @param nb The index of the band/neighboring block.
+ * @param a Pointer to input array 1 (Size: ng).
+ * @param b Pointer to input array 2 (Size: ng).
+ * @param c Pointer to input array 3 (Size: ng).
+ * @param d Pointer to output array (Size: ng) - modified in-place.
+ * 
+ * @note Using __restrict__ is critical here to allow the compiler to use SIMD.
+ */
+void PGrid::cccc_pack_Sum(const int nb, 
+                          const double *__restrict__ a, 
+                          const double *__restrict__ b, 
+                          const double *__restrict__ c, 
+                          double *__restrict__ d)
 {
-   int ng = 2*(nida[nb] + nidb[nb]);
+   // Use size_t to prevent overflow on very large grids
+   const size_t ng = 2 * (static_cast<size_t>(nida[nb]) + static_cast<size_t>(nidb[nb]));
 
-   for (auto i=0; i<ng; ++i)
-      d[i] = (a[i] + b[i] + c[i]);
+   // The compiler can now use AVX instructions to process multiple elements per cycle
+   for (size_t i = 0; i < ng; ++i)
+   {
+      d[i] = a[i] + b[i] + c[i];
+   }
 }
 
+
 /********************************
  *                              *
  *     PGrid:c_pack_addzeros    *

Nwpw/nwpwlib/psp_library/Psp1d_Hamann.cpp

@@ -1488,7 +1488,8 @@ void Psp1d_Hamann::vpp2_generate_stress_ray(Parallel *myparall,
            const double a = -sn[i]/(q*q) + rho[i]*cs[i]/q;
            f[i] = a * wp[i + indx[n + 0*5]*nrho] * vp[i + 0*nrho];
          }
-         dvnl_ray[0*lmaxnray + k1 + indx[n + 0*5]*nray] = P0 * util_simpson(nrho, f, drho) / q;
+         //dvnl_ray[0*lmaxnray + k1 + indx[n + 0*5]*nray] = P0 * util_simpson(nrho, f, drho) / q;
+         dvnl_ray[0*lmaxnray + k1 + indx[n + 0*5]*nray] = P0 * util_simpson(nrho, f, drho);
 
          /* B-term for l=0 is typically zero/unused; keep if your formulation requires it.
             Here we leave dvnl_ray[1*...] as already zeroed. */
@@ -1892,6 +1893,11 @@ void Psp1d_Hamann::vpp2_generate_stress_spline(PGrid *mygrid, int nray, double *
             const double DD = util_splint(G_ray, &(dvnlDD[ch*nray]), &(dvnlDD_splineray[ch*nray]), nray, nx, q);
 
             emit(D, DD, 1.0, 0.0, 0.0, 0.0);
+
+            //--lcount;
+            //dvnl[k + (0 + 3*lcount)*npack1] = D*ux;
+            //dvnl[k + (1 + 3*lcount)*npack1] = D*uy;
+            //dvnl[k + (2 + 3*lcount)*npack1] = D*uz;
          }
     }
   }

Nwpw/pspw/lib/psp/Pseudopotential.cpp

@@ -1544,6 +1544,9 @@ static void vpp_generate(PGrid *mygrid, char *pspname, char *fname, char *commen
          mygrid->t_read(5, tmp2, -1);
          mygrid->t_pack(1, tmp2);
          mygrid->tt_pack_copy(1, tmp2, &prj[(i+3*l) * mygrid->npack(1)]);
+
+         double sum = mygrid->tt_pack_dot(1, tmp2,tmp2);
+         std::cout << "read dvnl, i=" << i << " l=" << l << " dvnl|dvnl=" << sum << std::endl;
       }
    }
    if (*semicore) 
@@ -1744,6 +1747,10 @@ static void vpp2_write(PGrid *mygrid, char *fname2, char *comment, int psp_type,
          for (auto i=0; i<3; ++i)
          {
              mygrid->tt_pack_copy(1, &prj[(i + 3*l)*mygrid->npack(1)], tmp2);
+
+         double sum = mygrid->tt_pack_dot(1, tmp2,tmp2);
+         std::cout << "write dvnl, i=" << i << " l=" << l << " dvnl|dvnl=" << sum << std::endl;
+
              mygrid->t_unpack(1, tmp2);
              mygrid->t_write_buffer(6, tmp2, 0);
           }
@@ -1853,6 +1860,14 @@ static void vpp2_generate(// input
       psp1d.vpp2_generate_stress_spline(mygrid, nray, G_ray, dvl_ray, dvnl_ray,
                                         rho_sc_k_ray, dvl, *dvnl, *dncore);
 
+         double sum = mygrid->tt_pack_dot(1, *dvnl,*dvnl);
+         std::cout << "generate dvnl, l=" << 0 << " dvnl|dvnl=" << sum << std::endl;
+         double sum1 = mygrid->tt_pack_dot(1, *dvnl+(mygrid->npack(1)),*dvnl+(mygrid->npack(1)) );
+         std::cout << "generate dvnl, l=" << 1 << " dvnl|dvnl=" << sum1 << std::endl;
+         double sum2 = mygrid->tt_pack_dot(1, *dvnl+2*(mygrid->npack(1)),*dvnl+2*(mygrid->npack(1)) );
+         std::cout << "generate dvnl, l=" << 2 << " dvnl|dvnl=" << sum2 << std::endl;
+         double sum3 = mygrid->tt_pack_dot(1, *dvnl+3*(mygrid->npack(1)),*dvnl+3*(mygrid->npack(1)) );
+         std::cout << "generate dvnl, l=" << 3 << " dvnl|dvnl=" << sum3 << std::endl;
 
       /* deallocate ray formatted grids */
       delete[] rho_sc_k_ray;
@@ -3851,7 +3866,7 @@ void Pseudopotential::v_nonlocal_euv(const double *psi, double *stress, double *
 
 
          //**** calculate dF^(lm)_I/dhus ****
-         for (auto l=0; l < nprj[ia]; ++l)
+         for (auto l=0; l<nprj[ia]; ++l)
          {
             bool sd_function = !(l_projector[ia][l] & 1);
 
@@ -3861,21 +3876,33 @@ void Pseudopotential::v_nonlocal_euv(const double *psi, double *stress, double *
                auto dvnlprj = dvnl[ia] + s*npack1 + l*3*npack1;
                auto prj = prjtmp + (l*npack2) + u*nprj[ia]*npack2 + s*3*nprj[ia]*npack2;
 
+            double sum =  mypneb->tt_pack_dot(1, dvnlprj, dvnlprj);
+            std::cout << "ii=" << ii << " s=" << s << " l=" << l << " sum=" << sum << std::endl;
+
                mypneb->ttt_pack_Mul(1, dvnlprj, g_segments[u], tmp2);
 
+            double sum2=  mypneb->tt_pack_dot(1, tmp2, tmp2);
+            std::cout << "ii=" << ii << " s=" << s << " u=" << u << " l=" << l << " sum2=" << sum2 << std::endl;
+
                if (sd_function)
                   mypneb->tcc_pack_Mul(1, tmp2, exi, prj);
                else
                   mypneb->tcc_pack_iMul(1, tmp2, exi, prj);
 
-               //mypneb->cc_pack_ndot(1,nn psi,prj, sw3+u*nn+s*3*nn);
+               mypneb->cc_pack_ndot(1,nn,const_cast<double*>(psi),prj, sw3+u*nn+s*3*nn);
             }
-         }
-         mypneb->cc_pack_inprjdot(1, nn, 9*nprj[ia], const_cast<double*>(psi), prjtmp, sw3);
+
+         std::cout << "ii=" << ii << " sw2= ";
+         for (auto k = 0; k<nn*nprj[ia]; ++k)
+            std::cout << sw2[k] << " " ;
+         std::cout << std::endl;
+
+            std::cout << "ii=" << ii << " l=" << l << "  sw3= ";
+            for (auto k = 0; k<9*nn; ++k)
+               std::cout << sw3[k] << " " ;
+            std::cout << std::endl;
 
 
-         for (auto l=0; l<nprj[ia]; ++l)
-         {
             for (size_t i=0; i<nn; ++i)
             {
                // Pre-calculate the part that only depends on band and l
@@ -3886,8 +3913,14 @@ void Pseudopotential::v_nonlocal_euv(const double *psi, double *stress, double *
                {
                   //Bus[u+3*s] -= weight * sw3[i + u*nn + s*3*nn];
                   Bus[s+3*u] -= weight * sw3[i + u*nn + s*3*nn];
+                  std::cout << "i=" << i << "u=" << u << "s=" <<  s << " weight=" << weight << " sw3=" << sw3[i+u*nn+s*3*nn] << " bb=" << Bus[s+3*u] << std::endl;
                }
             }
+
+            std::cout << "ii=" << ii << " l=" << l << "  Bus= ";
+            for (auto k = 0; k<9; ++k)
+               std::cout << Bus[k] << " " ;
+            std::cout << std::endl;
          }