Форум программистов, компьютерный форум CyberForum.ru

С++ для начинающих

Войти
Регистрация
Восстановить пароль
 
Рейтинг: Рейтинг темы: голосов - 17, средняя оценка - 4.76
nullpointer
46 / 46 / 5
Регистрация: 30.03.2009
Сообщений: 521
#1

Решение СЛАУ большой размерности методом сопряженных градиентов - C++

19.12.2012, 17:28. Просмотров 2394. Ответов 3

Всем првиет! Возникла проблемка с методом сопряженных градиентов. Если задавать самому значения матрицы и правой части, то все решается отлично. Сделал rand'омное заполнение. В итоге получаю, что программа работает только с матрицами размера 2*2 и то не со всеми, даже 3*3 она уже не считает правильно, в результате выводит ответ в виде -1.#IND. А мне нужно, чтобы она работала с матрицами большого порядка, мне надо узнать время выполнения алгоритма. В чем можнт быть проблема? Я даже double заменил на long double и все равно не помогает.
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#include<stdio.h> 
#include <time.h>
#include <conio.h>
#include <iostream>
using namespace std;
 
#define M 2 //размерность системы
 
#define E 0.1 //точность вычислений
 
 
long double Xk[M], Zk[M];
long double Rk[M], Sz[M], alpha, beta, mf;
long double Spr, Spr1, Spz;
long double A[M][M];//={2.0, -1.0, 0.0, -1.0, 5.0, 3.0, 0.0, 3.0, 1.0};
long double F[M] ;//= {1.0, 7.0, 4.0};
 
int main(void){
    int  i, j, kl = 1;
 
    srand(time(0));
    for (i = 0; i < M; i++)
    {
        for (j = 0; j < M; j++){
     A[i][j] = rand() % 2 - rand() % 1; // заполняем массив случайными значениями в диапазоне от -1 до 1 включительно
     cout << A[i][j] << " "; // печать элементов одномерного массива
        }
        cout << endl;
    }
 
    for (j = 0; j < M; j++){
     F[j] = rand() % 2 - rand() % 2; // заполняем массив случайными значениями в диапазоне от -1 до 1 включительно
     cout << F[j] << " "; // печать элементов одномерного массива 
        }
 
    clock_t time;
    time = clock();//начальное время
 
/* Вычисляем сумму квадратов элементов вектора F*/  
    for(mf = 0,i = 0; i < M; i++){
        mf += F[i] * F[i];
      }
 
 
/* Задаем начальное приближение корней. В Хk хранятся значения корней
 * к-й итерации. */
    for(i = 0; i < M; i++){
        Xk[i] = 0.2;
    } 
 
/* Задаем начальное значение r0 и z0. */
    for(i = 0; i < M; i++){
        for(Sz[i]=0,j = 0; j < M; j++)
            Sz[i] += A[i][j] * Xk[j];
      Rk[i] = F[i] - Sz[i];
      Zk[i] = Rk[i];  
      }
    
    do{ 
  /* Вычисляем числитель и знаменатель для коэффициента
   * alpha = (rk-1,rk-1)/(Azk-1,zk-1) */
        Spz = 0;
        Spr = 0;
        for(i = 0; i < M; i++){ 
            for(Sz[i] = 0, j = 0; j < M; j++){
                Sz[i] += A[i][j] * Zk[j];
            }
          Spz += Sz[i] * Zk[i];
          Spr += Rk[i] * Rk[i];
          }
        alpha = Spr / Spz;             /*  alpha    */
    
 
/* Вычисляем вектор решения: xk = xk-1+ alpha * zk-1, 
    вектор невязки: rk = rk-1 - alpha * A * zk-1 и числитель для betaa равный (rk,rk) */
        Spr1 = 0;
        for(i = 0; i < M; i++){
            Xk[i] += alpha * Zk[i];
            Rk[i] -= alpha * Sz[i];
            Spr1 += Rk[i] * Rk[i];
            cout << "Iter #" << kl;
            cout << " " << "X[" << i << "] = " << Xk[i] << endl;
        }
        cout << endl;
        kl++;
 
/* Вычисляем  beta  */
        beta = Spr1 / Spr; 
 
/* Вычисляем вектор спуска: zk = rk+ beta * zk-1 */
        for(i = 0; i < M; i++)
            Zk[i] = Rk[i] + beta * Zk[i];       
      }
/* Проверяем условие выхода из итерационного цикла  */
    while(Spr1 / mf > E * E);
 
    time = clock() - time;//время выполнения
    cout << "Time result = " << (double)time << "ms" << endl;
    system("pause");
    return 0;
  }
Добавлено через 59 минут
UPD.Понял почему не считает! Матрица же симметричная и положительно определенная должна быть!!! Теперь возникает другой вопрос: как рандомно заполнить симметричную матрицу?

Добавлено через 1 час 1 минуту
UPD.Все. Во всем разобрался
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
19.12.2012, 17:28     Решение СЛАУ большой размерности методом сопряженных градиентов
Посмотрите здесь:

Создать СЛАУ размерности 1000х1000 в С++, чтобы система имела только одно решение - C++
Здравствуйте, как можно создать СЛАУ большой размерности в виде матрицы в c++, чтобы система имела только одно решение?

Решение СЛАУ методом вращения - C++
Доброго времени суток, товарищи. Имеется задание: дано интегральное уравнение: U(x) + I ( K(x,s)U(s)ds= f(x) ( I - определённый...

Решение СЛАУ методом прогонки - C++
Добрый день,нужно решить СЛАУ методом прогонки,вот СЛАУ и ее решение. коэфф перед х1считается как за а1. перед х2считается...

Решение СЛАУ методом отражений - C++
Добрый вечер :) Было две темы &quot;Решение СЛАУ методом отражений&quot;, но нет реализации) У меня есть код , в котором реализовано два метода...

Решение СЛАУ методом Зейделя - C++
Методом Зейделя решить с точностью до 0.001 заданную систему уравнений 3.3*x1+2.1*x2+2.8*x3=0.8 4.1*x1+3.7*x2+4.8*x3=5.7 ...

Решение СЛАУ методом Гаусса - C++
Всем привет))) Друзья, сколько я уже лазию по инету и по этому форуму, я никак не могу найти правильно работающей программки. Может кто...

После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
аия
0 / 0 / 0
Регистрация: 25.12.2010
Сообщений: 4
02.05.2013, 17:39     Решение СЛАУ большой размерности методом сопряженных градиентов #2
}{@k, не могли бы вы скинуть своё решение, пожалуйста) очень нужно
nullpointer
46 / 46 / 5
Регистрация: 30.03.2009
Сообщений: 521
05.05.2013, 21:15  [ТС]     Решение СЛАУ большой размерности методом сопряженных градиентов #3
аия, держите, только там еще и метод Гаусса реализован для сравнения скорости вычисления, но его несложно удалить
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
#include <stdio.h> 
#include <ctime>
#include <conio.h>
#include <iostream>
using namespace std;
 
#define M 50 //размерность системы
#define E 0.01 //точность вычислений
 
float Xk[M], Zk[M];
float Rk[M], Sz[M], alpha, beta, mf;
float Spr, Spr1, Spz;
float A[M][M];
float F[M] ;
 
float conj_grad(float A[M][M], float F[M]){
    int i,j,kl=1;
    double max_iter = 100000;
/* Вычисляем сумму квадратов элементов вектора F*/  
    for(mf = 0,i = 0; i < M; i++){
        mf += F[i] * F[i];
      }
 
 
/* Задаем начальное приближение корней. В Хk хранятся значения корней
 * к-й итерации. */
    for(i = 0; i < M; i++){
        Xk[i] = 0.2;
    } 
 
/* Задаем начальное значение r0 и z0. */
    for(i = 0; i < M; i++){
        for(Sz[i]=0,j = 0; j < M; j++)
            Sz[i] += A[i][j] * Xk[j];
      Rk[i] = F[i] - Sz[i];
      Zk[i] = Rk[i];  
      }
    int Iteration = 0;
    do{ 
        Iteration++;
  /* Вычисляем числитель и знаменатель для коэффициента
   * alpha = (rk-1,rk-1)/(Azk-1,zk-1) */
        Spz = 0;
        Spr = 0;
        for(i = 0; i < M; i++){ 
            for(Sz[i] = 0, j = 0; j < M; j++){
                Sz[i] += A[i][j] * Zk[j];
            }
          Spz += Sz[i] * Zk[i];
          Spr += Rk[i] * Rk[i];
          }
        alpha = Spr / Spz;             /*  alpha    */
    
 
/* Вычисляем вектор решения: xk = xk-1+ alpha * zk-1, 
    вектор невязки: rk = rk-1 - alpha * A * zk-1 и числитель для betaa равный (rk,rk) */
        Spr1 = 0;
        for(i = 0; i < M; i++){
            Xk[i] += alpha * Zk[i];
            Rk[i] -= alpha * Sz[i];
            Spr1 += Rk[i] * Rk[i];
            //cout << "Iter #" << kl;
            //cout << " " << "X[" << i << "] = " << Xk[i] << endl;
        }
        //cout << endl;
        kl++;
        cout << "Rko-vo iter" << kl << endl;
 
/* Вычисляем  beta  */
        beta = Spr1 / Spr; 
 
/* Вычисляем вектор спуска: zk = rk+ beta * zk-1 */
        for(i = 0; i < M; i++)
            Zk[i] = Rk[i] + beta * Zk[i];       
      }
/* Проверяем условие выхода из итерационного цикла  */
    while(Spr1 / mf > E * E && Iteration < max_iter);
 
    return 0;
}
 
void GlavElem(int k, float mas[][M + 1], int otv[])
{
  int i, j, i_max = k, j_max = k;
  float temp;
  //Ищем максимальный по модулю элемент
  for (i = k; i < M; i++)
  {
      for (j = k; j < M; j++)
      {
          if (fabs(mas[i_max][j_max]) < fabs(mas[i][j]))
          {
              i_max = i;
              j_max = j;
          }
      }
  }
  //Переставляем строки
  for (j = k; j < M + 1; j++)
  {
      temp = mas[k][j];
      mas[k][j] = mas[i_max][j];
      mas[i_max][j] = temp;
  }
  //Переставляем столбцы
  for (i = 0; i < M; i++)
  {
      temp = mas[i] [k];
      mas[i] [k] = mas[i] [j_max];
      mas[i] [j_max] = temp;
  }
  //Учитываем изменение порядка корней
  i = otv[k];
  otv[k] = otv[j_max];
  otv[j_max] = i;
}
 
float gauss(float A[M][M], float F[M]){
    int i, j, k; int otv[M];
    float x[M], mas[M][M + 1];
 
    for (i = 0; i < M; i++){
        for (j = 0; j < M; j++){
        mas[i][j] = A[i][j];
        }
        mas[i][M] = F[i];
    }
    
    for (i = 0; i < M; i++)
  {
      otv[i] = i;
  }
    for (k = 0; k < M; k++)
  { 
      //На какой позиции должен стоять главный элемент
      GlavElem(k, mas, otv); //Установка главного элемента
          
 
    for (j = M; j >= k; j--)
    {
        mas[k][j] /= mas[k][k];
    }
 
    for (i = k + 1; i < M; i++)
    {
        for (j = M; j >= k; j--)
        {
            mas[i][j] -= mas[k][j] * mas[i][k];
        }
    }
  }
  //Обратный ход
  for (i = 0; i < M; i++)
  {
      x[i] = mas[i][M];
  }
  for (i = M - 2; i >= 0; i--)
  {
      for (j = i + 1; j < M; j++)
      {
          x[i] -= x[j] * mas[i][j];
      }
  }
 
    return 0;
}
 
int main(void){
    int  i, j, kl = 1;
 
    srand(time(0));
    
    for (i = 0; i < M; i++){
        for (j = 0; j < M; j++){
            A[i][i] = rand() % 2;
            if ((i != j) && (i < j)){
                A[i][j] = rand() % 2 - rand() % 2;//заполняем массив случайными значениями
            }
                A[j][i] = A[i][j];
        }
    }
 
    for (j = 0; j < M; j++){
     F[j] = rand() % 2 - rand() % 2; // заполняем массив случайными значениями
    }
 
    clock_t time_conj_grad; 
    clock_t time_gauss;
    
    time_conj_grad = clock();//начальное время  
    conj_grad(A,F); 
    time_conj_grad = clock() - time_conj_grad;//время выполнения
    cout << "Running time the conjugate gradient algorithm = " << (double)time_conj_grad << "ms" << endl;
    
    cout << endl <<endl;
 
    time_gauss = clock();//начальное время
    gauss(A,F);
    time_gauss = clock() - time_gauss;//время выполнения
    cout << "Running time Gauss = " << (double)time_gauss << "ms" << endl << endl;
    system("pause");
    return 0;
  }
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
06.05.2013, 13:11     Решение СЛАУ большой размерности методом сопряженных градиентов
Еще ссылки по теме:

Решение слау методом релаксации - C++
Доброго времени суток. В качестве исходных данных имеются слау большой размерности (koeff__100.7z) и решение этой слау (result__100.7z). ...

Решение СЛАУ методом отражений - C++
Всем привет. Задали писать курсач. Нужно реализовать метод отражения. Предусмотреть ввод числа уравнений, матрицы коэффициентов и вектора...

Решение СЛАУ методом Гаусса - C++
помогите разобраться!!ВЫдает 85 ошибок!!!! #include &lt;stdio.h&gt; #include &lt;conio.h&gt; #include &lt;tchar.h&gt; #include &lt;iostream&gt; #include...

Решение СЛАУ методом Якоби - C++
Решить СЛАУ методом Якоби. Вывести значения решения, график зависимости нормы невязки от номера итерации и его значение, при котором...

Решение СЛАУ методом Гаусса - C++
У меня программа для решения слау методом Гаусса с выбором главного элемента по столбцам. Что-то мне не нравятся ответы, гляньте может в...

Решение СЛАУ методом Крамера - C++
Доброго времени суток! Пишу курсовую по методу Крамера. Вроде ничего сложного, но программированием я занимаюсь не часто, скачал с нета...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
аия
0 / 0 / 0
Регистрация: 25.12.2010
Сообщений: 4
06.05.2013, 13:11     Решение СЛАУ большой размерности методом сопряженных градиентов #4
}{@k, спасибо!
Yandex
Объявления
06.05.2013, 13:11     Решение СЛАУ большой размерности методом сопряженных градиентов
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Рейтинг@Mail.ru