Precision of Lapack function DGEES

[fatalme]

Recently I met a problem both in Lapack and Matlab.
When I use DGEES() in Lapack or schur() in Matlab to schur_factorize a
matrix A, i.e. to make it in this form A=U*T*U'.
When I got U and T, I try to verify whether it satisfy
((Temp=A-U*T*U')==0).
The result is very disappointing.
some elements in matrix temp are several magnitude larger than elements in matrix A at corresponding positions.
What's the reason of this problem?
How to solve it?

A is 15x15 matrix ranging very wide.

0.2316733D-11 0.0000000D+00 0.1060792D-16 -0.1060792D-16 -0.1364627D-21 0.1060805D-16 0.0000000D+00 0.0000000D+00 -0.1060792D-16 0.0000000D+00 0.1364627D-21 -0.1364627D-21 -0.1878819D-42 0.0000000D+00 0.0000000D+00
0.4163336D-16 0.0000000D+00 0.1588187D-21 -0.1588187D-21 -0.1615587D-26 0.1588187D-21 0.0000000D+00 0.0000000D+00 -0.1588187D-21 0.0000000D+00 0.1615587D-26 -0.1615587D-26 -0.6842278D-48 0.0000000D+00 0.0000000D+00
0.3641793D-07 0.0000000D+00 0.1690943D-12 -0.1690943D-12 0.2605464D-19 0.1690942D-12 0.0000000D+00 0.0000000D+00 -0.1690943D-12 0.0000000D+00 -0.2605464D-19 0.2605464D-19 -0.1248483D-38 0.0000000D+00 0.0000000D+00
-0.4397242D+12 0.0000000D+00 0.1892094D-10 0.2208738D+07 -0.6233284D+00 0.2190697D+07 -0.6233284D+00 0.0000000D+00 -0.2190696D+07 0.0000000D+00 0.1804122D+05 0.2605464D-19 -0.1804122D+05 0.0000000D+00 0.0000000D+00
-0.2410358D+13 0.0000000D+00 0.7539664D+02 -0.7539664D+02 0.4717222D+07 0.3046213D-11 -0.7539664D+02 0.0000000D+00 0.2575040D-14 0.0000000D+00 0.4717146D+07 -0.4717146D+07 0.1901244D-40 0.0000000D+00 0.0000000D+00
0.2423797D-08 0.0000000D+00 0.1372268D-13 -0.1372268D-13 -0.1388481D-18 0.1372282D-13 0.0000000D+00 0.0000000D+00 -0.1372268D-13 0.0000000D+00 0.1388481D-18 -0.1388481D-18 0.1901243D-40 0.0000000D+00 0.0000000D+00
-0.8204680D+11 0.0000000D+00 0.1045757D+05 -0.8283219D+06 0.2992852D-09 -0.2575064D-14 0.8283219D+06 -0.1045757D+05 0.2575038D-14 0.8410451D-09 0.8178643D+06 -0.8178643D+06 -0.8576857D-39 0.0000000D+00 0.0000000D+00
-0.1220541D-08 0.0000000D+00 -0.2575040D-14 0.2575040D-14 0.2605464D-19 -0.2575066D-14 -0.2056157D-29 0.2056157D-29 0.2575040D-14 0.0000000D+00 -0.2605464D-19 0.2605464D-19 0.1901242D-40 0.0000000D+00 0.0000000D+00
-0.4189072D+02 0.0000000D+00 0.5157582D-06 -0.5157583D-06 0.2605464D-19 0.2846071D-09 -0.1722371D-03 0.0000000D+00 0.1727523D-03 -0.5154675D-06 0.1722371D-03 -0.1722371D-03 -0.6160259D-11 0.0000000D+00 0.0000000D+00
0.5884200D-07 0.0000000D+00 -0.2575040D-14 0.2575040D-14 0.2605464D-19 -0.2575066D-14 -0.5843693D-12 0.0000000D+00 0.2575040D-14 0.2921846D-12 -0.2605464D-19 0.2605464D-19 0.1901244D-40 0.0000000D+00 0.0000000D+00
-0.1125002D-08 0.0000000D+00 -0.2575039D-14 0.2575039D-14 -0.6079415D-19 -0.2574979D-14 0.0000000D+00 0.0000000D+00 0.2575039D-14 0.0000000D+00 0.6079415D-19 -0.6079415D-19 -0.4436235D-40 0.0000000D+00 0.0000000D+00
-0.1351098D-08 0.0000000D+00 -0.2575039D-14 0.2575039D-14 0.2605464D-19 -0.2575066D-14 0.0000000D+00 0.0000000D+00 0.2575039D-14 0.0000000D+00 -0.2605464D-19 0.2605464D-19 0.1901243D-40 0.0000000D+00 0.0000000D+00
-0.1986585D+13 0.0000000D+00 0.9579757D-04 -0.1097734D+02 0.2605464D-19 0.1444103D+08 0.0000000D+00 0.0000000D+00 -0.1444103D+08 0.0000000D+00 -0.1444105D+08 0.2605464D-19 0.1444105D+08 -0.9579757D-04 0.0000000D+00
-0.1321972D-08 0.0000000D+00 -0.2575039D-14 0.2575039D-14 0.2605464D-19 -0.2575066D-14 0.0000000D+00 0.0000000D+00 0.2575039D-14 0.0000000D+00 -0.2605464D-19 0.2605464D-19 -0.4983222D-29 0.4983222D-29 0.0000000D+00
-0.1122824D-08 0.0000000D+00 -0.2575040D-14 0.2575040D-14 0.2605464D-19 -0.2575066D-14 0.0000000D+00 0.0000000D+00 0.2575040D-14 0.0000000D+00 -0.2605464D-19 0.2605464D-19 0.1901244D-40 0.0000000D+00 0.0000000D+00

[sven]

Just to check, did you remember to zero the lower triangular part of T? (When doing that remember that T may contain 2 by 2 diagonal blocks.)

Best wishes.

Sven Hammarling

[fatalme]

DGEES generates a upper triangle matrix T by default, and U*U'=I.
This triangle matrix T also include some 2 by 2 diagonal blocks.
By the way, I use DGEES from Campq Visual Fortran library CXML.

Many thanks for your reply.

Just to check, did you remember to zero the lower triangular part of T? (When doing that remember that T may contain 2 by 2 diagonal blocks.)

Best wishes.

Sven Hammarling

[sven]

The Schur factorization is normwise backward stable, not componentwise wise, so you cannot expect the small elements of A and UTU' to agree. But norm(A - UTU') should be small relative to norm (A).

Sven.

[bassabo]

hello,

I have the same problem.

I take a very simple example A =[1,3;2,4].

The results are : U=[0.447214,-0.894427;0.894427,0.447214] and T= [1,3;0,4]

U*U' =Id and T is upper triangle matrix but A -U*T*U' =[ - 1.1999982, 3.5999997 ; 5.5999998, 1.1999981].

Do anyone have an explanation ??

[sven]

I'm afraid that you must have a bug in your calling program. I get

U = ( -0.9094 -0.4160 ), T = ( -0.3723 1.0000 )
( 0.4160 -0.9094 ) ( 0 5.3723 )

and A - U*T*U' = 1.0E-15*( 0 0 ).
( 0.2220 0 )

With best wishes,

Sven.