Метод наискорейшего градиентного спуска на C

Как ни странно, тоже оооочень простой метод.



Теория
Дано:
Ax=b, где A- симметричная положительно определённая матрица
Найти:
X(x1...xm)
Алгоритм:
Прост и очевиден. Всего у нас имеется две формулы:
Xn+1=Xn-αn(AXn-b);
αn=(AXn-b, AXn-b) / (A*(AXn -b), AXn-b); //запятая - скалярное умножение векторов
Всё что нам нужно делать - находить альфа и Xn+1 до тех пор, пока сумма разностей векторов Xn-1 и Xn(Векторная норма) > Eps.

Практика
Да кому эта теория нужна! Намного интереснее писать программки!
Функции, которые нам понадобятся:
void matrix_m_vector(float **m, float *v, float *res)
{
for (int i=0; i<_n; i++)
{
  res[i] = 0;
  for (int j=0; j<_n; j++){
    res[i]+=m[i][j]*v[j];
  }
}
}

void vector_m_vector(float *v1, float *v2, float *res)
{
for (int i=0; i<_n; i++)
{
        *res+=v1[i]*v2[i];
}
}

void vector_minus_vector(float *v1, float *v2, float *res)
{
for (int i=0; i<_n; i++)
{
    res[i]=v1[i]-v2[i];
}
}
Названия у них вполне очевидные. Работают не совсем очевидно. Во-первых удобнее было бы юзать return а не строку параметров, а если совсем красиво - юзать ООП. Ну это опять же...По мере наличия желания. Главное - рабочая реализация алгоритма.
Сначала нужно выбрать начальные значения вектора X. Ну что тут думать-то? 0.1 и всего делов...

while(1)
{
cout << c << endl;
matrix_m_vector(ma, X, result);
vector_minus_vector(result, vector, result2); //result2 хранит результат AXn-B
for (int i=0; i<_n; i++) cout << "X: " <<X[i]<< " ";
matrix_m_vector(A, result2, result3);
float num, den;
vector_m_vector(result2, result2, &num);
vector_m_vector(result3, result2, &den);
alpha=num/den;  //альфа найдена по формуле
cout << "\n alpha: " << c << alpha <<endl;
//считаем следующий X, опять же по формуле...
for (int j=0; j<_n; j++) result2[j]*=alpha;
vector_minus_vector(X, result2, Xn1);
float eps=0;
for (int k=0; k<_n; k++)
{
    eps+=fabs(Xn1[k]-X[k]); //критерий окончания
}
if (eps<0.0001) break;
memcpy(X, Xn1, sizeof(Xn1)*_n);
c++;
}


Ну вот собственно и всё. Вообще если грамотно переписать этот цикл - кода получится намнооого меньше (:. Кому нечем заняться - вперёд, постим оптимизации в комментах!