@jimmyinhmb There's a lot to consider here. I'll have to leave parallelization of it until another day, hope that's ok.

Your elementwise approach is mostly getting bogged down by memory management, I think. Consider the statement,

        G1[Mid..Last] := G1[Mid..Last] + exp~( c +~ H[Mid..Last] ):

I count at least five overt Vector workspaces being created there. Two to extract H[Mid..Last] and G1[Mid..Last] on the RHS. One more to do elementwise adddition of `c`, and another to do elementwise exponentiation. And then another to do the final addition (prior to recopying into G1).

And the second method below cuts that down to just two overt Vector workspaces created, and gets some improvement. Further improvement comes from getting rid of all iterated temp workspace creation, and acting purely inplace. (A re-usable workspace or two, created just once outside the loop can really help performance.) It seems natural to conclude that the improvements are largely due to avoiding repeated Vector creation and collection as garbage: ie. avoidable memory management. And, indeed, this aspect of the performance improvement seems to match the reduction in kernelopts(bytesused) which is one measure of the amount of memory management.

It might be work noting that simply putting the whole job into a single procedure that acts inplace on G, and merely Compiling that, is almost fastest. And yet it is one of the easiest to code.

The following output is from a run with Maple 16.01 on a Core2 Duo running Windows XP. I saw some variation with regard to the relative performance of methods 3-to-5 below when run on Windows 7 or Linux 64bit.

Download maptry3.mw

@jimmyinhmb There's a lot to consider here. I'll have to leave parallelization of it until another day, hope that's ok.

Your elementwise approach is mostly getting bogged down by memory management, I think. Consider the statement,

        G1[Mid..Last] := G1[Mid..Last] + exp~( c +~ H[Mid..Last] ):

I count at least five overt Vector workspaces being created there. Two to extract H[Mid..Last] and G1[Mid..Last] on the RHS. One more to do elementwise adddition of `c`, and another to do elementwise exponentiation. And then another to do the final addition (prior to recopying into G1).

And the second method below cuts that down to just two overt Vector workspaces created, and gets some improvement. Further improvement comes from getting rid of all iterated temp workspace creation, and acting purely inplace. (A re-usable workspace or two, created just once outside the loop can really help performance.) It seems natural to conclude that the improvements are largely due to avoiding repeated Vector creation and collection as garbage: ie. avoidable memory management. And, indeed, this aspect of the performance improvement seems to match the reduction in kernelopts(bytesused) which is one measure of the amount of memory management.

It might be work noting that simply putting the whole job into a single procedure that acts inplace on G, and merely Compiling that, is almost fastest. And yet it is one of the easiest to code.

The following output is from a run with Maple 16.01 on a Core2 Duo running Windows XP. I saw some variation with regard to the relative performance of methods 3-to-5 below when run on Windows 7 or Linux 64bit.

Download maptry3.mw

@Alejandro Jakubi By "computing with identifiers" I meant just using them only as names, not parsing them. Clearly typeset names which look the same would also have to be the same (internally), or else mixing them arithmetically and otherwise would not be so useful.

Thanks, about the brackets.

By trying a few tricks I got some hints that Typesetting:-Typeset may not get called at all when one does right-click conversion to Atomic Identifier. That makes getting the exact same structures problematic. It'd be easier to use the same basic mechanism, than to fix up the results from something else.

@AliKhan Is this what you mean?

A:=LinearAlgebra:-RandomMatrix(3,outputoptions=[shape=diagonal]);

                                  [27   0    0]
                                  [           ]
                             A := [ 0  -4    0]
                                  [           ]
                                  [ 0   0  -74]

B:=Matrix(3,3):

for i from 1 to 3 do B[i,i]:=A; end do:

B;

                [[27   0    0]                              ]
                [[           ]                              ]
                [[ 0  -4    0]        0              0      ]
                [[           ]                              ]
                [[ 0   0  -74]                              ]
                [                                           ]
                [               [27   0    0]               ]
                [               [           ]               ]
                [      0        [ 0  -4    0]        0      ]
                [               [           ]               ]
                [               [ 0   0  -74]               ]
                [                                           ]
                [                              [27   0    0]]
                [                              [           ]]
                [      0              0        [ 0  -4    0]]
                [                              [           ]]
                [                              [ 0   0  -74]]

B[2,2];

                                [27   0    0]
                                [           ]
                                [ 0  -4    0]
                                [           ]
                                [ 0   0  -74]

Most of regular LinearAlgebra will not be able to handle B, although arithmetic and some simpler things might work. (You might be able to teach LinearAlgebra[Generic] to do more...)

@AliKhan Is this what you mean?

A:=LinearAlgebra:-RandomMatrix(3,outputoptions=[shape=diagonal]);

                                  [27   0    0]
                                  [           ]
                             A := [ 0  -4    0]
                                  [           ]
                                  [ 0   0  -74]

B:=Matrix(3,3):

for i from 1 to 3 do B[i,i]:=A; end do:

B;

                [[27   0    0]                              ]
                [[           ]                              ]
                [[ 0  -4    0]        0              0      ]
                [[           ]                              ]
                [[ 0   0  -74]                              ]
                [                                           ]
                [               [27   0    0]               ]
                [               [           ]               ]
                [      0        [ 0  -4    0]        0      ]
                [               [           ]               ]
                [               [ 0   0  -74]               ]
                [                                           ]
                [                              [27   0    0]]
                [                              [           ]]
                [      0              0        [ 0  -4    0]]
                [                              [           ]]
                [                              [ 0   0  -74]]

B[2,2];

                                [27   0    0]
                                [           ]
                                [ 0  -4    0]
                                [           ]
                                [ 0   0  -74]

Most of regular LinearAlgebra will not be able to handle B, although arithmetic and some simpler things might work. (You might be able to teach LinearAlgebra[Generic] to do more...)

Apparently I have trouble remembering my own advice.

A better double loop in the Compilable BlockAdd procedure above would walk entries by column (and not by row) in the innermost loop because my example uses Fortran_order Matrices by default. Ie,

   for j from sxj to exj do   
      for i from sxi to exi do
         x[i,j]:=x[i,j]+y[i+syi-sxi,j+syj-sxj];
      end do;
   end do;

With that change, and by reducing some of the indexing integer-arithmetic, the Compiled version gets about 30% faster (but varying with platform). Still not as fast as the block-copy & daxpy. But it also might be threaded, after splitting the action by column (which I have not tried).

Apparently I have trouble remembering my own advice.

A better double loop in the Compilable BlockAdd procedure above would walk entries by column (and not by row) in the innermost loop because my example uses Fortran_order Matrices by default. Ie,

   for j from sxj to exj do   
      for i from sxi to exi do
         x[i,j]:=x[i,j]+y[i+syi-sxi,j+syj-sxj];
      end do;
   end do;

With that change, and by reducing some of the indexing integer-arithmetic, the Compiled version gets about 30% faster (but varying with platform). Still not as fast as the block-copy & daxpy. But it also might be threaded, after splitting the action by column (which I have not tried).

@jimmyinhmb It looks as if the block-copying & daxpy apporach is faster than the (serial) Compiled.

First, though, I apologize for pointing you at some wrong lines of code above. In Maple 15, the call out to f06ecf/daxpy would be more like this line in the Library,

showstat((LinearAlgebra::LA_External)::MatrixAdd,52);
       ...
  52     HWCall[2](lengthA,evalf(c2),localB,0,1,localA,0,1)
       ...

Here's a comparison, done on 64bit Maple on Windows 7,

restart: #randomize():

(m,n):=10^5,10: maxiter:=10^2:

M1:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):
M2:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):

DAXPY:=proc(num,c,x,offsetx,incx,y,offsety,incy)
   # http://www.netlib.org/blas/daxpy.f
   # http://en.wikipedia.org/wiki/SAXPY
   local extlib, extdaxpy;
   extlib:=ExternalCalling:-ExternalLibraryName("linalg");
   extdaxpy:=ExternalCalling:-DefineExternal('hw_f06ecf',extlib);
   extdaxpy(num,evalf(c),x,offsetx,incx,y,offsety,incy);
   NULL;
end proc:

(tm,tn):=floor(m/2),n:
T1:=Matrix(tm,tn,datatype=float[8]):
T2:=Matrix(tm,tn,datatype=float[8]):

st,str:=time(),time[real]():
for iter from 1 to maxiter do
  ArrayTools:-BlockCopy(M1,m-tm,m,tm,tn,T1,m-tm);
  ArrayTools:-BlockCopy(M2,m-tm,m,tm,tn,T2,m-tm);
  DAXPY(tm*tn,1.0,T2,0,1,T1,0,1);
  ArrayTools:-BlockCopy(T1,0,1,1,tm*tn,M1,m-tm,m,tm,tn);
end do:
time()-st,time[real]()-str;

                          0.873, 0.874

restart: #randomize():

(m,n):=10^5,10: maxiter:=10^2:

M1:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):
M2:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):

BlockAdd:=proc(x::Matrix(datatype=float[8]),
               sxi::integer[4],exi::integer[4],
               sxj::integer[4],exj::integer[4],
               y::Matrix(datatype=float[8]),
               syi::integer[4],syj::integer[4])
   local i::integer[4], j::integer[4];
   for i from 0 to exi-sxi do
      for j from 0 to exj-sxj do
         x[i+sxi,j+sxj]:=x[i+sxi,j+sxj]+y[i+syi,j+syj];
      end do;
   end do;
   NULL;
end proc:

try
  cBlockAdd:=Compiler:-Compile(BlockAdd):
  compiled:=true:
catch:
end try:

if compiled then

st,str:=time(),time[real]():
for iter from 1 to maxiter do
   cBlockAdd(M1,floor(m/2)+1,m,1,n,M2,floor(m/2)+1,1);
end do:
print(time()-st,time[real]()-str);

end if:

                          2.746, 2.737

Interestingly, the MKL did not appear to use more than a single core, when running daxpy, on this 4-physical core (8-hyperthreaded) i7. That seems borne out by the total-cpu time vs real-time ration, as well as by not exceeding 12% Usage in Windows Task Manager. But MKL does use more cores for other tasks, like Matrix-Matrix multiplication. Maybe Intel has not yet threaded all level-2 BLAS. I don't know.

Hopefully, the block-copy approach will be fast enough.

acer

@jimmyinhmb It looks as if the block-copying & daxpy apporach is faster than the (serial) Compiled.

First, though, I apologize for pointing you at some wrong lines of code above. In Maple 15, the call out to f06ecf/daxpy would be more like this line in the Library,

showstat((LinearAlgebra::LA_External)::MatrixAdd,52);
       ...
  52     HWCall[2](lengthA,evalf(c2),localB,0,1,localA,0,1)
       ...

Here's a comparison, done on 64bit Maple on Windows 7,

restart: #randomize():

(m,n):=10^5,10: maxiter:=10^2:

M1:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):
M2:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):

DAXPY:=proc(num,c,x,offsetx,incx,y,offsety,incy)
   # http://www.netlib.org/blas/daxpy.f
   # http://en.wikipedia.org/wiki/SAXPY
   local extlib, extdaxpy;
   extlib:=ExternalCalling:-ExternalLibraryName("linalg");
   extdaxpy:=ExternalCalling:-DefineExternal('hw_f06ecf',extlib);
   extdaxpy(num,evalf(c),x,offsetx,incx,y,offsety,incy);
   NULL;
end proc:

(tm,tn):=floor(m/2),n:
T1:=Matrix(tm,tn,datatype=float[8]):
T2:=Matrix(tm,tn,datatype=float[8]):

st,str:=time(),time[real]():
for iter from 1 to maxiter do
  ArrayTools:-BlockCopy(M1,m-tm,m,tm,tn,T1,m-tm);
  ArrayTools:-BlockCopy(M2,m-tm,m,tm,tn,T2,m-tm);
  DAXPY(tm*tn,1.0,T2,0,1,T1,0,1);
  ArrayTools:-BlockCopy(T1,0,1,1,tm*tn,M1,m-tm,m,tm,tn);
end do:
time()-st,time[real]()-str;

                          0.873, 0.874

restart: #randomize():

(m,n):=10^5,10: maxiter:=10^2:

M1:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):
M2:=LinearAlgebra:-RandomMatrix(m,n,outputoptions=[datatype=float[8]]):

BlockAdd:=proc(x::Matrix(datatype=float[8]),
               sxi::integer[4],exi::integer[4],
               sxj::integer[4],exj::integer[4],
               y::Matrix(datatype=float[8]),
               syi::integer[4],syj::integer[4])
   local i::integer[4], j::integer[4];
   for i from 0 to exi-sxi do
      for j from 0 to exj-sxj do
         x[i+sxi,j+sxj]:=x[i+sxi,j+sxj]+y[i+syi,j+syj];
      end do;
   end do;
   NULL;
end proc:

try
  cBlockAdd:=Compiler:-Compile(BlockAdd):
  compiled:=true:
catch:
end try:

if compiled then

st,str:=time(),time[real]():
for iter from 1 to maxiter do
   cBlockAdd(M1,floor(m/2)+1,m,1,n,M2,floor(m/2)+1,1);
end do:
print(time()-st,time[real]()-str);

end if:

                          2.746, 2.737

Interestingly, the MKL did not appear to use more than a single core, when running daxpy, on this 4-physical core (8-hyperthreaded) i7. That seems borne out by the total-cpu time vs real-time ration, as well as by not exceeding 12% Usage in Windows Task Manager. But MKL does use more cores for other tasks, like Matrix-Matrix multiplication. Maybe Intel has not yet threaded all level-2 BLAS. I don't know.

Hopefully, the block-copy approach will be fast enough.

acer

@Markiyan Hirnyk The matter of whether to extrapolate the surface or to introduce new cutting planes, so as to obtain a closed region, is relatively trivial to the task of constructing a procedure that can determine the interpolated surface.

The second question asked was, "how to define whether the arbitrary point (x,y,z) is inside this surface or not?" Another interpretation is that he wishes to determine whether a selected point lies approximately on the surface (as opposed to meaning within some region bounded in part by the surface).

Only minor changes to the worksheet I uploaded are necessary, to address one of these interpretations versus the other. The key part is the procedure which interpolates the surface, between the data points.

He can test that func(n,mx)=my approximately, to test whether point (mx,my,n) is on the surface. Or he can bound a region (by extrapolation or introducing planes, or whatever) and test using an inequality. That's still up to him, and is relatively easy compared to procedurally generating the approximating surface.

For example, one could use extrapolation=undefined instead of extrapolation=true in the call to ArrayInterpolation within the defn of `S`. After which a plot of the points and surface (at grid=[100,100]) can look like this:

The procedure `func` may be more useful than the mere plot returned by `surfdata`, for answering his question 2).

I personally found it interesting because of the method of mapping to a regular grid, which I had not done before by can imagine using again. So I learned something.

acer

@Markiyan Hirnyk The matter of whether to extrapolate the surface or to introduce new cutting planes, so as to obtain a closed region, is relatively trivial to the task of constructing a procedure that can determine the interpolated surface.

The second question asked was, "how to define whether the arbitrary point (x,y,z) is inside this surface or not?" Another interpretation is that he wishes to determine whether a selected point lies approximately on the surface (as opposed to meaning within some region bounded in part by the surface).

Only minor changes to the worksheet I uploaded are necessary, to address one of these interpretations versus the other. The key part is the procedure which interpolates the surface, between the data points.

He can test that func(n,mx)=my approximately, to test whether point (mx,my,n) is on the surface. Or he can bound a region (by extrapolation or introducing planes, or whatever) and test using an inequality. That's still up to him, and is relatively easy compared to procedurally generating the approximating surface.

For example, one could use extrapolation=undefined instead of extrapolation=true in the call to ArrayInterpolation within the defn of `S`. After which a plot of the points and surface (at grid=[100,100]) can look like this:

The procedure `func` may be more useful than the mere plot returned by `surfdata`, for answering his question 2).

I personally found it interesting because of the method of mapping to a regular grid, which I had not done before by can imagine using again. So I learned something.

acer

Here's a worksheet which runs through the idea of using a procedure which interpolates S=(N,Mx)->My, where another scaling procedure Xadj=N->scalingfactor is first applied to the second argument of S. These are combined into a procedure `func`=(b,a)->S(b,Xadj(b)*a) which can be plotted as a surface.

It extrapolates from N=1285..1750 (which leaves something to be desired) and from N=0..64 (which is ok), in order to get a region for which the Asker wanted to test arbitrary (mx,my,n) points as to their inclusion.

LC_for_MP2.mw

Mapleprimes failed to inline that worksheet's contents into this Comment. But here is the final plot, which is produced using `plot3d` and procedure `func`, rather than `surfdata`.

acer

Here's a worksheet which runs through the idea of using a procedure which interpolates S=(N,Mx)->My, where another scaling procedure Xadj=N->scalingfactor is first applied to the second argument of S. These are combined into a procedure `func`=(b,a)->S(b,Xadj(b)*a) which can be plotted as a surface.

It extrapolates from N=1285..1750 (which leaves something to be desired) and from N=0..64 (which is ok), in order to get a region for which the Asker wanted to test arbitrary (mx,my,n) points as to their inclusion.

LC_for_MP2.mw

Mapleprimes failed to inline that worksheet's contents into this Comment. But here is the final plot, which is produced using `plot3d` and procedure `func`, rather than `surfdata`.

acer

@Markiyan Hirnyk Of course I cannot explain with 100% surety what the Asker wants there, and it is not covered by what information has been given so far.

I imagine that it might mean that he wants to extrapolate the surface, or else possibly that he wants a hard drop to the Mx-N and My-N planes as another boundary at that N-value.

There is a similar ambiguity near My=0.

But I figure that accomodating either of those two interpretations would be pretty straightforward provided that a general scheme to handle all the other points in the octant has been obtained.

acer

@Markiyan Hirnyk Of course I cannot explain with 100% surety what the Asker wants there, and it is not covered by what information has been given so far.

I imagine that it might mean that he wants to extrapolate the surface, or else possibly that he wants a hard drop to the Mx-N and My-N planes as another boundary at that N-value.

There is a similar ambiguity near My=0.

But I figure that accomodating either of those two interpretations would be pretty straightforward provided that a general scheme to handle all the other points in the octant has been obtained.

acer