Должен ли размер массива WORK в процедурах xGEHRD и xHSEQR быть равен при вычислении собственных значений?

Question

Я вычисляю собственные значения плотной несимметричной матрицы А. Для этой цели я использую xGEHRD и xHSEQR Лакпатные подпрограммы, чтобы сначала вычислить верхнюю гессенбергскую форму А, а затем вычислить только собственные значения полученной матрицы.

Для обеих подпрограмм требуется параметр 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. Если бы значения были одинаковыми, я мог бы значительно упростить мой код (нет необходимости в realloc и только один вызов для определения значения LWORK, меньше проверки ошибок ...).

Мой вопрос: для той же матрицы и тех же значений ILO и IHI, могу ли я предположить, что значения будут равны обеим процедурам?

Answer 1

sgehrd.f Глядя, кажется, что оптимальный размер WORK для рутинного sgehrd() является N*NB, где NB является

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

где NBMAX=64. Следовательно, оптимальный LWORK зависит от N, ILO и IHI.

Глядя на shseqr.f, вычисление оптимальной длины LWORK является более сложным: процедура slaqr0() называется ... Но документация в файле состояний:

Если LWORK = -1, то SHSEQR делает запрос рабочей области. В этом случае SHSEQR проверяет входные параметры и оценки оптимального размера рабочего пространства для заданных значений N, МОТ и 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); 

} 

Compile по gcc main.c -o main -llapack -lblas

Мой выхода:

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