<![CDATA[
#include <numeric.h>
int main() {
int m = 3;
array double complex a[3][3] = { complex(-149,1), -50, -154,
537, 180, 546,
-27, -9, -25};
int status;
array double complex h[m][m];
array double complex p[m][m];
status = hessdecomp(a, h);
if (status == 0) {
printf("h =\n%5.2f\n", h);
}
else
printf("Hessenberg matrix calculation error.\n");
status = hessdecomp(a, h, p);
if (status == 0) {
printf("h =\n%5.2f\n", h);
printf("p =\n%5.2f\n", p);
printf("transpose(p)*p =\n%5.2f\n", transpose(p)*p);
printf("conj(transpose(p))*a*p - h =\n%5.2f\n", conj(transpose(p)*a*p - h));
}
else
printf("Hessenberg matrix calculation error.\n");
}
]]>
<![CDATA[
h =
complex(-149.00, 1.00) complex(42.20, 0.00) complex(-156.32, 0.00)
complex(-537.68, 0.00) complex(152.55, 0.00) complex(-554.93, 0.00)
complex( 0.00, 0.00) complex( 0.07, 0.00) complex( 2.45, 0.00)
h =
complex(-149.00, 1.00) complex(42.20, 0.00) complex(-156.32, 0.00)
complex(-537.68, 0.00) complex(152.55, 0.00) complex(-554.93, 0.00)
complex( 0.00, 0.00) complex( 0.07, 0.00) complex( 2.45, 0.00)
p =
complex( 1.00, 0.00) complex( 0.00, 0.00) complex( 0.00, 0.00)
complex( 0.00, 0.00) complex(-1.00, 0.00) complex( 0.05, 0.00)
complex( 0.00, 0.00) complex( 0.05, 0.00) complex( 1.00, 0.00)
transpose(p)*p =
complex( 1.00, 0.00) complex( 0.00, 0.00) complex( 0.00, 0.00)
complex( 0.00, 0.00) complex( 1.00, 0.00) complex( 0.00, 0.00)
complex( 0.00, 0.00) complex( 0.00, 0.00) complex( 1.00, 0.00)
conj(transpose(p))*a*p - h =
complex( 0.00,-0.00) complex( 0.00,-0.00) complex( 0.00,-0.00)
complex( 0.00,-0.00) complex( 0.00,-0.00) complex( 0.00,-0.00)
complex( 0.00,-0.00) complex( 0.00,-0.00) complex(-0.00,-0.00)
]]>