Dave Linder

## 517 Reputation

18 years, 190 days
Maplesoft
Software Architect

Dave Linder
Mathematical Software, Maplesoft

## CLAPACK zgeevx...

If you are using the LinearAlgebra package in Maple 10 then it's likely that the eigenvalues of your complex Matrix are being computed by the zgeevx function in the CLAPACK library. The above will be true if the Matrix entries are all purely numeric and there is some floating-point value present in them. You can set, infolevel[LinearAlgebra]:=1: to see the name of the external function being used during the LinearAlgebra computation. If Digits < evalhf(Digits) and UseHardwareFloats is not equal to 'false' then a double-precision CLAPACK function will be used. That happens with the default precision settings. For extended precision behaviour, see the help-page ?UseHardwareFloats . LinearAlgebra routines such as Eigenvectors() obey that flag. If the hardware double precision version is used then the external function name will be prefixed by hw_, while in the software multiprecision case it will be prefixed by sw_ . An example might look like this, > restart: > with(LinearAlgebra): > infolevel[LinearAlgebra]:=1: > M := RandomMatrix(3,outputoptions=[datatype=complex]): > Eigenvectors(M): Eigenvectors: "calling external function" Eigenvectors: "CLAPACK" hw_zgeevx_ > Digits := 20: > Eigenvectors(M): Eigenvectors: "calling external function" Eigenvectors: "CLAPACK" sw_zgeevx_ If you are using the older linalg package, and a matrix rather than a Matrix, then the computation will likely be done by Maple Library routines called by `evalf/Eigenvals`. I believe that these (interpreted, uncompiled) Maple Library routines are based on translations of the older ALGOL instances of EISPACK routines. See Numer Math. 13 293-304, B.N. Parlett and C.Reinsch, for details. I would recommend using LinearAlgebra instead of linalg here. Dave Linder Mathematical Software, Maplesoft

## no options on rtable_elems()...

I didn't read the help-page carefully enough. There are options available for rtable_num_elems(), but not for rtable_elems(). It could be useful to have some similar options, where they make sense, for rtable_elems() as well. I suppose that 'All', 'Stored', 'NonZero', and 'NonZeroStored' could be of use. It seems to me that the current behaviour of rtable_elems() is similar to that of rtable_num_elems() when called with the NonZeroStored option. > A := Matrix(2,2,shape=symmetric,[[0,3],[3,0]]): > rtable_elems(A); {(1, 2) = 3} > rtable_num_elems(A,NonZeroStored); 1 Dave Linder Mathematical Software, Maplesoft

## hmmmm...

> A := Array(1..3,1..4,(i,j)->i-j): > map([lhs], select(t -> type(rhs(t),nonnegint), rtable_elems(A,'All'))); {[2, 1], [3, 1], [3, 2]} Where are the indices for the entries equal to zero? Maybe I'm misinterpreting the 'All' option to rtable_elems(). Dave Linder Mathematical Software, Maplesoft

## support@maplesoft.com...

Yes, emailing support@maplesoft.com is the very best way to report a bug, issue, or make a feature request. It is best for precisely the reason that Jacques has stated, that it will get into the database and the queue. There are some people from Maplesoft who read almost every single post here of late. I strive to read it all. And when I see an issue that deserves more attention, I ask a `bugkeeper' to enter it, or I enter it myself. I've also seen items submitted by some readers here to Maplesoft's beta tester's forum, from which point the issues also get recorded. One potential way for an item to get overlooked is for no reader here, other than the poster, to view a posting as being an issue that needs addressing. So if a poster is in doubt about whether anyone would recognize an issue, then by all means send it to support@maplesoft.com and be as explicit as possible. Dave Linder Mathematical Software, Maplesoft

A clear example of using DEtools,matrixDE might make a useful Task template, for under "Differential Equations" in the Tasks section of the help-browser. Dave Linder Mathematical Software, Maplesoft

## almost apropos...

Perhaps the attendees of ISSAC 2007 might find suggestions about local fare of some interest. ps. That Guelph pub also serves a nice maple bison chili. Dave Linder Mathematical Software, Maplesoft

## The Lion Brewery...

One of the interesting social spots in town is The Huether Hotel. It is the home of the Lion Brewery Restaurant, which serves beer produced on the premises. Their English Ale is good. Their Adlys Ale is very good. Their lagers are popular too. There was a retail store attached to the building, which also sold the beer, but I don't know whether the outlet is still open. Dave Linder Mathematical Software, Maplesoft

## On Linux, that should be...

On Linux, that should be getenv(LD_LIBRARY_PATH) instead of getenv(PATH), to see the path that the runtime linker would use. It might have worked, regardless, only because there may have been no gmp shared library in a subdirectory of bin.X86_64_LINUX . Dave Linder Mathematical Software, Maplesoft

## binary search on "ordered" univariate pi...

I'm not sure whether it is coincidence, but just yesterday I was discussing fast binary search on an "ordered" univariate piecewise. The other developer was explaining his ideas on implementing this. So there is hope, for your 3). I can try to look at BSplineCurve, but as always my to-do list is long. I do agree with you, that interpolating and plotting large data sets is important. Dave Linder Mathematical Software, Maplesoft

## Maple-NAG Connector...

Hi Axel, The Application Center URL that you mentioned relates to the Maple-NAG Connector, which is a separate add-on product. Your question seems more general. Maple 11 has a functionality enhancement illustrated below, which relates to a particular use of the NAG Library: the univariate quadrature functions may now work with integrands which are not currently evalhf'able. Here's an example. The new functionality is accessible in Maple 11 by specifying the method.
```maple11> restart:
maple11> infolevel[`evalf/int`]:=1:
maple11> evalf(Int(x^2*BesselJ(0,x),x=1..10,method=_d01ajc));
Control:   Entering NAGInt
d01ajc:   result=.535387806202787586
0.5353878062
```
This functionality is possible only up to the limits of hardware double-precision.
```maple11> restart:
maple11> Digits:=trunc(evalhf(Digits)); # I'm using Linux
Digits := 15

maple11> infolevel[`evalf/int`]:=1:
maple11> evalf(Int(x^2*BesselJ(0,x),x=1..10,epsilon=1.0e-12,method=_d01ajc));
Control:   Entering NAGInt
d01ajc:   result=.535387806202795247
0.535387806202795
```
In Maple 10, specifying method=_d01ajc for a non-evalhf'able integrand would not work. Ie,
```maple10> restart:
maple10> infolevel[`evalf/int`]:=1:
maple10> evalf(Int(x^2*BesselJ(0,x),x=1..10,method=_d01ajc));
Control:   Entering NAGInt
Error, (in evalf/int) remember tables are not supported in evalhf
```
In both Maple 11 and Maple 10, when the method is not specified then the `evalf/Int` control mechanism will attempt first to use evalhf to evaluate the integrand, and if that fails may select another method according to its own scheme. For the given example, it uses the _CCquad method and produces the result.
```> evalf(Int(x^2*BesselJ(0,x),x=1..10));
Control:   Entering NAGInt
evalf/int/control:   NAG failed   result = result
evalf/int/CreateProc:   Trying easyproc
evalf/int/CreateProc:   Trying makeproc
evalf/int/ccquad:   n = 2 integral estimate = -36.98473043720
n = 6 integral estimate = .5957900650800
evalf/int/ccquad:   n = 18 integral estimate = .5353878062050
error = .5058362401405e-8
evalf/int/ccquad:   n = 54 integral estimate = .5353878062166
error = .5353878062166e-11
From ccquad, result = .5353878062166   integrand evals = 55
error = .5353878062166e-11
0.5353878062
```
Dave Linder Mathematical Software, Maplesoft

## approximate Jordan form...

How does linalg[jordan] work for the floating-point form of this Matrix (taken from the first reference below) whose eigenvalues are all identically zero. A := Matrix(5,5, [[-9,11,-21,63,-252], [70,-69,141,-421,1684], [-575,575,-1149,3451,-13801], [3891,-3891,7782,-23345,93365], [1024,-1024,2048,-6144,24572]]); LinearAlgebra[EigenConditionNumbers](Matrix(A,datatype=float),output=[values,conditionvalues]); linalg[jordan](map(evalf,convert(A,matrix))); Section 10.8 of http://www.mathworks.com/moler/eigs.pdf. http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/f1-7191.html http://www.math.rwth-aachen.de/mapleAnswers/html/340.html http://mathforum.org/kb/message.jspa?messageID=5507069&tstart=15 Section 72.8, Examples 5, of http://www.apmaths.uwo.ca/~djeffrey/Offprints/C5106_C072.pdf Best regards, Dave Linder Mathematical Software, Maplesoft

## Thanks Joe, yet again...

Thank you Joe, for pointing out the correction with respect to loading from an archive, which was in the submitter's original post. It's odd that I forgot that, since I've recently been using the method in ModuleLoad myself. :) Dave Linder Mathematical Software, Maplesoft

## Thanks Joe, yet again...

Thank you Joe, for pointing out the correction with respect to loading from an archive, which was in the submitter's original post. It's odd that I forgot that, since I've recently been using the method in ModuleLoad myself. :) Dave Linder Mathematical Software, Maplesoft

## Eigenvector schemes...

You mentioned Matrices with complex floats. I gather that you're interested in Matrices with floating-point entries (and all numeric). If you set infolevel[LinearAlgebra]:=1 then you should see displayed some mention of the external routine that gets used. Here's part of the breakdown in Maple 10, where by "hermitian" below I mean that the Matrix was created with shape=hermitian, float, symmetric : CLAPACK dspevd float, no shape : CLAPACK dgeevx complex, hermitian : CLAPACK zhpevd complex, no shape : CLAPACK zgeevx You may have noticed that zhpevd orders the (necessarily) real eigenvalues. Yet zgeevx does not. Here is some code which (I believe) will order the eigenvalues and eigenvectors by descending value of the argument of the eigenvalues. with(LinearAlgebra): A := RandomMatrix(50,outputoptions=[datatype=float]): B := RandomMatrix(50,outputoptions=[datatype=float]): C := A+B*I: evls,evcs := Eigenvectors(C): temp := Vector(Dimension(evls),(i)->[evls[i],i]): temp := sort(temp,(a,b)->argument(a)>=argument(b)): new_evls := Vector(Dimension(evls),(i)->temp[i],datatype=complex(float)): new_evcs := Matrix(Dimension(evls),(i,j)->evcs[i,temp[j]],datatype=complex(float)): Norm( C.new_evcs - new_evcs.DiagonalMatrix(new_evls) ); seq(argument(new_evls[i]),i=1..Dimension(evls)); You could modify the sort() call above, to order them in some other way. Note that it does not further order eigenvalues which might happen to have the same argument. The code could be simpler, were it to only order the eigenvalues. The purpose of the temp Vector is to store the ordering index information for use in ordering the eigenvectors. There are doubtless more efficient ways to go about it, but that might appear more complicated. Dave Linder Mathematical Software, Maplesoft
﻿