What is supposed to be a upper quasi-triangular matrix has non-zero
elements below the first sub-diagonal. This is demonstrated in the
attached file. Computer and OS information are also in the file.
#ifdef X
System: Host: debbie Kernel: 4.19.0-16-amd64 x86_64 bits: 64 compiler: gcc v: 8.3.0
Desktop: Cinnamon 4.8.6 wm: muffin dm: LightDM Distro: LMDE 4 Debbie
base: Debian 10.2 buster
Machine: Type: Desktop Mobo: ASUSTeK model: P5Q SE PLUS v: Rev 1.xx serial: <filter>
BIOS: American Megatrends v: 2202 date: 06/23/2009
#endif
// compile:
// g++ schur.cpp -I/usr/local/include -L/usr/local/lib -lgsl -lgslcblas -lm -o schur
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
// m : dimension of a square matrix
// a : on input, a linear array to be stored in row order,
// on return, the quasi-triangular form of the input
// v : on return, orthognal transformation matrix
void schur_decomposition(const int m, double* a, double* v)
{
gsl_eigen_nonsymm_workspace* w = gsl_eigen_nonsymm_alloc(m);
gsl_vector_complex* vec = gsl_vector_complex_alloc(m);
gsl_eigen_nonsymm_params(1, 0, w); // calc T without balancing
gsl_matrix_view A = gsl_matrix_view_array(a, m, m);
gsl_matrix_view V = gsl_matrix_view_array(v, m, m);
gsl_eigen_nonsymm_Z(&A.matrix, vec, &V.matrix, w);
double* p = a;
double* q = v;
for(int i = 0; i<m; i++)
for(int j = 0; j<m; j++)
{
*p++ = gsl_matrix_get(&A.matrix, i, j);
*q++ = gsl_matrix_get(&V.matrix, i, j);
}
gsl_eigen_nonsymm_free(w);
gsl_vector_complex_free(vec);
}
void print(FILE* f, double* x, int m, const char* title = 0)
{
if(!f) return;
if(title) fprintf(f, "%s\n", title);
double* p = x;
for(int i = 0; i<m; i++)
{
for(int j = 0; j< m; j++, p++)
fprintf(f, "%12.6g", (fabs(*p)>0.000001) ? *p : 0.0);
fprintf(f, "\n");
}
}
int main ()
{
const int m = 5;
double a[] = {
-0.2, 0.1, 0, 0.1, 0,
0, -0.2, 0.2, 0, 0,
0.3, 0, -0.3, 0, 0,
0, 0, 0, -0.2, 0.2,
0.3, 0, 0, 0, -0.3};
double* v = new double[m*m];
print(stdout, a, m, "\nOriginal matrix");
schur_decomposition(m, a, v);
print(stdout, a, m, "\nSchur form");
print(stdout, v, m, "\nTransformation matrix");
return 0;
}