diff --git a/README.md b/README.md
index f1ea851..ecc5366 100644
--- a/README.md
+++ b/README.md
@@ -13,29 +13,34 @@ Pkg.clone("https://github.com/kibaekkim/Dsp.jl")
 DSP can read and solve model from JuMP:
 
 ```julia
-using Dsp, JuMP
+using JuMP, Dsp, MPI
+
+# Comment out this line if you want to run in serial
+MPI.Init()
 
 xi = [[7,7] [11,11] [13,13]]
 
-m = Model()
+# create JuMP.Model with number of blocks
+m = Model(3)
+
 @variable(m, 0 <= x[i=1:2] <= 5, Int)
 @objective(m, Min, -1.5 * x[1] - 4 * x[2])
 for s = 1:3
-    blk = Model()
+    # create a JuMP.Model block linked to m with id s and probability 1/3
+    blk = Model(m, s, 1/3)
     @variable(blk, y[j=1:4], Bin)
     @objective(blk, Min, -16 * y[1] + 19 * y[2] + 23 * y[3] + 28 * y[4])
     @constraint(blk, 2 * y[1] + 3 * y[2] + 4 * y[3] + 5 * y[4] <= xi[1,s] - x[1])
     @constraint(blk, 6 * y[1] + y[2] + 3 * y[3] + 2 * y[4] <= xi[2,s] - x[2])
-
-    # Model m has block blk indexed by s with objective function weight of 1/3.
-    @block(m, blk, s, 1/3)
-
 end
 
 solve_types = [:Dual, :Benders, :Extensive]
 solve(m, solve_type = solve_types[1], param = "myparam.txt")
 
 getobjectivevalue(m)
+
+# Comment out this line if you want to run in serial
+MPI.Finalize()
 ```
 
 or, it can also read and solve model from SMPS files:
diff --git a/examples/farmer.jl b/examples/farmer.jl
index 119106d..770eda9 100644
--- a/examples/farmer.jl
+++ b/examples/farmer.jl
@@ -20,15 +20,14 @@ Yield    = [3.0 3.6 24.0;
 Minreq   = [200 240 0]     # minimum crop requirement
 
 # JuMP model
-# m = Model(solver = DspSolver(method = :DD, param = "myparams.txt"))
-m = Model()
+m = Model(NS)
 
 @variable(m, x[i=CROPS] >= 0, Int)
 @objective(m, Min, sum{Cost[i] * x[i], i=CROPS})
 @constraint(m, const_budget, sum{x[i], i=CROPS} <= Budget)
 
 for s in 1:NS
-    blk = Model()
+    blk = Model(m, s, probability[s])
 
     @variable(blk, y[j=PURCH] >= 0)
     @variable(blk, w[k=SELL] >= 0)
@@ -38,7 +37,4 @@ for s in 1:NS
     @constraint(blk, const_minreq[j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
     @constraint(blk, const_minreq_beets, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
     @constraint(blk, const_aux, w[3] <= 6000)
-
-    # add blk to m as a block
-    @block(m, blk, s, probability[s])
 end
diff --git a/examples/farmer_mpi.jl b/examples/farmer_mpi.jl
new file mode 100644
index 0000000..ea2ebe4
--- /dev/null
+++ b/examples/farmer_mpi.jl
@@ -0,0 +1,46 @@
+# Kibaek Kim - ANL MCS 2016
+# Farmer example from Birge and Louveaux book.
+#
+# To use MPI, add MPI.Init() and MPI.Finalize() around the scritp lines.
+
+using MPI, JuMP, Dsp
+
+MPI.Init()
+
+NS = 3;                        # number of scenarios
+probability = [1/3, 1/3, 1/3]; # probability
+
+CROPS  = 1:3 # set of crops (wheat, corn and sugar beets, resp.)
+PURCH  = 1:2 # set of crops to purchase (wheat and corn, resp.)
+SELL   = 1:4 # set of crops to sell (wheat, corn, sugar beets under 6K and those over 6K)
+
+Cost     = [150 230 260]   # cost of planting crops
+Budget   = 500             # budget capacity
+Purchase = [238 210];      # purchase price
+Sell     = [170 150 36 10] # selling price
+Yield    = [3.0 3.6 24.0;
+            2.5 3.0 20.0;
+            2.0 2.4 16.0]
+Minreq   = [200 240 0]     # minimum crop requirement
+
+# JuMP model
+m = Model(NS)
+
+@variable(m, x[i=CROPS] >= 0, Int)
+@objective(m, Min, sum{Cost[i] * x[i], i=CROPS})
+@constraint(m, const_budget, sum{x[i], i=CROPS} <= Budget)
+
+for s in 1:NS
+    blk = Model(m, s, probability[s])
+
+    @variable(blk, y[j=PURCH] >= 0)
+    @variable(blk, w[k=SELL] >= 0)
+
+    @objective(blk, Min, sum{Purchase[j] * y[j], j=PURCH} - sum{Sell[k] * w[k], k=SELL})
+
+    @constraint(blk, const_minreq[j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
+    @constraint(blk, const_minreq_beets, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
+    @constraint(blk, const_aux, w[3] <= 6000)
+end
+
+MPI.Finalize()
\ No newline at end of file
diff --git a/src/Dsp.jl b/src/Dsp.jl
index a7685e1..225f8fd 100644
--- a/src/Dsp.jl
+++ b/src/Dsp.jl
@@ -1,4 +1,4 @@
-__precompile__()
+# __precompile__()
 
 module Dsp
 
@@ -8,21 +8,60 @@ include("block.jl")
 using Dsp.DspCInterface
 
 import JuMP
-export readSmps, getblocksolution, getdualobjval
+export
+    readSmps, 
+    getblocksolution, 
+    optimize, 
+    getprimobjval, 
+    getdualobjval, 
+    getprimvalue,
+    getdualvalue,
+    getsolutiontime, 
+    blockids
 
 # DspModel placeholder
 model = DspModel()
 
 ###############################################################################
-# JuMP.solvehook
+# Override JuMP.Model
 ###############################################################################
 
-# This function is hooked by JuMP (see block.jl)
-function dsp_solve(m::JuMP.Model; suppress_warnings = false, options...)
-
+# This is for the master problem.
+function JuMP.Model(nblocks::Integer)
+    # construct model
+    m = JuMP.Model()
     # free DspModel
     DspCInterface.freeModel(Dsp.model)
+    # set block ids
+    DspCInterface.setBlockIds(Dsp.model, nblocks)
+    # set extension
+    m.ext[:DspBlocks] = BlockStructure(nothing, Dict{Int,JuMP.Model}(), Dict{Int,Float64}())
+    # set solvehook
+    JuMP.setsolvehook(m, Dsp.dsp_solve)
+    # return model
+    m
+end
+
+# This is for the subproblems.
+function JuMP.Model(parent::JuMP.Model, id::Integer, weight::Float64)
+    # construct model
+    m = JuMP.Model()
+    # set extension
+    m.ext[:DspBlocks] = BlockStructure(nothing, Dict{Int,JuMP.Model}(), Dict{Int,Float64}())
+    # set block
+    m.ext[:DspBlocks].parent = parent
+    setindex!(parent.ext[:DspBlocks].children, m, id)
+    setindex!(parent.ext[:DspBlocks].weight, weight, id)
+    # return model
+    m
+end
+
+###############################################################################
+# JuMP.solvehook
+###############################################################################
 
+# This function is hooked by JuMP (see block.jl)
+function dsp_solve(m::JuMP.Model; suppress_warnings = false, options...)
     # parse options
     for (optname, optval) in options
         if optname == :param
@@ -60,37 +99,11 @@ function dsp_solve(m::JuMP.Model; suppress_warnings = false, options...)
     end
 
     if !(stat == :Infeasible || stat == :Unbounded)
-        try
-            getDspSolution(suppress_warnings)
-
-            # Don't corrupt the answers if one of the above two calls fails
-            m.objVal = Dsp.model.primVal
-            # parse solution to each block
-            n_start = 1
-            n_end = m.numCols
-            m.colVal = copy(Dsp.model.colVal[n_start:n_end])
-            n_start += m.numCols
-            if haskey(m.ext, :DspBlocks) == true
-                blocks = m.ext[:DspBlocks].children
-                for s in keys(blocks)
-                    n_end += blocks[s].numCols
-                    blocks[s].colVal = Dsp.model.colVal[n_start:n_end]
-                    n_start += blocks[s].numCols
-                end
-            end
-        end
-    end
-
-    # maximization?
-    if m.objSense == :Max
-        m.objVal *= -1
-        dsp.primVal *= -1
-        dsp.dualVal *= -1
+        getDspSolution(m)
     end
 
     # Return the solve status
     stat
-
 end
 
 ###############################################################################
@@ -130,7 +143,7 @@ function optimize(;suppress_warnings = false, options...)
     end
 
     if !(stat == :Infeasible || stat == :Unbounded)
-        getDspSolution(suppress_warnings)
+        getDspSolution()
     end
 
     # Return the solve status
@@ -140,10 +153,8 @@ JuMP.solve(;suppress_warnings = false, options...) = optimize(;suppress_warnings
 
 # Read model from SMPS files
 function readSmps(filename::AbstractString)
-
     # free DspModel
     DspCInterface.freeModel(Dsp.model)
-
     # read Smps file
     DspCInterface.readSmps(Dsp.model, filename)
 end
@@ -170,8 +181,14 @@ end
 getprimobjval() = Dsp.model.primVal
 # get dual objective value
 getdualobjval() = Dsp.model.dualVal
+# get primal value
+getprimvalue() = Dsp.model.colVal
 # get dual value
-JuMP.getdual() = Dsp.model.rowVal
+getdualvalue() = Dsp.model.rowVal
+# get solution time
+getsolutiontime() = DspCInterface.getWallTime(Dsp.model)
+# get block ids
+blockids() = Dsp.model.block_ids
 
 function parseStatusCode(statcode::Integer)
     stat = :NotSolved
@@ -199,15 +216,45 @@ function parseStatusCode(statcode::Integer)
     stat
 end
 
-function getDspSolution(suppress_warnings)
+function getDspSolution(m::JuMP.Model = nothing)
     Dsp.model.primVal = DspCInterface.getPrimalBound(Dsp.model)
     Dsp.model.dualVal = DspCInterface.getDualBound(Dsp.model)
 
     if Dsp.model.solve_type == :Dual
-        suppress_warnings || warn("Primal solution is not available in dual decomposition")
+        Dsp.model.rowVal = DspCInterface.getDualSolution(Dsp.model)
     else
         Dsp.model.colVal = DspCInterface.getSolution(Dsp.model)
+        if m != nothing
+            # parse solution to each block
+            n_start = 1
+            n_end = m.numCols
+            m.colVal = Dsp.model.colVal[n_start:n_end]
+            n_start += m.numCols
+            if haskey(m.ext, :DspBlocks) == true
+                numBlockCols = DspCInterface.getNumBlockCols(Dsp.model, m)
+                blocks = m.ext[:DspBlocks].children
+                for i in 1:Dsp.model.nblocks
+                    n_end += numBlockCols[i]
+                    if haskey(blocks, i)
+                        # @show b
+                        # @show n_start
+                        # @show n_end
+                        blocks[i].colVal = Dsp.model.colVal[n_start:n_end]
+                    end
+                    n_start += numBlockCols[i]
+                end
+            end
+        end
+    end
+    if m != nothing
+        m.objVal = Dsp.model.primVal
+        # maximization?
+        if m.objSense == :Max
+            m.objVal *= -1
+            dsp.primVal *= -1
+            dsp.dualVal *= -1
+        end
     end
 end
 
-end # module
\ No newline at end of file
+end # module
diff --git a/src/DspCInterface.jl b/src/DspCInterface.jl
index 1d9b337..6eba704 100644
--- a/src/DspCInterface.jl
+++ b/src/DspCInterface.jl
@@ -24,6 +24,9 @@ end
 type DspModel
     p::Ptr{Void}
 
+    # Number of blocks
+    nblocks::Int
+
     # solve_type should be one of these:
     # :Dual
     # :Benders
@@ -37,8 +40,21 @@ type DspModel
     colVal::Vector{Float64}
     rowVal::Vector{Float64}
 
+    # MPI settings
+    comm
+    comm_size::Int
+    comm_rank::Int
+
+    # Array of block ids:
+    # The size of array is not necessarily same as nblocks,
+    # as block ids may be distributed to multiple processors.
+    block_ids::Vector{Integer}
+
     function DspModel()
+        # assign Dsp pointer
         p = @dsp_ccall("createEnv", Ptr{Void}, ())
+        # initialize variables
+        nblocks = 0
         solve_type = :Dual
         numRows = 0
         numCols = 0
@@ -46,44 +62,58 @@ type DspModel
         dualVal = NaN
         colVal = Vector{Float64}()
         rowVal = Vector{Float64}()
-        prob = new(p, solve_type, numRows, numCols, primVal, dualVal, colVal, rowVal)
-        finalizer(prob, freeDSP)
-        return prob
+        comm = nothing
+        comm_size = 1
+        comm_rank = 0
+        block_ids = Vector{Integer}()
+        # create DspModel
+        dsp = new(p, nblocks, solve_type, numRows, numCols, primVal, dualVal, colVal, rowVal, comm, comm_size, comm_rank, block_ids)
+        # with finalizer
+        finalizer(dsp, freeDSP)
+        # return DspModel
+        return dsp
     end
 end
 
-function freeDSP(prob::DspModel)
-    if prob.p == C_NULL
+function freeDSP(dsp::DspModel)
+    if dsp.p == C_NULL
         return
     else
-        @dsp_ccall("freeEnv", Void, (Ptr{Void},), prob.p)
-        prob.p = C_NULL
+        @dsp_ccall("freeEnv", Void, (Ptr{Void},), dsp.p)
+        dsp.p = C_NULL
     end
-    prob.solve_type = nothing
-    prob.numRows = 0
-    prob.numCols = 0
-    prob.primVal = NaN
-    prob.dualVal = NaN
-    prob.colVal = Vector{Float64}()
-    prob.rowVal = Vector{Float64}()
+    dsp.nblocks = 0
+    dsp.solve_type = nothing
+    dsp.numRows = 0
+    dsp.numCols = 0
+    dsp.primVal = NaN
+    dsp.dualVal = NaN
+    dsp.colVal = Vector{Float64}()
+    dsp.rowVal = Vector{Float64}()
     return
 end
 
-function freeModel(prob::DspModel)
-    check_problem(prob)
-    @dsp_ccall("freeTssModel", Void, (Ptr{Void},), prob.p)
+function freeModel(dsp::DspModel)
+    check_problem(dsp)
+    @dsp_ccall("freeModel", Void, (Ptr{Void},), dsp.p)
+    dsp.nblocks = 0
+    dsp.numRows = 0
+    dsp.numCols = 0
+    dsp.primVal = NaN
+    dsp.dualVal = NaN
+    dsp.colVal = Vector{Float64}()
+    dsp.rowVal = Vector{Float64}()
 end
 
-function check_problem(prob::DspModel)
-    if prob.p == C_NULL
+function check_problem(dsp::DspModel)
+    if dsp.p == C_NULL
         error("Invalid DspModel")
     end
-    return true
 end
 
-function readParamFile(prob::DspModel, param_file::AbstractString)
-    check_problem(prob)
-    @dsp_ccall("readParamFile", Void, (Ptr{Void}, Ptr{UInt8}), prob.p, param_file);
+function readParamFile(dsp::DspModel, param_file::AbstractString)
+    check_problem(dsp)
+    @dsp_ccall("readParamFile", Void, (Ptr{Void}, Ptr{UInt8}), dsp.p, param_file);
 end
 
 function prepConstrMatrix(m::JuMP.Model)
@@ -117,60 +147,136 @@ function prepConstrMatrix(m::JuMP.Model)
     return sparse(rind, cind, value, length(m.linconstr), blocks.parent.numCols + m.numCols)
 end
 
+###############################################################################
+# Block IDs
+###############################################################################
+
+function setBlockIds(dsp::DspModel, nblocks::Integer)
+    check_problem(dsp)
+    # set number of blocks
+    dsp.nblocks = nblocks
+    # set MPI settings
+    if isdefined(:MPI) && MPI.Initialized()
+        dsp.comm = MPI.COMM_WORLD
+        dsp.comm_size = MPI.Comm_size(dsp.comm)
+        dsp.comm_rank = MPI.Comm_rank(dsp.comm)
+    end
+    #@show dsp.nblocks
+    #@show dsp.comm
+    #@show dsp.comm_size
+    #@show dsp.comm_rank
+    # get block ids with MPI settings
+    dsp.block_ids = getBlockIds(dsp)
+    #@show dsp.block_ids
+    # send the block ids to Dsp
+    @dsp_ccall("setIntPtrParam", Void, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Cint}), 
+        dsp.p, "ARR_PROC_IDX", convert(Cint, length(dsp.block_ids)), convert(Vector{Cint}, dsp.block_ids - 1))
+end
+
+function getBlockIds(dsp::DspModel)
+    check_problem(dsp)
+    # processor info
+    mysize = dsp.comm_size
+    myrank = dsp.comm_rank
+    # empty block ids
+    proc_idx_set = Int[]
+    # DSP is further parallelized with mysize > dsp.nblocks.
+    modrank = myrank % dsp.nblocks
+    # If we have more than one processor, 
+    # do not assign a sub-block to the master.
+    if mysize > 1
+        if myrank == 0
+            return proc_idx_set
+        end
+        # exclude master
+        mysize -= 1;
+        modrank = (myrank-1) % dsp.nblocks
+    end
+    # assign sub-blocks in round-robin fashion
+    for s = modrank:mysize:(dsp.nblocks-1)
+        push!(proc_idx_set, s+1)
+    end
+    # return assigned block ids
+    return proc_idx_set
+end
+
+function getNumBlockCols(dsp::DspModel, m::JuMP.Model)
+    check_problem(dsp)
+    # subblocks
+    blocks = m.ext[:DspBlocks].children
+    # get number of block columns
+    numBlockCols = Dict{Int,Int}()
+    if dsp.comm_size > 1
+        num_proc_blocks = convert(Vector{Cint}, MPI.Allgather(length(blocks), dsp.comm))
+        #@show num_proc_blocks
+        #@show collect(keys(blocks))
+        block_ids = MPI.Allgatherv(collect(keys(blocks)), num_proc_blocks, dsp.comm)
+        #@show block_ids
+        ncols_to_send = Int[blocks[i].numCols for i in keys(blocks)]
+        #@show ncols_to_send
+        ncols = MPI.Allgatherv(ncols_to_send, num_proc_blocks, dsp.comm)
+        #@show ncols
+        for i in 1:dsp.nblocks
+            setindex!(numBlockCols, ncols[i], block_ids[i])
+        end
+    else
+        for b in blocks
+            setindex!(numBlockCols, b.second.numCols, b.first)
+        end
+    end
+    return numBlockCols
+end
 
 ###############################################################################
 # Load problems
 ###############################################################################
 
-function readSmps(prob::DspModel, filename::AbstractString, dedicatedMaster::Bool)
+function readSmps(dsp::DspModel, filename::AbstractString)
     # Check pointer to TssModel
-    check_problem(prob)
-    @dsp_ccall("readSmps", Void, (Ptr{Void}, Ptr{UInt8}), prob.p, convert(Vector{UInt8}, filename))
-    nscen = getNumScenarios(prob);
-    proc_idx_set = collect(1:nscen);
-    if isdefined(:MPI) == true
-        proc_idx_set = getProcIdxSet(nscen, dedicatedMaster);
-    end
-    setProcIdxSet(prob, proc_idx_set);  
+    check_problem(dsp)
+    # read smps files
+    @dsp_ccall("readSmps", Void, (Ptr{Void}, Ptr{UInt8}), dsp.p, convert(Vector{UInt8}, filename))
+    # set block Ids
+    setBlockIds(dsp, getNumScenarios(dsp))
 end
-readSmps(prob::DspModel, filename::AbstractString) = readSmps(prob, filename, false);
 
-function loadProblem(prob::DspModel, model::JuMP.Model, dedicatedMaster::Bool)
-    check_problem(prob)
+function loadProblem(dsp::DspModel, model::JuMP.Model, dedicatedMaster::Bool)
+    check_problem(dsp)
     if haskey(model.ext, :DspBlocks)
-        loadStochasticProblem(prob, model, dedicatedMaster);
+        loadStochasticProblem(dsp, model, dedicatedMaster)
     else
         warn("No block is defined.")
-        loadDeterministicProblem(prob, model);
+        loadDeterministicProblem(dsp, model)
     end
 end
-loadProblem(prob::DspModel, model::JuMP.Model) = loadProblem(prob, model, false);
+loadProblem(dsp::DspModel, model::JuMP.Model) = loadProblem(dsp, model, true);
 
-function loadStochasticProblem(prob::DspModel, model::JuMP.Model, dedicatedMaster::Bool)
+function loadStochasticProblem(dsp::DspModel, model::JuMP.Model, dedicatedMaster::Bool)
     
     # get DspBlocks
     blocks = model.ext[:DspBlocks]
     
-    nscen  = convert(Cint, length(blocks.children))
-    ncols1 = convert(Cint, model.numCols)
-    nrows1 = convert(Cint, length(model.linconstr))
+    nscen  = dsp.nblocks
+    ncols1 = model.numCols
+    nrows1 = length(model.linconstr)
     ncols2 = 0
     nrows2 = 0
-    proc_idx_set = collect(1:nscen);
-    if isdefined(:MPI) == true && MPI.Initialized() == true
-        proc_idx_set = getProcIdxSet(nscen, dedicatedMaster);
+    for s in values(blocks.children)
+        ncols2 = s.numCols
+        nrows2 = length(s.linconstr)
+        break
     end
-    setProcIdxSet(prob, proc_idx_set);
-    for s in 1:length(proc_idx_set)
-        ncols2 = convert(Cint, blocks.children[s].numCols)
-        nrows2 = convert(Cint, length(blocks.children[s].linconstr))
-        break;
+    
+    # set scenario indices for each MPI processor
+    if dsp.comm_size > 1
+        ncols2 = MPI.allreduce([ncols2], MPI.MAX, dsp.comm)[1]
+        nrows2 = MPI.allreduce([nrows2], MPI.MAX, dsp.comm)[1]
     end
     
-    @dsp_ccall("setNumberOfScenarios", Void, (Ptr{Void}, Cint), prob.p, nscen)
+    @dsp_ccall("setNumberOfScenarios", Void, (Ptr{Void}, Cint), dsp.p, convert(Cint, nscen))
     @dsp_ccall("setDimensions", Void, 
         (Ptr{Void}, Cint, Cint, Cint, Cint), 
-        prob.p, ncols1, nrows1, ncols2, nrows2)
+        dsp.p, convert(Cint, ncols1), convert(Cint, nrows1), convert(Cint, ncols2), convert(Cint, nrows2))
     
     # get problem data
     start, index, value, clbd, cubd, ctype, obj, rlbd, rubd = getDataFormat(model)
@@ -178,23 +284,23 @@ function loadStochasticProblem(prob::DspModel, model::JuMP.Model, dedicatedMaste
     @dsp_ccall("loadFirstStage", Void, 
         (Ptr{Void}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, 
             Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
-            prob.p, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
-    
-    for s in 1:length(proc_idx_set)
-        # get model
-        blk = blocks.children[s]
-        probability = blocks.weight[s]
+            dsp.p, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
+
+    for id in dsp.block_ids
+        # model and probability
+        blk = blocks.children[id]
+        probability = blocks.weight[id]
         # get model data
         start, index, value, clbd, cubd, ctype, obj, rlbd, rubd = getDataFormat(blk)
         @dsp_ccall("loadSecondStage", Void, 
             (Ptr{Void}, Cint, Cdouble, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, 
                 Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), 
-            prob.p, proc_idx_set[s]-1, probability, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
+            dsp.p, id-1, probability, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
     end
     
 end
 
-function loadDeterministicProblem(prob::DspModel, model::JuMP.Model)
+function loadDeterministicProblem(dsp::DspModel, model::JuMP.Model)
     ncols = convert(Cint, model.numCols)
     nrows = convert(Cint, length(model.linconstr))
     start, index, value, clbd, cubd, ctype, obj, rlbd, rubd = getDataFormat(model)
@@ -202,7 +308,7 @@ function loadDeterministicProblem(prob::DspModel, model::JuMP.Model)
     @dsp_ccall("loadDeterministic", Void, 
         (Ptr{Void}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Cint, Cint, Cint, 
             Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
-            prob.p, start, index, value, numels, ncols, nrows, clbd, cubd, ctype, obj, rlbd, rubd)
+            dsp.p, start, index, value, numels, ncols, nrows, clbd, cubd, ctype, obj, rlbd, rubd)
 end
 
 
@@ -210,59 +316,52 @@ end
 # Get functions
 ###############################################################################
 
-function solveDe(prob::DspModel)
-    check_problem(prob)
-    return @dsp_ccall("solveDe", Void, (Ptr{Void},), prob.p)
-end
-
-function solveBd(prob::DspModel)
-    check_problem(prob)
-    return @dsp_ccall("solveBd", Void, (Ptr{Void},), prob.p)
+for func in [:freeSolver, 
+             :solveDe, 
+             :solveBd, 
+             :solveDd]
+    strfunc = string(func)
+    @eval begin
+        function $func(dsp::DspModel)
+            return @dsp_ccall($strfunc, Void, (Ptr{Void},), dsp.p)
+        end
+    end
 end
 
-function solveDd(prob::DspModel, comm)
-    check_problem(prob)
-    return @dsp_ccall("solveDd", Void, (Ptr{Void}, Cint), prob.p, convert(Cint, comm.val))
+for func in [:solveBdMpi, :solveDdMpi]
+    strfunc = string(func)
+    @eval begin
+        function $func(dsp::DspModel, comm)
+            return @dsp_ccall($strfunc, Void, (Ptr{Void}, Cint), dsp.p, convert(Cint, comm.val))
+        end
+    end
 end
 
-function solveBdMpi(prob::DspModel, comm)
-    check_problem(prob)
-    return @dsp_ccall("solveBdMpi", Void, (Ptr{Void}, Cint, Cint), prob.p, 1, convert(Cint, comm.val))
+function solve(dsp::DspModel)
+    check_problem(dsp)
+    if dsp.comm_size == 1
+        if dsp.solve_type == :Dual
+            solveDd(dsp);
+        elseif dsp.solve_type == :Benders
+            solveBd(dsp);
+        elseif dsp.solve_type == :Extensive
+            solveDe(dsp);
+        end
+    elseif dsp.comm_size > 1
+        if dsp.solve_type == :Dual
+            solveDdMpi(dsp, dsp.comm);
+        elseif dsp.solve_type == :Benders
+            solveBdMpi(dsp, dsp.comm);
+        elseif dsp.solve_type == :Extensive
+            solveDe(dsp);
+        end
+    end
 end
 
 ###############################################################################
 # Get functions
 ###############################################################################
 
-function getProcIdxSet(numScens::Integer, dedicatedMaster::Bool)
-    mysize = 1;
-    myrank = 0;
-    if isdefined(:MPI) == true && MPI.Initialized() == true
-        comm = MPI.COMM_WORLD
-        mysize = MPI.Comm_size(comm)
-        myrank = MPI.Comm_rank(comm)
-    end
-    # Round-and-Robin
-    proc_idx_set = Int[];
-    # DSP is further parallelized with mysize > numScens.
-    modrank = myrank % numScens;
-    
-    if dedicatedMaster == true
-        if myrank == 0
-            return proc_idx_set;
-        end
-        # exclude master
-        mysize -= 1;
-        modrank = (myrank-1) % numScens;
-    end
-    
-    for s = modrank:mysize:(numScens-1)
-        push!(proc_idx_set, s+1);
-    end
-    return proc_idx_set;
-end
-getProcIdxSet(numScens::Integer) = getProcIdxSet(numScens,false);
-
 function getDataFormat(model::JuMP.Model)
     # Get a column-wise sparse matrix
     mat = prepConstrMatrix(model)
@@ -304,53 +403,52 @@ for (func,rtn) in [(:getNumScenarios, Cint),
                    (:getSolutionStatus, Cint), 
                    (:getNumIterations, Cint), 
                    (:getNumNodes, Cint), 
-                   (:getSolutionTime, Cdouble), 
+                   (:getWallTime, Cdouble), 
                    (:getPrimalBound, Cdouble), 
                    (:getDualBound, Cdouble)]
     strfunc = string(func)
     @eval begin
-        function $func(prob::DspModel)
-            check_problem(prob)
-            return @dsp_ccall($strfunc, $rtn, (Ptr{Void},), prob.p)
+        function $func(dsp::DspModel)
+            check_problem(dsp)
+            return @dsp_ccall($strfunc, $rtn, (Ptr{Void},), dsp.p)
         end
     end
 end
 
-function getObjCoef(prob::DspModel)
-    check_problem(prob)
+function getObjCoef(dsp::DspModel)
+    check_problem(dsp)
     num = getTotalNumCols()
     obj = Array(Cdouble, num)
-    @dsp_ccall("getObjCoef", Void, (Ptr{Void}, Ptr{Cdouble}), prob.p, obj)
+    @dsp_ccall("getObjCoef", Void, (Ptr{Void}, Ptr{Cdouble}), dsp.p, obj)
     return obj
 end
 
-function getSolution(prob::DspModel, num::Integer)
+function getSolution(dsp::DspModel, num::Integer)
     sol = Array(Cdouble, num)
-    @dsp_ccall("getSolution", Void, (Ptr{Void}, Cint, Ptr{Cdouble}), prob.p, num, sol)
+    @dsp_ccall("getPrimalSolution", Void, (Ptr{Void}, Cint, Ptr{Cdouble}), dsp.p, num, sol)
     return sol
 end
-getSolution(prob::DspModel) = getSolution(prob, getTotalNumCols(prob))
+getSolution(dsp::DspModel) = getSolution(dsp, getTotalNumCols(dsp))
+
+function getDualSolution(dsp::DspModel, num::Integer)
+    sol = Array(Cdouble, num)
+    @dsp_ccall("getDualSolution", Void, (Ptr{Void}, Cint, Ptr{Cdouble}), dsp.p, num, sol)
+    return sol
+end
+getDualSolution(dsp::DspModel) = getDualSolution(dsp, getNumCouplingRows(dsp))
 
 ###############################################################################
 # Set functions
 ###############################################################################
 
-function setSolverType(prob::DspModel, solver)
-    check_problem(prob)
+function setSolverType(dsp::DspModel, solver)
+    check_problem(dsp)
     solver_types = [:DualDecomp, :Benders, :ExtensiveForm]
     if solver in solver_types
-        prob.solver = solver
+        dsp.solver = solver
     else
         warn("Solver type $solver is invalid.")
     end
 end
 
-function setProcIdxSet(prob::DspModel, scenarios::Array{Int,1})
-    check_problem(prob)
-    num = convert(Cint, length(scenarios));
-    scenarios = scenarios - 1;
-    @dsp_ccall("setIntPtrParam", Void, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Cint}), 
-        prob.p, "ARR_PROC_IDX", convert(Cint, num), convert(Vector{Cint}, scenarios))
-end
-
-end # end of module
\ No newline at end of file
+end # end of module
diff --git a/src/block.jl b/src/block.jl
index db973c7..6c49c18 100644
--- a/src/block.jl
+++ b/src/block.jl
@@ -1,28 +1,7 @@
 import JuMP
-export @block
 
 type BlockStructure
     parent
     children::Dict{Int,JuMP.Model}
     weight::Dict{Int,Float64}
-end
-
-macro block(self, child, id, weight)
-    return quote
-        if haskey($(esc(self)).ext, :DspBlocks) == false
-            $(esc(self)).ext[:DspBlocks] = BlockStructure(
-                nothing, Dict{Int,JuMP.Model}(), Dict{Int,Float64}())
-        end
-        if haskey($(esc(child)).ext, :DspBlocks) == false
-            $(esc(child)).ext[:DspBlocks] = BlockStructure(
-                $(esc(self)), Dict{Int,JuMP.Model}(), Dict{Int,Float64}())
-        else
-            $(esc(child)).ext[:DspBlocks].parent = $(esc(self))
-        end
-        $(esc(self)).ext[:DspBlocks].children[$(esc(id))] = $(esc(child))
-        $(esc(self)).ext[:DspBlocks].weight[$(esc(id))]   = $(esc(weight))
-        if $(esc(self)).ext[:DspBlocks].parent == nothing
-            JuMP.setsolvehook($(esc(self)), Dsp.dsp_solve)
-        end
-    end
 end
\ No newline at end of file