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

 

   Martix类声明参见：http://blog.csdn.net/lq_0402/archive/2010/11/13/6007719.aspx

*/

 

 

 

/* 两矩阵相加 */
Matrix& Matrix::operator +(const Matrix &genMatrix)
{
 if(&genMatrix==0)
  return *this;
 /* 如果两矩阵不能相加 则返回当前矩阵 */
 if(total_ln!=genMatrix.total_ln||total_col!=genMatrix.total_col)
 {
  cout<<"当前的两矩阵不能相加"<<endl;
  return *this;
 }
 /* 估算相加得到的矩阵最多的非零元素个数 不会超过两个矩阵的非零元素之和 */
 long totalElem = total_elem+genMatrix.total_elem-total_ln;
 Node tempData;
 tempData.GetSpace(totalElem+1); /* 辅助数组，存放临时的值 */
 long i=0;
 for(i=0; i<totalElem+1; i++)
 {
  tempData[i].value = 0;
  tempData[i].index = -1;
 }
 long *num = new long[total_ln]; /* 每行非零元素个数，得初始化为0，后面要通过它累加 */
 for(i=0; i<total_ln; i++) 
  num[i] = 0;

 /* 对角元素相加 以及total_ln位置元素相加 */
 long offset;
 for(offset=0; offset<=total_ln; offset++)
  tempData[offset].value = data[offset].value+genMatrix.data[offset].value;

 /* 实现非零非对角元素的相加 */
 /* 算法：
  1、针对每行实现相加，首先逐个的比较两个矩阵当前行的非零元素值的列标，将最小的列表值存入
  tempData中(两个矩阵的每行的列标值都是有序的，从小到大排列)，列标相同则相加再放入tempData
  中。
 */
 for(i=0; i<total_ln; i++)
 {
  /* */
  long j=data[i].index, k = genMatrix.data[i].index;
  for( ; j<data[i+1].index&&k<genMatrix.data[i+1].index; )
  {
   if(data[j].index<genMatrix.data[k].index)
   {
    tempData[offset].value = data[j].value;
    tempData[offset].index = data[j].index;
    if(tempData[offset].value!=0.0f){ /* 剔除相加之后小于精度要求的元素 */
     num[i]++;  /* 该行非零非对角元素个数加1 */
     offset++;
    }
    j++;
   }
   else if(data[j].index==genMatrix.data[k].index)
   {
    tempData[offset].value = data[j].value+genMatrix.data[k].value;
    tempData[offset].index = data[j].index;
    if(tempData[offset].value!=0.0f)
    {
     num[i]++;
     offset++;
    }
    j++;
    k++;
   }
   else
   {
    tempData[offset].value = genMatrix.data[k].value;
    tempData[offset].index = genMatrix.data[k].index;
    if(tempData[offset].value!=0.0f)
    {
     num[i]++;
     offset++;
    }
    k++;
   }
  }
  /* 上面的for循环结束时，有三种情况：
   1、j==data[i+1].index&&k==data[i+1].index 不用考虑
   2、j==data[i+1].index&&k<data[i+1].index 则说明当前矩阵的该行元素已添加完，而genMatrix
   矩阵的非零非对角元素还没有添加完，继续将其添加完
   3、j<data[i+1].index&&k==data[i+1].index 则说明genMatrix该行元素已添加完，而当前矩阵的
   该行的非零非对角元素还没有添加完，继续将其添加完
  */
  if(j==data[i+1].index&&k<genMatrix.data[i+1].index)
  {
   tempData[offset].value = genMatrix.data[k].value;
   tempData[offset].index = genMatrix.data[k].index;
   if(tempData[offset].value!=0.0f)
   {
    num[i]++;
    offset++;
   }
   k++;
  }else if(j<data[i+1].index&&k==genMatrix.data[i+1].index)
  {
   tempData[offset].value = data[j].value;
   tempData[offset].index = data[j].index;
   if(tempData[offset].value!=0)
   {
    num[i]++;
    offset++;
   }
   j++;
  }
 }

 /* 每行第一个非零元素的序列号 */
 tempData[0].index = total_ln+1;
 for(i=1; i<=total_ln; i++)
  tempData[i].index = tempData[i-1].index + num[i-1];

 /* 将tempData赋值给matrix.data */
 total_elem = offset-1;

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

/* 赋值运算，将创建一个与参数matrix一样的矩阵 */
Matrix& Matrix::operator =(const Matrix &matrix)
{
 if(&matrix == this)
  return *this;
 this->total_elem = matrix.total_elem;
 this->total_ln = matrix.total_ln;
 this->total_col = matrix.total_col;

 if (matrix.data.Length() == 0)
 {
  data.Delete();
  return *this;
 }

 try
 {
  this->data.GetSpace(matrix.total_elem+1);
 }
 catch(bad_alloc x)
 {
  cout<<"Memory alloc is fail!"<<endl;   //看看是不是输出alloc bad之类的提示
 }

 for(int i=0; i<matrix.total_elem+1; i++)
 {
  data[i].index = matrix.data[i].index;
  data[i].value = matrix.data[i].value;
 }
 return *this;
}

Matrix& Matrix::operator -(const Matrix& genMatrix)
{
 Matrix matrix = genMatrix;
 return (*this)+matrix*(-1);
}

ostream& operator <<(ostream& os, const Matrix& matrix)
{
 matrix.DisplayAsMatrix("", 1);

 return os;
}
/* 在矩阵中提取某一行作为一个矩阵 */
/* 算法：  提取矩阵中第ln行的元素作为一新矩阵(向量)
 1、创建一个1行total_col列的矩阵matrix
 2、替matrix.data开辟 ln行非零元素个数+1 个空间
 3、matrix的对角元素应该为data中ln行元素中下标为0的元素，将其添进matrix.data的第0位置
 4、将data中的元素添进matrix.data中，因为data中ln行的对角元素在data的前面，而其他的在后面，
 往matrix.data里填数据时，又得保持列标的有序，所以。该步骤得分为一下几种情况
  I、将data中的ln行的列标在ln之前的元素添进matrix.data中(得排除列标为0的元素,它为对角元素)
  II、将data中的ln行的列标为的ln的元素(对角元素)添进matrix.data中(得排除列标为0的元素,它为对角元素)
  III、将data中的ln行的列标在ln之后的元素添进matrix.data中(得排除列标为0的元素,它为对角元素)

*/
Matrix Matrix::GetVector(const long ln) const
{
 if(ln>=total_ln)
 {
  cout<<"所提取的行信息错误: "<<ln<<" >= "<<total_ln<<endl;
  exit(0);
 }

 Matrix matrix;
 matrix.total_ln = 1;
 matrix.total_col = total_col;

 long nonZero = 0;
 /* 当ln==0时，提取向量很容易 */
 if(ln==0)
 {
  matrix.total_elem = data[ln+1].index-data[ln].index+1;
  matrix.data.GetSpace(matrix.total_elem+1);
  long i=0;
  for(i=0; i<matrix.total_elem+1; i++)
  {
   matrix.data[i].value = 0;
   matrix.data[i].index =-1;
  }
  matrix.data[0].value = data[0].value;
  matrix.data[0].index = 2;
  matrix.data[1].value = 0;
  matrix.data[1].index = matrix.data[0].index+matrix.total_elem-matrix.total_ln;
  for(i=2; i<=matrix.total_elem; i++)
  {
   matrix.data[i].value = data[i-2+total_ln+1].value;
   matrix.data[i].index = data[i-2+total_ln+1].index;
  }
  return matrix;
 }

 // 当ln不等于0时，data[ln].value也算是非零非对角元素，假若他不为0
 // 当矩阵中存在下标[ln,0]不为0的元素时，其为对角元素，不是非零非对角元素
 nonZero = data[ln+1].index-data[ln].index; // 非零非对角元素个数
 long pos = QuickLocate(ln,0);
 if(data[ln].value!=0) // 假如data中ln行的对角元素不为0，则提取之后nonZero++;
  nonZero++;
 if(pos!=-1)
  nonZero--;

 matrix.total_elem = nonZero+1;  // 总的非零元素 等于 非零非对角 + 对角
 matrix.data.GetSpace(matrix.total_elem+1);
 long i=0;
 for(i=0; i<matrix.total_elem+1; i++)
 {
  matrix.data[i].value = 0;
  matrix.data[i].index = 0;
 }

 if(pos!=-1)
  matrix.data[0].value = data[pos].value;
 else
  matrix.data[0].value = 0;
 matrix.data[0].index = matrix.total_ln+1;

 matrix.data[1].value = 0;
 matrix.data[1].index = matrix.data[0].index + nonZero;

 long offset = matrix.total_ln+1;
 /* 将data中ln行的列标在ln位置前的元素置于matrix.data中 */
 for(i=data[ln].index; i<data[ln+1].index&&data[i].index<ln; i++)
 {
  if(data[i].index!=0)
  {
   matrix.data[offset].value = data[i].value;
   matrix.data[offset].index = data[i].index;
   offset++;
  }
 }
 /* 将data中ln行的列标为ln位置的元素添加到matrix.data中, 假如它不为0 */
 if(data[ln].value!=0)
 {
  matrix.data[offset].value = data[ln].value;
  matrix.data[offset].index = ln;
  offset++;
 }
 /* 将data中ln行的列标在ln位置后的置于matrix.data中 */
 for( ; i<data[ln+1].index; i++)
 {
  matrix.data[offset].value = data[i].value;
  matrix.data[offset].index = data[i].index;
  offset++;
 }

 return matrix;
}

/* 在矩阵中提取col列作为一个矩阵 */
Matrix Matrix::GetColVector(const long col) const
{
 Matrix matrix;
 if(col>=this->total_col)
 {
  cout<<"所提取的列错误 col>="<<this->total_col<<endl;
  return matrix;
 }

 matrix.total_col = 1;
 matrix.total_ln = this->total_ln;
 matrix.data.GetSpace(this->total_ln+this->total_ln);
 matrix.data[this->total_ln].value = 0;

 long i=0, j=0, offset = this->total_ln+1;
 for(i=0; i<this->total_ln; i++)
 {
  if(col==0)
  {
   if(i==0)
   {
    matrix.data[i].value = this->data[0].value;
    matrix.data[i].index = this->total_ln+1;
    matrix.data[i+1].index = this->total_ln+1;
   }
   else
   {
    long nonDZero = 0;
    long pos = this->QuickLocate(i,col);
    if(pos!=-1)
    {
     matrix.data[offset].value = this->data[pos].value;
     matrix.data[offset].index = 0;
     offset++;
     nonDZero = 1;
    }
    matrix.data[i].value = 0;
    matrix.data[i+1].index = matrix.data[i].index + nonDZero;
   }
  }
  else
  {
   if(i==0)
   {
    matrix.data[i].value = 0;
    long pos = this->QuickLocate(i,col);
    if(pos!=-1)
    {
     matrix.data[i].value = this->data[pos].value;
    }
    matrix.data[i].index = this->total_ln+1;
    matrix.data[i+1].index = this->total_ln+1;
   }
   else if(i==col)
   {
    if(this->data[i].value!=0)
    {
     matrix.data[offset].value = this->data[i].value;
     matrix.data[offset].index = 0;
     offset++;
     matrix.data[i+1].index = matrix.data[i].index + 1;
    }
    else
    {
     matrix.data[i+1].index = matrix.data[i].index;
    }
    matrix.data[i].value = 0;
   }
   else
   {
    long nonDZero = 0;
    long pos = this->QuickLocate(i,col);
    if(pos!=-1)
    {
     matrix.data[offset].value = this->data[pos].value;
     matrix.data[offset].index = 0;
     offset++;
     nonDZero = 1;
    }
    matrix.data[i].value = 0;
    matrix.data[i+1].index = matrix.data[i].index + nonDZero;
   }
  }
 }
 matrix.total_elem = offset-1;
 return matrix;
}

/* 求矩阵的范数(模) */
Double Matrix::GetModulus() const
{
 Double sum = 0;
 
 for(long i=0; i<=total_elem; i++)
 {
  if(i!=total_ln)
   sum += data[i].value*data[i].value;
 }
 sum = sqrt(sum.GetData());

 return sum;
} 

/* 返回两矩阵的点积 */
Double Matrix::GetDotMetrix(const Matrix& matrix)const
{
 Double dot = 0;
 if (matrix.data.Length() == 0)
 {
  return dot;
 }
 long i=0;
 for(i=0; i<total_ln; i++)
  dot += data[i].value*matrix.data[i].value;
 for(i=0; i<total_ln; i++)
 {
  int j=data[i].index, k=matrix.data[i].index;
  for(; j<data[i+1].index&&k<matrix.data[i+1].index; )
  {
   if(data[j].index==matrix.data[k].index)
   {
    dot += data[j].value*matrix.data[k].value;
    k++;
    j++;
   }
   else if(data[j].index<matrix.data[k].index)
    j++;
   else
    k++;
  }
 }
// return sqrt(dot);
 return dot;
}

/* 向矩阵中添加一行，参数vector是一个只有一行n列的矩阵 */
int Matrix::AddLine(const Matrix &vector)
{
 if(vector.total_ln!=1)
 {
  cout<<"新添加的行信息不正确"<<endl;
  return 0;
 }

 if(this->data.point == 0)
 {
  *this = vector;
  return 1;
 }
 /* 开辟一辅助Node数组，记录data和vector.data中的值 */
 Node tempData;
 tempData.GetSpace(total_elem+vector.total_elem+1+vector.total_ln);
 long i=0;
 for(i=0; i<total_elem+vector.total_elem+1; i++)
 {
  tempData[i].value = 0;
  tempData[i].index = -1;
 }
 /* 将data中前total_ln个(对角元素信息及每行pos值)值赋值给tempData */
 for(i=0; i<total_ln; i++)
 {
  tempData[i].value = data[i].value;
  tempData[i].index = data[i].index+1; /* 添加了一个对角元素，所以data中每行第一个非零非对角元素列标要加1 */
 }
 /* 添加新行的对角元素，及该行第一个非零非对角元素的序号 */
 long pos = vector.QuickLocate(0,total_ln);
 if(pos!=-1)
  tempData[i].value = vector.data[pos].value; /* i==total_ln */
 else
  tempData[i].value = 0;
 tempData[i].index = data[i].index+1;
 i++;
 tempData[i].value = 0;  /* i==total_ln+1 */
 i++;
 /* 将data中非零非对角元素信息复制到tempData中 */
 for( ; i<total_elem+1+1; i++)
 { /* 注意循环的次数 */
  tempData[i].value = data[i-1].value;
  tempData[i].index = data[i-1].index;
 }
 /* 将vector.data中的非零非对角元素信息复制到tempData中 */
 long j = vector.total_ln+1;
 if(vector.data[0].value!=0)
 {
  tempData[i].value = vector.data[0].value;
  tempData[i].index = 0;
  i++;
 } 
 for( ;j<vector.total_elem+1 ;j++)
 {
  if(vector.data[j].index!=total_ln)
  {
   if(vector.data[j].value!=0)
   {
    tempData[i].value = vector.data[j].value;
    tempData[i].index = vector.data[j].index;
    i++;
   }
  }
 }
 tempData[total_ln+1].index = i;//tempData[i-1].index+vector.data[1].index-vector.data[0].index;

 /* 更新总元素 和 总行数 信息 */
 total_elem = i-1;
 total_ln++;

 data = tempData;
 return 1;
}

/* 先矩阵中添加一个 列向量，成功返回1，否则返回0 */
int Matrix::AddCol(const Matrix& vector){
 if(&vector==0||&vector.data==0||vector.total_col!=1)
 {
  cout<<"添加的列向量信息不正确"<<endl;
  return 0;
 }
 if(this==0||(this->data.point!=0&&this->total_ln!=vector.total_ln) ){
  cout<<"当前矩阵信息不正确，有可能是当前矩阵为NULL，或者当前矩阵的行与添加的向量行数不一致"<<endl;
  return 0;
 }

 // 假如当前矩阵里没元素，则直接赋值为vector即可
 if(this->data.point==0)
 {
  this->data = vector.data;
  this->total_col = vector.total_col;
  this->total_elem = vector.total_elem;
  this->total_ln = vector.total_ln;
  return 1;
 }
 Node tempData;
 tempData.GetSpace(this->total_elem+vector.total_ln);
 tempData[this->total_ln].value = 0;

 // 假如当前矩阵的列数小于vector的总行数，则说明vector的第this->total_col个元素应该以对角元素插入this中
 // 否则以非零非对角元素插进来
 long i = 0, j = 0;
 long offset = this->total_ln+1, nonDZero = 0;

 if(this->total_col<vector.total_ln)
 {
  if(vector.data[0].value!=0)
  {
   tempData[0].index = this->data[0].index;
   tempData[0].value = this->data[0].value;
   tempData[1].index = this->data[1].index + 1;
   tempData[1].value = this->data[1].value;
   for(j=this->data[0].index; j<this->data[1].index; j++)
   {
    tempData[offset].value = this->data[j].value;
    tempData[offset].index = this->data[j].index;
    offset++;
   }
   tempData[offset].value = vector.data[0].value;
   tempData[offset].index = this->total_col;
   offset++;
  }
  else
  {
   tempData[0].index = this->data[0].index;
   tempData[0].value = this->data[0].value;
   tempData[1].index = this->data[1].index;
   tempData[1].value = this->data[1].value;
   for(j=this->data[0].index; j<this->data[1].index; j++)
   {
    tempData[offset].value = this->data[j].value;
    tempData[offset].index = this->data[j].index;
    offset++;
   }
  }
  for(i=1; i<vector.total_ln; i++)
  {
   if(i==this->total_col)
   {
    nonDZero = vector.data[i+1].index-vector.data[i].index;
    if(nonDZero!=0)
    {
     tempData[i].value = vector.data[vector.data[i].index].value; // ##
    }
    else
    {
     tempData[i].value = 0;           // ##
    }
    for(j=this->data[i].index; j<this->data[i+1].index; j++)
    {
     tempData[offset].value = this->data[j].value;
     tempData[offset].index = this->data[j].index;
     offset++;
    }     
    tempData[i+1].index = tempData[i].index + this->data[i+1].index - this->data[i].index;
   }
   else
   {
    nonDZero = vector.data[i+1].index-vector.data[i].index;
    for(j=this->data[i].index; j<this->data[i+1].index; j++)
    {
     tempData[offset].value = this->data[j].value;
     tempData[offset].index = this->data[j].index;
     offset++;
    }
    if(nonDZero!=0)
    {
     tempData[i].value = this->data[i].value;
     tempData[offset].value = vector.data[vector.data[i].index].value;
     tempData[offset].index = this->total_col;
     offset++;
    }
    else
    {
     tempData[i].value = this->data[i].value;
    }
    tempData[i+1].index = tempData[i].index + nonDZero + this->data[i+1].index - this->data[i].index;
   }
  } 
 }
 else
 {
  tempData[0].value = this->data[0].value;
  tempData[0].index = this->data[0].index;

  for(i=0; i<this->total_ln; i++)
  {
   nonDZero = vector.data[i+1].index - vector.data[i].index;
   tempData[i].value = this->data[i].value;
   for(j=this->data[i].index; j<this->data[i+1].index; j++)
   {
    tempData[offset].value = this->data[j].value;
    tempData[offset].index = this->data[j].index;
    offset++;
   }
   if(nonDZero!=0)
   {
    tempData[offset].value = vector.data[vector.data[i].index].value;
    tempData[offset].index = this->total_col;
    offset++;
   }
   else
   {
   }
   tempData[i+1].index = tempData[i].index + nonDZero + this->data[i+1].index - this->data[i].index;
  }
 }
 this->data = tempData;
 this->total_ln = this->total_ln;
 this->total_col = this->total_col + 1;
 this->total_elem = offset-1;

 return 1;
}