adding in new_pardiso solver, that should be fast

#ifndef __RANDOM_H__

#define __RANDOM_H__

typedef unsigned long int U_Int_32_t;

typedef unsigned long long int  U_Int_64_t;

class Random {

public:

    U_Int_64_t u,v,w;

    Random(U_Int_64_t j) : v(4101842887655102017LL), w(1) {

        u = j ^ v; int64();

        v = u; int64();

        w = v; int64();

    }

    inline U_Int_64_t int64() {

        u = u * 2862933555777941757LL + 7046029254386353087LL;

        v ^= v >> 17; v ^= v << 31; v ^= v >> 8;

        w = 4294957665U*(w & 0xffffffff) + (w >> 32);

        U_Int_64_t x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;

        return (x + v) ^ w;

    }

    inline double doub() { return 5.42101086242752217E-20 * int64(); }

    inline unsigned int int32() { return (unsigned int)int64(); }

};

#endif

mpiicpc -cxx=icpx -qmkl=parallel -O3 -std=c++17 main_mpipardiso.cpp -o i_mpipardiso -I ./ -I ~/

icpx -qmkl=parallel -O3 -std=c++17 main_pardiso.cpp -o i_pardiso -I ./ -I ~/

#include <cstdint>

#include <iostream>

#include <sstream>

#include <mpi.h>

#include <Random.h>

#include <lattice3d.h>

// #include <rfm.h>

#include <rfm_pardiso.h>

#include <algorithm>

#include <vector>

#include <array>

#include <list>

#include <Eigen/Core>

#include <Eigen/Sparse>

using namespace std;

using namespace Eigen;

template<typename T>

void get_from_command_line(int argc, char **argv, string name, T &a) {

    string value;

    for(int i=1; i<argc; ++i) {

        string argument(argv[i]);

        size_t found = argument.find(name);

        if(found!=string::npos) {

            size_t eqsgn = argument.find("=");

            value = string(argument,eqsgn+1,argument.size());

            break;

        }

    }

    stringstream vss(value);

    vss >> a;

}

double wei(double r, double b) {

            double rb = 1.0/b;

            double rg = 1.0/tgamma(1.0+rb);

            double rl = pow(log(1.0/r), rb);

            return  rg*rl;

        }

void export_sel(const char* filename, vector<uint32_t> &sel, vector<double> &ths) {

        uint32_t nedges0 = sel.size()/2;

        ofstream f;

        f.precision(17);

        f.open(filename);

        for(int e=0; e<nedges0; ++e) {

            f << sel[2*e] << " " << sel[2*e+1] << " " << ths[e] << "\n";

        }

        f.close();

}

int main(int argc, char **argv)

{

    MPI_Init(NULL,NULL);

    // Get the number of processes

    int comm_size = 1;

    MPI_Comm_size(MPI_COMM_WORLD, &comm_size);

    // Get the rank of the process

    int comm_rank = 1;

    MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);

    uint64_t seed = 1;

    uint32_t steps = 4;

    char type = 's';

    uint32_t nremove = 0;

    double weik = 4.0;

    uint32_t precr = 0;

    uint32_t niter = 1;

    uint32_t networks = 1;

    bool optional = true; //output data for individual realizations

    bool storeiv = true;

    get_from_command_line<>(argc,argv,"seed",seed);

    get_from_command_line<>(argc,argv,"steps",steps);

    get_from_command_line<>(argc,argv,"type",type);

    get_from_command_line<>(argc,argv,"nremove",nremove);

    get_from_command_line<double>(argc,argv,"weik",weik);

    get_from_command_line<>(argc,argv,"niter",niter);

    get_from_command_line<>(argc,argv,"networks",networks);

    get_from_command_line<>(argc,argv,"optional",optional);

    get_from_command_line<>(argc,argv,"storeiv",storeiv);

    get_from_command_line<>(argc,argv,"precr",precr);

    seed = seed*1000ULL + comm_rank;

    cout << "seed = " << seed << endl;

    Random ra(seed);

    int realizations = networks;

    //if(realizations>1) optional = false;

    for(int r=0; r<realizations; ++r) {

        Lattice lt(ra,steps,type,nremove,weik,precr);

        Lattice lb(ra,steps,'r',nremove,weik,0);

        Lattice l(ra,lb,lt);

        auto sel = l.sorted_edgelist<uint32_t>();

        vector<double> ths(l.edges0);

        for(auto& t: ths) {t=wei(ra.doub(), weik);}

        stringstream e_stream;

        e_stream << "edgelist_0_seed" << seed << "_" << r;

        export_sel(e_stream.str().c_str(),sel,ths);

        RFM rfm(l.N,l.NxNy,sel,ths);

        rfm.run(niter);

        if(storeiv) {

            stringstream namestream;

            namestream << "iv_seed" << seed << "_" << r;

            rfm.export_iv(namestream.str().c_str());

        }

        stringstream bs_stream;

        bs_stream << "brokenlist_sorted_seed" << seed << "_" << r;

        rfm.export_sbl(bs_stream.str().c_str());

    }

    MPI_Finalize();

}