@stanto As you mention, Maple considers the names to be possibly nonreal here, unless qualified under assumptions.

Compare what you are now asking with the results from,

simplify(e,radical) assuming a::real, Db::real;

and

simplify(e,radical) assuming real;

@stanto As you mention, Maple considers the names to be possibly nonreal here, unless qualified under assumptions.

Compare what you are now asking with the results from,

simplify(e,radical) assuming a::real, Db::real;

and

simplify(e,radical) assuming real;

@majd This is an interesting question, and sorry that I could not respond earlier because of vacation and the recent security lock-out of the site (my valid cookie was expired/clobbered). Presumably this is not something where pdsolve/numeric could be used instead of ArrayInterpolation. It might be possible to compute the finite differences (which fdiff does) oneself, and to apply this in a Compiled procedure after first using ArrayInterpolation to quickly generate the needed subgrids of evaluation points in a fast vectorized manner. Easier for 1D than for the 2D case. I'm afraid that I'll likely be too busy to give it a try, for some weeks.

To get a smoother (ie. not piecewise flat) 2nd derivative from ArrayInterpolation might mean using method=spline and degree=4 or higher. The higher and higher the order used for the spline then the worse might become the 2nd deriv approximations for points farther and farther from the boundary.

@majd This is an interesting question, and sorry that I could not respond earlier because of vacation and the recent security lock-out of the site (my valid cookie was expired/clobbered). Presumably this is not something where pdsolve/numeric could be used instead of ArrayInterpolation. It might be possible to compute the finite differences (which fdiff does) oneself, and to apply this in a Compiled procedure after first using ArrayInterpolation to quickly generate the needed subgrids of evaluation points in a fast vectorized manner. Easier for 1D than for the 2D case. I'm afraid that I'll likely be too busy to give it a try, for some weeks.

To get a smoother (ie. not piecewise flat) 2nd derivative from ArrayInterpolation might mean using method=spline and degree=4 or higher. The higher and higher the order used for the spline then the worse might become the 2nd deriv approximations for points farther and farther from the boundary.

@Markiyan Hirnyk The suggestion to use rhs(sol[1]) instead of eval(k[1],sol) is not so widely reliable, because in older versions of Maple (or in recent Maple if launched with --setsort=0) using rhs(sol[1]) may produce an incorrect answer.

The reason is that in some older Maple versions (or with --setsort=0) the elements of the set sol could be in a session dependent order. The k[2]=... could be sol[1], so sol[1] is not generally reliable as a means to extract the value in the equation k[1]=...

Eg, in Maple 11.02 (or Maple 16.0 launched with --setsort=0) I see,

> restart:
> k[2]: k[2]=666.0: k[1]=5.0:
> sol:={k[1]=5.0,k[2]=666.0};
                       sol := {k[2] = 666.0, k[1] = 5.0}

> rhs(sol[1]);
                                     666.0

> restart:
> sol:={k[1]=5.0,k[2]=666.0};
                       sol := {k[1] = 5.0, k[2] = 666.0}

> rhs(sol[1]);
                                      5.0

@Markiyan Hirnyk The suggestion to use rhs(sol[1]) instead of eval(k[1],sol) is not so widely reliable, because in older versions of Maple (or in recent Maple if launched with --setsort=0) using rhs(sol[1]) may produce an incorrect answer.

The reason is that in some older Maple versions (or with --setsort=0) the elements of the set sol could be in a session dependent order. The k[2]=... could be sol[1], so sol[1] is not generally reliable as a means to extract the value in the equation k[1]=...

Eg, in Maple 11.02 (or Maple 16.0 launched with --setsort=0) I see,

> restart:
> k[2]: k[2]=666.0: k[1]=5.0:
> sol:={k[1]=5.0,k[2]=666.0};
                       sol := {k[2] = 666.0, k[1] = 5.0}

> rhs(sol[1]);
                                     666.0

> restart:
> sol:={k[1]=5.0,k[2]=666.0};
                       sol := {k[1] = 5.0, k[2] = 666.0}

> rhs(sol[1]);
                                      5.0

@dj_gssst You are using sol[1] and sol[2] inappropriately.

Look at what and are.

sol[1];

         k[1] = 0.001637935477698018732276851530308612938004131156280124580147547622858728899732152932156471141157968415

sol[2];

         k[2] = 0.8738880449255668655493761487129800277295707224229803653538821365150932040118518225052273427259289080

Those are both equations, not names or scalar values.

If you want to use the name k[1] then just use it, not sol[1]. If you want to use the value in that equation returned by fsolve then use eval(k[1],sol), and not sol[1].

It's somewhat irritating to wait for this site to load all the image of 2D Math that you've inlined into your posts. It would be quicker and nicer (in future) if you could please just type in your comments here, along with the link to your uploaded worksheets heavy with 2D Math.

@dj_gssst You are using sol[1] and sol[2] inappropriately.

Look at what and are.

sol[1];

         k[1] = 0.001637935477698018732276851530308612938004131156280124580147547622858728899732152932156471141157968415

sol[2];

         k[2] = 0.8738880449255668655493761487129800277295707224229803653538821365150932040118518225052273427259289080

Those are both equations, not names or scalar values.

If you want to use the name k[1] then just use it, not sol[1]. If you want to use the value in that equation returned by fsolve then use eval(k[1],sol), and not sol[1].

It's somewhat irritating to wait for this site to load all the image of 2D Math that you've inlined into your posts. It would be quicker and nicer (in future) if you could please just type in your comments here, along with the link to your uploaded worksheets heavy with 2D Math.

@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).