Could it be that B hasn't actually been assigned a matrix or Matrix?

acer

If you are using Maple 12, then the command Explore(%) in the code above should pop-up a box. It's a tool to explore expressions or plots with parameters in them.

acer

If you are using Maple 12, then the command Explore(%) in the code above should pop-up a box. It's a tool to explore expressions or plots with parameters in them.

acer

Here's a possible way to visualize it.

with( LinearAlgebra ):
A := <<1|1>,<4|1>>;
B := MatrixExponential(A,t);
C := Vector(2,[c1,c2]);
M := B.C;

It seems that you want a parameterized plot, with x being M[1] and y being M[2]. So you could try something like this.

When the Explore set-up window pops up, change the lower and upper values for both c1 and c2 to -3.0 and 3.0 say. Check the boxes to skip t and view.

'plot'([M[1], M[2],t=-10..10],view=[-10..10,-10..10]):
Explore(%);

You could also make that 'plot' call with t=a..b, to get sliders which controlled the length. Set the slider range for a from -2.0 to 0.0 and for b from 0.0 to 2.0 say.

acer

Here's a possible way to visualize it.

with( LinearAlgebra ):
A := <<1|1>,<4|1>>;
B := MatrixExponential(A,t);
C := Vector(2,[c1,c2]);
M := B.C;

It seems that you want a parameterized plot, with x being M[1] and y being M[2]. So you could try something like this.

When the Explore set-up window pops up, change the lower and upper values for both c1 and c2 to -3.0 and 3.0 say. Check the boxes to skip t and view.

'plot'([M[1], M[2],t=-10..10],view=[-10..10,-10..10]):
Explore(%);

You could also make that 'plot' call with t=a..b, to get sliders which controlled the length. Set the slider range for a from -2.0 to 0.0 and for b from 0.0 to 2.0 say.

acer

Just to be more clear, I'm not asking about powering float Matrices, per se.

Other plausible techniques are known. Eg. since "almost all" float Matrices are diagonalizable, the product (P^(-1).D.P)^n telescopes down to P^(-1).D^n.P and  only  the scalar  eigenvalues  would be powered.

I'm wondering more about the prospects for optimization following enhancement for iterative specialization and, perhaps more difficult, handling the in-place semantic.

acer

I thoroughly enjoyed reading the thesis, and I look forward to poking at the code.

I have just one or two quite specific questions, which I'd like to ask now if I may. In subsection 7.2 there is an example on binary powering (Listings 7.7 and 7.8). In the specialized code of Listing 7.8 there are more local variables than are strictly necessary. Ok, so locals are not too expensive, in themselves. But what about optimizing for the least number of such locals? Is that something that might fall under the possible future work of "Iterative Specialization" as mentioned in subsection 8.3.2? I ask because I'm thinking of the case of powering a Matrix rather than an scalar expression. In the hardware datatype case, the entries of Matrices are never garbage, but their container Matrices are. So the more locals that get produced (and assigned to, rather than copied to in-place), the more expensive it is.

I'll try to explain by example. Consider the binary powering of hardware (any, actually) Matrices as done in Maple 9.5. It is akin to Listing 7.7, and it produces whole new Matrices galore. And it's really slow, on account of all the ensuing garbage collection.

    |\^/|     Maple 9.5 (IBM INTEL LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2004
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> interface(verboseproc=3):
> eval(`rtable/Power/Power`);
proc(M, n)
local P;
option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`;
    if n = 1 then M
    elif n = 2 then LinearAlgebra:-MatrixMatrixMultiply(M, M)
    else
        P := `rtable/Power/Power`(M, iquo(n, 2));
        P := LinearAlgebra:-MatrixMatrixMultiply(P, P);
        if n mod 2 = 1 then
            LinearAlgebra:-MatrixMatrixMultiply(P, M, 'inplace' = true)
        end if;
        P
    end if
end proc

Now look at the current implementation. I won't show the whole listing, but it can be described easily enough.

interface(verboseproc=3);
eval(`rtable/Power/float`);

It has distinct "specializations" for powers 1,2,3,4, and 5. And each of those minimizes the number of allocated hardware Matrices. Matrices are never assigned to locals except when first created. No other instances of assigning a Matrix to a variable occurs, and all updates are done by in-place copying. And then it has a version for powers 6 and higher, which follows a pattern (common for any p>6) and thus uses a common scheme minimized for the total number of allocated Matrices.

What I'm wondering is, is there any hope of generating programs such as that using MapleMix techniques (present, or with any future idea from section 8)? I would like to create such optimized programs, but I don't want to have to work it out by hand.

acer

Just a minor point: that exact syntax has been available for many releases now.

What's new in Maple 12 (I think) is this,

> A := Array(1..10,1..10,1..10, (i,j,k)->i+j/k ):

> A[1,1,1..];
               [                                            11]
               [2, 3/2, 4/3, 5/4, 6/5, 7/6, 8/7, 9/8, 10/9, --]
               [                                            10]

> A[1,1,4..]:

> Vector[row](A[1,1,4..]);
                      [                               11]
                      [5/4, 6/5, 7/6, 8/7, 9/8, 10/9, --]
                      [                               10]

acer

Just a minor point: that exact syntax has been available for many releases now.

What's new in Maple 12 (I think) is this,

> A := Array(1..10,1..10,1..10, (i,j,k)->i+j/k ):

> A[1,1,1..];
               [                                            11]
               [2, 3/2, 4/3, 5/4, 6/5, 7/6, 8/7, 9/8, 10/9, --]
               [                                            10]

> A[1,1,4..]:

> Vector[row](A[1,1,4..]);
                      [                               11]
                      [5/4, 6/5, 7/6, 8/7, 9/8, 10/9, --]
                      [                               10]

acer

The try..catch mechanism allows one to control the program flow, even in the presence of an error.

For example,

> p := proc(x)
> local z;
> try
>   if 0 < x then error "whoops";
>   else z := x^2;
>   end if;
> catch "whoops":
>   print("found a whoops");
>   z := sin(z);
> end try;
> end proc:

> p(-4);
                                      16
 
> p(4);
                               "found a whoops"
 
                                    sin(z)

So, simply catch the particular error message which LPSolve issues when it's infeasible.

See the ?try help-page.

acer

The try..catch mechanism allows one to control the program flow, even in the presence of an error.

For example,

> p := proc(x)
> local z;
> try
>   if 0 < x then error "whoops";
>   else z := x^2;
>   end if;
> catch "whoops":
>   print("found a whoops");
>   z := sin(z);
> end try;
> end proc:

> p(-4);
                                      16
 
> p(4);
                               "found a whoops"
 
                                    sin(z)

So, simply catch the particular error message which LPSolve issues when it's infeasible.

See the ?try help-page.

acer

Am I just missing it, or is there no reference or link to ?Rounding from either the ?evalf or ?evalf,details help-pages? Also, maybe the ?Rounding page could benefit from an example such as,

> evalf[5](Pi);
                                    3.1416

> Rounding:=0:
> forget(evalf):
> evalf[5](Pi);
                                    3.1415

Yes, as Alec points out iquo (or irem) is usual. See also,

interface(verboseproc=3):
eval(`convert/base`);

acer

Am I just missing it, or is there no reference or link to ?Rounding from either the ?evalf or ?evalf,details help-pages? Also, maybe the ?Rounding page could benefit from an example such as,

> evalf[5](Pi);
                                    3.1416

> Rounding:=0:
> forget(evalf):
> evalf[5](Pi);
                                    3.1415

Yes, as Alec points out iquo (or irem) is usual. See also,

interface(verboseproc=3):
eval(`convert/base`);

acer

That's a good point, about the timing function. A more convenient version of Maple's time() command, which returns both the result and the elapsed time, would be useful.

As far as the speed goes, does anyone know how many pre-computed digits of Pi Mathematica has stored (hard-coded)? It's even more interesting if one knows that both systems had to compute the answer in full.

acer

That's a good point, about the timing function. A more convenient version of Maple's time() command, which returns both the result and the elapsed time, would be useful.

As far as the speed goes, does anyone know how many pre-computed digits of Pi Mathematica has stored (hard-coded)? It's even more interesting if one knows that both systems had to compute the answer in full.

acer