walnutcy

详细讲解矩阵求逆的快速算法 转载

0
阅读(12808)

矩阵求逆在3D程序中很常见,主要应用于求Billboard矩阵。按照定义的计算方法乘法运算,严重影响了性能。在需要大量Billboard矩阵运算时,矩阵求逆的优化能极大提高性能。这里要介绍的矩阵求逆算法称为全选主元高斯-约旦法。
  高斯-约旦法(全选主元)求逆的步骤如下:
  首先,对于 k 从 0 到 n - 1 作如下几步:
  从第 k 行、第 k 列开始的右下角子阵中选取绝对值最大的元素,并记住次元素所在的行号和列号,在通过行交换和列交换将它交换到主元素位置上。这一步称为全选主元。
  m(k, k) = 1 / m(k, k)
  m(k, j) = m(k, j) * m(k, k),j = 0, 1, ..., n-1;j != k
  m(i, j) = m(i, j) - m(i, k) * m(k, j),i, j = 0, 1, ..., n-1;i, j != k
  m(i, k) = -m(i, k) * m(k, k),i = 0, 1, ..., n-1;i != k
  最后,根据在全选主元过程中所记录的行、列交换的信息进行恢复,恢复的原则如下:在全选主元过程中,先交换的行(列)后进行恢复;原来的行(列)交换用列(行)交换来恢复。
  实现(4阶矩阵)
float Inverse(CLAYMATRIX& mOut, const CLAYMATRIX& rhs)
{
  CLAYMATRIX m(rhs);
  DWORD is[4];
  DWORD js[4];
  float fDet = 1.0f;
  int f = 1;
  for (int k = 0; k
  {
        // 第一步,全选主元
       float fMax = 0.0f;
       for (DWORD i = k; i
       {
           for (DWORD j = k; j
           {
                const float f = Abs(m(i, j));
                if (f fMax)
               {
                     fMax= f;
                     is[k]= i;
                     js[k]= j;
              }
         }
  }
  if (Abs(fMax)   return 0;
  if (is[k] != k)
  {
        f = -f;
        swap(m(k, 0), m(is[k], 0));
        swap(m(k, 1), m(is[k], 1));
        swap(m(k, 2), m(is[k], 2));
        swap(m(k, 3), m(is[k], 3));
  }
  if (js[k] != k)
  {
        f = -f;
        swap(m(0, k), m(0, js[k]));
        swap(m(1, k), m(1, js[k]));
        swap(m(2, k), m(2, js[k]));
        swap(m(3, k), m(3, js[k]));
  }
  // 计算行列值
  fDet *= m(k, k);
  // 计算逆矩阵
  // 第二步
  m(k, k) = 1.0f / m(k, k);
  // 第三步
  for (DWORD j = 0; j
  {
       if (j != k)
            m(k, j) *= m(k, k);
  }
  // 第四步
  for (DWORD i = 0; i
  {
       if (i != k)
      {
           for(j = 0; j
           {
                  if (j != k)
                      m(i, j) = m(i, j) - m(i, k) * m(k, j);
           }
      }
   }
  // 第五步
  for (i = 0; i
  {
       if (i != k)
           m(i, k) *= -m(k, k);
  }
  }
  for(k = 3; k = 0; k --)
  {
      if (js[k] != k)
     {
          swap(m(k, 0), m(js[k], 0));
          swap(m(k, 1), m(js[k], 1));
          swap(m(k, 2), m(js[k], 2));
          swap(m(k, 3), m(js[k], 3));
    }
   if (is[k] != k)
   {
          swap(m(0, k), m(0, is[k]));
          swap(m(1, k), m(1, is[k]));
          swap(m(2, k), m(2, is[k]));
          swap(m(3, k), m(3, is[k]));
   }
  }
  mOut = m;
  return fDet * f;
}