RTKLIB学习笔记——rtkcmn.c（矩阵运算二）

 LAPACK

#ifdef LAPACK /* with LAPACK/BLAS or MKL */

LAPACK包含了求解科学与工程计算中最常见的数值线性代数问题，如求解线性方程组、线性最小二乘问题、特征值问题和奇异值问题等。

• 乘法运算

extern void matmul(const char *tr, int n, int k, int m, double alpha,
                   const double *A, const double *B, double beta, double *C)
{
    int lda=tr[0]=='T'?m:n,ldb=tr[1]=='T'?k:m;
    
    dgemm_((char *)tr,(char *)tr+1,&n,&k,&m,&alpha,(double *)A,&lda,(double *)B,
           &ldb,&beta,C,&n);
}

C(n*k)=alpha*A(n*m)B(m*k) +beta*C

*tr : N 正常 T转置

int lda=tr[0]=='T'?m:n,  对于A  取行数n    转置 取列数m

ldb=tr[1]=='T'?k:m;       对于B   取行数m   转置 去列数k

dgemm_为LAPACK包含的函数

• 逆矩阵

extern int matinv(double *A, int n)
{
    double *work;
    int info,lwork=n*16,*ipiv=imat(n,1);
    
    work=mat(lwork,1);
    dgetrf_(&n,&n,A,&n,ipiv,&info);
    if (!info) dgetri_(&n,A,&n,ipiv,work,&lwork,&info);
    free(ipiv); free(work);
    return info;
}

 

• 线性方程求解

extern int solve(const char *tr, const double *A, const double *Y, int n,
                 int m, double *X)
{
    double *B=mat(n,n);
    int info,*ipiv=imat(n,1);
    
    matcpy(B,A,n,n);
    matcpy(X,Y,n,m);
    dgetrf_(&n,&n,B,&n,ipiv,&info);
    if (!info) dgetrs_((char *)tr,&n,&m,B,&n,ipiv,X,&n,&info);
    free(ipiv); free(B); 
    return info;
}

ELSE  无LAPACK

• 矩阵乘法

extern void matmul(const char *tr, int n, int k, int m, double alpha,
                   const double *A, const double *B, double beta, double *C)

 

C(n*k)=alpha*A(n*m)B(m*k) +beta*C

*tr : N 正常 T转置

tr[0]	tr[1]	f	d
N	N	1	d+=A[i+x*n]*B[x+j*m]
N	F	2	d+=A[i+x*n]*B[x+j*k]
F	N	3	d+=A[i+x*m]*B[x+j*m]
F	F	4	d+=A[i+x*m]*B[x+j*k]

• LU分解

static int ludcmp(double *A, int n, int *indx, double *d)

LU decomposition 

• LU回代

static void lubksb(const double *A, int n, const int *indx, double *b)

LU back-substitution 

• 逆矩阵

extern int matinv(double *A, int n)

inverse of matrix 

• 线性方程求解

extern int solve(const char *tr, const double *A, const double *Y, int n,
                 int m, double *X)

X=A\Y  X=A'\Y 

A(n*n)  Y(n*m)  输出X