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