基本的GMRES基于cpp实现
理论部分和matlab版本可以查看我的其他文章:
以下为代码,头文件里面有cuda的东西。无视即可。
int main()
{
//time start
clock_t startTime = clock();
//load
FILE *fp;
static double data[100];
if (fp = fopen("data.bin", "rb"))
{
for (int i = 0; i<100; i++)
{
fread(&data[i], sizeof(double), 100, fp);
}
}
//data
double *A, *x0, *b;
A = data;
x0 = new double[10];
b = new double[10];
for (int i = 0; i < 10; i++)
{
x0[i] = 0;
b[i] = 0;
}
x0[0] = 1;
//GMRES
double *y0, *r0, *beta;
y0 = new double[10];
r0 = new double[10];
for (int i = 0; i < 10; i++)
{
y0[i] = 0;
r0[i] = 0;
}
for (int i = 0; i < 10; i++)
{
for (int j = 0; j < 10; j++)
{
y0[i] += A[i * 10 + j] * x0[j];
}
r0[i] = b[i] - y0[i];
}
beta = new double[1];
beta[0] = 0;
double *v = new double[10];
for (int i = 0; i < 10; i++)
{
v[i] = 0;
}
//V1,单位化
double temp = 0;
for (int j = 0; j < 10; j++)
{
temp = temp + r0[j] * r0[j];
}
beta[0] = sqrt(temp);
for (int j = 0; j < 10; j++)
{
if (beta[0] == 0)
{
v[j] = 0;
}
else
{
v[j] = r0[j] / beta[0];
}
}
//data
double *omiga, *V, *H;
omiga = new double[10];//vector
for (int j = 0; j < 10; j++)
{
omiga[j] = 0;
}
V = new double[100];//matrix
for (int j = 0; j < 100; j++)
{
V[j] = 0;
}
H = new double[110];//matrix
for (int j = 0; j < 110; j++)
{
H[j] = 0;
}
//V(:,1) = v1
for (int i = 0; i < 10; i++)
{
V[10 * i] = v[i];
}
//Do
double resVec[10],resVal;
for (int j = 0; j < 10; j++)
{
//Compute omiga(j) := A*v(j)
for (int i = 0; i < 10; i++)
{
for (int k = 0; k < 10; k++)
{
omiga[i] += A[i * 10 + k] * V[k * 10 + j];
}
}
//h(i,j) and res
for (int i = 0; i <= j; i++)
{
//h(i,j)
for (int k = 0; k < 10; k++)
{
H[i * 10 + j] += omiga[k] * V[10 * k + i];
}
//omiga(i,j)
for (int k = 0; k < 10; k++)
{
omiga[k] = omiga[k] - H[i * 10 + j] * V[10 * k + i];
}
}
//h(j+1,j)
for (int i = 0; i < 10; i++)
{
H[10 * (j + 1) + j] += omiga[i] * omiga[i];
}
H[10 * (j + 1) + j] = sqrt(H[10 * (j + 1) + j]);
//if
if (H[10 * (j + 1) + j] < 1e-11)//10-10
{
break;
printf("done , residual is %d\n", H[10 * (j + 1) + j]);
}
v(j+1)
//for (int i = 0; i <10; i++)
//{
// v[i] = omiga[i] / H[10 * (j + 1) + j];
//}
//V(:,j+1)
for (int i = 0; i <10; i++)
{
if (j < 9)
{
V[10 * i + (j + 1)] = omiga[i] / H[10 * (j + 1) + j];
}
else
{
resVec[i] = omiga[i];
}
}
}
//visualizion
//data
printf("data:\n");
for (int i = 0; i < 10; i++)
{
for (int j = 0; j < 10; j++)
{
cout << data[10*i+j] << " ";
}
cout << endl;
}
printf("\n");
//A
printf("A:\n");
for (int i = 0; i < 10; i++)
{
for (int j = 0; j < 10; j++)
{
cout << A[10 * i + j] << " ";
}
cout << endl;
}
printf("\n");
//x0
printf("x0:\n");
for (int i = 0; i < 10; i++)
{
cout << x0[i] << " ";
}
printf("\n");
//b
printf("b:\n");
for (int i = 0; i < 10; i++)
{
cout << b[i] << " ";
}
printf("\n");
//test
printf("v:\n");
for (int i = 0; i < 10; i++)
{
cout << v[i] << " ";
}
printf("\n");
//V
printf("V:\n");
for (int i = 0; i < 10; i++)
{
for (int j = 0; j < 10; j++)
{
cout << V[10 * i + j] << " ";
}
cout << endl;
cout << "next row" << endl;
}
printf("\n");
//H
printf("H:\n");
for (int i = 0; i < 11; i++)
{
for (int j = 0; j < 10; j++)
{
cout << H[10 * i + j] << " ";
}
cout << endl;
cout << "next row" << endl;
}
printf("\n");
//res vector
printf("res vector:\n");
for (int i = 0; i < 10; i++)
{
cout << resVec[i] << " ";
}
printf("\n");
//res Val
resVal = 0;
for (int i = 0; i < 10; i++)
{
resVal += resVec[i] * resVec[i];
}
cout << "残差值为"<<sqrt(resVal) <<endl;;
//window
clock_t endTime = clock();
cout << "整个程序用时:" << double(endTime - startTime) / CLOCKS_PER_SEC << "s" << endl;
system("PAUSE");
return 0;
}
其中data是10*10的一个矩阵,然后new出来的空间都没有delete掉,用的时候可以注意一下