2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 метод Хаусхолдера+QL-алгоритм
Сообщение11.02.2009, 18:24 


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

Код:
   
#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) {
  int l,k,j,i;
  float scale,hh,h,g,f;
  /* Проход по стадиям процесса редукции */
  for(i=n;i>=2;i--) {
    l=i-1; h=scale=0.;
    if(l>1) {
      /* вычислить шкалу */
      for(k=1;k<=l;k++) scale += fabs(a[i][k]);
      /* малая величина шкалы -> пропустить преобразование */
      if(scale==0.) e[i]=a[i][l];
      else {
       /* отмасштабировать строку и вычислить s2 в h */
        for(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(r);
        else g-=fabs(r);
        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

 Профиль  
                  
 
 Re: метод Хаусхолдера+QL-алгоритм
Сообщение03.02.2011, 18:12 


03/02/11
1
в 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];

 Профиль  
                  
 
 Re: метод Хаусхолдера+QL-алгоритм
Сообщение22.08.2011, 06:02 


22/08/11
13
пишу, подобную программу, но на c#. помогите код отладить. ошибку выдает, не могу найти ее
код: [ скачать ] [ спрятать ]
Используется синтаксис C#
using System;
using System.IO;
namespace prog
{

    class program
    {
        static float[,] Input(out float n)
        {
            StreamReader file = new StreamReader("matrix.txt");
            string s = file.ReadToEnd();
            file.Close();
            string[] строка = s.Split('\n');
            string[] столбцы = строка[0].Split(' ');
            float[,] a = new float[строка.Length, столбцы.Length];
            float t = 0;
            n = 0;
            for (int i = 0; i < строка.Length; i++)
            {
                столбцы = строка[i].Split(' ');
                for (int j = 0; j < столбцы.Length; j++)
                {
                    t = float.Parse(столбцы[j]);
                    a[i, j] = t;
                }
            }
            return a;
        }
        static void Print(float[,] a)
        {
            for (int i = 0; i < a.GetLength(0); ++i, Console.WriteLine())
                for (int j = 0; j < a.GetLength(1); ++j)
                    Console.Write("{0} ", a[i, j]);
        }
        static void Main()
        {
            int n,i;
            float p;

            float[,] myArray = Input(out p);
            Console.WriteLine("Исходный массив:");
            Print(myArray);
            n = myArray.GetLength(0);
            float[] d = new float[n+ 1];
            float[] e = new float[n + 1];
            for (i = 0; i < n + 1; i++)
                d[i]=e[i]=0;
           
            tred2(myArray, n, d, e);
            Print(myArray);
            tqli(d,e,n,myArray);

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

                            //a[j][i]=a[i][j]/h;
                             /* сформировать элемент Au (в g) */
                            g=0f;
                            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]=0f;
             /* Все содержание цикла, кроме одной инструкции, можно опустить, если не
     требуются собственные вектора */

            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]=1;
                for(j=1;j<=1;j++)
                    a[j,i]=a[i,j]=0f;
            }
   
        }
        static void tqli(float []d, float []e, int n, float [,]z)
        {
            int MAXITER =30;
            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]=0f;
            /* главный цикл идет по строкам матрицы */
            for(l=1;l<=n;l++)
            {
                /* обнуляем счетчик итераций для этой строки */
                iter=0;
                /* цикл проводится, пока минор 2х2 в левом верхнем углу начиная со строки l
               не станет диагональным */

            do {
              /* найти малый поддиагональный элемент, дабы расщепить матрицу */
              for(m=l;m<=n-1;m++)
              {
                dd=Math.Abs(d[m])+Math.Abs(d[m+1]);
                if((float)(Math.Abs(e[m]+dd))==dd) break;
              }
              /* операции проводятся, если верхний левый угол 2х2 минора еще не диагональный */
              if(m!=l) {
                /* увеличить счетчик итераций и посмотреть, не слишком ли много.  */
                if(++iter>=MAXITER)
                    Console.WriteLine("Too many iterations in program");
                /* сформировать сдвиг */
                g=(d[l+1]-d[l])/(2f*e[l]); r=(float)hypot(1f,g);
                /* здесь d_m - k_s */
                if(g>=0f) g+=Math.Abs(r);
                else g-=Math.Abs(r);
                g=d[m]-d[l]+e[l]/g;
                /* инициализация s,c,p */
                s=c=1f; p=0f;
                /* плоская ротация оригинального QL алгоритма, сопровождаемая ротациями
                   Гивенса для восстановления трехдиагональной формы */

                for(i=m-1;i>=l;i--)
                {
                  f=s*e[i]; b=c*e[i];
                  e[i+1]=r=(float)hypot(f,g);
                  /* что делать при малом или нулевом знаменателе */
                  if(r==0f) {d[i+1]-=p; e[m]=0f; break;}
                  /* основные действия на ротации */
                  s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2f*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==0f && i>=l)
                    continue;
                /* новые значения на диагонали и "под ней" */
                d[l]-=p; e[l]=g; e[m]=0f;
              }
            } while(m!=l);
          }
        }
        static double hypot(double x, double y)
        {
            double a = Math.Sqrt(x * x + y * y);
           
            return a;
        }
   

    }
}
 

 Профиль  
                  
 
 Re: метод Хаусхолдера+QL-алгоритм
Сообщение22.08.2011, 06:18 


06/04/09
156
Воронеж
какую ошибку?

 Профиль  
                  
 
 Re: метод Хаусхолдера+QL-алгоритм
Сообщение23.08.2011, 08:54 


22/08/11
13
в конце tred2
Используется синтаксис C#
/* эта инструкция остается */
                d[i]=a[i,i];
 

индекс находится за границами массива

 Профиль  
                  
 
 Re: метод Хаусхолдера+QL-алгоритм
Сообщение31.08.2011, 04:24 


22/08/11
13
хелп, я уже код пересмотрел не найду никак ошибку

 Профиль  
                  
 
 Re: метод Хаусхолдера+QL-алгоритм
Сообщение31.08.2011, 19:16 
Заслуженный участник


09/08/09
3438
С.Петербург
Насколько я понял, у Вас в tred2 массив a имеет размерность [n, n] (т. е. "последний" элемент - a[n-1, n-1]), поэтому при i=n a[i, i] ссылается за границы массива (индексы в c# начинаются с нуля).

 Профиль  
                  
 
 Re: метод Хаусхолдера+QL-алгоритм
Сообщение19.10.2015, 17:47 
Аватара пользователя


19/10/15
1
:lol: Посмотрел код, полагаю, что в цикле после комментария
Код:
/* Все содержание цикла, кроме одной инструкции, можно опустить, если не
    требуются собственные вектора */

надо просто вместо "<="

Код:
for (i = 1; i <= n; i++)
            {


написать просто "<":

Код:
for (i = 1; i < n; i++)
            {


Так программа, как минимум, запускается без ошибок.


Однако, есть такая проблема, что при этом на любой входящей матрице что tred2, что tqli выдаёт на выходе всегда ту же самую матрицу, что и на входе, лишь с парой отличий: последний (самый правый нижний элемент) всегда будет один, и третий элемент второй строки в любом случае зануляется.

К автору темы в связи с этим вопрос по программе. Возможно, за четыре года удалось добиться прогресса? Потому как самому нужно тоже вычисление собственных чисел и векторов на C#, а эта реализация - ну почти последний шанс. :-)

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

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



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

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


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

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