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

/*
* Copyright (c) 2009 湖南师范大学数计院 一心飞翔项目组
* All Right Reserved
*
* 文件名：matrix.cpp 定义Point、Node、Matrix类的各个方法
* 摘 要：定义矩阵类，包括矩阵的相关信息和方法
*
* 作 者：刘 庆
* 修改日期：2009年7月19日21:15:12
*
*/

/* 该类方法均基于稀疏线性系统的MCRF存储结构，MCRF相关介绍参见：http://blog.csdn.net/lq_0402/archive/2010/11/11/6003822.aspx */

/* 两个矩阵相乘 */
Matrix& Matrix::operator *(const Matrix &genMatrix)
{
if(&genMatrix==0)
return *this;
/* 如果两矩阵不能相乘，则返回*this */
if(total_col!=genMatrix.total_ln)
{
cout<<"两矩阵不能相乘"<<endl;
return *this;
}
/* 创建一个新Matrix对象matrix (ab*bc==>ad) */
total_ln = total_ln;
total_col = genMatrix.total_col;
total_elem = total_ln*total_col;
Node tempData;
tempData.GetSpace(total_elem+total_ln);
long i=0;
for( i=0; i<total_elem+total_ln; i++)
{
tempData[i].index = -1;
tempData[i].value = 0;
}

/* this 和 genMatrix 相乘，结果存储在matrix中 */
/* 最里层循环就是this的第i行的所有元素与genMatrix的所有列(j)的元素一一相乘累加 */
long offset = total_ln+1; /* tempData存储元素的偏移位置 */
Double sum = 0; /* 临时变量，存储累加的结果 */
long *num = new long[total_ln]; /* 每行非零非对角元素个数 */
for(i=0; i<total_ln; i++) /* num[]应初始化为0，后面用到它自加 */
num[i] = 0;
// long DiagonalZero = total_ln; /* 在目前对角线上的有total_ln个元素 */

for(i=0; i<total_ln; i++)
{ /* 控制this的行 */
for(long j=0; j<genMatrix.total_col; j++)
{ /* 控制genMatrix的列 */
/* 控制this(genMatrix)行的特定行(列)的元素 */
/* this矩阵第i行非零非对角元素与genMatrix第j列的对应元素相乘累加 */
sum = 0;
long pos = -1;
for(long k=data[i].index; k<data[i+1].index; k++)
{
/* this当前元素的坐标为(i,data[k].index)，所以genMatrix元素的坐标为(data[k].index,j) */
pos = genMatrix.QuickLocate(data[k].index, j);
if(pos!=-1)
{ /* pos值为-1，则(data[k].index,j)的元素为0，否则为序列号的值 */
sum += data[k].value*genMatrix.data[pos].value;
}
}
/* this矩阵的Aii与genMatrix矩阵的第j列的第i个元素相乘并累加*/
if(data[i].value!=0.0f)
{
pos = genMatrix.QuickLocate(i,j);
if(pos!=-1)
sum += data[i].value*genMatrix.data[pos].value;
}
/* i==j，则说明得到的matrix矩阵的元素为对角元素
i!=j&&sum!=0，则说明得到的matrix矩阵的元素不为对角元素，并且最终结果不为0，为0不存储 */
if((sum!=0.0f)&&i!=j)
{
tempData[offset].value = sum;
tempData[offset].index = j;
num[i]++; /* matrix第i行的非零非对角元素加1 */
offset++; /* 偏移位置也加1 */
}else if(i==j)
{ /* i、j相等则为对角元素 */
tempData[i].value = sum;
}
}
}
/* 第一行的第一个非零非对角元素值为 总行数加1 */
tempData[0].index = total_ln+1;
for(i=1; i<total_ln+1; i++) /* 公式 */
tempData[i].index = tempData[i-1].index + num[i-1];

/* 将tempData中的数据复制到data中 */

data.GetSpace(offset);
total_elem = offset-1;
for(i=0; i<=total_elem; i++)
{
data[i].value = tempData[i].value;
data[i].index = tempData[i].index;
}
if(num!=NULL)
{
delete[] num;
num = NULL;
}
return *this;
}

/* 矩阵的数乘 */
Matrix& Matrix::operator *(const Double& gen)
{
long offset = total_ln+1;
long *num = new long[total_ln];
long i=0;
for(i=0; i<total_ln; i++)
num[i] = 0;

for(i=0; i<total_ln; i++)
{
data[i].value *= gen;
if(data[i].value==0)
data[i].value = 0;
}
for(i=0; i<total_ln; i++)
{
for(long j=data[i].index; j<data[i+1].index; j++)
{
data[j].value *= gen;
if(data[j].value!=0.0f)
{
data[offset].value = data[j].value;
data[offset].index = data[j].index;
num[i]++;
offset++;
}
}
}
for(i=1; i<total_ln+1; i++)
data[i].index = data[i-1].index + num[i-1];
total_elem = offset-1;
delete[] num;
num = NULL;

return *this;
}
Matrix& Matrix::operator /(const Double& gen)
{
if(gen==0.0f)
{
cout<<"除数不能为0"<<endl;
return *this;
}else
{
double f = 1/gen.GetData();
Double gen2(f);
return (*this)*gen2;
}
}
/* 删除ln,col位置的元素
如果ln,col位置为0，则直接返回
如果ln,col位置存在，假设为对角元素则将其置为0，否则删除该值
*/
Matrix& Matrix::DeleteElement(const long ln, const long col)
{
long pos = (*this).QuickLocate(ln, col);
/* 如果ln,col位置的元素为0，则直接返回 */
if(pos==-1)
return *this;

/* 如果ln,col位置的元素存在，且为对角元素，则将其置为0 */
if(ln==col)
{
data[ln].value = 0;
return *this;
}
/* 程序进入以下部分说明，ln,col位置存在且不为对角元素 */
long i=0;
// Display("this");
/* 将ln行之后的每行第一个非零非对角元素的位置向前挪一个位置 */
for(i = ln+1; i<total_ln+1; i++)
{
data[i].index--;
}
/* 将pos后的元素向前挪一个位置，以达到删除该元素的目的 */
for(i = pos; i<total_elem; i++)
{
data[i].value = data[i+1].value;
data[i].index = data[i+1].index;
}
total_elem--; /* 总非零元素个数减一 */
return *this;
}

/* 增加ln和col位置的元素 等价与PosSetValue(ln, col, value);
如果ln,col位置为0，且value大于门槛值，则将其添进data中
如果ln,col位置不为0，假设value大于门槛值，则将其改为vlaue，否则删除该位置元素
*/
Matrix& Matrix::AddElement(const long ln, const long col, const Double& value)
{
PosSetValue(ln, col, value);
return *this;
}
/* 第ln,col位置的元素乘以gen
假如相乘之后该数小于 门槛值，则得将其视为0元素，在data中删除该元素
*/
Matrix& Matrix::PosMultiply(const long ln, const long col, const Double& gen)
{
long pos = this->QuickLocate(ln, col);
/* 该位置元素为0，则返回 */
if(pos==-1)
return *this;

/* 该位置元素不为0，则乘以gen */
data[pos].value *= gen;
/* 如果该位置对角元素，则返回，否则判断其是否大于门槛值 */
if(pos<total_ln)
return *this;

/* 如果该位置元素大于门槛值，则返回，否则删除该元素 */
if(data[pos].value!=0.0f)
return *this;

/* 删除该元素 */
DeleteElement(ln, col);
return *this;
}
/* 替ln,col位置赋值
如果data中无数据
如果ln==col，则 data.length = ln + 1;
如果ln!=col,则data.legnth = ln + 2;
return
如果ln,col位置的元素存在,则将其赋值为value,如果value==0.且ln!=col,则删除该元素,否则赋值为0
如果ln,col位置的元素不存在
假如value==0,
假如col>=total_col,total_col = col+1;
假如ln>=total_ln,重新开辟空间将data中的值赋值给tempData中,然后用tempData覆盖data
(tempData.length = total_elem+ln+1-total_ln)
return
假如value!=0,则将其添在适当的位置.
假如col>=total_col,total_col = col+1;
假如ln>=total_ln,重新开辟空间tempData,将data中的值复制到temData中,同时添进value,覆盖data
(tempData.length = total_elem + (ln+1-total_ln) + 1.
return

(ln,col必定不为对角元素)重新开辟空间tempData,将data中的值复制到data中,以及添进value值,覆盖data
(tempData.length = total_elem + 1;
return
*/
Matrix& Matrix::PosSetValue(const long ln, const long col, const Double& value)
{
if(data.point==0)
{
if(ln==col)
{
total_ln = ln + 1;
total_col = col + 1;
total_elem = total_ln;
data.GetSpace(total_elem+1);
long i = 0;
for( ; i<total_ln-1; i++)
{
data[i].value = 0;
data[i].index = total_ln+1;
}
if(value!=0.0f) // 假如value小于门槛值，则视为0
data[i].value = value;
else
data[i].value = 0;
data[i].index = total_ln+1;
i++;
data[i].value = 0;
data[i].index = total_ln+1;

return *this;
}
else
{
total_ln = ln + 1;
total_col = col + 1;
total_elem = total_ln + 1;
data.GetSpace(total_elem+1);
long i=0;
for( ; i<total_ln; i++)
{
data[i].value = 0;
data[i].index = total_ln+1;
}
data[i].value = 0;
data[i].index = total_ln+2;
i++;
if(value!=0.0f)
data[i].value = value;
else
data[i].value = 0;
data[i].index = col;

return *this;
}
}
else
{ // data.point != 0
// 如果ln,col位置的元素存在,则将其赋值为value,如果value==0.且ln!=col,则删除该元素,否则赋值为0
long pos = (*this).QuickLocate(ln,col);
if(pos!=-1)
{
if(ln==col)
{
if(value!=0.0f)
data[pos].value = value;
else
data[pos].value = 0;
}
else
{
if(value!=0.0f)
data[pos].value = value;
else
DeleteElement(ln,col);
}
return *this;
}
else
{ // ln,col元素不存在
/*假如value==0,
假如col>=total_col,total_col = col+1;
假如ln>=total_ln,重新开辟空间将data中的值赋值给tempData中,然后用tempData覆盖data
(tempData.length = total_elem+ln+1-total_ln)
return
*/
if(value==0)
{
if(col>=total_col)
total_col = col+1;
if(ln>=total_ln)
{ // 包含ln==col，将其赋值为0
if(data.Length()>=total_elem+1+ln+1-total_ln)
{
long i = total_elem;
// 将total_elem~total_ln+1元素后移(ln+1-total_ln)位
for( ; i>total_ln; i--)
{
data[i+(ln+1-total_ln)].value = data[i].value;
data[i+(ln+1-total_ln)].index = data[i].index;
}
// 添加新的空闲位置元素，i==total_ln
data[i+(ln+1-total_ln)].value = 0;
data[i+(ln+1-total_ln)].index = data[i].index + (ln+1-total_ln);
i--;
// 添加新添加的行对角元素
long j = 0;
for( ; j<(ln+1-total_ln); j++)
{
data[i+(ln+1-total_ln)-j].value = 0;
data[i+(ln+1-total_ln)-j].index = data[total_ln].index+(ln+1-total_ln);
}
// 修改data中total_ln~0的index值
for( ; i>=0; i--)
{
data[i].index += (ln+1-total_ln);
}
}
else
{
Node tempData;
tempData.GetSpace(total_elem+1+ln+1-total_ln+total_ln); // 多开辟total_ln个空间
long i=0;
for( ; i<total_ln; i++)
{ // 0~total_ln-1 对角元素
tempData[i].value = data[i].value;
tempData[i].index = data[i].index+ln+1-total_ln; // #####
}
for( ; i<ln+1; i++)
{ // 包括total_ln~ln对角元素 total_ln行之后无非零非对元素
tempData[i].value = 0;
tempData[i].index = data[total_ln].index+(ln+1-total_ln);
}
tempData[i].value = 0; // 空闲位置
tempData[i].index = data[total_ln].index+(ln+1-total_ln);
i++;
for( ; i<total_elem+1+ln+1-total_ln; i++ )
{ // 复制后面的非零非对角元素
tempData[i].value = data[i-ln-1+total_ln].value;
tempData[i].index = data[i-ln-1+total_ln].index;
}
data = tempData;
}
total_elem = total_elem + (ln+1-total_ln);
total_ln = ln + 1;
return *this;
}
}
else
{
/*假如value!=0,则将其添在适当的位置.
假如col>=total_col,total_col = col+1;
假如ln>=total_ln,重新开辟空间tempData,将data中的值复制到temData中,同时添进value,覆盖data
[tempData.length = total_elem + (ln+1-total_ln) + 1]
return

重新开辟空间tempData,将data中的值复制到data中,以及添进value值,覆盖data
(tempData.length = total_elem + 1;
return
*/
if(col>=total_col)
total_col = col + 1;
if(ln>=total_ln)
{
if(data.Length()>=total_elem+1+(ln+1-total_ln)+1)
{
if(ln==col)
{ // 添加的元素为超过行标总数的对角元素
long i = total_elem;
// 将total_elem~total_ln+1元素后移(ln+1-total_ln)位
for( ; i>total_ln; i--)
{
data[i+(ln+1-total_ln)].value = data[i].value;
data[i+(ln+1-total_ln)].index = data[i].index;
}
// 添加新的空闲位置元素，i==total_ln
data[i+(ln+1-total_ln)].value = 0;
data[i+(ln+1-total_ln)].index = data[total_ln].index + (ln+1-total_ln);
i--;
// 添加新添加的行对角元素
data[i+(ln+1-total_ln)].value = value;
data[i+(ln+1-total_ln)].index = data[total_ln].index + (ln+1-total_ln);
i--;
long j = 1;
for( ; j<(ln+1-total_ln); j++){
data[i+(ln+1-total_ln)-j].value = 0;
data[i+(ln+1-total_ln)-j].index = data[total_ln].index+(ln+1-total_ln);
}
// 修改data中total_ln~0的index值
for( ; i>=0; i--){
data[i].index += (ln+1-total_ln);
}
total_elem = total_elem + (ln+1-total_ln);
total_ln = ln + 1;
return *this;
}
else
{ // 添加的元素为超过行标总数的非对角元素
long i=total_elem;
// 添加新元素至末尾
data[i+(ln+1-total_ln+1)].value = value;
data[i+(ln+1-total_ln+1)].index = col;
// 将 total_elem~total_ln+1元素后移
for( ; i>total_ln; i--){
data[i+(ln+1-total_ln)].value = data[i].value;
data[i+(ln+1-total_ln)].index = data[i].index;
}
// 添加新的空闲位置元素，i==total_ln
data[i+(ln+1-total_ln)].value = 0;
data[i+(ln+1-total_ln)].index = data[total_ln].index + (ln+1-total_ln)+1;
i--;
// 添加新添加的行对角元素
long j = 0;
for( ; j<(ln+1-total_ln); j++){
data[i+(ln+1-total_ln)-j].value = 0;
data[i+(ln+1-total_ln)-j].index = data[total_ln].index+(ln+1-total_ln);
}
// 修改data中total_ln~0的index值
for( ; i>=0; i--){
data[i].index += (ln+1-total_ln);
}
total_elem = total_elem + (ln+1-total_ln+1);
total_ln = ln + 1;
return *this;
}
}
else
{
Node tempData;
tempData.GetSpace(total_elem+1+(ln+1-total_ln)+1+total_ln); // 多开辟total_ln个空间
long i = 0;
// 0~total_ln-1
for( ; i<total_ln; i++)
{
tempData[i].value = data[i].value;
tempData[i].index = data[i].index + (ln+1-total_ln);
}
// total_ln ~ ln-1
for( ; i<ln; i++)
{
tempData[i].value = 0;
tempData[i].index = data[total_ln].index + (ln+1-total_ln);
}
// ln ~
if(ln==col)
{
// ln
tempData[i].value = value;
tempData[i].index = data[total_ln].index + (ln+1-total_ln);
i++;
// 空闲位置
tempData[i].value = 0;
tempData[i].index = data[total_ln].index + (ln+1-total_ln);
i++;
for( ; i<total_elem+1+(ln+1-total_ln); i++)
{ // =====
tempData[i].value = data[i-(ln+1-total_ln)].value;
tempData[i].index = data[i-(ln+1-total_ln)].index;
}
total_elem = total_elem + (ln+1-total_ln);
total_ln = ln + 1;
data = tempData;
return *this;
}
else
{ // ln != col
// ln
tempData[i].value = 0;
tempData[i].index = data[total_ln].index + (ln+1-total_ln);
i++;
// 空闲位置
tempData[i].value = 0;
tempData[i].index = data[total_ln].index + (ln+1-total_ln)+1;
i++;
for( ; i<total_elem+1+(ln+1-total_ln); i++)
{ // =====
tempData[i].value = data[i-(ln+1-total_ln)].value;
tempData[i].index = data[i-(ln+1-total_ln)].index;
}
tempData[i].value = value;
tempData[i].index = col;
total_elem = total_elem + (ln+1-total_ln) + 1;
total_ln = ln + 1;
data = tempData;

return *this;
}
}
}
else
{ // ln < total_ln 有可能为对角元素，col>total_col&&ln==col情况
if(ln==col)
{
// 为对角元素，则直接修改其值即可
data[ln].value = value;
}
else
{
if(data.Length()>=total_elem+1+1)
{
long i = total_elem;
// 第ln+1行及以后的非零非对角元素后移一位
for( ; i>=data[ln+1].index; i--)
{
data[i+1].value = data[i].value;
data[i+1].index = data[i].index;
}
// 第ln行列标大于col的元素后移一位
for( ; i>=data[ln].index; i--)
{
if(data[i].index>col)
{
data[i+1].value = data[i].value;
data[i+1].index = data[i].index;
}else
break;
}
// 插入新元素
data[i+1].value = value;
data[i+1].index = col;
i = ln+1;
for( ; i<total_ln+1; i++)
{
data[i].index++;
}
total_elem++;

return *this;
}
else
{
Node tempData;
tempData.GetSpace(total_elem+1+1);
long i=0;
for( ; i<ln+1; i++)
{
tempData[i].value = data[i].value;
tempData[i].index = data[i].index;
}
// ln行之后的每行第一个非零非对角元素的序号加1
for( ; i<total_ln+1; i++)
{
tempData[i].value = data[i].value;
tempData[i].index = data[i].index + 1;
}
// 复制ln行之前的元素
for( ; i<data[ln].index; i++)
{
tempData[i].value = data[i].value;
tempData[i].index = data[i].index;
}
// 复制ln行的元素
// 假如ln有元素，则执行for循环，假如没有，则会跳过次for循环，而添加的元素必须添加在这个位置
long tempI = 0;
for( ; i<data[ln+1].index; i++)
{
if(data[i].index<col)
{
tempData[i].value = data[i].value;
tempData[i].index = data[i].index;
}
else
{ //
tempData[i].value = value;
tempData[i].index = col;
tempI = 1;
break;
}
}
if(tempI == 0)
{
tempData[i].value = value;
tempData[i].index = col;
i++;
}
for( ; i<total_elem+1+1; i++)
{
tempData[i].value = data[i-1].value;
tempData[i].index = data[i-1].index;
}
total_elem = total_elem + 1;
data = tempData;

return *this;
}
}
}
}
}
}
return *this;
}