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

 

/* 行初等变换 第base行乘以gene, 没有过滤掉乘以gene之后小于门槛值的元素 */
void Matrix::LnTransformMultiply(const long base,const Double& gene)
{
 long offset=data[base].index;
 /* 该行的非零非对角元素乘以gene */
 long i=0;
 for(i=data[base].index; i<data[base+1].index; i++)
 {
  data[i].value *= gene;
   data[offset].value = data[i].value;
   data[offset].index = data[i].index;
   offset++;
 }
 if(offset!=i)
 {
  for( ; i<total_elem+1; i++)
  {
   data[offset].value = data[i].value;
   data[offset].index = data[i].index;
   offset++;
  }
  total_elem = offset-1;
  for(long j=base+1; j<total_ln+1; j++)
   data[j].index -= i-offset;
 }

 /* 该行对角元素乘以gene */
 data[base].value *= gene;
}

/* 行初等变换 第base行加上第ln行  */
/* 算法：
 1、开辟一个辅助数组tempData存放一行的数据，其长度为total_col
 2、将base行的数据根据下标放于tempData中。
 3、将ln行的数据根据下标，加到tempData相应下标位置。
 4、开辟一个辅助数组tempDataAll，存放base和ln相加之后的整个新data的值
 5、将tempDataAll覆盖原有的data
*/
void Matrix::LnTransfromAdd(const long base,const long ln)
{
 Double *tempData = new Double[(total_col>total_ln?total_col:total_ln)+10];
 long nonZero = 0; // 非零非对角元素个素 */
 long i=0;
 for( i=0; i<total_col+1; i++)
  tempData[i] = 0;
 /* 将base行的元素置于tempData中 */
 /* 处理非零非对角元素 */
 for(i=data[base].index; i<data[base+1].index; i++ )
 {
  if(data[i].value!=0.0f)
  {
   tempData[data[i].index] = data[i].value;
   nonZero++;
  }
 }
 /* 处理对角元素 */
 tempData[base] = data[base].value;

 /* 将ln行的元素加到tempData相应位置 */
 /* 处理非零非对角元素 */
 for(i=data[ln].index; i<data[ln+1].index; i++)
 {
  // 先判断其值是否为0
  if(tempData[data[i].index]==0&&data[i].index!=base)
   nonZero++;
  tempData[data[i].index] += data[i].value;
  if(tempData[data[i].index]==0&&data[i].index!=base)
   nonZero--;
 }
 /* 处理对角元素 */
 if(tempData[ln]==0)
  nonZero++;
 tempData[ln] += data[ln].value;
 if(tempData[ln]==0)
  nonZero--;

 /* 将数组tempData中的数复制到data中的相应位置，涉及到三个步骤：
  1、新增一辅助data空间 tempDataAll
  2、将data和tempData中的数据复制到tempDataAll中
  3、用tempDataAll覆盖data
 */
 long overNum = nonZero-(data[base+1].index-data[base].index);
 total_elem = total_elem+overNum;
 Node tempDataAll;
 tempDataAll.GetSpace(total_elem+1);
 /* 将辅助数组初始化 */
 for(i=0; i<total_elem+1; i++)
 {
  tempDataAll[i].value = 0;
  tempDataAll[i].index = -1;
 }

 /* 前total_ln+1个位置 */
 for(i=0; i<total_ln+1; i++)
 {
  tempDataAll[i].value = data[i].value;
  if(i<=base)
   tempDataAll[i].index = data[i].index;
  else
  {
   tempDataAll[i].index = data[i].index+overNum;
  }
 }
 /* tota_ln+2 到 data[base].index-1 个位置 */
 for( ; i<data[base].index; i++){
  tempDataAll[i].value = data[i].value;
  tempDataAll[i].index = data[i].index;
 }
 long j=0;
 /* 新的 base 行的数据 */
 for( ; j<total_col; j++)
 {
  if(j==base){
   tempDataAll[base].value = tempData[j];
  }
  else if(tempData[j]!=0.0f)
  {
   tempDataAll[i].value = tempData[j];
   tempDataAll[i].index = j;
   i++;
  }
 }
 /* data[base+1].index 到最后 */
 j=data[base+1].index;
 for( ; i<total_elem+1; )
 {
  tempDataAll[i].value = data[j].value;
  tempDataAll[i].index = data[j].index;
  i++;
  j++;
 }

 data = tempDataAll;
 delete[] tempData;
}

/* 矩阵转置 考虑nm矩阵，当n>m、n=m、n<m三种情况，不管哪种情况，对角元素个数为小的那个 */
Matrix& Matrix::Transpose()
{
 if(this->data.point==0)
  return *this;

 long over = 0;
 over = total_col - total_ln;
 Node tempData;
 tempData.GetSpace(total_elem+1+over); /* 缓存器 */
 long *numCol = new long[total_col];   /* 记录每列非零非对角元素个数 */
 long *countCol = new long[total_col];
 long i=0;

 for(i=0; i<total_col; i++)
 {
  numCol[i] = 0;
  countCol[i] = 0;
 }
 /* 将辅助数组初始化 */
 for(i=0; i<total_elem+1+over; i++)
 {
  tempData[i].value = 0;
  tempData[i].index = -1;
 }

 long smallNum = total_ln<total_col?total_ln:total_col;
 /* tempData前面的对角元素赋值，以及计算每列的非零元素个数 */
 for(i=0; i<total_elem+1; i++)
 {
  if(i<smallNum&&i!=total_ln)
  {
   if(data[i].value!=0.0f)
    tempData[i].value = data[i].value;
   else
    tempData[i].value = 0;
  }
  else if(i>=smallNum&&i!=total_ln)
  {
   if(i<total_ln)
   {
    if(data[i].value!=0.0f)
    {
     numCol[i]++;
    }
   }
   else
   {
    numCol[data[i].index]++;
   }
  }
 }
 tempData[total_col].value = data[total_ln].value;
 tempData[0].index = total_col+1;

 /* 计算每行第一个非零非对角元素的位置 */
 for(i=1; i<=total_col; i++)
  tempData[i].index = tempData[i-1].index + numCol[i-1];

 /* smallNum前的元素已经被赋值了。 */
 /* 替smallNum及以后的元素赋值，index从total_col+1开始 */
 long j = 1;
 for(i=total_ln+1; i<=total_elem; i++)
 {
  /* 当前元素的列标 */
  long tempCol = data[i].index;
  /* 当前元素转置之后的这列第一个非零元素的位置 */
  long temp = 0;
  temp = tempData[tempCol].index;  // tempCol行第一个非零非对角元素的position
  tempData[temp+countCol[tempCol]].value = data[i].value;
  /* 得到当前元素的行标 */
  while(j<=total_ln)
  {
   /* 假如当前元素的位置小于第j行第一个非零非对角元素元素的序列号 */
   /* 那么这一行必定在第j-1行 【第一次符合条件便退出循环】*/
   /* data 中的下一元素的行标必定大于或等于当前元素的行标 */
   if(i<data[j].index)
   {
    tempData[temp+countCol[tempCol]].index = j-1;
    /* 这一列已添加一个元素，该列计数器自加1 */
    countCol[tempCol]++;
    break;
   }
   j++;
  }
 }
 total_elem += over;
 long temp = total_ln;
 total_ln = total_col;
 total_col = temp;
 /* 用缓存器里的元素值覆盖data中的值，以实现转置 */

 data = tempData;
 delete[] numCol;
 numCol = NULL;
 delete[] countCol;
 countCol = NULL;

 return *this;
}

/* 单位化一个矩阵 */
void Matrix::ToBeIdentityMatrix()
{
 total_elem = total_ln;
 data.GetSpace(total_elem+1);

 for(long i=0; i<=total_ln;i++)
 {
  if(i<total_ln)
  {
   data[i].value = 1;
   data[i].index = total_ln+1;
  }
  else if(i==total_ln)
  {
   data[i].value = 0;
   data[i].index = total_ln+1;
  }
 }
}
/* 将当前矩阵化为0矩阵 */
void Matrix::ToBeZeroMatrix()
{
 total_elem = total_ln;
 for(long i=0; i<total_elem+1; i++)
 {
  data[i].value = 0.0;
  data[i].index = total_ln+1;
 }
}
// 求逆矩阵 返回当前矩阵的逆矩阵
/* 思路：因为矩阵可逆，所以，任何一行或者任意列，都至少存在一个非零元素。
 1、假若A00不为零(若为零，则查找第0列不为零的行m，将其加上第m行，使A00不为零)，将第0行所元素都除以
  A00，将A00化为1，然后针对任意第0列不为零的行m，将其加上-Am0*A0j(j=0,1,2,total_col-1)，即将Am0
  化为0。
 2、假若A11不为零，....
 */

void Matrix::ContraryMatrix()
{
 Matrix matrix = *this;
 this->total_col = matrix.total_col;
 this->total_ln = matrix.total_ln;
 this->total_elem = this->total_ln+1;
 this->ToBeIdentityMatrix();

 for(long i=0; i<matrix.total_ln; i++)
 {
  /* 第i行元素都除以 Aii,将Aii 化为1 */
  /* 假如Aii==0，那么，得找到一个元素Aji!=0的行j，将其元素加到第i行，使Aii!=0 */
  /* Note: 假如找到的所有Aji!=0中的都有j<i，则说明此矩阵不满秩，因为对这些Aji，都可以通过
  Ajj(第j列只有Ajj==1,其他已为0)将其转化为0，以导致i列全为0。
  */
  if(matrix.data[i].value==0)
  {
   /* 在第i列找不为0的元素Aji，如果j<i，则继续查找。直到j>i，或者自动退出循环。*/
   /* 此处代码有冗余 变量ok可以省掉 */
   long ok=0, j = 0;
   for( j=0; j<total_ln; j++)
   {
    long pos = -1;
    if(j!=i)
     pos = matrix.QuickLocate(j,i);
    if(pos!=-1)
    {
     if(j<i)
      ok = j;
     else if(j>i)
      break;
    }
   }
   if(j==total_ln&&ok==0)
   {
    cout<<"此矩阵不满秩!"<<endl;
    return;
   }
   else if(ok==0)
    matrix.LnTransfromAdd(i, j);
   else
   {
    cout<<"此矩阵不满秩!"<<endl;
    return;
   }
  }
  /* 将matrix和indentityMatrix的第i行所有元素除以Aii */
  long j = 0;
  for(j=matrix.data[i].index; j<matrix.data[i+1].index; j++)
   matrix.data[j].value = matrix.data[j].value/matrix.data[i].value;

  for(j = data[i].index; j<data[i+1].index; j++)
   data[j].value = data[j].value/matrix.data[i].value;
  data[i].value /= matrix.data[i].value;
  matrix.data[i].value /= matrix.data[i].value;

  /* 通过行初等变换，将第i列所有元素(除了Aii)都化为0
     此时i是定值，将除了第i行以外的行都减去 Aji*Aii (j=1,2,3,..total_ln且j!=i) */
  for(j=0; j<matrix.total_ln; j++)
  {
   /* 将第i行乘以Aji */
   if(j!=i)
   {
    /* 查找Aji，如果为非零，则通过行变换Aji+(-1)*Aji*Aii变成0 */
    long pos = matrix.QuickLocate(j,i);
    /* pos == -1 说明此元素为0 */
    if(pos!=-1)
    {
     /* 将第i行临时存储，包括右端向量 */
     /* 对矩阵matrix进行操作 */
     long k;
     long temp = matrix.data[i+1].index-matrix.data[i].index;
     Double *tempLn = new Double[temp+2];
     tempLn[0] = matrix.data[i].value;
     for( k=1; k<=temp; k++)
      tempLn[k] = matrix.data[matrix.data[i].index+k-1].value;

     /* 对 单位矩阵进行操作 */
     long tempIden = data[i+1].index-data[i].index;
     Double *tempLnIden = new Double[tempIden+2];
     tempLnIden[0] = data[i].value;
     for(k=1; k<tempIden+1; k++)
      tempLnIden[k] = data[data[i].index+k-1].value;

     matrix.LnTransformMultiply(i, matrix.data[pos].value*(-1));
     LnTransformMultiply(i,matrix.data[pos].value*(-1));

     /* 将第i行加到第j行 */
     matrix.LnTransfromAdd(j, i);
     LnTransfromAdd(j,i);

     /* 将第i行复原 */
     matrix.data[i].value = tempLn[0];
     for(k=1; k<=temp; k++)
      matrix.data[matrix.data[i].index+k-1].value = tempLn[k];
     
     data[i].value = tempLnIden[0];
     for(k=1; k<tempIden+1;k++)
      data[data[i].index+k-1].value = tempLnIden[k];

     delete[] tempLn;
     tempLn = NULL;
     delete[] tempLnIden;
     tempLnIden = NULL;
    }
   }
  }
//  this->DisplayAsMatrix("E", 1);
  matrix.DisplayAsMatrix("matrix", 1);
 }
}

/* 求矩阵的秩 */
long Matrix::GetMatrixRank() const
{
 long rank = 0;
 Matrix AA = *this;

 long i = 0, j=0;
 for( i=0; i< AA.total_col; i++)
 {
  if( i>=(AA.total_col<AA.total_ln ? AA.total_col:AA.total_ln) )
  {
   return AA.total_col<AA.total_ln ? AA.total_col:AA.total_ln;
  }
  // 假如i行的对角元素为0，则需要通过行变换将其变为不等于0，方法如下：
  //  在第i列(j,i)[j>i]中查找非零的元素，假若存在，则将j行的元素加到第i行，假若不存在，则说明第i列可用之前列线性表出
  if(AA.data[i].value==0)
  {
   for(j=i+1; j<AA.total_ln; j++)
   {
    long pos = AA.QuickLocate(j, i);
    if(pos!=-1)
    {
     AA.LnTransfromAdd(i, j);
     break;
    }
   }
   if(j==AA.total_ln)
   {
    i++;
    break;
   }
  }
  j  = i+1;
  for( ; j<AA.total_ln; j++)
  {
   long pos = AA.QuickLocate(j, i);
   // 假如j,i位置的元素为0，则该行不需要变化
   if(pos!=-1)
   {
    // 假如j,i位置的元素不为0，则第j行的元素都必须加上(j,i)/(i,i)*(i,k)[k>=i]，将(j,i)化为0
    // 完成上述步骤，分为两步：
    //   1、首先将第i行的数据备份
    //   2、将第i行乘以(j,i)/(i,i) 注意：在此程序中，某行乘以一个数，没有过滤掉小于门槛值的元素
    //  3、将第j行加上第i行
    //  4、将第i行复原

    // 保存第i行的数据，以便后续步骤中恢复i-1行数据
    long length = AA.data[i+1].index-AA.data[i].index+1;
    Double *tempValue = new Double[length];
    tempValue[0] = AA.data[i].value;
    long kk = 0;
    for(kk=AA.data[i].index; kk<AA.data[i+1].index; kk++)
    {
     tempValue[kk-AA.data[i].index+1] = AA.data[kk].value;
    }

    Double gene = AA.data[pos].value/AA.data[i].value*(-1);
    AA.LnTransformMultiply(i,gene);
    AA.LnTransfromAdd(j,i);

    AA.data[i].value = tempValue[0];
    for(kk=AA.data[i].index; kk<AA.data[i+1].index; kk++)
    {
     AA.data[kk].value = tempValue[kk-AA.data[i].index+1];
    }

    delete[] tempValue;
   }
  }
  rank++;
 }
 return rank;
}

/* 求增广矩阵的秩 */
long Matrix::GetExtendMatrixRank(const Matrix& b) const
{
 if(&b==0||&b.data==0||b.total_col!=1||b.total_ln!=this->total_ln)
 {
  cout<<"右端向量不合法"<<endl;
  return 0;
 }

 long rank = 0;
 Matrix AA = *this;
 Matrix bb = b;

 long i = 0, j=0;
 for( i=0; i< AA.total_col; i++)
 {
  if( i>=(AA.total_col<AA.total_ln ? AA.total_col:AA.total_ln) )
  {
   return AA.total_col<AA.total_ln ? AA.total_col:AA.total_ln;
  }
  // 假如i行的对角元素为0，则需要通过行变换将其变为不等于0，方法如下：
  //  在第i列(j,i)[j>i]中查找非零的元素，假若存在，则将j行的元素加到第i行，假若不存在，则说明第i列可用之前列线性表出
  if(AA.data[i].value==0)
  {
   for(j=i+1; j<AA.total_ln; j++)
   {
    long pos = AA.QuickLocate(j, i);
    if(pos!=-1)
    {
     AA.LnTransfromAdd(i, j);
     bb.LnTransfromAdd(i, j);
     break;
    }
   }
   if(j==AA.total_ln)
   {
    i++;
    break;
   }
  }
  j  = i+1;
  for( ; j<AA.total_ln; j++)
  {
   long pos = AA.QuickLocate(j, i);
   // 假如j,i位置的元素为0，则该行不需要变化
   if(pos!=-1)
   {
    // 假如j,i位置的元素不为0，则第j行的元素都必须加上(j,i)/(i,i)*(i,k)[k>=i]，将(j,i)化为0
    // 完成上述步骤，分为两步：
    //   1、首先将第i行的数据备份
    //   2、将第i行乘以(j,i)/(i,i) 注意：在此程序中，某行乘以一个数，没有过滤掉小于门槛值的元素
    //  3、将第j行加上第i行
    //  4、将第i行复原

    // 保存第i行的数据，以便后续步骤中恢复i-1行数据
    long length = AA.data[i+1].index-AA.data[i].index+1;
    Double *tempValue = new Double[length+1];
    tempValue[0] = AA.data[i].value;
    long kk = 0;
    for(kk=AA.data[i].index; kk<AA.data[i+1].index; kk++)
    {
     tempValue[kk-AA.data[i].index+1] = AA.data[kk].value;
    }
    long bPos = bb.QuickLocate(i,0);
    if(bPos!=-1)
    {
     tempValue[kk-AA.data[i].index+1] = bb.data[bPos].value;
    }else
     tempValue[kk-AA.data[i].index+1] = 0;

    Double gene = AA.data[pos].value/AA.data[i].value*(-1);
    AA.LnTransformMultiply(i,gene);
    AA.LnTransfromAdd(j,i);
    bb.LnTransformMultiply(i, gene);
    bb.LnTransfromAdd(j,i);

    AA.data[i].value = tempValue[0];
    for(kk=AA.data[i].index; kk<AA.data[i+1].index; kk++)
    {
     AA.data[kk].value = tempValue[kk-AA.data[i].index+1];
    }
    if(bPos!=-1)
     bb.data[bPos].value = tempValue[kk-AA.data[i].index+1];

    delete[] tempValue;
   }
  }
  rank++;
 }
 // 判断增广矩阵的秩 和 非增光矩阵的秩是否一致
 for(i=rank; i<bb.total_ln; i++)
 {
  long pos = bb.QuickLocate(i,0);
  if(pos!=-1)
   rank++;
 }

 return rank;
}

/* 判断矩阵和其增广矩阵的秩是否一致，一致返回1，不一致返回0 */
int Matrix::IsEMRequalMR(const Matrix& b) const
{
 if(&b==0||&b.data==0||b.total_col!=1||b.total_ln!=this->total_ln)
 {
  cout<<"右端向量不合法"<<endl;
  return 0;
 }

 long rank = 0;
 Matrix AA = *this;
 Matrix bb = b;

 long i = 0, j=0;
 for( i=0; i< AA.total_col; i++)
 {
  if( i>=(AA.total_col<AA.total_ln ? AA.total_col:AA.total_ln) )
  {
   return AA.total_col<AA.total_ln ? AA.total_col:AA.total_ln;
  }
  // 假如i行的对角元素为0，则需要通过行变换将其变为不等于0，方法如下：
  //  在第i列(j,i)[j>i]中查找非零的元素，假若存在，则将j行的元素加到第i行，假若不存在，则说明第i列可用之前列线性表出
  if(AA.data[i].value==0)
  {
   for(j=i+1; j<AA.total_ln; j++)
   {
    long pos = AA.QuickLocate(j, i);
    if(pos!=-1)
    {
     AA.LnTransfromAdd(i, j);
     bb.LnTransfromAdd(i, j);
     break;
    }
   }
   if(j==AA.total_ln)
   {
    i++;
    break;
   }
  }
  j  = i+1;
  for( ; j<AA.total_ln; j++)
  {
   long pos = AA.QuickLocate(j, i);
   // 假如j,i位置的元素为0，则该行不需要变化
   if(pos!=-1)
   {
    // 假如j,i位置的元素不为0，则第j行的元素都必须加上(j,i)/(i,i)*(i,k)[k>=i]，将(j,i)化为0
    // 完成上述步骤，分为两步：
    //   1、首先将第i行的数据备份
    //   2、将第i行乘以(j,i)/(i,i) 注意：在此程序中，某行乘以一个数，没有过滤掉小于门槛值的元素
    //  3、将第j行加上第i行
    //  4、将第i行复原

    // 保存第i行的数据，以便后续步骤中恢复i-1行数据
    long length = AA.data[i+1].index-AA.data[i].index+1;
    Double *tempValue = new Double[length+1];
    tempValue[0] = AA.data[i].value;
    long kk = 0;
    for(kk=AA.data[i].index; kk<AA.data[i+1].index; kk++)
    {
     tempValue[kk-AA.data[i].index+1] = AA.data[kk].value;
    }
    long bPos = bb.QuickLocate(i,0);
    if(bPos!=-1)
    {
     tempValue[kk-AA.data[i].index+1] = bb.data[bPos].value;
    }
    else
     tempValue[kk-AA.data[i].index+1] = 0;

    Double gene = AA.data[pos].value/AA.data[i].value*(-1);
    AA.LnTransformMultiply(i,gene);
    AA.LnTransfromAdd(j,i);
    bb.LnTransformMultiply(i, gene);
    bb.LnTransfromAdd(j,i);

    AA.data[i].value = tempValue[0];
    for(kk=AA.data[i].index; kk<AA.data[i+1].index; kk++)
    {
     AA.data[kk].value = tempValue[kk-AA.data[i].index+1];
    }
    if(bPos!=-1)
     bb.data[bPos].value = tempValue[kk-AA.data[i].index+1];

    delete[] tempValue;
   }
  }
  rank++;
 }

 // 判断增广矩阵的秩 和 非增光矩阵的秩是否一致
 for(i=rank; i<bb.total_ln; i++)
 {
  long pos = bb.QuickLocate(i,0);
  if(pos!=-1)
   return 0;
 }

 return 1;
}