Two comments:
1) It's pretty easy to use MPI from Julia. For some use cases it may make
more sense than Julia's built-in approach to parallelism, especially if
you're already comfortable with MPI. Though if you're well served by pmap,
that's much simpler.
2) It looks like the subproblem you're solving is a 1d optimization. I
wouldn't be surprised if you could get a significant speedup by tweaking
the optimization algorithm (e.g., using a derivative-based method) and
making the function evaluations more efficient (e.g., avoiding temporary
allocations).
On Sunday, November 23, 2014 11:02:20 PM UTC-5, Gabriel Mihalache wrote:
>
> I made some progress on getting the code to run correctly in parallel
> while using all the cores of one machine, by using pmap and running some
> prep code on all the workers:
> The next steps are to get this running on multiple machines (using
> addprocs_ssh?) and then benchmarking. Here's the state of my solution so
> far, for anyone interested/stuck at this too:
>
> The SLURM batch file:
>
> #!/bin/bash
> #SBATCH -J sgmVfi
>
> #SBATCH -e elog_sgm
> #SBATCH -o olog_sgm
>
> #SBATCH -t 0:10:00
>
> #SBATCH --nodes=1
> #SBATCH --ntasks-per-node=1
> #SBATCH --cpus-per-task=24
>
> #SBATCH -p standard
> #SBATCH --mail-type=fail
>
> srun julia -q -p 24 -L loader.jl sgm.jl
>
>
> the "loader" file which defines constants for all workers:
>
> using MAT;
> using Grid;
> using NLopt;
>
> const alpha = 0.25;
> const uBeta = 0.90;
> const delta = 1.00;
>
> const kGridMin = 0.01;
> const kGridMax = 0.50;
> const kGridSize = 100;
> const kGrid = linrange(kGridMin, kGridMax, kGridSize);
>
> const vErrTol = 1e-8;
> const maxIter = 200;
>
> file = matopen("shock.mat")
> yData = read(file, "y");
> yGridSize = size(yData, 1);
> const y = linrange(yData[1], yData[end], yGridSize);
> const yPi = read(file, "yPi")
> close(file)
>
> gdp(yIx, kIx) = exp(y[yIx]) * kGrid[kIx].^alpha;
>
>
> and finally, the code running on the master worker/driver:
>
> function sgmVfi()
> V0 = zeros(yGridSize, kGridSize);
> V1 = zeros(size(V0));
> kPr = zeros(size(V0));
>
> iter = 0;
> err = 1.0;
> while iter < maxIter && err > vErrTol
> tic();
> iter += 1;
> println("Iteration ", iter);
> V0 = copy(V1);
> V0f = CoordInterpGrid( (y, kGrid), V0, BCnearest, InterpQuadratic);
>
> function workAt(iix)
> yIx = iix[1];
> kIx = iix[2];
> opt = Opt(:LN_BOBYQA, 1)
> lower_bounds!(opt, kGridMin)
> upper_bounds!(opt, min(kGridMax, gdp(yIx, kIx) + (1-delta) * kGrid[kIx] -
> 1e-5))
>
> myObj(kk) = log(gdp(yIx, kIx) - kk + (1-delta) * kGrid[kIx]) + uBeta *
> sum( [ yPi[yIx,yPrIx] * V0f[y[yPrIx], kk] for yPrIx = 1:yGridSize ]);
>
> max_objective!(opt, (x::Vector, grad::Vector) -> myObj(x[1]) );
> (maxVal, maxK, ret) = optimize!(opt, [kGridMin] );
> return maxVal, maxK[1];
> end
>
> tmp = pmap(workAt, [[yIx, kIx] for yIx = 1:yGridSize, kIx = 1:kGridSize]);
> tmp = reshape(tmp, size(V1));
> for yIx = 1:yGridSize, kIx = 1:kGridSize
> V1[yIx, kIx] = tmp[yIx, kIx][1];
> kPr[yIx, kIx] = tmp[yIx, kIx][2];
> end
>
> err = maximum(abs(V1[:] - V0[:]));
> toc()
> println("Err = ", err)
> end
>
> return kGrid, V1, kPr;
> end
>
> kGrid, V, kPr = sgmVfi()
>
> file = matopen("results.mat", "w")
> write(file, "kPr", kPr)
> write(file, "V", V)
> write(file, "y", [y])
> write(file, "yPi", yPi)
> write(file, "kGrid", [kGrid])
> close(file)
>
>
>
>
>
