VPF::Метод Хаусхолдера+QL-алгоритм - Форум программистов

allsolovey

Дата 11.2.2009, 22:26 (ссылка)

(нет голосов)

Загрузка ...

Empty

Профиль
Группа: Участник
Сообщений: 29
Регистрация: 3.1.2008

Репутация: нет
Всего: нет

Решил реализовать Метод Хаусхолдера+QL-алгоритм (на С++ в билдере), но есть проблемы. Собственные значения считает правильно а вот векторы нет. Два метода и две процедуры соответственно - это tred2 и tqli. Может я не так понял как заполняются поддиагональные элементы из вектора на входе tqli в z[][] (то ли по стб то ли по диагоналям то ли по стр).Помогите плиз разобраться smile

. Вот исходник :

Код


#include <math.h>
#include <iostream>
#include <fstream>
#include "nrutil.h"
using namespace std;


/* максимальное число итераций */
#define MAXITER 30

/* Здесь определяются некоторые утилиты типа выделения памяти */


/* Редукция Хаусхолдера действительной симметричной матрицы a[1...n][1...n].
На выходе a заменяется ортогональной матрицей трансформации q.
d[1...n] возвращает диагональ трехдиагональной матрицы.
e[1...n] возвращает внедиагональные элементы, причем e[1]=0. */

void tred2(float **a, int n, float *d, float *e) {
    float scale,hh,h,g,f;
    /* Проход по стадиям процесса редукции */
    for(int i = n; i >= 2; i--)
    {
        int l = i - 1;
        float h = 0.;
        float scale = 0;
        if(l>1)
        {
            /* вычислить шкалу */
            for(int k = 1; k <= l ; ++k)
                scale += fabs(a[i][k]);
            /* малая величина шкалы -> пропустить преобразование */
            if(scale==0.)
                e[i] = a[i][l];
            else
            {
                /* отмасштабировать строку и вычислить s2 в h */
                for(int k = 1; k<=l; ++k)
                {
                    a[i][k] / =scale;
                    h += a[i][k] * a[i][k];
                }
                /* вычислить вектор u */
                f=a[i][l];
                g=(f >= 0. ? -sqrt(h) : sqrt(h));
                e[i] = scale*g;
                h -= f*g;
                /* записать u на место i-го ряда a */
                a[i][l]=f-g;
                /* вычисление u/h, Au, p, K */
                f=0.;
                for(j=1;j<=l;j++)
                {
                    a[j][i]=a[i][j]/h;
                    /* сформировать элемент Au (в g) */
                    g=0.;
                    for(k=1;k<=j;k++)
                        g += a[j][k]*a[i][k];
                    for(k=j+1;k<=l;k++)
                        g += a[k][j]*a[i][k];
                    /* загрузить элемент p во временно неиспользуемую область e */
                    e[j]=g/h;
                    /* подготовка к формированию K */
                    f += e[j]*a[i][j];
                }
                /* Сформировать K */
                hh=f/(h+h);
                for(j=1;j<=l;j++)
                {
                    /* Сформировать q и поместить на место p (в e) */
                    f = a[i][j];
                    e[j] = g =e [j]-hh*f;
                    /* Трансформировать матрицу a */
                    for(k=1;k<=j;k++)
                        a[j][k] -= (f*e[k]+g*a[i][k]);
                }
            }
        }
        else
        {
            e[i]=a[i][l];
        }
        d[i]=h;
    }

    d[1]=0.;
    e[1]=0.;
    for(i=1;i<=n;i++)
    {
        l=i-1;
        /* этот блок будет пропущен при i=1 */
        if(d[i]!=0.)
        {
            for(j=1;j<=l;j++)
            {
                g=0.;
                /* формируем PQ, используя u и u/H */
                for(k=1;k<=l;k++)
                    g += a[i][k]*a[k][j];

                for(k=1;k<=l;k++)
                    a[k][j] -= g*a[k][i];
            }
        }
        d[i]=a[i][i];
        /* ряд и колонка матрицы a преобразуются к единичной, для след. итерации */
        a[i][i]=0.;
        for(j=1;j<=l;j++)
            a[j][i]=a[i][j]=0.;
    }
}



/* QL-алгоритм с неявными сдвигами для определения собственных значений (и собственных
векторов) действительной, симметричной, трехдиагональной матрицы. Эта матрица может
быть предварительно получена с помощью программы tred2. На входе d[1...n] содержит
диагональ исходной матрицы, на выходе - собственные значения. На входе e[1...n]
содержит поддиагональные элементы, начиная с e[2]. На выходе массив e разрушается.
При необходимости поиска только собственных значений в программе следует
закомментировать или удалить инструкции, необходимые только для поиска собственных
векторов. Если требуются собственные вектора трехдиагональной матрицы, массив
z[1...n][1...n] необходимо инициализировать на входе единичной матрицей. Если
требуются собственные вектора матрицы, сведенной к трехдиагональному виду с помощью
программы tred2, в массив z требуется загрузить соответствующий выход tred2. В
обоих случаях на выходе массив z возвращает матрицу собственных векторов, расположенных
по столбцам.
*/

/* максимальное число итераций */
#define MAXITER 30

void tqli(float *d, float *e, int n, float **z) {
    int m,l,iter,i,k;
    float s,r,p,g,f,dd,c,b;
    /* удобнее будет перенумеровать элементы e */
    for(i=2;i<=n;i++) e[i-1]=e[i];
    e[n]=0.;
    /* главный цикл идет по строкам матрицы */
    for(l=1;l<=n;l++) {
        /* обнуляем счетчик итераций для этой строки */
        iter=0;
        /* цикл проводится, пока минор 2х2 в левом верхнем углу начиная со строки l
        не станет диагональным */
        do {
            /* найти малый поддиагональный элемент, дабы расщепить матрицу */
            for(m=l;m<=n-1;m++) {
                dd=fabs(d[m])+fabs(d[m+1]);
                if((float)(fabs(e[m]+dd)==dd)) break;
            }
            /* операции проводятся, если верхний левый угол 2х2 минора еще не диагональный */
            if(m!=l) {
                /* увеличить счетчик итераций и посмотреть, не слишком ли много. Функция
                nerror завершает программу с диагностикой ошибки. */
                if(++iter>=MAXITER) {cout<<"Error!!!"; return;};
                /* сформировать сдвиг */
                g=(d[l+1]-d[l])/(2.*e[l]); r=hypot(1.,g);
                /* здесь d_m - k_s */
                if(g>=0.) g+=fabs®;
                else g-=fabs®;
                g=d[m]-d[l]+e[l]/g;
                /* инициализация s,c,p */
                s=c=1.; p=0.;
                /* плоская ротация оригинального QL алгоритма, сопровождаемая ротациями
                Гивенса для восстановления трехдиагональной формы */
                for(i=m-1;i>=l;i--) {
                    f=s*e[i]; b=c*e[i];
                    e[i+1]=r=hypot(f,g);
                    /* что делать при малом или нулевом знаменателе */
                    if(r==0.)
                    {
                        d[i+1]-=p;
                        e[m]=0.;
                        break;
                    }
                    /* основные действия на ротации */
                    s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2.*c*b; d[i+1]=g+(p=s*r); g=c*r-b;
                    /* Содержимое следующего ниже цикла необходимо опустить, если
                    не требуются значения собственных векторов */
                    for(k=1;k<=n;k++) {
                        f=z[k][i+1]; z[k][i+1]=s*z[k][i]+c*f; z[k][i]=c*z[k][i]-s*f;
                    }
                }
                /* безусловный переход к новой итерации при нулевом знаменателе и недоведенной
                до конца последовательности ротаций */
                if(r==0. && i>=l) continue;
                /* новые значения на диагонали и "под ней" */
                d[l]-=p; e[l]=g; e[m]=0.;
            }
        } while(m!=l);
    }
}

main()
{
    const int n = 3;
    int i, j;
    float z;
    float **a;
    a=new float*[n+1];
    for(i=0;i<n+1;i++) a[i]=new float[n+1];

    ifstream f("chisla.dat");
    for(i=1; i<n+1; i++)
        for(j=1; j<n+1; j++) {f>>z; a[i][j]=z;}
        f.close();

        float d[n+1],e[n+1];
        for(i=0; i<n+1; i++) d[i]=e[i]=0;
        tred2(a, n, d, e);

        printf("Diag. 3-h diag-y matritsy: ");
        for(i=1; i<n+1; i++) printf("%f ",d[i]);
        printf("\n");

        printf("Vnediag-e elementy: ");
        for(i=1; i<n+1; i++) printf("%f ",e[i]);
        printf("\n");

        // обнуление
        for(i=0;i<n+1;i++)
            for(j=0;j<n+1;j++)a[i][j]=0;

        //заполнение диагональных
        for(i=1;i<n+1;i++)a[i][i]=d[i];

        //заполнение поддиагональных
        for(i=1;i<n;i++)a[i+1][i]=e[i+1];

        tqli(d, e, n, a);

        printf("\n\n\n---------------\n");
        printf("D[]=");
        for(i=1; i<n+1; i++) printf("%f ",d[i]);
        printf("\n");

        printf("A[][]=\n");
        for(i=1;i<n+1;i++)
        {
            for(j=1;j<n+1;j++)printf("%f ", a[i][j]);
            printf("\n");
        }

        for(i=0; i<n+1; i++) delete a[i];
        delete a;
        scanf("sdsd");
}

тут nrutil.h

Это сообщение отредактировал(а) allsolovey - 13.2.2009, 19:16

Присоединённый файл ( Кол-во скачиваний: 14 )

nrutil.h___chisla.dat.rar 0,89 Kb

Dmi3ev

Дата 12.2.2009, 19:44 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 1698
Регистрация: 28.11.2007

Репутация: 16
Всего: 41

Могу дать лишь один совет, сначала почитать немного про язык программирования С++, а потом уже заниматься подобными вещами...
Совсем недавно ты не мог заполнить матрицу, я тебе тогда помог, а уже...

P.S. Конечно, читают книги, мануалы, статьи, ... по программированию только полные лохи, я согласен. Тебе это ни к чему. Но все же попробуй, может, поможет...
А еще этот раздел лучше поместить в С++ Общие вопросы, тк здесь непосредственно к Билдеру мало что относится...
А если не понимаешь алгоритм, то программу не стоит начинать писать, пока не разберешься что к чему...
А ты уже вон сколько наконтролцелил smile

И всё зря

Есть еще и раздел Алгоритмы на этом форуме...

Это сообщение отредактировал(а) Dmi3ev - 12.2.2009, 19:45

--------------------

allsolovey

Дата 12.2.2009, 22:28 (ссылка)

(нет голосов)

Загрузка ...

Empty

Профиль
Группа: Участник
Сообщений: 29
Регистрация: 3.1.2008

Репутация: нет
Всего: нет

Эти 2 функции (tred2 и tqli ) считаются стандартными в них ошибок нет.Ошибки в main(). Времени просто катастрофически нет для разбора самих методов.Разбираться придется потом все равно мне показывать завтра надо че нить науч. руководителю. smile

А теорию я буду рассказывать через 5 дней. Меня тут интересует вопрос правильно ли я инициализировал матрицу и все ли правильно с входными и выходными данными.

Dmi3ev

Дата 12.2.2009, 23:27 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 1698
Регистрация: 28.11.2007

Репутация: 16
Всего: 41

1) зачем так делать?

Код


const int n = 3;
float **a;
a=new float*[n+1];
for(i=0;i<n+1;i++) 
 a[i]=new float[n+1];

если размер заранее известен 4х4. зачем объявлять n=3? Тоже неясно...
просто:

Код


const int n = 4;
float a[n][n];

было 5 строчек, стало 2, да еще и простых))) волшебство
2) зачем объявлять это тут?

Код


int i, j;

эти переменные лучше объявлять в цикле непосредственно, когда используешь...
3) main имеет тип возвращаемого значения...

Код


int main()
{
//...
return 0;
}

4)а почему так?

Код


for(i=1; i<n+1; i++)
for(j=1; j<n+1; j++) {f>>z; a[i][j]=z;}

сразу нельзя в массив?
5) есть и попроще способ присвоить всем элементам 0, чем

Код


for(i=0; i<n+1; i++) d[i]=e[i]=0;

7) а это что?

Код


scanf("sdsd");

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

Это сообщение отредактировал(а) Dmi3ev - 12.2.2009, 23:39

--------------------

slavikzh

Дата 3.2.2011, 18:01 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 1
Регистрация: 25.11.2008

Репутация: нет
Всего: нет

Цитата

Эти 2 функции (tred2 и tqli ) считаются стандартными в них ошибок нет.

в tred2 вконце стоит:

Код

    /* ряд и колонка матрицы a преобразуются к единичной, для след. итерации */
    a[i][i]=0.;

а надо:

Код

    /* ряд и колонка матрицы a преобразуются к единичной, для след.  итерации */
    a[i][i]=1;

и еще надо вырезать кусок из main():

Код


  printf("Diag. 3-h diag-y matritsy: ");
  for(i=1; i<n+1; i++) printf("%f ",d[i]);
  printf("\n");

  printf("Vnediag-e elementy: ");
  for(i=1; i<n+1; i++) printf("%f ",e[i]);
  printf("\n");

  // обнуление
  for(i=0;i<n+1;i++)
  for(j=0;j<n+1;j++)a[i][j]=0;

  //заполнение диагональных
  for(i=1;i<n+1;i++)a[i][i]=d[i];

  //заполнение поддиагональных
  for(i=1;i<n;i++)a[i+1][i]=e[i+1];

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| C++ Builder \| Следующая тема »