@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

Looking down on the 3D point-plot in the negative-My direction, the view projected onto the Mx-N plane looks like this:

plots[display](seq(pl[j], j = 1 .. 20), orientation = [-90, -90, 0])

A quick check of the SOL data indicates that the points are a set of 20 groups, where the N value is the same for each point in any of the groups. In other words, they are arranged in a grid that is uniform in the N direction. If a mapping in the Mx direction were devised so that -- for each group -- the newMx values agreed then one would have a nice regular grid in the newMx-N plane. That would allow very easy construction of an interpolating procedure F (using ArrayInterpolation) which could return the My "height" for any given point in the newMx-N plane.

Once that mechanism is at hand, then it's easy to see whether any arbitary point in the octant lies above (in the My sense) of below that interpolated surface. Given point [mx,my,n] then use the mapping to get [newmx,my,n], then use the interpolating function F to return the surface height sy = F(newmx,n). That either greater than (above) `my`, or not.

Now consider construction of the mapping which first scales in the Mx direction. It appears to be a function only of N-value, since the spacing in the Mx direction look uniform (up to 7 digits or so).

> seq(SOL[1, 1][j+1]-SOL[1, 1][j], j = 1 .. 9);

          -4.96317706, -4.96317707, -4.96317706, -4.96317707, -4.96317707, 
          -4.96317706, -4.963177068, -4.963177066, -4.963177066

It should thus be possible to construct such an Mx-scaling mapping by interpolating w.r.t the N values of the 20 sets.

Got to go now...

acer

Looking down on the 3D point-plot in the negative-My direction, the view projected onto the Mx-N plane looks like this:

plots[display](seq(pl[j], j = 1 .. 20), orientation = [-90, -90, 0])

A quick check of the SOL data indicates that the points are a set of 20 groups, where the N value is the same for each point in any of the groups. In other words, they are arranged in a grid that is uniform in the N direction. If a mapping in the Mx direction were devised so that -- for each group -- the newMx values agreed then one would have a nice regular grid in the newMx-N plane. That would allow very easy construction of an interpolating procedure F (using ArrayInterpolation) which could return the My "height" for any given point in the newMx-N plane.

Once that mechanism is at hand, then it's easy to see whether any arbitary point in the octant lies above (in the My sense) of below that interpolated surface. Given point [mx,my,n] then use the mapping to get [newmx,my,n], then use the interpolating function F to return the surface height sy = F(newmx,n). That either greater than (above) `my`, or not.

Now consider construction of the mapping which first scales in the Mx direction. It appears to be a function only of N-value, since the spacing in the Mx direction look uniform (up to 7 digits or so).

> seq(SOL[1, 1][j+1]-SOL[1, 1][j], j = 1 .. 9);

          -4.96317706, -4.96317707, -4.96317706, -4.96317707, -4.96317707, 
          -4.96317706, -4.963177068, -4.963177066, -4.963177066

It should thus be possible to construct such an Mx-scaling mapping by interpolating w.r.t the N values of the 20 sets.

Got to go now...

acer