Метод сопряженных градиентов

[nemoart]

Доброго времени суток, не могу разобраться в алгоритмом сопряженных градиентов. Не могу понять что делаю не так...
Для решения системы
$ \left\{ \begin{array}{l} 3x_0-x_1 = 3\\ -x_0+3x_1 = 7 \end{array} \right. $
При первом шаге получаю что приближенное значение $x=(1,318, 3,076)$... а потом значение начинает убывать....
Вот код на С++:

#include <conio.h>
#include <stdio.h>
#include <tchar.h>
#include <math.h>
#define n 2

void gaxpy (float A[n][n],float x[n], float b[n])
{
   for (int i=0; i<n; i++)
      for (int j=0; j<n; j++)
      {
         b[i]+=A[i][j]*x[j];
      }
}

float dot (float x[n], float y[n])
{
   float s;
   for (int i=0; i<n; i++)
      s+=x[i]*y[i];
   return s;
}

void saxpy (float alf,float x[n],float y[n],float z[n])
{
   for (int i=0; i<n; i++)
      z[i]=x[i]+alf*y[i];
}

void print_vec (float b[n])
{
   for (int i=0; i<n; i++)
   {
      printf(" %3.1f",b[i]);
   }
}

void main(void)
{
  float A[n][n];
  float b[n], x[n], r[n], p[n], v[n];
  b[0]=3; b[1]=7;
  A[0][0]=3; A[0][1]=-1; A[1][0]=-1; A[1][1]=3;
  for (int i=0; i<n; i++) x[i]=0;
  float alfa,alfa_,la,al;
  float eps=0.001;
  for (int i=0; i<n; i++) p[i]=x[i];
  goxpy(A,p,v);
  saxpy(-1,b,v,p);
  for (int i=0; i<n; i++) r[i]=p[i];
  alfa=dot(r,r);
  int k;
  k=0;
  while (k<=1000)
  {
     gaxpy(A,p,v);
     la=alfa/dot(p,v);
     saxpy(la,x,p,x);
     saxpy(-la,r,v,r);
     alfa_=dot(r,r);
     al=alfa_/alfa;
     if (dot(r,r)<eps) break;
     saxpy(al,r,p,p);
     alfa=alfa_;
     k++;
  }
  print_vec(x);
  getch();
}

Помогите разобраться...

[nemoart]

Ну или хотя бы подскажите достойную литературу с алгоритмом, уже готовым для программирования... без теоретических выкладок. Буду очень благодарен.

[petywilo]

Ищи книжку Yousef Saad "Iterative methods for sparse linear systems" - имхо лучшая книга по решению СЛАУ. Алгоритмы есть на все общеупотребительный методы.

[nestoklon]

http://www.nr.com/
Думаю, в интернетах можно найти в украденном виде. Я находил =) .
Там про сопряжённые градиенты чуть ли не глава целая была.
UPD А, я не украденную находил. Есть старая версия. Которой с точки зрения описания алгоритмов более чем достаточно.

[nemoart]

petywilo в сообщении #308774 писал(а):

Ищи книжку Yousef Saad "Iterative methods for sparse linear systems" - имхо лучшая книга по решению СЛАУ. Алгоритмы есть на все общеупотребительный методы.

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

Используется синтаксис C

 float dot (float x[n], float y[n])
{
   float s;
   s=0; //<-----!!!!!!!!!!!!!!!!
   for (int i=0; i<n; i++)
      s+=x[i]*y[i];
   return s;
}
 

Сумму то надо занулять :) А вроде третий курс