#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;
  }