ODE数学运算

初始化

将实数数组都填成0或者某个值,其中n是元素个数。

void dSetZero (dReal *a, int n);
void dSetValue (dReal *a, int n, dReal value);

向量点积

返回两个n维向量的点积，如果n<=0返回0

dReal dDot (const dReal *a, const dReal *b, int n);

矩阵乘法

rows\columns	1	2	n
1	1,1	1,2	1,n
2	2,1	2,2	2,n
m	m,1	m,2	m,n

上面描述一个m*n矩阵，ODE内部一个矩阵是一整块线性内存。第一行在开始，然后一段用来对齐。

以4的倍数对齐边界(dPAD(n))，然后是第二行，一直到第m行。这种存储叫做row format。

column format指的是存储方式不同，比如同样上面的m*n矩阵，第一行变成第一列。第二行变成第二列。

其他对齐保存不变，这时候是(dPAD(m))。

- dMultiply0 A = B * C (A:p*r B:p*q C:q*r)

- dMultiply1 A = B’ * C (A:p*r B:p*q C:q*r)

- dMultiply2 A = B * C’ (A:p*r B:p*q C:q*r)

上诉运算中加’的标示一个存储格式是column format的矩阵。

B’就是B的转置矩阵,这样增加两个函数可以减少转置矩阵的操作。这样效率略微提高点。

void dMultiply0 (dReal *A, const dReal *B, const dReal *C, int p,int q,int r);
void dMultiply1 (dReal *A, const dReal *B, const dReal *C, int p,int q,int r);
void dMultiply2 (dReal *A, const dReal *B, const dReal *C, int p,int q,int r);

Cholesky分解

A是一个n*n对称矩阵(symmetric matrix)，A=L*L’。函数通过A计算L。其中L’是L的转置矩阵。

计算的结果通过A返回，其中低三角是L，上三角保持原来的A不变。

函数返回1正确进行分解，0分解失败矩阵不是一个正定矩阵(not positive definite)

int dFactorCholesky (dReal *A, int n);

测试程序用来检测dFactorCholesky的正确性。下面是伪代码：

dMakeRandomMatrix(A,n,n)
B=A*A' (A是n*n矩阵)
dFactorCholesky(B) //这将对B做分解，同时存入B
dClearUpperTriangle(B,n) //将矩阵的上三角清零
B*B' == A*A'

**注意**B和A并不相等。

Cholesky解

解x: A=L*L’ , A*x = b, L是一个n*n矩阵, b是一个n*1矩阵。仅考虑L的低三角部分。

最后将解x写入到b中。

void dSolveCholesky (const dReal *L, dReal *b, int n);

测试代码的伪代码：

dMakeRandomMatrix(A,n,n)
L=A*A'
A=L
dFactorCholesky(L)
dClearUpperTriangle(L)
x=b
dSolveCholesky(L,x)
A*x == b

dInvertPDMatrix

计算n*n正定矩阵A的逆矩阵，将结构放入到Ainv中。如果成功返回1(A确定是一个正定矩阵)

否则返回0(不是一个正定矩阵not Positive Definite Matrix)

int dInvertPDMatrix (const dReal *A, dReal *Ainv, int n);

测试伪代码:

dMakeRandomMatrix(A,n,n)
Ainv=A*A' //确保有逆矩阵
dInvertPDMatrix(A,Ainv,n)
L=A*Ainv
L==I //矩阵乘它的逆矩阵应该等于单位矩阵

检测正定矩阵

检测看A是否是一个正定矩阵，返回1/0(yes/no)。

正定意思是任何x有x’*A*x > 0，实际执行dFactorCholesky分解(decomposition)，如果通过就被认为是正定矩阵。

A的内容不会被改变。

int dIsPositiveDefinite (const dReal *A, int n);

dFactorLDLT

分解一个矩阵A到L*D*L’模式，这里L是一个低三角矩阵并且对角线都为1。

D是对角线。A是一个row格式的n*n矩阵，L被写入到A的低三角区。(对角线1不写入)

D的对角线元素的倒数写入到d中。

void dFactorLDLT (dReal *A, dReal *d, int n, int nskip);

测试伪代码:

dMakeRandomMatrix(B,n,n)
L=B*B' //确保L可以被分解？
A=L
dFactorLDLT(L,d,n,dPAD(n))
dClearUpperTriangle(L,n)
diagonal_one(L)//对角线都设置成1
DL=L/d //一个矩阵DL,每一行有DL[i][j] = L[i][j]/d[j]
//上式，等价于L*D,其中D是一个对角线上是d。
//这里如果D是对角矩阵,等式成立DL' = (DL)'
A==L*DL'

dSolveL1

解方程L*x=b,这里L是一个对角线都为1的n*n低三角矩阵，x,b是一个n*1向量。

x解出来后覆盖b。

void dSolveL1 (const dReal *L, dReal *b, int n, int nskip);

dSolveL1T

解L’*x=b,这里L是一个对角线都为1的n*n低三角矩阵，x,b是一个n*1向量。

x解出来后覆盖b。

void dSolveL1T (const dReal *L, dReal *b, int n, int nskip);

dVectorScale

a[i] *= d[i]

void dVectorScale (dReal *a, const dReal *d, int n);

dSolveLDLT

给出L’，L是一个对角线都为1的n*n低三角矩阵，d’ 是一个n*1向量，它是矩阵D的对角线

向量的倒数。就是将d’的倒数填入一个D的对角线。

解L*D*L’*x=b， x,b是一个n*1向量。

x解出来后覆盖b。

void dSolveLDLT (const dReal *L, const dReal *d, dReal *b, int n, int nskip);

测试伪代码：

dMakeRandomMatrix(B,n,n) //创建一个n*n随机矩阵
A=B*B' //确保L可以被分解？
L=A
dMakeRandomMatrix(b,n,1)
dFactorLDLT(L,d)
dSolveLDLT(L,d,x)
A*x == d

dLDLTAddTL

给出一个n*n矩阵A的L*D*L’分解，返回跟新的分解L2*D2*L2’

L是一个n*n对角线为1的低三角矩阵，仅仅L的低三角部分被使用。

d是对角矩阵D,对角上的元素的倒数组成的向量n*1。

结果被写入到L,除非L的左行和d[0]没有被修改。

void dLDLTAddTL (dReal *L, dReal *d, const dReal *a, int n, int nskip);

测试伪代码：

dMakeRandomMatrix(B,n,n)
A=B*B'
L=A
dFactorLDLT(L,d)
a[i]=-A[i*nskip] //i=0;i<n,将A的第一列的负值赋予a向量
a[0]+=1
dLDLTAddTL(L,d,a)
L[i*nskip]=0 //i=1;i<nL的第一列去除第一元素等于0
d[0]=1
dClearUpperTriangle(L)
diagonal_one(L) //将L的对角线置为1
dSetZero(DL)
DL[i*n+j]=L[i*n+j]/d[j]

A[i*n]=A[i]=0 //i=1;i<n
A[0] = 1;
A == L*DL'

dLDLTRemove

void dLDLTRemove (dReal **A, const int *p, dReal *L, dReal *d,
int n1, int n2, int r, int nskip);

dRemoveRowCol

void dRemoveRowCol (dReal *A, int n, int nskip, int r);

dPAD

取一个4的倍数的边界，如dPAD(2)等于4，dPAD(3)等于4，dPAD(4)等于4，dPAD(5)等于8

#define dPAD(a) (((a) > 1) ? ((((a)-1)|3)+1) : (a))