我在计算稠密非对称矩阵a的特征值。为此,我使用xGEHRDxHSEQRLapack例程,首先计算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

我做了一些测试,得到了工作数组维数的最优值。如果值相同,我可以简化代码很多(不需要realloc,只需要一个调用来确定LWORK的值,更少的错误检查…)。
我的问题是,对于ILO和IHI的相同矩阵和相同值,我能假设两个例程的值将相等吗?

最佳答案

sgehrd.f来看,常规WORK的最佳sgehrd()大小是N*NB,其中NB

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

其中NBMAX=64。因此,最佳LWORK取决于NILOIHI
shseqr.f来看,计算最佳长度LWORK更为复杂:例程slaqr0()被称为。。。但文件中的文件说:
如果LWORK=-1,则SHSEQR执行工作区查询。
在这种情况下,SHSEQR检查输入参数和
估计给定的最佳工作空间大小
N值、ILO值和IHI值。估计数已退回
工作中(1)。没有与LWORK相关的错误消息是
由XERBLA发行。H和Z都不可访问。
WORKsgehrd()的最佳长度可能不同。下面是一个例子:
#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);

}

编译人shseqr()
我的输出是:
sgehrd,length of optimal work : 320
shseqr,length of optimal work : 10

10-07 23:51