A bug in the schur decomposition function?

[david allen via Gsl-discuss]

I am using g++ on linux mint and getting a wrong result for the schur 
decomposition function. I have attached a C++ program that produces 
wrong results and octave program that produces correct results.
// To compile
// g++ -Wall schur-test.cpp -lgsl -lgslcblas -lm -o schur-test
#include <stdio.h>
#include <string.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int gslprint(FILE* f, gsl_matrix* A, const char* title)
  {
  if(!f) return 1;
  if(title) fprintf(f, "%s\n", title);
  int m = (int) A->size1;
  int n = (int) A->size2;
  int k = 0;
  for (int i = 0; i < m; i++)
    {
    for (int j = 0; j < n; j++)
      fprintf(f, "%12.4lf", gsl_matrix_get(A, i, j));
    fprintf(f, "\n");
    }  
  return 0;
  }

int main ()
  {
  int m = 3;
  double data[] = {
                  -0.5000,   0.1000,   0.4000,
                   0.3000,  -0.3000,        0,
                        0,        0,        0};

  gsl_eigen_nonsymm_workspace*    w = gsl_eigen_nonsymm_alloc (m);
  gsl_matrix*                     M = gsl_matrix_alloc (m, m);
  gsl_matrix*                     Z = gsl_matrix_alloc (m, m);
  gsl_vector_complex*           vec = gsl_vector_complex_alloc (m);

  memcpy(M->data, data, m*m*sizeof(double));
  gslprint(stdout, M, "Original matrix");
  gsl_eigen_nonsymm_Z (M, vec, Z, w);
  gslprint(stdout, M, "Schur matrix");
  gslprint(stdout, Z, "Transformation matrix");
  
  gsl_eigen_nonsymm_free(w);
  gsl_vector_complex_free(vec);
  gsl_matrix_free (M);
  gsl_matrix_free (Z);
  return 0;
  }
  
A = [-0.5 0.1 0.4;
      0.3 -0.3 0;
      0    0   0];
disp("A matrix");   
disp(A);   
[T, S] = schur(A);
disp("Schur matrix");
disp(S);
disp("Transformation matrix");
disp(T);