고유 값 계산시 xGEHRD 및 xHSEQR 루틴의 WORK 배열 크기가 동일해야합니까?

Question

밀도가 높은 비대칭 행렬 A의 고유 값을 계산하고 있습니다. 그 목적을 위해 먼저 xGEHRD 및 xHSEQR Lapack 루틴을 사용하여 먼저 헤센 버그 형식의 A를 계산 한 다음 얻은 행렬의 고유 값을 계산합니다.

두 루틴 모두 LWORK 매개 변수가 필요하며 두 매개 변수 모두 최적 값을 계산하는 메커니즘을 제공합니다. 이 매개 변수는 버퍼 기술의 내부 차단과 관련이 있다고 생각하지만 어떻게 결정되는지는 알 수 없습니다. 내가 WORK 배열의 차원에 대해 항상 같은 최적의 값을 취득, 일부 테스트를 수행 한

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 

: 같은 워크 플로가 있어야한다 LWORK의 최적 값을 얻기 위해 쿼리 메커니즘을 사용

. 값이 같으면 코드를 훨씬 단순하게 만들 수 있습니다 (realloc이 필요없고 LWORK의 값을 결정할 때 하나의 호출만으로 오류 검사가 줄어 듭니다 ...).

내 질문은 동일한 행렬과 ILO와 IHI의 동일한 값에 대해 입니다. 그 값이 두 루틴과 동일하다고 가정 할 수 있습니까? NB이

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

NBMAX=64입니다

Answer 1

sgehrd.f 보면, 일상적인 sgehrd()에 대한 WORK의 최적 크기가 N*NB 것 같다. 따라서 최적의 LWORK은 N, ILO 및 IHI에 따라 다릅니다.

shseqr.f 보면은, 최적의 길이 LWORK의 계산은 더 복잡하다 : 루틴 slaqr0()이라고합니다 ...하지만 파일 상태에서 문서 : LWORK = -1, 다음 SHSEQR가하는

경우 작업 영역 쿼리 이 경우 SHSEQR은 추정치 에 대해 주어진 값 N, ILO 및 IHI에 대한 최적 작업 영역 크기를 확인합니다. 예상치는 작업 (1)에서 으로 반환됩니다. LWORK 관련 오류 메시지는 XERBLA가 발행 한 입니다. H도 Z도 액세스하지 않습니다.

WORK의 최적 길이는 sgehrd() 및 shseqr() 다를 수있다.

#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