Dear all, I have now spent several days baffled by this issue so please, any help and comment is appreciated.
The code below generates NaN values when ran on two processors, that is, one processors sends garbage to MPI while the other resolves perfectly. I dont know what is going on but apparently it could be down to pointers sent to MPI? Can anyone help?

Code: Select all

                #include <iostream>
                #include <Eigen/Dense>
                #include <mpi.h>


                //using Eigen::MatrixXd;
                using namespace Eigen;
                using namespace std;

                template <class a_type> class steepest_descent
                {
                public:
                VectorXd xmin;
                double tol;
                int iter;

                //initialise
                  steepest_descent(const double toll = 0.00000003, const int interr = 100) : tol(toll), iter(interr) {}
                // I want to have a matrix A multiply it by R
                //define loopp and epsilon later
                  a_type minimise(a_type& A, a_type& v, a_type& x_0, int rankk, int nprocss);
                };

                //define methods
                template <class a_type> a_type steepest_descent<a_type>::minimise(a_type& A, a_type& v, a_type& x_0, int rankk, int nprocss)
                {


                // we want to submit each column as a start vector x_0 to each processor
                int rank = rankk;
                int nprocs = nprocss;
                int count = x_0.cols()/nprocs;


                // we use the rank of this process to work out which
                // iterations to perform.
                int start = rank * count;
                int stop = start + count;

                for (int ii = start; ii<stop; ii++) {
                    VectorXd x0;
                    x0 = x_0.col(ii);
                    VectorXd r(v.rows());
                    double alpha;
                    VectorXd xx;
                    //Vector to hold minimum values and signal
                    VectorXd sig(x0.rows()+1);

                    for (int i = start; i<iter; i++) {


                    if (i ==0) {
                        r = v - A*x0;
                        VectorXd ar = A*r;
                        alpha = r.dot(r)/r.dot(ar);
                        xx = x0 + alpha*r;
                        r = r - alpha*ar;
                    }
                    else {
                        VectorXd ar = A*r;
                        alpha = r.dot(r)/r.dot(ar);
                        xx = x0 + alpha*r;
                        r = r - alpha*ar;
                    }


                        VectorXd err = xx - x0;
                        double eff = err.cwiseAbs().maxCoeff();
                        if (eff  > tol) {
                            cout << x0(0,0) << ',' << xx(0,0) << ',' << x0(1) << ',' << xx(1) << ',' << eff << '\n';
                            x0 = xx;
             // sig.segment(0,x0.rows()) = x0;
             //sig(sig.rows()-1) = 0;
             for (int it = 0; it < sig.rows(); it++){
               if (it < (sig.rows()-1)){
            sig(it) = x0(it);
               }
               else{
            sig(it) = 0;
               }
             }
                            MPI::COMM_WORLD.Send(&sig, sig.rows(), MPI_DOUBLE, 0, 0 );
                            if (rank ==0){
                                for (int ip = 0; ip<nprocs; ip++){
                                  MPI::COMM_WORLD.Recv(&sig, sig.rows(), MPI_DOUBLE, ip, 0 );
                                  if(sig(sig.rows()-1) == 1){
                                    cout <<"Process of rank "<< ip << " found mininum at " << sig(0) << ',' << sig(1) << '\n';
                                      MPI_Abort(MPI_COMM_WORLD,911);
                                    }

                                }
                            }

                        }
                        else {
                            xmin = xx;
                            cout << xx(0,0) << ',' << xx(1,0) <<'\n';
                            //sig.segment(0,x0.rows()) = xx;
                            //sig(sig.rows()-1) = 1;
             for (int it = 0; it < sig.rows(); it++){
               if (it < (sig.rows()-1)){
            sig(it) = x0(it);
               }
               else{
            sig(it) = 1;
               }
             }
                            MPI::COMM_WORLD.Send(&sig, sig.rows(), MPI_DOUBLE, 0, 0 );
                            if (rank ==0){
                                for (int ip = 0; ip<nprocs; ip++){
                                  MPI::COMM_WORLD.Recv(sig.data(), sig.rows(), MPI_DOUBLE, ip, 0 );
                                     if(sig(sig.rows()-1) == 1){
                                       cout <<"Process of rank "<< ip << " found a mininum at " << sig(0) << ',' << sig(1) << '\n';
                                       MPI_Abort(MPI_COMM_WORLD,911);
                                     }

                                }
                            }


                        }

                    }


                }

                };

                int main(int argc, char **argv){
                MPI::Init(argc, argv);
                // get the number of processes, and the id of this process
                int rank = MPI::COMM_WORLD.Get_rank();
                int nprocs = MPI::COMM_WORLD.Get_size();


                MatrixXd com(2,1);
                com(0,0) = rank;
                com(1,0) = nprocs;
                MatrixXd A(2,2);
                A(0,0) = 3;
                A(1,0) = 2;
                A(0,1) = 2;
                A(1,1) = 6;
                MatrixXd v(2,1);
                v(0,0) = 2;
                v(1,0) = -8;
                MatrixXd x(2,2);
                x(0,0) = -2;
                x(1,0) = -2;
                x(0,1) = 4;
                x(1,1) = 3;

                steepest_descent <MatrixXd> steepest;
                steepest.minimise(A,v,x,rank,nprocs);

                MPI::Finalize();

                return 0;

                }


I am no expert at all, so please be patient with me.
I should also make clear that the version below of the code works absolutely fine, but sending eigen type vectors (or maybe any vectors?) to MPI fails spectacularly. This code is fine

Code: Select all

          #include <iostream>
          #include <Eigen/Dense>
          #include <mpi.h>


          //using Eigen::MatrixXd;
          using namespace Eigen;
          using namespace std;

          template <class a_type> class steepest_descent
          {
         public:
         VectorXd xmin;
         double tol;
         int iter;

         //initialise
            steepest_descent(const double toll = 0.00000003, const int interr = 100) : tol(toll), iter(interr) {}
         // I want to have a matrix A multiply it by R
         //define loopp and epsilon later
            a_type minimise(a_type& A, a_type& v, a_type& x_0, int rankk, int nprocss);
          };
          //define methods
template <class a_type> a_type steepest_descent<a_type>::minimise(a_type& A, a_type& v, a_type& x_0, int rankk, int nprocss)
          {


         // we want to submit each column as a start vector x_0 to each processor
         int rank = rankk;
         int nprocs = nprocss;
         int count = x_0.cols()/nprocs;


         // we use the rank of this process to work out which
         // iterations to perform.
         int start = rank * count;
         int stop = start + count;

         for (int ii = start; ii<stop; ii++) {
                 int n_0 = 0;
                 VectorXd x0;
            x0 = x_0.col(ii);
            VectorXd r(v.rows());
            double alpha;
            VectorXd xx;
            //Vector to hold minimum values and signal

            for (int i = start; i<iter; i++) {


            if (i ==0) {
                r = v - A*x0;
                VectorXd ar = A*r;
                alpha = r.dot(r)/r.dot(ar);
                xx = x0 + alpha*r;
                r = r - alpha*ar;
            }
            else {
               VectorXd ar = A*r;
               alpha = r.dot(r)/r.dot(ar);
               xx = x0 + alpha*r;
               r = r - alpha*ar;
            }


               VectorXd err = xx - x0;
               double eff = err.cwiseAbs().maxCoeff();
               if (eff  > tol) {
                  cout << x0(0,0) << ',' << xx(0,0) << ',' << x0(1) << ',' << xx(1) << ',' << eff << '\n';
                  x0 = xx;
                  MPI::COMM_WORLD.Send(&n_0, 1, MPI_INT, 0, 0 );
                  if (rank ==0){
                     for (int ip = 0; ip<nprocs; ip++){
                        MPI::COMM_WORLD.Recv(&n_0, 1, MPI_INT, ip, 0 );
                        if(n_0 == 1){
                          cout <<"Process of rank "<< ip << " caused an abort.\n";
                          MPI_Abort(MPI_COMM_WORLD,911);
                        }

                     }
                  }

               }
               else {
                   xmin = x0;
                    cout << x0(0,0) << ',' << x0(1,0) <<'\n';
                       n_0 = 1;
                  MPI::COMM_WORLD.Send( &n_0, 1, MPI_INT, 0, 0 );
                  if (rank ==0){
                     for (int ip = 0; ip<nprocs; ip++){
                       MPI::COMM_WORLD.Recv( &n_0, 1, MPI_INT, ip, 0 );
                         if(n_0 == 1){
                           cout <<"Process of rank "<< ip << " caused an abort.\n";
                           MPI_Abort(MPI_COMM_WORLD,911);
                         }

                     }
                  }


               }

            }


         }

          };

          int main(int argc, char **argv){
         MPI::Init(argc, argv);
         // get the number of processes, and the id of this process
         int rank = MPI::COMM_WORLD.Get_rank();
         int nprocs = MPI::COMM_WORLD.Get_size();


         MatrixXd com(2,1);
         com(0,0) = rank;
         com(1,0) = nprocs;
         MatrixXd A(2,2);
         A(0,0) = 3;
         A(1,0) = 2;
         A(0,1) = 2;
         A(1,1) = 6;
         MatrixXd v(2,1);
         v(0,0) = 2;
         v(1,0) = -8;
         MatrixXd x(2,2);
         x(0,0) = -2;
         x(1,0) = -2;
         x(0,1) = 4;
         x(1,1) = 3;

         steepest_descent <MatrixXd> steepest;
         steepest.minimise(A,v,x,rank,nprocs);

         MPI::Finalize();

         return 0;

          }


I know there is a 'sig.data()' instead of the reference '&sig' somewhere there, so feel free to change it to that when playing with the code, but you still get the NaN values.

Apparently this is the issue, does anyone know how to fix it

https://lists.sdsc.edu/pipermail/npaci- ... 38886.html


Please, also note that the start vector [4,3] should descend faster to the true answer (minimum), which is why I know for sure there is a problem. If I run the code by sending integers to MPI, the start vector [4,3] always causes the abort.