@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

Someone asked me what makes the final color red, since changing the color=red to color=green in the line of code which creates the initial plot (and assigns it to `p`) does not affect the final result of a point-plot with a red color gradient.

The color=red in the initial call,

    p:=plots:-pointplot(xy,color=red,symbolsize=4):

is only there so that the PLOT structure has a COLOR section which can be replaced. I probably ought to have chosen some other color and mentioned that the particular choice -- at this stage -- is irrelevent as far as the goal goes.

The actual coloring of the result was determined by the originally posted bit of code,

   (i)->`if`(irem(i,3)=1,xy[(i+2)/3,2],0)

in the Array constructor call which creates Array `c`. And n the ArrayTools:-Copy example, the final color was determined by choosing a stride which only copied every third entry, starting with the first (ie. the "red" parts of each triple).

This post is just an exposition, to show that the float[8] Array can be used to set a coloring scheme. It allows for all kinds of schemes, where rgb or hue can change gradually in any direction, or even along a conceptualized curve in space.

If someone were to turn the process into a robust procedure, which could be applied after the fact to any pointplot PLOT strucure, then it would test for the presence of a COLOR section. If it already existed then it would be replaced, and if it did not already exist then it would get added.

In much older Maple releases the location data itself of a 2D plot could not be supplied as  float[8] rtables. Then, for some release(s), such rtables (eg. Matrix, Array) were used internally but not acceptable as arguments to the plot commands. That's around the time frame at which Robert Israel noticed that a PLOT/PLOT3D structure could be manually formed with such rtables, and his fast plotters in his Maple Advisor database ensued. And in recent Maple versions many plot commands have been changed to accept their data argument as such float[8] rtables. But that's just the point location data. It's a natural progression that the color-specification data -- which has become float[8] internally -- would finally be acceptable as input as well. That's my opinion, anyway.

acer

@Jimmy First thing: please don't repost this as a duplicate question. We've all seen it here, and it won't help to ask it over again.

When I asked you earlier whether variable `i` was supposed to be dependent or independent I had in mind that it might be considered as independent. (I realize that you answered the opposite.)

That is, suppose you took your equation i = i0*(exp((1000*(v-i*rs))/(25.9*n))-1) and got from that form i0*(exp((1000*(v-i*rs))/(25.9*n))-1) - i = 0.

Then suppose that the dependent variable (we don't even have to name it...) has value zero for every measured pair of values of `v` and `i` the independent variables.

In this view, the independent data can be taken as the m-by-2 Matrix formed by concatenating the column Vectors `V` and `idata`. In Maple syntax, something like,

dim := LinearAlgebra:-Dimension(V);

NonlinearFit( i0*(exp((1000*(v-i*rs))/(25.9*n))-1) - i, 
              < V | idata >,
              Vector(dim),
              [v,i],
              parameterranges=[...]
              output=[...] );

But I then ran into the problem that your stated parameter-range for `rs` could not be right. For your "shortened" version of the equation the parameter rs might have to be quite different than 1..500.

And since the variables and parameters can vary by a factor of about 1e6 then raising the working precision (in Maple, the environment setting `Digits`) to about 20 or 30 seems like not a bad idea.

So, how about using the full equation, with this approach? By full equation I mean use all of the terms of the right-hand-side equation that you showed at the very beginning of this post.

There may have some more constant parameter there, whose value you haven't yet given us, like R_SH, n_R, and I_SR. Perhaps by using the proper full formula the discrepancy with your stated range for parameter `rs` would be resolved.

@Jimmy First thing: please don't repost this as a duplicate question. We've all seen it here, and it won't help to ask it over again.

When I asked you earlier whether variable `i` was supposed to be dependent or independent I had in mind that it might be considered as independent. (I realize that you answered the opposite.)

That is, suppose you took your equation i = i0*(exp((1000*(v-i*rs))/(25.9*n))-1) and got from that form i0*(exp((1000*(v-i*rs))/(25.9*n))-1) - i = 0.

Then suppose that the dependent variable (we don't even have to name it...) has value zero for every measured pair of values of `v` and `i` the independent variables.

In this view, the independent data can be taken as the m-by-2 Matrix formed by concatenating the column Vectors `V` and `idata`. In Maple syntax, something like,

dim := LinearAlgebra:-Dimension(V);

NonlinearFit( i0*(exp((1000*(v-i*rs))/(25.9*n))-1) - i, 
              < V | idata >,
              Vector(dim),
              [v,i],
              parameterranges=[...]
              output=[...] );

But I then ran into the problem that your stated parameter-range for `rs` could not be right. For your "shortened" version of the equation the parameter rs might have to be quite different than 1..500.

And since the variables and parameters can vary by a factor of about 1e6 then raising the working precision (in Maple, the environment setting `Digits`) to about 20 or 30 seems like not a bad idea.

So, how about using the full equation, with this approach? By full equation I mean use all of the terms of the right-hand-side equation that you showed at the very beginning of this post.

There may have some more constant parameter there, whose value you haven't yet given us, like R_SH, n_R, and I_SR. Perhaps by using the proper full formula the discrepancy with your stated range for parameter `rs` would be resolved.

To be honest I did not give optimizing this as much consideration as I might have.

Partly that is because it should be utterly clear by this date in time that a float[8] Array will soundly beat a list (or sequence) of floats at these sizes, both in terms of speed and memory.

But there is so much more the the float[8] Array brings with it. It allows the color data to be held in a truly mutable data structure, which means it can be used and re-used with minimal penalty of memory managment. (...animation comes to mind.) There are a wealth of Library and kernel routines which have themselves been optimized for handling float[8]. And it will always take proportionally less memory allocation. (The only single thing that could take less memory while offering the color depth is float[4], and there are less routines optimized speedwise for that.)

Let's compare the two previous suggestions (both faster than my original) with one more faster way, which is to use ArrayTools:-Copy.

I have suppressed display of the plot, in the timing runs below, because I found that repeatedly rerunning the whole set of code would produce somewhat random timings. In general, however, I thought that I saw timings like what follows below, but on average with a constant addition of the time it takes for the GUI to display the plot. Hence, timings are without the plots being rendered on the worksheet.  (32bit Windows below. In 64bit the memory size of Herclau's method rises to abotu 47.0 MB, and of the other to about 3.3 MB.)

restart:
N:=25000:
xy:=LinearAlgebra:-RandomMatrix(N,2,generator=0.0..1.0,
                                outputoptions=[datatype=float[8]]):
str:=time[real]():
plots:-pointplot(xy,
                    color=[seq(RGB(xy[i, 2], 0., 0.), i = 1 .. N)],
                    symbolsize=4):
time[real]()-str;

                             0.327

kernelopts(bytesalloc);

                            18150148

restart:
N:=25000:
xy:=LinearAlgebra:-RandomMatrix(N,2,generator=0.0..1.0,
                                outputoptions=[datatype=float[8]]):
str:=time[real]():
f:=proc(A::Array(datatype=float[8],order=C_order),
        xy::Matrix(datatype=float[8]), n)
 local i;
 for i from 1 to n*3 do
  A[i]:=`if`(irem(i,3)=1,xy[(i+2.0)/3.0,2],0);
 end do:
 NULL:
end proc:
p:=plots:-pointplot(xy,color=red,symbolsize=4):
c:=Array(1..N*3, datatype=float[8],order=C_order):
evalhf(f(c,xy,N)):
subsindets(p,specfunc(anything,COLOUR),z->'COLOUR'('RGB',c)):
time[real]()-str;

                             0.046

kernelopts(bytesalloc);

                            2293340

restart:
N:=25000:
xy:=LinearAlgebra:-RandomMatrix(N,2,generator=0.0..1.0,
                                outputoptions=[datatype=float[8]]):
str:=time[real]():
p:=plots:-pointplot(xy,color=red,symbolsize=4):
c:=Array(1..N*3, datatype=float[8],order=C_order):
ArrayTools:-Copy(N,xy,N,1,c,0,3);
subsindets(p,specfunc(anything,COLOUR),z->'COLOUR'('RGB',c)):
time[real]()-str;

                             0.016

kernelopts(bytesalloc);

                            2293340

acer

@PatrickT The 2D Math in the uploaded worksheet (mirrored in 1D by Patrick) does not look right for initial conditions that specify derivatives evaluated at particular points.

Eg, diff(t1(0), y) = -(diff(t2(0),y)) becomes 0=0.

What's likely intended is D(t1)(0) = - D(t2)(0) or some equivalent using eval and diff.

But when I try such revisions, Maple 16 issues an error message, stating that it cannot convert to an explicit 1st order system. Maybe someone else could check that.