Parallelizing For Loops - Issue

Hi there, I'm attempting to parallelize a nested for loop as seen below that is running linear programs via a Gurobi solver.

for ir = 1:n-1
    
    Aeq = AugS; beq = zeros(m,1);
    
    newRow = zeros(1,n);
    newRow(ir) = 1;
    AeqNew = [Aeq; newRow];
    beqNew = [beq; 1];
    
    Aeqlb = [modelUpd.lb; -1];
    Aeqlb(idxCSource) = CsourceLB;
    Aeqlb(idxObj) = BiomassLB;
    Aequb = modelUpd.ub;
    
    Aineq = [diag(diag(ones(length(modelUpd.rxns)))) -Aequb(1:end);
        -diag(diag(ones(length(modelUpd.rxns)))) Aeqlb(1:end-1);
        zeros(1,length(modelUpd.rxns)) Aeqlb(end)];
    
    bineq = zeros(size(Aineq,1),1);
    
    
    for jr = 1:1:n-1
        
        %Build Gurobi model;        
        modelG.A = sparse([Aineq; AeqNew]);
        modelG.rhs = [bineq; beqNew];
        modelG.lb = -1000*ones(n,1);
        modelG.ub = 1000*ones(n,1);
        modelG.sense = [repmat('<',1,size(Aineq,1)) repmat('=',1,size(AeqNew,1))];
        
        %To find the minimum flux values.  
        f = zeros(1,n);
        f(jr) = 1;        
        modelG.modelsense = 'min';
        modelG.obj = f;
        
        params.outputflag = 0;
        result = gurobi(modelG,params);
        
        if strcmp(result.status, 'OPTIMAL')
            fluxMin1(ir,jr) = result.objval;
        else
            fluxMin1(ir,jr) = "N/A";
        end
        
        
        %To find the maximum flux values.
        f(jr) = 1;
        modelG.modelsense = 'max';
        modelG.obj = f;
        
        params.outputflag = 0;
        result = gurobi(modelG,params); 
        
        if strcmp(result.status, 'OPTIMAL')
            fluxMax1(ir,jr) = result.objval;
        else
            fluxMax1(ir,jr) = "N/A";
        end
        
        f(jr) = 0;
    end
end

Unfortunately, I keep getting the error: 'The PARFOR look cannot run due to the way variable 'modelG' is used.'

However, I cannot fix this variable as I'd seen similarly done as the modelG variable is a structure containing sparse matrices, modelsense options, etc.

Any help with somehow parallelizing this code would be great as it's currently running 8.4 million linear programs which is taking 0.5 secs each -- so I'm looking at a code that is running 3.5 days. :'D

Thanks!

Answer 1

[Matt J]

I would use parfor just to run the n gurobi optimizations. After you've collected all the results in a struct array Results(jr), you can post-process them.

        modelG.A = sparse([Aineq; AeqNew]);
        modelG.rhs = [bineq; beqNew];
        modelG.lb = -1000*ones(n,1);
        modelG.ub = 1000*ones(n,1);
        modelG.sense = [repmat('<',1,size(Aineq,1)) repmat('=',1,size(AeqNew,1))];
        modelG.modelsense = 'min'; 
        modelG.obj=zeros(1,n);

        params.outputflag = 0;      
        
        parfor jr=1:n
 
         mdl=modelG;
         mdl.obj(jr) = 1;     
         s = gurobi(mdl,params);

         Results(jr).status=s.status;
         Results(jr).optval=s.optval;
         
        end

[Walter Roberson]

Structuring it this way to avoid re-assigning all of those fields is a good idea.

Answer 2

[Walter Roberson]

It looks like it is not treating modelG as a local variable. That implies that there is a reference to modelG either before the parfor or else after the parfor.

I suggest that at the beginning of the for j loop (which is what I presume you are turning into parfor) that you add

modelG = struct();

That will force parfor j to treat modelG as local to the parfor iteration.

Answer 3

[Sam Marshalik]

You may be able to put the contents of the parfor-loop into a separate function and calling that function from that parfor-loop. It resolves the error message, but I have not had a chance to run it.