0

我正在尝试使用 LAPACK 重现 schur 函数，但在获得与 MATLAB 匹配的结果方面面临困难，尽管 MATLAB 在内部使用 LAPACK 作为其 schur。我尝试使用 LAPACK_dgehrd 和 LAPACK_dhseqr，但是 T 的最后一列与 MATLAB 不匹配。然后我尝试使用 dgees，但是矩阵 T 中特征值的顺序与 MATLAB 给出的顺序不同，这反过来又导致 T 矩阵的其余部分不同。我已经插入了用于 LAPACK 的 C 代码。dgehrd 和 dhesqr 部分已被注释。还附上了一张图片，使用三种方式获得了一个小的 4x4 矩阵的结果：

  1. MATLAB
  2. dgehrd&dhesqr
  3. 吉斯

我使用的 MATLAB 函数非常简单：[Q,T]=schur(A);

请注意，当我使用 dgehrd 和 dhesqr 时，输出甚至不满足 schur 的 U T U'=A 条件；而 dgees 和 MATLAB 输出满足这一点。

您能否指导我使用 dgehrd&dhesqr 做错了什么，导致 T 的最后一列不正确？我曾尝试在 LAPACK_dgehrd 之前使用 dgebal，在 LAPACK_dhseqr 之前使用 LAPACK_dorghr，但这也无济于事。此外，我尝试在 C 代码末尾使用 lapack_logical 函数对 dgees 使用排序“S”选项，将 sdim 设置为 nrows-1 和；但这导致了崩溃，你能指导我哪里出错了吗？

来自 MATLAB 和两个 LAPACK 包的 schur 比较

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

0 回答 0