CLAPACK, ZHETRD приведение эрмитовой к трехдиагональной.

[Firth]

Доброго дня.

Вот уже который день бьюсь над этой функцией ZHETRD, которая на вход получает эрмитову матрицу, а на выходе выдает трехдиагональную.

Раздобыл документацию с примером в PDF формате:
http://www.nag.co.uk/numeric/FL/nagdoc_ ... F08FSF.pdf

Написал небольшой майн для теста:

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include "f2c.h"
#include "blaswrap.h"
#include "clapack.h"

int main()
{//
     doublecomplex H[4][4];          //Исходная матрица и та, что получается после разложения
     doublecomplex WORK1[256]; //Рабочая область
     doublecomplex TAU[3];          //Элементы унитарной Q
     doublereal D[4];                    // Выход диагональных элементов
     doublereal E[3];                    //Выход задиагональных элементов
     char UPLO;                         //Тип тип матрицы
     integer LDA;                       // Leading dimension
     int N;                               //Размер H
     int LWORK1;                     //Размер рабочей области
     int INFO;                           //Информация о выходе
     int i;
     int j;
     time_t start, end;

     UPLO='U'; //Матрица верхнетреугольная
     N=4;       // Размерности 4х4
     LWORK1=256; // рабочая область 16*16
     LDA=N;


    TAU[0].r=0; TAU[0].i=0;
    TAU[1].r=0; TAU[1].i=0;
    TAU[2].r=0; TAU[2].i=0;

     E[0]=0.0;
     E[1]=0.0;
     E[2]=0.0;

     D[0]=0.0;
     D[1]=0.0;
     D[2]=0.0;
     D[3]=0.0;

// Задаю матрицу

     H[0][0].r=-2.28; H[0][1].r=1.78; H[0][2].r=2.26; H[0][3].r=-0.12;
     H[1][0].r=1.78;  H[1][1].r=-1.12;H[1][2].r=0.01; H[1][3].r=-1.07;
     H[2][0].r=2.26;  H[2][1].r=0.01; H[2][2].r=-0.37;H[2][3].r=2.31;
     H[3][0].r=-0.12; H[3][1].r=-1.07;H[3][2].r=2.31; H[3][3].r=-0.73;


     H[0][0].i=0.0;   H[0][1].i=-2.03; H[0][2].i=0.1;  H[0][3].i=2.53;
     H[1][0].i=2.03;  H[1][1].i=0.0;   H[1][2].i=0.43; H[1][3].i=0.86;
     H[2][0].i=-0.1;  H[2][1].i=-0.43; H[2][2].i=0.0;  H[2][3].i=-0.92;
     H[3][0].i=-2.53; H[3][1].i=-0.86; H[3][2].i=0.92; H[3][3].i=0.0;

//Контрольный вывод
    for (i=0;i<N;i++){
        for(j=0;j<N;j++){
        printf("%8.2f%",H[i][j].r);
        printf("%5.2f",H[i][j].i);
        printf("(i)");
        }
        printf("\n");
    }
    printf("\n");

//вызов функции
start = time(NULL);
   zhetrd_(&UPLO,&N,&H,&LDA,&D,&E,&TAU,&WORK1,&LWORK1,&INFO);
end = time(NULL);

//вывод результата
    for (i=0;i<N;i++){
        for(j=0;j<N;j++){
        printf("%8.2f%",H[i][j].r);
        printf("%5.2f",H[i][j].i);
        printf("(i)");
        }
        printf("\n");
    }
    printf("\n");



printf("\nThe interval was %.2f seconds\n", difftime(end, start));

    return 0;
}

Такое ощущение, как будто, я не так подаю исходную матрицу или другие параметры. Очень смущает вот эта строчка в описании комментариях к ZHETRD:

Цитата:

On entry, the Hermitian matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced.

Остается впечатление, что можно передавать только треугольную матрицу, я пробовал по всякому, но результат существенно не меняется.

[Firth]

Эмм... я не знаю как вы это сделали, но все заработало. Спасибо большое.