OpenCV中LU分解实现矩阵求逆invert(DECOMP_LU)-定点化

布满荆棘的人生 2022-08-21 14:51 160阅读 0赞

基于LU分解的矩阵求逆定点化版本，由于需要频繁移位，因此定点比浮点还耗时

浮点版本参考上一篇博客：[http://blog.csdn.net/xidianzhimeng/article/details/51284206][http_blog.csdn.net_xidianzhimeng_article_details_51284206]

话不过多，直接上代码，有不明白的地方留言

typedef char S8;  
typedef short S16;  
typedef int S32;  
  
  
typedef unsigned char U8;  
typedef unsigned short U16;  
typedef unsigned int U32;  
  
  
typedef int CPLX16;  
typedef float F32;  
  
  
typedef long long S64;  
typedef unsigned long long U64;  
  
  
typedef long long S40;  
typedef unsigned long U40;  
  
  
typedef int BOOL;

\#define INVERT\_EXPAND1 10  
\#define INVERT\_EXPAND1\_N (1<<INVERT\_EXPAND1)  
\#define INVERT\_EXPAND2 20  
\#define INVERT\_EXPAND2\_N (1<<INVERT\_EXPAND2)  
\#define INVERT\_EXPAND3 30  
\#define INVERT\_EXPAND3\_N (1<<INVERT\_EXPAND3)  
\#define INVERT\_EXPAND4 40  
  
  
\#define MAX\_SD\_ROW 200  
\#define MAX\_SD\_COL 200

S64 matrixT\[MAX\_SD\_ROW\]\[MAX\_SD\_COL\];  
S64 matrixB\[MAX\_SD\_ROW\]\[MAX\_SD\_COL\];

// A(symmetry) : diagonal element is 18+8=26bit , other element is 18bit  
BOOL MatrixInvLUFixed(S32 A\[MAX\_SD\_ROW\]\[MAX\_SD\_COL\], const U32 n)  
\{  
S32 ii;  
U32 i, j, k;  
F32 d;  
S64 dd, s, alpha;  
  
  
for (i = 0; i < n; i++)  
\{

// 由于已知输入为对称矩阵并且对角线上元素最大，所以省去选择主元  
if (\_abs(matrixT\[i\]\[i\] >> 32) == 0 && \_abs(matrixT\[i\]\[i\] & 0xffffffff) == 0)  
return FALSE;  
  
d = \_rcpsp((F32)matrixT\[i\]\[i\]);  
dd = (S64)(d \* INVERT\_EXPAND3\_N);  
  
for (j = i + 1; j < n; j++)  
\{  
alpha = matrixT\[j\]\[i\] \* dd;  
for (k = i + 1; k < n; k++)  
\{  
matrixT\[j\]\[k\] = matrixT\[j\]\[k\] - (alpha\*matrixT\[i\]\[k\] >> INVERT\_EXPAND3);  
matrixB\[j\]\[k\] = matrixB\[j\]\[k\] - (alpha\*matrixB\[i\]\[k\] >> INVERT\_EXPAND3);  
\}  
for (k = 0; k <= i; k++)  
matrixB\[j\]\[k\] = matrixB\[j\]\[k\] - (alpha\*matrixB\[i\]\[k\] >> INVERT\_EXPAND3);  
\}  
  
  
matrixT\[i\]\[i\] = dd;  
\}  
  
  
// A occupy 10bit  
for (ii = n - 1; ii >= 0; ii--)  
for (j = 0; j < n; j++)  
\{  
s = matrixB\[ii\]\[j\] << INVERT\_EXPAND3;  
for (k = ii + 1; k < n; k++)  
s -= (matrixT\[ii\]\[k\] \* matrixB\[k\]\[j\]);  
matrixB\[ii\]\[j\] = (s >> INVERT\_EXPAND3) \* matrixT\[ii\]\[ii\];  
A\[ii\]\[j\] = (S32)(matrixB\[ii\]\[j\] >> INVERT\_EXPAND2);  
\}  
  
  
return TRUE;  
\}

[http_blog.csdn.net_xidianzhimeng_article_details_51284206]: http://blog.csdn.net/xidianzhimeng/article/details/51284206