【刘庆源码共享】稀疏线性系统求解算法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 */

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <String>
#include <math.h>
#include <assert.h>
#include "matrix.h"
#include "MyException.h"

using namespace std;

/* Node类的 方法定义 */
Node::Node()
{
point = 0;
length = 0;
}
/* 拷贝构造函数，新创建一个与参数完全一样的Node对象 */
Node::Node(const Node& node)
{
if (point != NULL)
{
Delete();
}
GetSpace(node.length);
for(long i=0; i<node.length; i++)
{
point[i].value = node.point[i].value;
point[i].index = node.point[i].index;
}
length = node.length;
}
/* 析构函数，释放所占空间 */
Node::~Node()
{
Delete();
}
/* 释放point指针所占空间 */
void Node::Delete()
{
if(point!=NULL)
{
TOTAL_SPACE -= Length()*sizeof(Point);
delete[] point;
point = NULL;
length = 0;
}
}
/* 重载[]，判断参数指定的下标是否越界，越界则抛出异常 */
inline Point& Node::operator[](const long index) const
{
if(length<=0)
{
cout<<"还未申请空间！"<<endl;
throw NullPointException("Null Point Exception");
}
else if(index>=length||index<0)
{
cout<<"下标越界啦!"<<" index="<<index<<"不属于[0, "<<length<<")"<<endl;
throw OutOfRangeException("Out Of Range Exception!");
}
return point[index];
}
/* 重载=，返回一个与参数完全一样的新Node对象 */
inline Node& Node::operator=(const Node& node){
if(this==&node)
return *this;

if (point != NULL)
{
Delete();
}
GetSpace(node.length);
for(long i=0; i<node.length; i++)
{
point[i].value = node.point[i].value;
point[i].index = node.point[i].index;
}
length = node.length;
return *this;
}
/* 申请空间，替point指针申请空间 */
void Node::GetSpace(const long len)
{
if (point != NULL)
{
Delete();
}
point = new Point[len];

assert(point != NULL); // 判断内存分配是否成功

length = len;
// 累加 申请的空间
TOTAL_SPACE += Length()*sizeof(Point);
if(TOTAL_SPACE>MAX_TOTAL_SPACE)
MAX_TOTAL_SPACE = TOTAL_SPACE;
}
/* 返回point指针的长度 */
long Node::Length() const
{
return length;
}

/* Matrix 类的方法定义 */
Matrix::Matrix()
{
InitZero();
}
/* 生成一个与matrix一样的Matrix */
Matrix::Matrix(const Matrix &matrix)
{
InitZero();
if (matrix.data.Length() != NULL)
{
data.GetSpace(matrix.total_elem+1);
for(long i=0; i<matrix.total_elem+1; i++)
{
data[i].value = matrix.data[i].value;
data[i].index = matrix.data[i].index;
}
total_ln = matrix.total_ln;
total_col = matrix.total_col;
total_elem = matrix.total_elem;
}
}
/* 析构函数 清除所有空间 */
Matrix::~Matrix()
{
total_elem = 0;
total_ln = 0;
total_col = 0;
}
/* 将矩阵信息置为0或者NULL */
void Matrix::InitZero()
{
total_elem = 0;
total_ln = 0;
total_col = 0;
}
/* 下面是total_elem、total_ln、total_col的set/get方法 */
void Matrix::SetTotal_elem(long totalElem)
{
total_elem = totalElem;
}
long Matrix::GetTotal_elem()
{
return total_elem;
}
void Matrix::SetTotal_ln(long totalLn)
{
total_ln = totalLn;
}
long Matrix::GetTotal_ln()
{
return total_ln;
}
void Matrix::SetTotal_col(long totalCol)
{
total_col = totalCol;
}
long Matrix::GetTotal_col()
{
return total_col;
}

/* 显示矩阵信息 参数为提示信息 */
void Matrix::Display(char *str) const
{
if(data.point==0){
cout<<"矩阵没有任何信息！"<<endl;
return;
}
if(str!=0)
cout<<"/t"<<str<<endl;
cout<<"total_elem: "<<total_elem<<endl;
cout<<"toatl_ln: "<<total_ln<<endl;
cout<<"total_col: "<<total_col<<endl;
cout<<"Quote:";
long i = 0;
for( i=0; i<total_elem+1; i++)
{
cout.width(6);
cout<<right<<i;
}
cout<<"/nValue:";
for(i=0; i<total_elem+1; i++)
{
cout.precision(2);
cout.width(6);
cout<<right<<fixed<<data[i].value;
}
cout<<"/nIndex:";
for(i=0; i<total_elem+1; i++)
{
cout.width(6);
cout<<right<<data[i].index;
}
cout<<endl;
}
/* 以矩阵的形式输出矩阵 zero == 1 输出0，否则不输出 */
void Matrix::DisplayAsMatrix(char *str, int zero) const
{
if(data.point==0)
{
cout<<"矩阵没有任何信息！"<<endl;
return;
}
if(str!=0)
cout<<"/t"<<str<<endl;

cout<<"Quote:";
long i = 0;
for( ; i<total_col; i++)
{
cout<<" ";
cout.width(2);
cout<<i;
cout<<" ";
}
cout<<"/n/n";
for( i=0; i<total_ln; i++)
{
long start = data[i].index;
long end = data[i+1].index;
long currentIndex = 0;
long endIndex = total_col;
cout<<" ";
cout.width(2);
cout<<left<<i;
cout<<" ";
for( ; currentIndex<endIndex; currentIndex++)
{
cout.width(6);
if(currentIndex==i){
if(data[i].value!=0.0f)
{
cout.precision(2);
cout<<right<<fixed<<data[i].value;
}
else if(zero==1)
cout<<right<<0;
else
cout<<right<<"";
continue;
}
if(start<end){
if(data[start].index == currentIndex)
{
if(start<end)
{
cout.precision(2);
cout<<right<<fixed<<data[start].value;
start++;
}
}
else if(zero==1)
cout<<right<<0;
else
cout<<right<<"";
}
else if(zero==1)
cout<<right<<0;
else
cout<<right<<"";
}
cout<<endl;
}
}
/* 将矩阵信息输出到参数指定的文件中 */
void Matrix::OutFile(char* fileName)const
{
ofstream outfile(fileName);
if(!outfile)
throw FileNotFoundException("File Not Found Exception!");

outfile.width(6);
outfile<<total_elem;
outfile.width(6);
outfile<<total_ln;
outfile.width(6);
outfile<<total_col<<"/n";
long sum = 0;
for(long i=0; i<total_ln; i++)
{
if(i<total_col){
outfile.width(6);
outfile<<i;
outfile.width(6);
outfile<<i;
outfile.width(12);
outfile.precision(9);
outfile<<fixed<<data[i].value;
sum++;
}
if(sum%3==0)
{
outfile.width(6);
outfile<<"/n";
}
for(long j=data[i].index; j<data[i+1].index; j++)
{
outfile.width(6);
outfile<<i;
outfile.width(6);
outfile<<data[j].index;
outfile.width(12);
outfile.precision(9);
outfile<<fixed<<data[j].value;
sum++;
if(sum%3==0)
{
outfile.width(6);
outfile<<"/n";
}
}
}
outfile.close();
}
/* 得到一n*1值全为1的矩阵 */
Matrix Matrix::GetE(const long ln)
{
Matrix m;
m.total_col = 1;
m.total_ln = ln;
m.total_elem = ln*2-1;

m.data.GetSpace(m.total_elem+1);
m.data[0].value = 0;
m.data[0].index = m.total_ln+1;
m.data[1].value = 0;
m.data[1].index = m.total_ln+1;
if(ln<2)
return m;

long i;
for(i=2; i<m.total_ln+1; i++)
{
m.data[i].value = 0;
m.data[i].index = m.data[i-1].index + 1;
}
for(; i<m.total_elem+1; i++)
{
m.data[i].value = 0;
m.data[i].index = 0;
}
return m;
}
/* 得到一个total_ln行的单位列矩阵，其中第j行的值为1，其余的都为0 */
Matrix Matrix::GetE(const long total_ln, const long j)
{
Matrix m;
m.total_ln = total_ln;
m.total_col = 1;
if(j==0)
{
m.total_elem = m.total_ln;
m.data.GetSpace(m.total_elem+1);
m.data[0].value = 1;
m.data[0].index = m.total_ln+1;
for(long i=1; i<m.total_elem+1; i++)
{
m.data[i].index = m.data[i-1].index;
m.data[i].value = 0;
}
}
else
{
m.total_elem = m.total_ln+1;
m.data.GetSpace(m.total_elem+1);
long i=0;
for(i=0; i<m.total_ln+1; i++)
{
m.data[i].value = 0;
if(i<=j)
m.data[i].index = m.total_ln+1;
else
m.data[i].index = m.total_ln+2;
}
m.data[i].value = 1;
m.data[i].index = 0;
}
return m;
}
/* 用折半查找定位元素 返回元素的位置，返回-1说明不存在该元素 */
long Matrix::QuickLocate(const long ln, const long col)const
{
if(data.point==0)
{
cout<<"矩阵中不存在任何元素"<<endl;
return -1;
}
if(ln>=total_ln)
return -1; // 行标越界
if(col>=total_col)
return -1; // 列标越界
long endPos = data[ln+1].index-1;
long startPos = data[ln].index;

if(ln==col)
return ln;
while(startPos<=endPos)
{
long pos = (startPos + endPos)/2;
if(data[pos].index == col)
return pos;
else if(data[pos].index<col)
startPos = pos+1;
else
endPos = pos-1;
}
return -1; // 有元素 但是无此元素
}

/* 用快速排序法将ln、data[].value 和 data[].index 分区，
将小于 关键码 的元素前置，大于关键码的后置，关键码在中间 */
long Matrix::QSPartitionLn(long *ln, long low, long high)
{
if(ln==0)
{
cout<<"待分区数组为空，无法排序!"<<endl;
throw QuickSortException("Quick Sort Exception!");
// return -1;
}
if(data.Length()==0){
cout<<"矩阵中无数据，无法实现数组分区!"<<endl;
throw QuickSortException("Quick Sort Exception!");
}

long tempLn = ln[low]; /* 行标 关键码，此时low位置可以被覆盖 */
Double tempValue = data[total_ln+1+low].value; /* value随 ln 变动 */
long tempIndex = data[total_ln+1+low].index; /* index随 ln 变动 */

/* 因为是将关键码置于中间，所以当low==high时，就退出while
然后将关键码的值置于low(high)处 */
while(low<high)
{
/* 在高端搜索比关键码小的值，找到或者low==hign退出此while */
while(low<high&&ln[high]>=tempLn)
high--;
/* 上面退出的while循环保证：
Case One：假如有比关键码小的数，则此时的high代表此元素的位置。
下面的语句就是实现将大于关键码的数前置，以覆盖low位置的数
Case Two：假如没有比关键码小的数，则此时low==high，说明关键码就是最小值
此时，下面的while循环也将不会进入，因为low==high，
两个位置的元素覆盖无意义，因为low==high，可以添加限制条件，不等才覆盖 */
ln[low] = ln[high]; /* high 位置可以被覆盖 */
data[total_ln+1+low].value = data[total_ln+1+high].value;
data[total_ln+1+low].index = data[total_ln+1+high].index;

/* 在低端搜索比关键码大的值，找到或者low==hign退出此循环 */
/* 此while的第次循环必将不用，因为必定有ln[low]<=ln[high]，上面循环所致。 */
while(low<high&&ln[low]<=tempLn)
low++;
/* 上面的此while，如果退出循环，可以保证：
Case One：假如有比关键码大的数，则此时的low代表此元素的位置
下面的语句就是实现将小于关键码的数后置，一覆盖high位置的数，
为什么可以覆盖high位置的数，因为上面已将high位置的数覆盖了low位置
Case two：假如没有比关键码大的数，则此时low==high，说明关键码就是最大值
下面的语句无意义，因为 low==high， 可以添加限制条件，不等才覆盖 */
ln[high] = ln[low];
data[total_ln+1+high].value = data[total_ln+1+low].value;
data[total_ln+1+high].index = data[total_ln+1+low].index;
}

/* 此时的low或者high(low==high)必定是中间位置，前面的必将小于或等于关键码，
后面的必将大于或等于关键码，所以此位置放关键码 */
ln[low] = tempLn;
data[total_ln+1+low].value = tempValue;
data[total_ln+1+low].index = tempIndex;

/* 返回 支点位置*/
return low;
}

/* 实现 递归*/
void Matrix::QuickSortLn(long *ln, const long low, const long high)
{
if(ln==0)
return ;
/* 只有当里面至少有2个元素时才有必要继续递归。 */
if(low<high)
{
/* 得到支点位置，pivot位置前的必将小于或等于它的值 */
/* pivot位置后的元素必将大于或等于它的值 */
long pivot = QSPartitionLn(ln, low, high);
/* 将pivot前面的作为一个分区空间，再将其分为两个区，直到只有一个元素 */
QuickSortLn(ln, low, pivot-1);
/* 将pivot后面的作为一个分区空间，再将其分为两个区，直到只有一个元素 */
QuickSortLn(ln, pivot+1, high);
}
}

/* 用快速排序法 实现将列标分为两个区间
/* 算法和 QSPartitionLn一样 */
long Matrix::QSPartitionCol(long low, long high)
{
if(data.Length()==0)
{
cout<<"矩阵中无数据，无法实现列标分区!"<<endl;
throw QuickSortException("Quick Sort Exception!");
}
long tempIndex = data[low].index;
Double tempValue = data[low].value;

while(low<high)
{
while(low<high&&data[high].index>=tempIndex)
high--;
data[low].index = data[high].index;
data[low].value = data[high].value;

while(low<high&&data[low].index<=tempIndex)
low++;
data[high].index = data[low].index;
data[high].value = data[low].value;
}
data[low].index = tempIndex;
data[low].value = tempValue;

return low;
}

/* 实现递归 算法 同QuickSortLn 完全一样 */
void Matrix::QuickSortCol(const long low,const long high){
if(low<high)
{
long pivot = QSPartitionCol(low, high);
QuickSortCol(low, pivot-1);
QuickSortCol(pivot+1, high);
}
}

/* 初始化矩阵信息 需在Evaluate()函数中替total_elem等数据成员赋值 */
void Matrix::Init(Double value)
{
if (data.point)
{
data.Delete();
}
data.GetSpace(total_elem+1);

for(long i=0; i<total_elem+1; i++)
{
data[i].index = -1;
data[i].value = value;
}
}