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