Гаусс для коллекции

@taras atavin · Регистрация: 24.11.2009

C++

double abs(double x)
{
 if (x<0)
 {
  return -x;
 }
 return x;
}
void Gauss(double *a, double *b, double *&x, unsigned short int n) // а - матрица коэффициентов, строки подряд, b - свободные члены, x - решение.
{
 unsigned hort int i,j,k,t;
 double kof,s;
 if (x) // Если какой-то вектор решений уже дан
 {
  delete [] x; // Сотрём его нафиг!
 }
 for (i=n-1; i>0; --i) // Цикл по уравнениям - вычитаемым
 {
  for (t=i, j=i-1; j>=0; --j) Ищем строку с максимальным в i-том столбце коэффициентом.
  {
   if (abs(a[i+t*n])<abs(a[i+j*n])
   {
    t=j;
   }
  }
  if (a[i+t*n]==0.0)
  {
   return;
  }
  if (t!=i) // Если она не i-тая
  {
   for (k=n-1; k>=0; --k) // Меняем её с i-той
   {
    a[k+t*n]&=a[k+i*n];
    a[k+i*n]&=a[k+t*n];
    a[k+t*n]&=a[k+i*n];
   }
  }
  for (j=i-1; j>=0; --j) // Цикл пл уравнениям-уменьшаемым
  {
   kof=a[i+j*n]/a[i+i*n]; // kof=a[j][i]/a[i][i];
   for (a[i+j*n]=0.0, b[j]-=b[i]*kof, k=i-1; k>=0; --k) // Цикл по столбцам
   {
    a[k+j*n]-=a[k+i*n]*kof; // a[k][j]-=a[k][i]*kof
   }
  }
 }
 x=new double [n];
 if (x)
 {
  for (i=0; i<n; ++i)
  {
   for (s=0.0, j=i-1, j>=0; --j)
   {
    s+=a[j+i*n]*x[j]; //s+=a[j][i]*x[j]
   }
   x[i]=(b[i]-s)/a[i+i*n]; // x[i]=(b[i]-s)/a[i][i];
  }
 }
}

Добавлено через 16 часов 59 минут
Решенная здесь задача.
Дано: матрица a коэффициентов системы уравнений в в виде линерализованного массива (a[0][0], a[1][0], a[2][0]... a[n-1][0], a[0][1], a[1][1], a[2][1]...a[n-1][1], a[0][2], a[1][2], a[2][2]...a[n-1][2]...a[0][n-1], a[1][n-1], a[2][n-1]...a[n-1][n-1]), вектор-столбец b свободных членов той же системы уравнений и число уравнений. Матрица квадратная, число элементов вектора равно числу строк матрицы. Система эквивалентна матричному уравнению a*x=b. Требуется найти: вектор x решений системы уравнений. Если система не имеет однозначного решения, вернуть пустой указатель на массив x.

@vet · 16.03.2010, 05:15

Сообщение от taras atavin

for (k=n-1; k>=0; --k) // Меняем её с i-той
{
a[k+t*n]&=a[k+i*n];
a[k+i*n]&=a[k+t*n];
a[k+t*n]&=a[k+i*n];
}

taras atavin, а каким компилятором пользуешься?
У меня в этом блоке и Dev-C++ и VS 2008 пишут одно и тоже: что нельзя использовать тип double с операцией поразрядного присваивания &= ...?

@taras atavin · 16.03.2010, 05:56 **[ТС]**

Для инта почему-то работало, вот я и заменил. Придётся заменить назад: в оргинале был обмен через temp.

C++

double abs(double x)
{
 if (x<0)
 {
  return -x;
 }
 return x;
}
void Gauss(double *a, double *b, double *&x, unsigned short int n) // а - матрица коэффициентов, строки подряд, b - свободные члены, x - решение.
{
 unsigned hort int i,j,k,t;
 double kof,s;
 double temp;
 if (x) // Если какой-то вектор решений уже дан
 {
  delete [] x; // Сотрём его нафиг!
 }
 for (i=n-1; i>0; --i) // Цикл по уравнениям - вычитаемым
 {
  for (t=i, j=i-1; j>=0; --j) Ищем строку с максимальным в i-том столбце коэффициентом.
  {
   if (abs(a[i+t*n])<abs(a[i+j*n])
   {
    t=j;
   }
  }
  if (a[i+t*n]==0.0)
  {
   return;
  }
  if (t!=i) // Если она не i-тая
  {
   for (k=n-1; k>=0; --k) // Меняем её с i-той
   {
    temp=a[k+t*n];
    a[k+t*n]=a[k+i*n];
    a[k+i*n]=temp;
   }
  }
  for (j=i-1; j>=0; --j) // Цикл пл уравнениям-уменьшаемым
  {
   kof=a[i+j*n]/a[i+i*n]; // kof=a[j][i]/a[i][i];
   for (a[i+j*n]=0.0, b[j]-=b[i]*kof, k=i-1; k>=0; --k) // Цикл по столбцам
   {
    a[k+j*n]-=a[k+i*n]*kof; // a[k][j]-=a[k][i]*kof
   }
  }
 }
 x=new double [n];
 if (x)
 {
  for (i=0; i<n; ++i)
  {
   for (s=0.0, j=i-1, j>=0; --j)
   {
    s+=a[j+i*n]*x[j]; //s+=a[j][i]*x[j]
   }
   x[i]=(b[i]-s)/a[i+i*n]; // x[i]=(b[i]-s)/a[i][i];
  }
 }
}

@vet · 16.03.2010, 05:58

Сообщение от taras atavin

Для инта почему-то работало, вот я и заменил.

Все верно поразрядные либо с int либо,что чаще, с unsigned int используются

@vet · 22.03.2010, 14:45

taras atavin, код немного был с помарками и немного глючил =))), я вроде исправил,но щас другая проблема, стал корни считать как - то в разнобой, вот пример(но код оч хороший, нули на глав диаг и те "пережевывает", что не во всех реализациях Гаусса бывает):

C++

#include <iostream>
#include <cmath>
#include <ctime>
 
using namespace std;
 
void Gauss(double *a, double *b, double *x,int n) // à - ìàòðèöà êîýôôèöèåíòîâ, ñòðîêè ïîäðÿä, b - ñâîáîäíûå ÷ëåíû, x - ðåøåíèå.
{
 int i,j,k,t;
 double kof,s;
 double temp;
 for (i=n-1; i>0; --i) // Öèêë ïî óðàâíåíèÿì - âû÷èòàåìûì
 {
  for (t=i, j=i-1; j>=0; --j) //Èùåì ñòðîêó ñ ìàêñèìàëüíûì â i-òîì ñòîëáöå êîýôôèöèåíòîì.
  {
   if (fabs(a[i*n+t])<fabs(a[i*n+j]))
   {
    t=j;
   }
  }
  if (a[i*n+t]==0.0)
  {
   return;
  }
  if (t!=i) // Åñëè îíà íå i-òàÿ
  {
   for (k=n-1; k>=0; --k) // Ìåíÿåì å¸ ñ i-òîé
   {
    temp=a[k*n+t];
    a[k*n+t]=a[k*n+i];
    a[k*n+i]=temp;
   }
  }
  for (j=i-1; j>=0; --j) // Öèêë ïo óðàâíåíèÿì-óìåíüøàåìûì
  {
   kof=a[i+j*n]/a[i*n+i]; // kof=a[j][i]/a[i][i];
   for (a[i+j*n]=0.0, b[j]-=b[i]*kof, k=i-1; k>=0; --k) // Öèêë ïî ñòîëáöàì
   {
    a[k+j*n]-=a[k+i*n]*kof; // a[k][j]-=a[k][i]*kof
   }
  }
 }
 x=new double [n];
 if (x)
 {
  for (i=0; i<n; ++i)
  {
   for (s=0.0, j=i-1; j>=0; --j)
   {
    s+=a[j+i*n]*x[j]; //s+=a[j][i]*x[j]
   }
   x[i]=(b[i]-s)/a[i*n+i]; // x[i]=(b[i]-s)/a[i][i];
   //cout << x[i] << "\n";
  }
 }
 cout << "\n";
 for (i=0; i<n; ++i)
cout << "x " << x[i] << endl;
}
 
 
main()
{
srand((unsigned)time(NULL));
      int n = 6;
     double gauss2[6] ={1,5,3,7,6,3};
      double M[36]={0, 0, 1,0, 0, 0,0, 0, 0, -1, 0, 0,0, 8, 0, 0, 0, 4,0, 0, 0, 0, 6,0, 0, 6, 0, 0, 0, 0,1, 0, 0, 0, 0, 0 };
        
      double *y = new double[n];
 
    /*  for(int i=0; i<n; i++)
        {
            for(int j=0; j<n; j++)
                {
                    //cin >> M[i+j*n];
                     M[i+j*n] = rand()/100;
//cout << "gauss1["<< i << "][" << j << "]=";
//cin >> gauss1[i][j];
                }
        }
*/
for(int i=0; i<n; i++)
        {
            for(int j=0; j<n; j++)
                {
                     cout << M[i*n+j] << "  ";
                }
                cout << gauss2[i] << endl;
        }   
      Gauss(M,gauss2,y,n);
       cin.get();
}

Ответ такой: 1, -5, -1.25, 1.17, 1, 3
А должен быть: 3, 1, 1, -5, 1.17, -1.25 (подстановкой в ур-ния проверял)
Весь день седня долбаюсь никак не хочет в нормальном порядке корни искать =(((
Может есть идеи по этому поводу???

Добавлено через 23 часа 56 минут
Может есть мысли по этому поводу???

26.03.2010, 18:31

код обратной матрицы бы

@taras atavin · 17.06.2010, 10:08 **[ТС]**

там ещё

Сообщение от taras atavin

x[i]=(b[i]-s)/a[i+i*n]; // x[i]=(b[i]-s)/a[i][i];

надо заменить на

C++

   if (a[i+i*n]<>0.0)
   {
    x[i]=(b[i]-s)/a[i+i*n]; // x[i]=(b[i]-s)/a[i][i];
   }
   else
   {
    delete [] x; 
    x=NULL;
    return;
   }

@vet · 17.06.2010, 14:55

Сообщение от taras atavin

C++

1
2
3

   if (a[i+i*n]<>0.0)
 
   }

А что такое <> ? в С/С++ нет такого условного оператора

@taras atavin · 18.06.2010, 11:07 **[ТС]**

Сообщение от vet

Ответ такой: 1, -5, -1.25, 1.17, 1, 3
А должен быть: 3, 1, 1, -5, 1.17, -1.25 (подстановкой в ур-ния проверял)

Я немножко пропустил, что вместе со строками a должны точно также переставляться и элементы b:

C++

for (k=n-1; k>=0; --k) // Меняем её с i-той
{
 temp=a[k+t*n];
 a[k+t*n]=a[i+t*n];
 a[i+t*n]=temp;
}
temp=b[t];
b[t]=b[i];
b[i]=temp;

@easybudda · 02.10.2010, 22:23

Сообщение от taras atavin

double abs(double x)

а чем стандартная fabs() не угодила?

@Vergil · 12.10.2010, 08:27

Выкладываю свой не консольный вариант, сделанный по алгоритму преподавателя. На особую гениальность не претендует, но работает. Метод Гаусса с выбором наибольшего элемента столбца. Возможно не правильная проверка на линейность. Но у меня приняли. Кажется, здесь есть схожий, а, возможно, и такой же вариант(самый первый пост). Если есть ошибки в программе, прошу мне написать. Сделано на Borland Developer Studio(C++ Builder).

16.10.2010, 23:58

поделюсь своими трудами...
код "чуть" короче
имеются недоработки (было лень пришивать лишние проверки)

@Крюгер · 10.11.2013, 09:38

Сообщение от vet

А что такое <> ? в С/С++ нет такого условного оператора

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

@CEBEP · 14.04.2014, 14:30

Честно говоря офигел... Искал вчера несколько часов решение СЛАУ общего вида. Все либо зависали либо были для частных случаев. В итоге решил довести до ума то, что валялось здесь. Я выяснил (методом внимательного всматривания), что правильным порядок иксов в ответе нарушается из-за изменения порядка строк в A (кстати, так и не понял, почему при этом в b ничего менять не надо и ответ правильный). Быстро дописал костыль, проверил, всё идеально.
Работает, естественно, не слишком быстро но для размерности меньше 1000 сойдёт... Решает меньше секунды...

C++

#include <stack>
void EDM::Gauss(double** A, double *b, double *x, int n)
{
    std::stack<std::pair<size_t,size_t> > swaps;
    int i, j, k, t;
    double kof, s;
    for (i = n - 1; i > 0; --i) {
        for (t=i, j=i-1; j >= 0; --j) {//Поиск максимума в столбце
            if (fabs(A[i][t]) < fabs(A[i][j])) {
                t = j;
            }
        }
        if (A[i][t] == 0.0) {
            return;
        }
        if (t != i) {
            for (k = n - 1; k >= 0; --k) {
                std::swap(A[k][t],A[k][i]);
            }
            swaps.push(std::pair<size_t, size_t>(i,t));
        }
        for (j=i-1; j>=0; --j) {
            kof = A[j][i] / A[i][i];
            A[j][i] = 0.0;
            b[j] -= b[i] * kof;
            for (k = i - 1; k >= 0; --k) {
                A[j][k] -= A[i][k] * kof;
            }
        }
    }
    for (i = 0; i < n; ++i)
    {
        s = 0.0;
        for (j = i - 1; j >= 0; --j) {
            s += A[i][j] * x[j];
        }
        x[i] = (b[i] - s) / A[i][i];
    }
    while (!swaps.empty()) {
        std::swap(x[swaps.top().first], x[swaps.top().second]);
        swaps.pop();
    }
}

Добавлено через 3 минуты
СЛАУ размерностью 5184 решалась 228 секунд... ну да ладно, сойдёт для сельской местности...

@vet 182 / 183 / 55 Регистрация: 08.04.2009 Сообщений: 1,309
	16.03.2010, 05:15	2
	Сообщение от taras atavin for (k=n-1; k>=0; --k) // Меняем её с i-той { a[k+tn]&=a[k+in]; a[k+in]&=a[k+tn]; a[k+tn]&=a[k+in]; } taras atavin, а каким компилятором пользуешься? У меня в этом блоке и Dev-C++ и VS 2008 пишут одно и тоже: что нельзя использовать тип double с операцией поразрядного присваивания &= ...? 1

@vet 182 / 183 / 55 Регистрация: 08.04.2009 Сообщений: 1,309
	16.03.2010, 05:58	4
	Сообщение от taras atavin Для инта почему-то работало, вот я и заменил. Все верно поразрядные либо с int либо,что чаще, с unsigned int используются 0

sergeyq
	26.03.2010, 18:31	6
	код обратной матрицы бы

@easybudda Модератор 11655 / 7167 / 1702 Регистрация: 25.07.2009 Сообщений: 13,132
	02.10.2010, 22:23	10
	Сообщение от taras atavin double abs(double x) а чем стандартная fabs() не угодила? 0

@Крюгер 40 / 60 / 8 Регистрация: 16.11.2012 Сообщений: 459 Записей в блоге: 7
	10.11.2013, 09:38	13
	Сообщение от vet А что такое <> ? в С/С++ нет такого условного оператора а разве непонятно что в коде это значит что там число отличное от нуля должно быть 0