【刘庆源码共享】稀疏线性系统求解算法MGMRES(m) 之 MGMRES类定义（C++）

/*
* Copyright (c) 2009 湖南师范大学数计院 一心飞翔项目组
* All Right Reserved
*
* 文件名：mgmres.cpp
* 摘 要：MGMRES类的函数定义
*
* 作 者：刘 庆
* 完成日期：2009年4月28日
*
*/
#include "stdafx.h"

#include <math.h>
#include <iostream>
#include <time.h>
#include "mgmres.h"
#include "EvalFile.h"

using namespace std;

MGMRES::MGMRES()
{
max_outer_itr = 1;
resultOk = 0;
}
MGMRES::~MGMRES()
{
if(resultOk!=NULL)
delete[] resultOk;
resultOk = NULL;
}

/* 替算法所涉及的矩阵赋初始值，所有数据在fileName的文件中 */
void MGMRES::SetMatrixsValue(const char* fileName)
{
EvalFile file = EvalFile(fileName);
file.SetWithSub(A); /* 替系数矩阵赋值 */
file.SetWithoutSub(b); /* 替右端向量赋值 n*1 */

if(A.GetTotal_col()<A.GetTotal_ln())
{
A.SetTotal_col(A.GetTotal_ln());
b.PosSetValue(A.GetTotal_col()-1, 0, 0);
}
if(A.GetTotal_col()>A.GetTotal_ln())
{
A.PosSetValue(A.GetTotal_col()-1,0,0);
b.PosSetValue(A.GetTotal_col()-1,0,0);
}
x = Matrix::GetE(A.GetTotal_col()); /* 将x的各个值初始化为0 */

// A.DisplayAsMatrix("A", 1);
}

/* A*B、A+B运算中A都会发生变化，B不会 */
Matrix& MGMRES::MGMRESCode()
{
Double mu=0, residue=0, h=0;
Double *c = new Double[A.GetTotal_ln()+1], *s = new Double[A.GetTotal_ln()+1], *g = new Double[A.GetTotal_ln()+1];
max_inner_itr = A.GetTotal_ln()>A.GetTotal_col() ? A.GetTotal_ln():A.GetTotal_col();
long i=0, j=0, k=0, outer_itr=0;

do
{
Matrix r, v, w, y, A_backup, vv, Vm, H;
A_backup = A;
A_backup = A_backup*x;
r = (A_backup-b)*(-1); /* 计算 r = b - Ax */

residue = r.GetModulus(); /* residue = |r| */

v = r/residue; /* v0 = r */
w = v;

for(i=0; i<A.GetTotal_ln()+1; i++)
{
g[i] = 0;
c[i] = 0;
s[i] = 0;
}
g[0] = residue;
long k_backup = 0;
for(k=0; k<max_inner_itr; k++)
{
k_backup = k;
Vm.AddCol(w);

for(j=0; j<k+1; j++)
{
v = Vm.GetColVector(j);
A_backup = A;
w = A_backup*v;
w = w/w.GetModulus();

A_backup = A;
vv = A_backup*w;

for (i = 0; i <= j; i++)
{
v = Vm.GetColVector(i);
h = w.GetDotMetrix(v);
H.AddElement(i,j, h);

vv = vv - v*h;
}
h = vv.GetModulus();
v = vv/h;
H.AddElement(j+1, j, h);
}

/* for (j = k+1; j < max_inner_itr; j++)
{
A_backup = A;
vv = A_backup*w;
for (i = 0; i <= k; i++)
{
v = Vm.GetColVector(i);
h = w.GetDotMetrix(v);
H.AddElement(i,j, h);

vv = vv - v*h;
}
h = vv.GetModulus();
v = vv/h;
H.AddElement(j+1, j, h);
}
*/
if(k>0)
{
for(j=0; j<k; j++){
Double yk = 0, yk1 = 0;
long pos = H.QuickLocate(j,k);
if(pos!=-1)
yk = H.data[pos].value; // y[k]
pos = H.QuickLocate(j+1, k);
if(pos!=-1)
yk1 = H.data[pos].value; // y[k+1]

Double g1 = c[j]*yk-s[j]*yk1; // hjk = c[j]*y[k]-s[j]*y[k+1]
Double g2 = s[j]*yk+c[j]*yk1; // hj+1,k = s[j]*y[k]-c[j]*y[k+1]
H.PosSetValue(j,k,g1);
H.PosSetValue(j+1,k,g2);
}
}

Double hkk = 0, hk1k = 0;
long posk = H.QuickLocate(k,k);
long posk1 = H.QuickLocate(k+1,k);
if(posk!=-1)
hkk = H.data[posk].value;
if(posk1!=-1)
hk1k = H.data[posk1].value;
mu = sqrt(hkk*hkk+hk1k*hk1k);
c[k] = hkk/mu;
s[k] = hk1k/mu*(-1);
hkk = mu;
H.PosSetValue(k,k,hkk);
H.PosSetValue(k+1, k, 0);

Double g1 = c[k]*g[k]-s[k]*g[k+1];
Double g2 = s[k]*g[k]+c[k]*g[k+1];
g[k] = g1;
g[k+1] = g2;

residue = fabs(g[k+1].GetData());
cout<<"第 "<<outer_itr+1<<" 次外循环/t第 "<<k+1<<" 次内循环/t误差为"<<residue<<endl;
if(residue==0)
break;
}

cout<<endl;
k = k_backup;
h = H.data[k].value;
Double yk;
if(h!=0)
yk = g[k]/h;
else
yk = 0;
y.PosSetValue(k,0,yk);

for(i=k-1; i>=0; i--)
{
yk = g[i];
for(j=i+1; j<k+1; j++)
{
long pos = H.QuickLocate(i,j);
if(pos!=-1)
h = H.data[pos].value;
Double yj = 0;
pos = y.QuickLocate(j,0);
if(pos!=-1)
yj = y.data[pos].value;
yk = yk - h*yj;
}
yk = yk/H.data[i].value;
y.PosSetValue(i,0,yk);
}
Matrix V_backup = Vm;
x = x + V_backup*y;

A_backup = A;
A_backup = (A_backup*x-b)*(-1);
residue = A_backup.GetModulus();

cout<<"第 "<<outer_itr+1<<" 次循环/t误差为"<<residue<<endl;

outer_itr++;
if(outer_itr>=max_outer_itr)
break;
}while( residue!=0 );

delete[] c;
delete[] s;
delete[] g;
c = NULL;
s = NULL;
g = NULL;

return x;
}

Matrix MGMRES::ValidateResult()
{
Matrix x_backup = x;
Matrix r = b-A*x_backup;
resultOk = new int[r.GetTotal_ln()];
for(int i=0; i<r.GetTotal_ln(); i++)
resultOk[i] = 0;

long sum=0;
for(int i=0; i<r.GetTotal_ln(); i++)
{
long pos=r.QuickLocate(i,0);
if(pos==-1)
{
resultOk[i] = 1;
}
else
{
if(r.data[pos].value==0)
resultOk[i] = 1;
}
}
cout<<"该组解已经满足以下方程："<<endl;
for(int i=0; i<r.GetTotal_ln(); i++)
{
if(resultOk[i]==1){
cout<<i+1<<" ";
if(++sum%10==0)
cout<<endl;
}
}
cout<<endl;
cout<<"该组解还不满足以下方程："<<endl;
for(int i=0; i<r.GetTotal_ln(); i++)
{
if(resultOk[i]==0){
cout<<i+1<<" ";
if(++sum%10==0)
cout<<endl;
}
}
cout<<endl;
return r;
}

int MGMRES::IsEMRequalMR() const
{
return A.IsEMRequalMR(b);
}

void MGMRES::PrintRunTime(clock_t start, clock_t finsh, char* str) const
{
double totalTime = (double)(finsh-start)/CLOCKS_PER_SEC;
cout<<str<<totalTime<<"s"<<endl;
}