LAPACK

Michael Monagan, monagan@inf.ethz.ch

LAPACK is a library of numerical software packages for linear algebra. On the cover of the LAPACK User's guide you will find this matrix

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> lapack := matrix(6,6, [ L, A, P, A, C, K,
L,-A, P,-A, C,-K,
L, A, P, A,-C,-K,
L,-A, P,-A,-C, K,
L, A,-P,-A, C, K,
L,-A,-P, A, C,-K ] );

Now this is clearly not a random matrix. But more than that, there is bound to be something interesting about this matrix. But what? Unfortunately, since LAPACK is a numerical only software library, you won't be able to find out what that is using LAPACK!! Clearly a symbolic system is helpful here. Let's try computing the determinant of this matrix

> det(lapack);

Bingo! So that's the secret. The determinant = -128 LAPACK. Is there anything else interesting about this matrix? Lets look at the inverse

> evalm(1/lapack);

> evalm(4*%);

It has a nice block structure similar to the orginal matrix. Let's look at the characteristic polynomial

> cp := charpoly(lapack,lambda);

> factor(cp);

> collect(%,lambda);

Probably there is something interesting here. But I was not able to express the eigenvalues in closed form. I looked instead at the product

> LAT := evalm(lapack &* transpose(lapack));

Bingo, the diagonal elements = LAPACK^2. Definitely something interesting there!

> cp := charpoly(LAT,lambda);

> cp := factor(cp);

So the singular values have a closed form. Given by the square roots of

> svals := solve(cp,lambda);

Let's play with the singular values a little. The first one is

> sqrt(svals[1],symbolic);

> nullspace(LAT-8*A^2);

The second one is

> svals[2];

> nullspace(LAT-svals[2]);

What is interesting here is that these are independent of L,P,C. Working with the cube roots is too messy. Better to work in an implicit form with RootOfs. So lets pick up the cubic factor again

> select( proc(f) degree(f,lambda)=3 end, cp );

> alias(alpha = RootOf(%,lambda));

> nullspace(LAT-alpha);

Nothing interesting here? But Maple could crank it out - it took quite a lot of time and memory though.