2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 CLAPACK, ZHETRD приведение эрмитовой к трехдиагональной.
Сообщение28.02.2017, 12:15 


06/06/11
60
Доброго дня.

Вот уже который день бьюсь над этой функцией 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.


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

 Профиль  
                  
 
 Re: CLAPACK, ZHETRD приведение эрмитовой к трехдиагональной.
Сообщение28.02.2017, 15:27 


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 2 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group