c - xGEHRD 和 xHSEQR 例程中 WORK 数组的维数在特征值计算中是否应该相等？

Question

我正在计算稠密非对称矩阵 A 的特征值。为此，我使用 xGEHRD和 xHSEQR Lapack 例程首先计算 A 的上 Hessenberg 形式，然后仅计算所获得矩阵的特征值。

这两个例程都需要参数 LWORK，并且都提供了一种计算其最佳值的机制。我相信这个参数与缓冲技术的内部阻塞有关，但我不知道它是如何确定的。

使用查询机制获取LWORK的最优值的工作流应该是这样的:

int LWORK = -1;
float* OPT_LWORK = (float*) malloc( sizeof(float));

sgehrd_ (..., OPT_LWORK ,&LWORK, ...) // query optimal value for sgehrd

LWORK = (int) OPT_WORK
float* WORK = (float*) malloc( (int) sizeof(float) * LWORK );

sgehrd_ (..., WORK ,&LWORK, ...) // calculate Hessenberg

int LWORK = -1;
shseqr_ (..., OPT_LWORK ,&LWORK, ...) // query optimal value for shseqr

LWORK = (int) OPT_WORK
float* WORK = // possibly realloc with the new LWORK value

shseqr_ (..., WORK ,&LWORK, ...) // finally obtain eigenvalues

我做了一些测试，始终获得 WORK 数组维度的相同最佳值。如果值相同，我可以大大简化我的代码(不需要重新分配，只需一次调用即可确定 LWORK 的值，减少错误检查...)。

我的问题是，对于相同的矩阵和相同的 ILO 和 IHI 值，我可以假设两个例程的值相等吗？

Answer 1

查看sgehrd.f , 例程 sgehrd() 的 WORK 的最佳大小似乎是 N*NB，其中 NB 是

NB = MIN( NBMAX, ILAENV( 1, 'SGEHRD', ' ', N, ILO, IHI, -1 ) )

其中 NBMAX=64。因此，最优的LWORK 取决于N、ILO 和IHI。

查看shseqr.f ，最佳长度的计算 LWORK 更复杂:例程 slaqr0() 被调用...但是文件中的文档指出:

If LWORK = -1, then SHSEQR does a workspace query. In this case, SHSEQR checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.

对于 sgehrd() 和 shseqr()，WORK 的最佳长度可能不同。这是一个例子:

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

extern void sgehrd_(int *n,int* ilo, int* ihi, float* a,int* lda,float* tau,float* work, int* lwork,int* info);

extern void shseqr_(char* job,char* compz,int *n,int* ilo, int* ihi,float* h,int* ldh,float* wr,float* wi,float* z,int* ldz,float* work, int* lwork,int* info);

int main()
{
    int n=10;
    int ilo=n;
    int ihi=n;
    float* a=malloc(sizeof(float)*n*n);
    int lda=n;
    float* tau=malloc(sizeof(float)*(n-1));
    int info;

    char job='S';
    char compz='I';
    float* h=malloc(sizeof(float)*n*n);
    int ldh=n;
    float* wr=malloc(sizeof(float)*(n));
    float* wi=malloc(sizeof(float)*(n));
    float* z=malloc(sizeof(float)*n*n);
    int ldz=n;

    int LWORK = -1;
    float* OPT_LWORK =(float*) malloc( sizeof(float));

    sgehrd_ (&n,&ilo, &ihi, a,&lda,tau,OPT_LWORK ,&LWORK,&info); // query optimal value for sgehrd
    if(info!=0){printf("sgehrd lwork=-1 : failed\n");}

    LWORK = (int) OPT_LWORK[0];
    printf("sgehrd,length of optimal work : %d \n",LWORK);

    float* WORK = (float*) malloc( (int) sizeof(float) * LWORK );

    sgehrd_ (&n,&ilo, &ihi, a,&lda,tau,WORK ,&LWORK,&info);// calculate Hessenberg
    if(info!=0){printf("sgehrd execution : failed\n");}

    LWORK = -1;
    shseqr_ (&job,&compz,&n,&ilo, &ihi, h,&ldh, wr, wi,z,&ldz, OPT_LWORK ,&LWORK, &info); // query optimal value for shseqr
    if(info!=0){printf("shgeqr lwork=-1 : failed\n");}

    LWORK = (int) OPT_LWORK[0];
    printf("shseqr,length of optimal work : %d \n",LWORK);

    WORK = realloc(WORK,(int) sizeof(float) * LWORK );// possibly realloc with the new LWORK value

    shseqr_ (&job,&compz,&n,&ilo, &ihi, h,&ldh, wr, wi,z,&ldz, WORK ,&LWORK, &info);  // finally obtain eigenvalues
    if(info!=0){printf("shgeqr execution : failed\n");}

    free(OPT_LWORK);
    free(WORK);

    free(a);
    free(tau);
    free(h);
    free(wr);
    free(wi);
    free(z);

}

通过gcc main.c -o main -llapack -lblas编译

我的输出是:

sgehrd,length of optimal work : 320 
shseqr,length of optimal work : 10