Форум программистов «Весельчак У»
  *
Добро пожаловать, Гость. Пожалуйста, войдите или зарегистрируйтесь.
Вам не пришло письмо с кодом активации?

  • Рекомендуем проверить настройки временной зоны в вашем профиле (страница "Внешний вид форума", пункт "Часовой пояс:").
  • У нас больше нет рассылок. Если вам приходят письма от наших бывших рассылок mail.ru и subscribe.ru, то знайте, что это не мы рассылаем.
   Начало  
Наши сайты
Помощь Поиск Календарь Почта Войти Регистрация  
 
Страниц: [1]   Вниз
  Печать  
Автор Тема: Входные и выходные данные  (Прочитано 12991 раз)
0 Пользователей и 1 Гость смотрят эту тему.
allsolovey
Гость
« : 12-02-2009 12:38 » 

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

Код:
#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 & chisla.dat.rar (0.89 Кб - загружено 885 раз.)
« Последнее редактирование: 13-02-2009 05:35 от LogRus » Записан
Алексей++
глобальный и пушистый
Глобальный модератор

ru
Offline Offline
Сообщений: 13


« Ответ #1 : 12-02-2009 17:31 » 

код редкий по ужасности на морду лица

Напиши алгоритм, никто не будет с таким кошмаром сидеть разбираться. Заодно, вдруг, и сам разберёшься Отлично
Записан

Антон (LogRus)
Глобальный модератор

ru
Offline Offline
Пол: Мужской
Внимание! Люблю сахар в кубиках!


WWW
« Ответ #2 : 13-02-2009 05:37 » 

поправил форматирование, но блин мало помогло, особенно плохо читаются if в одну строчку и несколько операций в одну строчку, опять же любовь к паскалевскому стилю угнетает
Записан

Странно всё это....
Алексей++
глобальный и пушистый
Глобальный модератор

ru
Offline Offline
Сообщений: 13


« Ответ #3 : 13-02-2009 05:58 » 

во первых, матрицу - сразу в отдельный класс засунуть. Если класс отлажен будет, то остальной код сократится раз в 10 и будет более читабельным и более рабочим Улыбаюсь


однострочные комментарии лучше писать не так

/* ывапывапвапрпатитмит */

а так

// ывапывапвапрпатитмит

- и красивее, и удобнее, есоли что, большой блок закомментировать через /**/
« Последнее редактирование: 13-02-2009 06:00 от Алексей1153++ » Записан

allsolovey
Гость
« Ответ #4 : 13-02-2009 16:12 » 

Спасибо за форматирование,только зашел отформатировал бы сам.Функции считаются стандартными (проверено) так что в них не должно быть ошибки,ошибка в main(). Только вот где...
Записан
Алексей++
глобальный и пушистый
Глобальный модератор

ru
Offline Offline
Сообщений: 13


« Ответ #5 : 13-02-2009 18:32 » 

allsolovey, давай, для начала напиши класс, инкапсулирующий матрицу N*N . После этого будем искать ошибки - это будет уже несложно
Записан

allsolovey
Гость
« Ответ #6 : 13-02-2009 21:09 » new

если бы умел написал бы))
Записан
Алексей++
глобальный и пушистый
Глобальный модератор

ru
Offline Offline
Сообщений: 13


« Ответ #7 : 13-02-2009 21:40 » 

в сыром и нетестированном виде, примерно так
Код:
class CMtx
{
  float* m_ppM;
  int m_nDim;
public:
  CMtx()
  {
    m_ppM=0;
    m_nDim=0;
  }

  CMtx(int N)
  {
    m_ppM=0;
    m_nDim=0;
    Create(int N);
  }

  bool Create(int N)
  {
    Destroy();

    if(!N)
    {
       return false;
    }

    m_ppM=new float[N*N];
    if(!m_ppM)
    {
       return false;
    }

    m_nDim=N;
    return true;
  }

  ~CMtx()
  {
     Destroy();
  }
 
  void Destroy()
  {
    if(m_ppM)
    {
       delete m_ppM;
    }
    m_ppM=0;
    m_nDim=0;
  }

};

осталось нужные методы и операторы прописать Улыбаюсь
Записан

Страниц: [1]   Вверх
  Печать  
 

Powered by SMF 1.1.21 | SMF © 2015, Simple Machines