Problem in Using Parallel Feature for a "For Loop"!

I have written a code for Newton Raphson which includes two "for loops" inside of each other. Basically, I have a 1000 nonlinear equations, which are in a function called "f_int". In the "for loop" I try to construct tangent stiffness of these 1000 equations by finite difference, which is as below:

f_int = internal_eqns1000(X);
for i=1:1000
    for j=1:nvar
        X(j)=X(j)+h;
        f_intp1 = internal_eqns1000(X);
        k_k(i,j)=(f_intp1(i)-f_int(i))/h;
        X(j)=X(j)-h;
    end
end

After the completion of tangent stiffness, another "X" will be created and the loop will run again until convergence. This for loop might be required to run more than 100 times, and each time it takes 30 minutes average. I would like to use parallel to decrease the computational timing, but I was not successful. The change I made is as below:

f_int = internal_eqns1000(X);
parfor i=1:1000
    for j=1:nvar
        X(j)=X(j)+h;
        f_intp1 = internal_eqns1000(X);
        k_k(i,j)=(f_intp1(i)-f_int(i))/h;
        X(j)=X(j)-h;
    end
end

But this does not work and gives me an error: "due to the way "X" is defined, parallel computation cannot be used here".

Please let me know how can I write the code such that it will be run with parallel feature. Thank you.

Answer 1

[OCDER]

X has to be sliced correctly. https://www.mathworks.com/help/distcomp/sliced-variable.html

But I didn't realize Newton Raphson iterations can be done in parallel. Or are you doing independent NR calc for many situations?

X = zeros(1000, nvar);
f_int = internal_eqns1000(X); 
parfor i=1:1000 
    for j=1:nvar 
        X(i,j)=X(i,j)+h; 
        f_intp1 = internal_eqns1000(X(i,j)); 
        k_k(i,j)=(f_intp1(i)-f_int(i))/h; 
        X(i,j)=X(i,j)-h; 
    end
end

[Maryam]

Thank you for your response. I am not using parallel for NR, I just want to use parallel for calculating the tangent stiffness (using finite difference).

[Maryam]

Using X(i,j) makes the parallel pool to run but it is not working. I have tried "X(i,j)", but this won't work and the reason is that X is an array from X(1) to X(1000); therefore, giving it two subscripts cause MATLAB to give an error of "Index exceeds matrix dimensions".

[OCDER]

Is X a constant? You'll have to make it clear it's a broadcast variable then as such:

X = zeros(1000, 1);
f_int = internal_eqns1000(Xt); 
parfor i=1:1000 
    Xt = X; %broadcast variable
    for j=1:nvar 
        Xt(j)=Xt(j)+h; 
        f_intp1 = internal_eqns1000(Xt(j)); 
        k_k(i,j)=(f_intp1(i)-f_int(i))/h;  %sliced variable
        Xt(j)=Xt(j)-h; 
    end
end

[Maryam]

X is constant, and it changes every time tangent stiffness calculated. My entire code is as below:

h = 10e-4; g = load('ext-load_1000'); Fext = (g.F_ext)'; g_1 = load('initial_1000'); q = (g_1.X)'; X=q;

while norm(f_residual) > errtol

    iter_NR = iter_NR+1;

    f_int = internal_eqns1000(X);
    f_residual = f_int'-Fext;

for i=1:nvar for j=1:nvar X(j)=X(j)+h; f_intp1 = internal_eqns1000(X); k_k(i,j)=(f_intp1(i)-f_int(i))/h; X(j)=X(j)-h; end end

b = k_k\(f_int'-Fext); X(i) = X(i) - b(i); f_residual = (f_int'-Fext); N_F = norm(f_residual); det_k_k = abs(det(k_k));

end

So basically, x will be calculated and assigned a new array each time after tangent stiffness calculation.

The new way you defined Xt does not work. It gives me dimensional error again. The internal_equations1000 contains 1000 nonlinear equations, all of which are function of array X (X(1) to X(1000)). The first two equations are as follow:

f_int(1) = -56.0d0*X(1) + 56.0d0*X(2) ... + 0.0d0;

f_int(2) = 56.0d0*X(2) + -12060.0d0*X(2) ... + 11299.0d0*X(3) ... + 563.0d0*X(102) ... + 0.0d0;

Thank you so much for your time and consideration.

[OCDER]

Can you do us a big favor and edit your code by pushing the "{} Code" to make your code readable? EXAMPLE:

for j=1:nvar end

After pushing {}Code -->

for j=1:nvar
end

[Maryam]

of course! Sorry for the mistake.

 h = 10e-4;
g = load('ext-load_1000');
Fext = (g.F_ext)';
g_1 = load('initial_1000');
q = (g_1.X)';
X=q;                                 

 while norm(f_residual) > errtol

    iter_NR = iter_NR+1;

    f_int = internal_eqns1000(X);
    f_residual = f_int'-Fext;

 for i=1:nvar
    for j=1:nvar
        X(j)=X(j)+h;
        f_intp1 = internal_eqns1000(X);
        k_k(i,j)=(f_intp1(i)-f_int(i))/h;
        X(j)=X(j)-h;
    end
 end

 b = k_k\(f_int'-Fext);
 for ji=1:nvar
    X(ji) = X(ji) - b(ji);
 end
 f_residual = (f_int'-Fext);
 N_F = norm(f_residual);
 det_k_k = abs(det(k_k));

 end

[Maryam]

I have used the last method you recommended, and I just changed it a bit and it is working now. Thank you. The changed I made is at "f_intp1 = internal_eqns1000(Xt)", in which Xt, should not have subscript.

 parfor i=1:nvar
    Xt=X;
    for j=1:nvar
        Xt(j)=Xt(j)+h;
        f_intp1 = internal_eqns1000(Xt);
        k_k(i,j)=(f_intp1(i)-f_int(i))/h;
        Xt(j)=Xt(j)-h;
    end
end

Thank you so much for your responses and your patience with me :)

[OCDER]

h = 10e-4;
g = load('ext-load_1000');
Fext = (g.F_ext)';
g_1 = load('initial_1000');
X = (g_1.X)';
f_residual = Inf; %What's the initial value of this?
nvar = length(X); %What's nvar?
while norm(f_residual) > errtol
    iter_NR = iter_NR+1;
    f_int = internal_eqns1000(X);
    %f_residual = f_int'-Fext;  %Unused?
    k_k = zeros(nvar); %Initialize
    parfor i=1:nvar
        for j=1:nvar
            %X(j)=X(j)+h; %So this is a temporary change in X?
            Xt = X;               %Maybe make the temporary var here as Xt + h
            Xt(j) = Xt(j) + h; 
            f_intp1 = internal_eqns1000(Xt);
            k_k(i,j)=(f_intp1(i)-f_int(i))/h;
            %X(j)=X(j)-h; %So this restore X back to original?
        end
    end
    b = k_k\(f_int'-Fext);
    for ji=1:nvar
        X(ji) = X(ji) - b(ji);
    end
    f_residual = (f_int'-Fext); %
    N_F = norm(f_residual);
    det_k_k = abs(det(k_k));
end

[Maryam]

Thank you so much. It is completely clear and fine. I just have a quick question about "number of workers in a parallel pool". It uses 4 workers for me automatically. How can I increase it to 1000 for example? Thanks again for your response.

[OCDER]

You can't get more workers than the number of cores in your computer. To change the number of workers to less than maximum, you could do this:

 NumCores = 2; %Want 2 cores
 delete(gcp('nocreate'));
 parpool(NumCores);

[Maryam]

Thank you!