#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#define N 4
#define epsilon ?0.00005


/*
矩陣A和向量b，相乘結(jié)果存在Ab里
*/
void matrixTimesVec(double A[N][N], double b[N], double Ab[N])
{
int i, j;


for (i = 0; i < N; i++)
{
Ab[i] = 0.0;
for (j = 0; j < N; j++)
{
Ab[i] = Ab[i] + A[i][j] * b[j];
}
}
}
/*
返回兩個向量的內(nèi)積
*/
double scalarProduct(double vec1[], double vec2[])
{
double s = 0;
int i;
for (i = 0; i < N; i++)
{
s = s + vec1[i] * vec2[i];
}
return s;
}
/*
向量vec1，vec2之和存在vec里
*/
void vecPlus(double vec1[], double vec2[], double vec[])
{
int i;
for (i = 0; i < N; i++)
{
vec[i] = vec1[i] + vec2[i];
}
}
/*
vec0向量每個元素乘以系數(shù)num，保存在vec中
*/
void numPlusVec(double num, double vec0[], double vec[])
{
int i;
for (i = 0; i < N; i++)
vec[i] = num * vec0[i];

}

int main()
{
//原問題：Ax=b
//求解二次泛函：fi(x)=1/2x^TAx-b^Tx
int i, j;
/*
初始化
*/
static double A[N][N] = { -16,0 ,0,0,
0,26,0,0,
0,0,13,0,
0,0,0,13 };//A(要求對稱 + 正定)
static double b[N] = { 2,21,34,10 };//b
static double x0[N] = { 0,0, 0, 0 };//初始解x0
double x[N], r[N], p[N], w[N], alpha, rho00, rho0, rho1, beta;
//打印
printf("\n要求解的示例方程組為：\n A ||| b \n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
printf("%f ", A[i][j]);
}
printf("||| %f\n", b[i]);
}
printf("初始解x0為：\n");
for (int i = 0; i < N; i++) {
printf("%f ", x0[i]);