How to vectorize inside user defined functions to avoid nested loops?

Hello everyone, I'm trying to vectorize some code so as to avoid making a nested loop. However, at the innermost loop I have to run a function that I've defined myself, and I can't seem to find a way to make that function "understand" its inputs properly. What I have right now, which works but is rather slow, is something like this:

for i = 1:n*m
  for j = 1:n*m
      I_ww(i,j) = InfluenceLine(XcWing(i,j),YcWing(i,j),ZcWing(i,j),XvWing,YvWing,ZvWing,XnWing(i,j),YnWing(i,j),ZnWing(i,j),AngleOfAttack);
  end
end

I'd like to rewrite that so it looks something like this:

I_ww(:,:) = InfluenceLine(XcWing(:,:),YcWing(:,:),ZcWing(:,:),XvWing,YvWing,ZvWing,XnWing(:,:),YnWing(:,:),ZnWing(:,:),AngleOfAttack);

I've realized, however, that with this syntax, the inputs which were supposed to be scalars, namely XcWing(i,j) and so forth, are now the entire matrices XcWing, which my function does not support. Is there a way to get the function to compute only one element of each of the input matrices at a time?

[jgg]

You can call arrayfun twice, but this isn't going to be faster than looping. Unfortunately, it's likely that if you want to vectorize this, you need to rewrite your function to take either vectors or matrices instead of scalars.

Matlab is really only faster when using its vector operations. Any way to recombine scalar operations isn't going to improve performance significantly.

[Stephen23]

Vectorization does not mean replacing my slow loops with some other code, it means writing functions that operate on entire arrays at once.

This means it is the function that is important, not the loops. However you do not tell us anything about the function, yet the function is the only important thing to consider. If you give us details about the function then we can give advice about vectorization.

[André Camões]

Allright, I think I understand. Anyhow, what the function does is actually call this function:

 function [Line] =   InfluenceLine(xc,yc,zc,Xv,Yv,Zv,xn,yn,zn,AngleOfAttack)
[n,m] = size(Xv);
Line = zeros(n-1,m-1);
%Calculates the influence of each vortex ring on the colocation point
%(xc,yc,zc) whose normal is (xn,yn,zn)
Line(:,:) = RingInfluenceCoefficient(xc,yc,zc,xn,yn,zn,Xv(1:end-1,1:end-1),Yv(1:end-1,1:end-1),Zv(1:end-1,1:end-1),...
    Xv(2:end,1:end-1),Yv(2:end,1:end-1),Zv(2:end,1:end-1),Xv(2:end,2:end),Yv(2:end,2:end),Zv(2:end,2:end),...
    Xv(1:end-1,2:end),Yv(1:end-1,2:end),Zv(1:end-1,2:end));
Line = transpose(line);
Line = reshape(Line,[1,numel(Line)]);
end

So that function alone has the only purpose to "avoid" even further loops (and to keep things organized, really). The inputs xc, yc and zc, as well as xn, yn and zn are supposed to be scalars, as they were when I was calling the function InfluenceLine in the form XcWing(i,j). Anyhow, we then have the function RingInfluenceCoefficient, which does the following:

 function [influence] = RingInfluenceCoefficient(xc,yc,zc,xn,yn,zn,x1,y1,z1,...
x2,y2,z2,x3,y3,z3,x4,y4,z4)
%Calculates the influence coefficient of a single vortex ring with vertices
%(xi,yi,zi), where i = 1:4 in that order, on the collocation point (xc,yc,zc)
%whose normal vector is (xn,yn,zn). This calculates the influence of both the
%vortices at the positive y region and their mirror images about the xz plane.
[u1,v1,w1] = RingVortexInducedVelocity(xc,yc,zc,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,1);
[u2,v2,w2] = RingVortexInducedVelocity(xc,yc,zc,x4,-y4,z4,x3,-y3,z3,x2,-y2,z2,x1,-y1,z1,1);
influence = dot([u1+u2,v1+v2,w1+w2],[xn,yn,zn]);
end

By this point, all the inputs should be scalar in the way I'm currently doing it, so there are no further loops inside any function from this point on. At the end of this function I calculate a dot-product, so maybe I'll need to implement some matrix-friendly version of that myself when I vectorize. Anyway, then there is the next function:

 function [u,v,w] = RingVortexInducedVelocity(x,y,z,x1,y1,z1,x2,y2,z2,...
x3,y3,z3,x4,y4,z4,Gamma)
%Calculates the velocity on an arbitrary point (x,y,z) in space induced
%by a vortex ring of corners (xi,yi,zi) with i = 1:4 in that order and
%of strength Gamma.
[u1,v1,w1] = VortexSegmentInducedVelocity(x,y,z,x1,y1,z1,x2,y2,z2,Gamma);
[u2,v2,w2] = VortexSegmentInducedVelocity(x,y,z,x2,y2,z2,x3,y3,z3,Gamma);
[u3,v3,w3] = VortexSegmentInducedVelocity(x,y,z,x3,y3,z3,x4,y4,z4,Gamma);
[u4,v4,w4] = VortexSegmentInducedVelocity(x,y,z,x4,y4,z4,x1,y1,z1,Gamma);
u = u1+u2+u3+u4;
v = v1+v2+v3+v4;
w = w1+w2+w3+w4;
end

Finally, we get to:

 function [u,v,w] = VortexSegmentInducedVelocity(x,y,z,x1,y1,z1,x2,y2,z2,Gamma)
%Calculates the velocity (u,v,w) induced on an arbitrary point (x,y,z)
%in space by a vortex segment from (x1,y1,z1) to (x2,y2,z2) of strength
%Gamma.
%First we define the vectors r:
r0 = [x2 - x1,y2 - y1,z2 - z1];
r1 = [x - x1,y - y1,z - z1];
r2 = [x - x2,y - y2,z - z2];
%The following if statement measures if the point of induced velocity
%is too close to the vortex, and if it is, in order to avoid values
%going to infinity, returns zero. Otherwise, the induced velocity is:
if norm(cross(r0,r1)) > 0.001*norm(r0)*norm(r1);
    q = (Gamma/(4*pi))*(cross(r1,r2)/(norm(cross(r1,r2))^2))*dot(r1/norm(r1) - r2/norm(r2), r0);
    u = q(1);
    v = q(2);
    w = q(3);
else
    u = 0;
    v = 0;
    w = 0;
end
end

By now, I'm using cross products, and it is precisely here that things started to go awry, as when I tried to change syntax from loops to vectorization, the inputs of cross weren't 1 by 3 arrays. Anyway, I'm sorry I didn't explain my problem to full detail before, but as you can see it turns to be rather lengthy, and I thought I could solve it with a simple implementation of a "vectorized" statement. Perhaps there is a way to vectorize the dot and cross functions? Or would I have to write those functions by hand?

P.S. Giving this question some further thought, I guess the most worrying part is not the dot or cross functions, but rather the if statement inside VortexSegmentInducedVelocity. Since I have to do a verification for each element of the matrices, would this render the vectorization approach useless? In other words, if I were to rewrite VortexLineInducedVelocity to support matrices as inputs, would it be possible to do an if statement that also supports matrices as inputs and then operates only on each element for which its requirements are met?

[Stephen23]

The usual way to resolve if-like decisions in vectorized code is to use indexing (logical indexing is the fastest), typically either:

• pick an array susbset using indices, apply operation to it, replace into a subset of the array. The rest of the array remain unchanged.
• apply an operator to the entire array, then use indices to replace a subset of elements with zero/NaN/...

Also note that both cross and dot can be applied to multidimensional arrays, so you can vectorize these operations too. Read the documentation to see how this works: I would highly recommend that you use the syntax using the optional dimension argument, like this example from the documentation:

>> A = cat(3,[1 1;1 1],[2 3;4 5],[6 7;8 9]);
>> B = cat(3,[2 2;2 2],[10 11;12 13],[14 15; 16 17]);
>> C = dot(A,B,3)
C =
 106   140
 178   220

"The result, C, contains four separate dot products. The first dot product, C(1,1) = 106, is equal to the dot product of A(1,1,:) with B(1,1,:)."

Both of these changes will mean going through your entire code and rewriting it all... but is seems like it might be possible. Good luck!

[André Camões]

Okay, I think I know how to do it now. Thank you very much, sir!

Answer 1

[Walter Roberson]

I_ww = arrayfun(@InfluenceLine, XcWing, YcWing, ZcWing, XvWing, YvWing, ZvWing, XnWing, YnWing, ZnWing, AngleOfAttack);

provided that all of the arguments are either scalars or arrays the same size as each other.

This will not necessarily be any faster than looping, and typically would be slower.

It would be far better to write the function to accept arrays.