Входные и выходные данные

Тема в разделе "Borland C++ Builder & Kylix", создана пользователем allsolovey, 12 фев 2009.

Статус темы:
Закрыта.
  1. allsolovey

    allsolovey Гость

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

    Код (Text):
    #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®;
    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
     

    Вложения:

  2. allsolovey

    allsolovey Гость

    Помогите КТО НИБУДЬ!!!!!!!! :wacko:
     
Статус темы:
Закрыта.

Поделиться этой страницей