Opencv 三对角线矩阵(Tridiagonal Matrix)解法之(Thomas Algorithm)
2024-09-27 03:17:30
1. 简介
三对角线矩阵(Tridiagonal Matrix),结构如公式(1)所示:
aixi−1+bixi+cixx+1=di(1)
其中a1=0,cn=0。写成矩阵形式如(2):
⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢b1a20c1b2a3c2b3⋱⋱⋱cn−1an0bn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢x1x2x3⋮xn⎤⎦⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢d1d2d3⋮dn⎤⎦⎥⎥⎥⎥⎥⎥⎥(2)
常用的解法为Thomas algorithm,又称为The Tridiagonal matrix algorithm(TDMA). 它是一种高斯消元法的解法。分为两个阶段:向前消元(Forward Elimination)和回代(Back Substitution)。
向前消元(Forward Elimination):
c′i=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪cibicibi−aic′i−1;i=1;i=2,3,…,n−1(3)d′i=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪dibidi−aid′i−1bi−aic′i−1;i=1;i=2,3,…,n.(4)回代(Back Substitution):
xn=d′nxi=d′i−c′ixi+1;i=n−1,n−2,…,1.(5)
2.代码
- 维基百科提供的C语言版本:
void solve_tridiagonal_in_place_destructive(float * restrict const x, const size_t X, const float * restrict const a, const float * restrict const b, float * restrict const c)
{
/*
solves Ax = v where A is a tridiagonal matrix consisting of vectors a, b, c
x - initially contains the input vector v, and returns the solution x. indexed from 0 to X - 1 inclusive
X - number of equations (length of vector x)
a - subdiagonal (means it is the diagonal below the main diagonal), indexed from 1 to X - 1 inclusive
b - the main diagonal, indexed from 0 to X - 1 inclusive
c - superdiagonal (means it is the diagonal above the main diagonal), indexed from 0 to X - 2 inclusive
Note: contents of input vector c will be modified, making this a one-time-use function (scratch space can be allocated instead for this purpose to make it reusable)
Note 2: We don't check for diagonal dominance, etc.; this is not guaranteed stable
*/
/* index variable is an unsigned integer of same size as pointer */
size_t ix;
c[0] = c[0] / b[0];
x[0] = x[0] / b[0];
/* loop from 1 to X - 1 inclusive, performing the forward sweep */
for (ix = 1; ix < X; ix++) {
const float m = 1.0f / (b[ix] - a[ix] * c[ix - 1]);
c[ix] = c[ix] * m;
x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
}
/* loop from X - 2 to 0 inclusive (safely testing loop condition for an unsigned integer), to perform the back substitution */
for (ix = X - 1; ix-- > 0; )
x[ix] = x[ix] - c[ix] * x[ix + 1];
}
- 本人基于Opencv的版本:
bool caltridiagonalMatrices(
cv::Mat_<double> &input_a,
cv::Mat_<double> &input_b,
cv::Mat_<double> &input_c,
cv::Mat_<double> &input_d,
cv::Mat_<double> &output_x )
{
/*
solves Ax = v where A is a tridiagonal matrix consisting of vectors input_a, input_b, input_c, and v is a vector consisting of input_d.
input_a - subdiagonal (means it is the diagonal below the main diagonal), indexed from 1 to X - 1 inclusive
input_b - the main diagonal, indexed from 0 to X - 1 inclusive
input_c - superdiagonal (means it is the diagonal above the main diagonal), indexed from 0 to X - 2 inclusive
input_d - the input vector v, indexed from 0 to X - 1 inclusive
output_x - returns the solution x. indexed from 0 to X - 1 inclusive
*/
/* the size of input_a is 1*n or n*1 */
int rows = input_a.rows;
int cols = input_a.cols;
if ( ( rows == 1 && cols > rows ) ||
(cols == 1 && rows > cols ) )
{
const int count = ( rows > cols ? rows : cols ) - 1;
output_x = cv::Mat_<double>::zeros(rows, cols);
cv::Mat_<double> cCopy, dCopy;
input_c.copyTo(cCopy);
input_d.copyTo(dCopy);
if ( input_b(0) != 0 )
{
cCopy(0) /= input_b(0);
dCopy(0) /= input_b(0);
}
else
{
return false;
}
for ( int i=1; i < count; i++ )
{
double temp = input_b(i) - input_a(i) * cCopy(i-1);
if ( temp == 0.0 )
{
return false;
}
cCopy(i) /= temp;
dCopy(i) = ( dCopy(i) - dCopy(i-1)*input_a(i) ) / temp;
}
output_x(count) = dCopy(count);
for ( int i=count-2; i > 0; i-- )
{
output_x(i) = dCopy(i) - cCopy(i)*output_x(i+1);
}
return true;
}
else
{
return false;
}
}
参考文献:https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
最新文章
- C语言中struct位域的定义和使用
- php 读取文件
- 获取B表数据添加到A表中作为一个下拉列表元素存在
- Scalaz(45)- concurrency :Task-函数式多线程编程核心配件
- Maven发布工程到公共库
- mysql锁机制总结
- C语言基础--数组及相关
- Fragemnt
- 提取数据表保存为XML文件
- 给destoon商城的列表中和首页添加购物车功能
- uva 10560 - Minimum Weight(数论)
- C++中实现回调机制的几种方式(一共三种方法,另加三种)
- 【Xamarin挖墙脚系列:Mono项目的图标为啥叫Mono】
- 【转载】Java重构示例【1】
- BZOJ 3211: 花神游历各国( 线段树 )
- CentOS下安装go语言编译环境
- ACM Super Jumping! Jumping! Jumping!
- 能够玩转BKY皮肤的 geek,有一半最后都成为了前端大师
- JavaScript复制文本探究
- php数组和对象转换函数