<![CDATA[
#include <numeric.h>
int main() {
int m = 3;
array double a[3][3] = { 0.8,0.2,0.1,
0.1,0.7,0.3,
0.1,0.1,0.6};
int status;
array double p[m][m], h[m][m];
status = hessdecomp(a, h);
if (status == 0) {
printf("h =\n%f\n", h);
}
else
printf(" Hessenberg matrix calculation error.\n");
status = hessdecomp(a, h, p);
if (status == 0) {
printf("h =\n%f\n", h);
printf("p =\n%f\n", p);
printf("transpose(p)*p =\n%f\n",transpose(p)*p );
printf("transpose(p)*a*p - h =\n%f\n", transpose(p)*a*p - h);
}
else
printf(" Hessenberg matrix calculation error.\n");
}
]]>
<![CDATA[
h =
0.800000 -0.212132 -0.070711
-0.141421 0.850000 -0.050000
0.000000 0.150000 0.450000
h =
0.800000 -0.212132 -0.070711
-0.141421 0.850000 -0.050000
0.000000 0.150000 0.450000
p =
1.000000 0.000000 0.000000
0.000000 -0.707107 -0.707107
0.000000 -0.707107 0.707107
transpose(p)*p =
1.000000 0.000000 0.000000
0.000000 1.000000 -0.000000
0.000000 -0.000000 1.000000
transpose(p)*a*p - h =
0.000000 -0.000000 -0.000000
0.000000 0.000000 -0.000000
0.000000 0.000000 0.000000
]]>