矩陣運算（C++）

寫了個矩陣運算的項目，把算法部分貼出來，大家多多指點啊。 不可多得哦

//矩陣加法

MATRIXDLL_API void Add(float MatrixA[255][255], int mA, int nA, float MatrixB[255][255], int

mB, int nB,float MatrixC[255][255],int& mC,int& nC)

{

 if(mA != mB || nA != nB) return;

 mC = mA;

 nC = nA;

 for(int i = 0; i < mA; i++)

  for(int j = 0; j < nA; j++)

   MatrixC[i][j] = MatrixA[i][j] + MatrixB[i][j];

}

//矩陣乘法

MATRIXDLL_API void Multiply(float MatrixA[255][255], int mA, int nA, float MatrixB[255]

[255], int mB, int nB,float MatrixC[255][255],int& mC,int& nC)

{

 if(nB != mA) return;

 mC = mB;

 nC = nA;

 for(int i = 0; i < nC; i++)

  for (int j = 0; j < mC; j++)

  {

   MatrixC[i][j] = 0;

   for(int p = 0; p < mA; p++) MatrixC[i][j] += MatrixA[i][p] * MatrixB[p][j];

  }

}

//求逆矩陣

MATRIXDLL_API void Inverse(float MatrixA[255][255], int mA, int nA, float MatrixC[255]

[255],int& mC,int& nC)

{

//Calculate the inverse matrix of MatrixA, use MatrixC to save the results. No return.

//Maybe the inverse matrix doesn't exist

 if(mA != nA)

 {

  mC = 0;

  nC = 0;

  return;

 }

 int rank = Rank(MatrixA, mA, nA);

 if (rank != nA)

 {

  mC = 0;

  nC = 0;

  return;

 }

 mC = mA;

 nC = nA;

//Gauss Erosion to upper triangular

 float tmp_MatrixC[255][255];

 int tmp_col = 2 * mA;

 for(int i = 0; i < mA; i++)

  for(int j = 0; j < tmp_col; j++)

  {

   if(j < mA) tmp_MatrixC[j][i] = MatrixA[j][i];

   else if(j - i == mA) tmp_MatrixC[i][j] = 1;

   else tmp_MatrixC[i][j] = 0;

  }

 for(int i = 0; i < nA; i++)

 {

  if (tmp_MatrixC[i][i] == 0)

  {

   int  p;

   for (p = i+1; p < nA; p++)

    if (tmp_MatrixC[p][i] != 0)

    {

     for(int m = 0; m < tmp_col; m++)

     {

      float temp;

      temp = tmp_MatrixC[i]

­­;

      tmp_MatrixC[i]

­ = tmp_MatrixC[p]

­;

      tmp_MatrixC[p] [M]­ = temp;

     }

     break;

    }

  }

  for (int j = i + 1; j < nA; j++)

  {

   float tmp = -(tmp_MatrixC[j][i] / tmp_MatrixC[i][i]);

   for (int q = i; q < tmp_col; q++)

    tmp_MatrixC[j][q] += tmp_MatrixC[i][q] * tmp;

  }

 }

// Change to E(n)

 for (int i = nA - 1; i >= 0; i--)

 {

  for(int j = tmp_col; j >= i; j--) tmp_MatrixC[i][j] /= tmp_MatrixC[i][i];

  for (int j = i-1; j >= 0; j--)

   for(int p = tmp_col; p >= i; p--) tmp_MatrixC[j][p] -= (tmp_MatrixC[j][i] * tmp_MatrixC[/M]

[i][p]);

 }

//Save the result

 for(int i = 0; i < mA; i++)

  for(int  j = 0; j < mA; j++)

   MatrixC[i][j] = tmp_MatrixC[i][j+mA];

}

//高斯消元法

MATRIXDLL_API void GaussErosion(float MatrixA[255][255], int mA, int nA, float MatrixC[255]

[255],int& mC,int& nC)

{

//Achieve GE by your code, conversion MatrixA to a Upper Triangular Matrix, use MatrixC to

save the Upper Triangular Matrix. No return.

 mC = mA;

 nC = nA;

 for(int i = 0; i < mA; i++)

  for(int j = 0; j < nA; j++)

   MatrixC[j][i] = MatrixA[j][i];

 int row = 0, col = 0;

 while(row < nA && col < mA)

 {

  int i = row;

  if (MatrixC[row][col] == 0)

  {

   int  p;

   for (p = row; p < nA; p++)

    if (MatrixC[p][col] != 0)

    {

     for(int m = 0; m < mA; m++)

     {

      double temp;

      temp = static_cast<double>(MatrixC[i]

­­);

      MatrixC[i]

­ = MatrixC[p]

­;

      MatrixC[p] [M]­ = static_cast<float>(temp);

     }

     break;

    }

   if(p == nA)

   {

    col++;

    continue;

   }

  }

  for (int j = row + 1; j < nA; j++)

  {

   double tmp = -(MatrixC[j][i] / MatrixC[i][col]);

   for (int q = col; q < mA; q++)

     MatrixC[j][q] += static_cast<float>(MatrixC[i][q] * tmp);

  }

  row++;

  col++;

 }

}

//求極大無關組

MATRIXDLL_API void MaximalIndependentGroup(float MatrixA[255][255], int mA, int nA, float [/M]

MatrixC[255][255],int& mC,int& nC)

{

//Find the maximal_independent_group of MatrixA, use MatrixC to save the

maximal_independent_group------column vectors.

//Better not to change the origin relative order, No return.

 nC = nA;

 mC = 0;

 double tmp_MatrixC[255][255];

 for(int i = 0; i < mA; i++)

  for(int j = 0; j < nA; j++)

   tmp_MatrixC[j][i] = static_cast<double>(MatrixA[j][i]);

 int row = 0, col = 0;

 while(row < nA && col < mA)

 {

  int i = row;

  if (tmp_MatrixC[row][col] == 0)

  {

   int  p;

   for (p = row; p < nA; p++)

    if (tmp_MatrixC[p][col] != 0)

    {

     for(int m = 0; m < mA; m++)

     {

      double temp;

      temp = tmp_MatrixC[i]

­­;

      tmp_MatrixC[i]

­ = tmp_MatrixC[p]

­;

      tmp_MatrixC[p]

­ = temp;

     }

     break;

    }

    if(p == nA)

    {

     col++;

     continue;

    }

  }

  for (int j = row + 1; j < nA; j++)

  {

   double tmp = -(tmp_MatrixC[j][i] / tmp_MatrixC[i][col]);

   for (int q = col; q < mA; q++)

    tmp_MatrixC[j][q] += tmp_MatrixC[i][q] * tmp;

  }

  row++;

  col++;

 }

 for (int i = 0; i < nA; i++)

 {

  int j;

  for (j = 0; fabs(tmp_MatrixC[i][j]) < MIN; j++)

   if(j == mA-1 && tmp_MatrixC[i][j] == 0) return;

  for(int p = 0; p < nA; p++)

   MatrixC[p][mC] = MatrixA[p][j];

  mC++;

 }

}

//求矩陣的秩

MATRIXDLL_API int Rank(float MatrixA[255][255], int mA, int nA)

{

//Calculate the rank of Matrix A, return the rank value.

 int mC = 0;

 double tmp_MatrixC[255][255];

 for(int i = 0; i < mA; i++)

  for(int j = 0; j < nA; j++)

   tmp_MatrixC[j][i] = static_cast<double>(MatrixA[j][i]);

 int row = 0, col = 0;

 while(row < nA && col < mA)

 {

  int i = row;

  if (tmp_MatrixC[row][col] == 0)

  {

   int  p;

   for (p = row; p < nA; p++)

    if (tmp_MatrixC[p][col] != 0)

    {

     for(int m = 0; m < mA; m++)

     {

      double temp;

      temp = tmp_MatrixC[i]

­;

      tmp_MatrixC[i]

­ = tmp_MatrixC[p]

­;

      tmp_MatrixC[p] [M]­ = temp;

     }

     break;

    }

    if(p == nA)

    {

     col++;

     continue;

    }

  }

  for (int j = row + 1; j < nA; j++)

  {

   double tmp = -(tmp_MatrixC[j][i] / tmp_MatrixC[i][col]);

   for (int q = col; q < mA; q++)

    tmp_MatrixC[j][q] += tmp_MatrixC[i][q] * tmp;

  }

  row++;

  col++;

 }

 for (int i = 0; i < nA; i++)

 {

  int j;

  for (j = 0; j < mA; j++)

   if(fabs(tmp_MatrixC[i][j]) > MIN)

   {

    mC++;

    break;

   }

 }

    return (mC);

}[/M]