Выводы Schur из LAPACK не совпадают с MATLAB

Question

Я пытаюсь воспроизвести функцию schur, используя LAPACK, но сталкиваюсь с трудностями при получении результатов, соответствующих MATLAB, несмотря на то, что MATLAB внутренне использует LAPACK для своего schur. Я попытался использовать LAPACK_dgehrd и LAPACK_dhseqr, однако последний столбец T не соответствует MATLAB. Затем я попытался использовать dgees, однако порядок собственных значений в матрице T отличается от порядка, заданного MATLAB, что, в свою очередь, приводит к тому, что остальная часть матрицы T отличается. Я вставил код C, который использовал для LAPACK. Часть dgehrd и dhesqr прокомментирована. Также приложено изображение с результатами для небольшой матрицы 4х4 с использованием трех способов:

1. MATLAB
2. dgehrd&dhesqr
3. dgees

Функция MATLAB, которую я использовал, была очень простой: [Q,T]=schur(A);

Обратите внимание, что когда я использую dgehrd и dhesqr, выходные данные даже не удовлетворяют условию UT U '= A schur; в то время как выходы dgees и MATLAB удовлетворяют этому.

Не могли бы вы подсказать мне, что я делаю неправильно с использованием dgehrd&dhesqr, из-за которого только последний столбец T был неверным? Я пытался использовать dgebal до LAPACK_dgehrd и LAPACK_dorghr до LAPACK_dhseqr, но это тоже не помогло. Кроме того, я попытался использовать опцию sort "S" с dgees, используя функцию lapack_logical в конце кода C, установив sdim в nrows-1 и; но это стало причиной аварии, не могли бы вы подсказать мне, где я ошибаюсь?

#include "lapacke.h"
#include <stdio.h>
#include <stdlib.h>
#ifndef free_NULL
    #define free_NULL(T) {free(T); T = NULL;}
#endif

int schur(double *val, int nrows, double *U, double *T)
{
    int i, j, lwork = -1;
    int info;
    int sdim = 0;//nrows-1;//
    int ilo = 1, ihi = nrows -1;
    double *tau = NULL;
    double wkopt;
    double *wr = NULL, *wi = NULL, *Z = NULL;
    double  *work = NULL;
    extern LAPACK_D_SELECT2          SELECT;
    lapack_logical    *bwork;
    FILE *file = NULL;  
    file = fopen("Schurpack.csv","w");
    wr = (double *)calloc(nrows,sizeof(double));
    wi = (double *)calloc(nrows,sizeof(double));
    Z  = (double *)calloc(nrows*nrows,sizeof(double));
    bwork  = (lapack_logical *)calloc(nrows,sizeof(lapack_logical));

    LAPACK_dgees ( "V", "N", SELECT, &nrows, val, &nrows, &sdim, wr, wi, Z, &nrows, &wkopt, &lwork, bwork, &info ); 
    lwork = (int)wkopt;
    work = (double*)malloc( lwork*sizeof(double) );
    LAPACK_dgees ( "V", "N", SELECT, &nrows, val, &nrows, &sdim, wr, wi, Z, &nrows, work, &lwork, bwork, &info ); 


    //tau = (double *)calloc(nrows-1,sizeof(double));
    //lwork = -1; 
    //LAPACK_dgehrd( &nrows, &ilo, &ihi, val, &nrows, tau, &wkopt, &lwork, &info);
    //
    //lwork = (int)wkopt;
    //work = (double *)calloc(lwork,sizeof(double));
    //LAPACK_dgehrd( &nrows, &ilo, &ihi, val, &nrows, tau, work, &lwork, &info);

    //lwork = nrows;
    //work = (double *)calloc(nrows,sizeof(double));
    //LAPACK_dhseqr ( "S", "I", &nrows, &ilo, &ihi, val, &nrows, wr, wi, Z, &nrows, work, &lwork, &info );   
    if( info > 0 ) 
    {
        printf( "The algorithm computing Schur Decompsition failed to converge.\n" );
        exit( 1 );
    }

    for( i = 0; i < nrows; i++)
    {
        for ( j = 0; j < nrows; j++)
        { 
            U[i*nrows+j] = Z[j*nrows+i];
            T[i*nrows+j] = val[j*nrows+i];
        }
    }

    fprintf(file, "U\n\n");
    for (j = 0; j < nrows; j++)
    {
        for (i = 0; i < nrows; i++)
            fprintf(file,"%f,",U[j*(nrows)+i]);
        fprintf(file,"\n");
    }


    fprintf(file,"T\n\n");
    for (j = 0; j < nrows; j++)
    {
        for (i = 0; i < nrows; i++)
            fprintf(file,"%f,",T[j*(nrows)+i]);
        fprintf(file,"\n");
    }
    fclose(file);
    file = NULL;


    if (tau)
        free_NULL(tau);
    if (wr)
        free_NULL(wr);
    if (wi)
        free_NULL(wi);
    if (Z)
        free_NULL(Z);
    if (work)
        free_NULL(work);

    return(0);
}


lapack_logical SELECT ( double ER, double EI)
{
    if ((ER==1) && (EI==0))
        return(0);
    else
        return(1);
}