Issues with MPI in PPPM-like computation

[my-tan]

I am making a code that calculates hydrodynamic interactions between particles with Rotne-Prager matrix (I will call it M here). The velocity of particle is calculated as U = M \cdot F + M^0.5 \cdot P, where U is velocity (unknown), F is force (known), and P is a random vector. The computation is split into two parts by using Ewald summation. The wave-part computation mostly follows the procedures of PPPMForceCompute. The simulation goes well when using single CPU, but fails with multiple CPUs by using MPI.

The issue most likely happens in the computation M^0.5 \cdot P. I have two major problems.

1. P is a random vector. I can generate a random "force" of each particle, and then follow the procedures assignParticles --> updateMeshes --> interpolateForces. To save a FFT time, or I can generate random numbers on half of the mesh, and the value of the other half is automatically set because of the conjugacy after taking FFT of a real-valued function. It is straightforward to do with one CPU, but with multiple CPUs, the mesh and its conjugate mesh are not on the same CPU. When I assign a random number to a grid on a CPU, how do I assign its conjugate value of a grid on another CPU?
2. The second issue is the computation of M^0.5 \cdot P in the real space. I use Lanczos algorithm to calculate the square root of a matrix. In Lanczos algorithm, I need to construct a tridiagonal symmetric matrix T that has alpha's on its main diagonal and beta's on its upper/lower diagonals, where alpha = v_{j} \cdot (M \cdot v_{j}) and beta = norm(v_{j}), v_{j} is calculated iteratively by doing matrix-vector multiplication with v_{0} = P. Since each CPU deals with a piece of M and P, I need to sum up alpha and beta from all CPUs by using MPI_Allreduce. But it does not seem to work. At some point, the simulation just gets stuck. I am not sure what goes wrong.

I copy and paste the Lanczos code here. The goal is to find "u", and "psi" is the random vector. "mobilityRealUF(vj, v)" does the matrix-vector multiplication.

`
void TwoStepRPY::brownianLanczos(Scalar * psi,
Scalar * iter_ff_v,
Scalar * iter_ff_vj,
Scalar * iter_ff_vjm1,
Scalar * iter_ff_V,
Scalar * iter_ff_uold,
Scalar * u)
{
unsigned int group_size = m_group->getNumMembers();
unsigned int numel = group_size * 6;

unsigned int m_in = m_m_lanczos_ff;
unsigned int mmax = m_mmax;

Scalar * v = iter_ff_v;
Scalar * V = iter_ff_V;
Scalar * vj = &V[numel];
Scalar * vjm1 = V;
Scalar * uold = iter_ff_uold;

double * alpha = new double[mmax];
double * alpha_save = new double[mmax];
double * beta = new double[mmax + 1];
double * beta_save = new double[mmax + 1];
double * W = new double[mmax * mmax];
double * W1 = new double[mmax];
double * Tm = new double[mmax];

// copy psi to v0
for (unsigned int i = 0; i < numel; i++)
    {
    vj[i] = psi[i];
    }
Scalar vnorm, psinorm;
vnorm = cblas_ddot(numel,
                   vj, 1,
                   vj, 1);

#ifdef ENABLE_MPI
if (m_pdata->getDomainDecomposition())
{
MPI_Barrier(m_exec_conf->getMPICommunicator());
MPI_Allreduce(MPI_IN_PLACE,
&vnorm,
1,
MPI_FLOAT,
MPI_SUM,
m_exec_conf->getMPICommunicator());
}
#endif
vnorm = sqrt(vnorm);
psinorm = vnorm;

unsigned int m = m_in - 1;
m = m < 1 ? 1 : m;

Scalar alpha_temp;
Scalar beta_temp = 0.0;

// first iteration
// scale v_{j-1} by 0 and v_{j} by 1 / vnorm
Scalar scale = 0.0;
cblas_dscal(numel,
            scale,
            vjm1, 1);
scale = 1.0 / vnorm;
cblas_dscal(numel,
            scale,
            vj, 1);

beta[0] = beta_temp;
for (unsigned int i = 0; i < numel; i++)
    {
    V[i] = vj[i];
    }

for (unsigned int j = 0; j < m; j++)
    {
    // v = M . v_{j} - beta_{j} * v_{j - 1}
    mobilityRealUF(vj,
                   v);
    scale = - beta_temp;
    cblas_daxpy(numel,
                scale,
                vjm1, 1,
                v, 1);

    // alpha_{j} = v_{j} \cdot v
    alpha_temp = cblas_ddot(numel,
                            vj, 1,
                            v, 1);

#ifdef ENABLE_MPI
if (m_pdata->getDomainDecomposition())
{
MPI_Barrier(m_exec_conf->getMPICommunicator());
MPI_Allreduce(MPI_IN_PLACE,
&alpha_temp,
1,
MPI_FLOAT,
MPI_SUM,
m_exec_conf->getMPICommunicator());
}
#endif
// store alpha_{j}
alpha[j] = alpha_temp;

    // v = v - alpha_{j} * v_{j}
    scale = - alpha_temp;
    cblas_daxpy(numel,
                scale,
                vj, 1,
                v, 1);

    // beta_{j + 1} = v \cdot v
    vnorm = cblas_ddot(numel,
                       v, 1,
                       v, 1);

#ifdef ENABLE_MPI
if (m_pdata->getDomainDecomposition())
{
MPI_Barrier(m_exec_conf->getMPICommunicator());
MPI_Allreduce(MPI_IN_PLACE,
&vnorm,
1,
MPI_FLOAT,
MPI_SUM,
m_exec_conf->getMPICommunicator());
}
#endif
vnorm = sqrt(vnorm);
beta_temp = vnorm;

    // store beta
    beta[j + 1] = beta_temp;
    
    if (beta_temp < 1e-8)
        {

#ifdef ENABLE_MPI
MPI_Barrier(m_exec_conf->getMPICommunicator());
#endif
m = j + 1;
break;
}

    // v = v / beta_{j + 1} 
    scale = 1.0 / beta_temp;
    cblas_dscal(numel,
                scale,
                v, 1);

    // store current basis vector V
    for (unsigned int i = 0; i < numel; i++)
        {
        V[(j + 2) * numel + i] = v[i];
        }
    
    // swap pointers
    vjm1 = &V[(j + 1) * numel];
    vj = &V[(j + 2) * numel];

#ifdef ENABLE_MPI
MPI_Barrier(m_exec_conf->getMPICommunicator());
#endif
} // for j

// save alpha and beta (will be overwritten by LAPACK)
for (unsigned int i = 0; i < m; i++)
    {
    alpha_save[i] = alpha[i];
    beta_save[i] = beta[i];
    }
beta_save[m] = beta[m];

// compute eigen-decomposition of tridiagonal matrix
int INFO = LAPACKE_dpteqr(LAPACK_ROW_MAJOR, 'I', m, alpha, &beta[1], W, m);
if (INFO != 0)
    {
    std::cout << "Eigenvalue decomposition #1 failed." << std::endl;
    std::cout << "INFO = " << INFO << std::endl;
    std::cout << "m = " << m << std::endl;
        std::cout << "alpha: " << std::endl;
        for (unsigned int k = 0; k < m; k++)
            {
            std::cout << alpha[k] << std::endl;
            }
        std::cout << std::endl;

        std::cout << "beta: " << std::endl;
        for (unsigned int k = 0; k < m; k++)
            {
            std::cout << beta[k] << std::endl;
            }
        std::cout << std::endl;
    exit(EXIT_FAILURE); 
    }
for (unsigned int i = 0; i < m; i++)
    {
    W1[i] = sqrt(alpha[i]) * W[i];
    }
double sum_temp;
for (unsigned int i = 0; i < m; i++)
    {
    sum_temp = 0.0;
    for (unsigned int k = 0; k < m; k++)
        {
        unsigned int idx = i * m + k;
        sum_temp += W[idx] * W1[k];
        } // for j

        Tm[i] = sum_temp;
    }

// multiply basis vectors by Tm to get x = A \cdot Tm
lanczosMapBack(V,
               Tm,
               m,
               group_size,
               6,
               u);
for (unsigned int i = 0; i < numel; i++)
    {
    uold[i] = u[i];
    }
// recover alpha and beta
for (unsigned int i = 0; i < m; i++)
    {
    alpha[i] = alpha_save[i];
    beta[i] = beta_save[i];
    }
beta[m] = beta_save[m];

std::cout << "Initial Lanczos iteration is done" << std::endl;
// keep adding to basis vectors untill step norm is small enough
Scalar stepnorm = 1.0;
unsigned int j;
while (stepnorm > m_error && m < mmax)
    {
    m++;
    j = m - 1;

    // v = M . v_{j} - beta_{j} * v_{j - 1}
    mobilityRealUF(vj,
                   v);
    scale = - beta_temp;
    cblas_daxpy(numel,
                scale,
                vjm1, 1,
                v, 1);

    // alpha_{j} = v_{j} . v
    alpha_temp = cblas_ddot(numel,
                            vj, 1,
                            v, 1);

#ifdef ENABLE_MPI
if (m_pdata->getDomainDecomposition())
{
MPI_Barrier(m_exec_conf->getMPICommunicator());
MPI_Allreduce(MPI_IN_PLACE,
&alpha_temp,
1,
MPI_FLOAT,
MPI_SUM,
m_exec_conf->getMPICommunicator());
}
#endif
// store alpha_{j}
alpha[j] = alpha_temp;

    // v = v - alpha_{j} * v_{j}
    scale = - alpha_temp;
    cblas_daxpy(numel,
                scale,
                vj, 1,
                v, 1);

    // beta_{j + 1} = v \cdot v
    vnorm = cblas_ddot(numel,
                       v, 1,
                       v, 1);

#ifdef ENABLE_MPI
if (m_pdata->getDomainDecomposition())
{
MPI_Barrier(m_exec_conf->getMPICommunicator());
MPI_Allreduce(MPI_IN_PLACE,
&vnorm,
1,
MPI_FLOAT,
MPI_SUM,
m_exec_conf->getMPICommunicator());
}
#endif
vnorm = sqrt(vnorm);
beta_temp = vnorm;

    // store beta
    beta[j + 1] = beta_temp;
    
    if (beta_temp < 1e-8)
        {

#ifdef ENABLE_MPI
MPI_Barrier(m_exec_conf->getMPICommunicator());
#endif
m = j + 1;
break;
}

    // v = v / beta_{j + 1} 
    scale = 1.0 / beta_temp;
    cblas_dscal(numel,
                scale,
                v, 1);

    // store current basis vector V
    for (unsigned int i = 0; i < numel; i++)
        {
        V[(j + 2) * numel + i] = v[i];
        }
    
    // swap pointers        
    vjm1 = &V[(j + 1) * numel];
    vj = &V[(j + 2) * numel];

    // save alpha and beta (will be overwritten by LAPACK)
    for (unsigned int i = 0; i < m; i++)
        {
        alpha_save[i] = alpha[i];
        beta_save[i] = beta[i];
        }
    beta_save[m] = beta[m];

    // compute eigen-decomposition of tridiagonal matrix
    int INFO = LAPACKE_dpteqr(LAPACK_ROW_MAJOR, 'I', m, alpha, &beta[1], W, m);
    if (INFO != 0)
        {
        std::cout << "Eigenvalue decomposition #2 failed." << std::endl;
        std::cout << "INFO = " << INFO << std::endl;
        std::cout << "m = " << m << std::endl;
        std::cout << "alpha: " << std::endl;
        for (unsigned int k = 0; k < m; k++)
            {
            std::cout << alpha[k] << std::endl;
            }
        std::cout << std::endl;

        std::cout << "beta: " << std::endl;
        for (unsigned int k = 0; k < m; k++)
            {
            std::cout << beta[k] << std::endl;
            }
        std::cout << std::endl;
        exit(EXIT_FAILURE); 
        }
    for (unsigned int i = 0; i < m; i++)
        {
        W1[i] = sqrt(alpha[i]) * W[i];
        }
    
    for (unsigned int i = 0; i < m; i++)
        {
        sum_temp = 0.0;
        for (unsigned int k = 0; k < m; k++)
            {
            unsigned int idx = i * m + k;
            sum_temp += W[idx] * W1[k];
            } // for j

            Tm[i] = sum_temp;
        }
    // multiply basis vectors by Tm to get x = A \cdot Tm
    lanczosMapBack(V,
                   Tm,
                   m,
                   group_size,
                   6,
                   u);

    // compute step norm
    scale = - 1.0;
    cblas_daxpy(numel,
                scale,
                u, 1,
                uold, 1);
    stepnorm = cblas_ddot(numel,
                          uold, 1,
                          uold, 1);

#ifdef ENABLE_MPI
if (m_pdata->getDomainDecomposition())
{
MPI_Barrier(m_exec_conf->getMPICommunicator());
MPI_Allreduce(MPI_IN_PLACE,
&stepnorm,
1,
MPI_FLOAT,
MPI_SUM,
m_exec_conf->getMPICommunicator());
}
#endif
stepnorm = sqrt(stepnorm);

    Scalar fnorm = cblas_ddot(numel,
                              u, 1,
                              u, 1);

#ifdef ENABLE_MPI
if (m_pdata->getDomainDecomposition())
{
MPI_Barrier(m_exec_conf->getMPICommunicator());
MPI_Allreduce(MPI_IN_PLACE,
&fnorm,
1,
MPI_FLOAT,
MPI_SUM,
m_exec_conf->getMPICommunicator());
}
#endif
fnorm = sqrt(fnorm);
stepnorm = stepnorm / fnorm;

    for (unsigned int i = 0; i < numel; i++)
        {
        uold[i] = u[i];
        }
    // recover alpha and beta
    for (unsigned int i = 0; i < m; i++)
        {
        alpha[i] = alpha_save[i];
        beta[i] = beta_save[i];
        }
        beta[m] = beta_save[m];

#ifdef ENABLE_MPI
MPI_Barrier(m_exec_conf->getMPICommunicator());
#endif
} // while

std::cout << "Lanczos iteration number = " << m << std::endl;

// scale by the norm of psi
cblas_dscal(numel,
            psinorm,
            u, 1); 

v = NULL;
vj = NULL;
vjm1 = NULL;
V = NULL;
uold = NULL;
delete [] alpha;
delete [] alpha_save;
delete [] beta;
delete [] beta_save;
delete [] W;
delete [] W1;
delete [] Tm;
} /* end of brownianLanczos */

`

[joaander]

When I assign a random number to a grid on a CPU, how do I assign its conjugate value of a grid on another CPU?

Use the same seed. See

hoomd-blue/hoomd/RandomNumbers.h

Lines 115 to 131 in daf7303

	class Seed
	{
	public:
	/** Construct a Seed from a class ID, timestep and user seed.
	The seed is 8 bytes. Construct this from an 1 byte class id, 2 byte seed, and the lower
	5 bytes of the timestep.
	When multiple class instances can instantiate RandomGenerator objects, include values in the
	Counter that are unique to each instance. Otherwise the separate instances will generate
	identical sequences of random numbers.
	Code inside HOOMD should name the class ID value in RNGIdentifiers.h and ensure that each
	class has a unique id. External plugins should use values 200 or larger.
	id seed1 seed0 timestep4 | timestep3 timestep2 timestep1 timestep0
	*/

But it does not seem to work. At some point, the simulation just gets stuck. I am not sure what goes wrong.

I have no idea.

[my-tan]

Thank you for the reply. I also have questions about the domain decomposition in FFT.

1. When the charges are assigned to the mesh, the index is off compared to one-processor situation (probably due to the ghost layer). I would appreciate it if you could explain how the assignment of charge is done in multiple-processor situation, or point me to the materials that I could learn.
2. In my case, I need to calculate the distance between a particle and its nearby mesh. I used "rx = (Scalar(i) + dx) * m_h.x", where "rx" the distance between the particle and a mesh, "i" and "dx" follow the variables in PPPMForceCompute::assignParticles(), and "m_h.x" is the distance between to adjacent mesh. Did I do this correctly?

[joaander]

I usually use a serial debugger to determine the cause of a MPI deadlock. OpenMPI documents the process here: https://docs.open-mpi.org/en/main/app-debug/serial-debug.html

@jglaser wrote the MPI implementation of PPPMForceCompute. He may be able to answer your specific questions about that code.