精华内容
下载资源
问答
  • C语言矩阵运算

    2019-05-02 22:19:55
    /*运行环境codeBlocks*/ /*C语言*/ #include #include #define M 3 #define N 3 #define m 3 #define n 3 #define P 3 #define Q 2
  • c语言矩阵运算程序

    2012-09-14 22:37:22
    c语言矩阵运算程序,简单的加减乘除运算,用c语言实现
  • C语言矩阵运算库大起底

    千次阅读 2016-11-30 17:11:04
    C语言矩阵运算

    GSL

    GNU Scientific Library自带的矩阵运算,据说速度一般。

    Blitz++

    Blitz++ 与 MTL 都是基于 C++ template 高效数值计算程序库,不过他们专注于不同的方向。
    Blitz++ 提供了一个 N 维( 1—10 )的 Array 类 , 这个 Array 类以 reference counting 技术实现,支持任意的存储序 (row-major 的 C-style 数组, column-major 的 Fortran-style 数组 ) ,数组的切割 (slicing), 子数组的提取 (subarray), 灵活的 Array 相关表达式处理。另外提供了可以产生不同分布的随机数 (F,Beta,Chi-Square ,正态,均匀分布等 ) 的类也是很有特色的。
    MTL 专注于线性代数相关的计算任务,如各种形式矩阵的生成 ( 对角,共轭,稀疏,对称等 ) ,相关的计算,变换,以及与一维向量的运算。
    两个程序库对于从 Matlab 导入导出数据都有不错的支持。

    MTL4

    类似Eigen和Armadillo,有开源版本。

    Eigen

    http://eigen.tuxfamily.org/dox/
    Eigen[1] 目前最新的版本是3.2.5,除了C++标准库以外,不需要任何其他的依赖包。Eigen使用的CMake建立配置文件和单元测试,并自动安装。如果使用Eigen库,只需包特定模块的的头文件即可。最新版本为Eigen 3.3。
    使用类似Matlab的方式操作矩阵,可以在这里查看官方的与Maltab的对应关系,个人感觉单纯讲和Matlab的对应的话,可能不如Armadillo对应的好,但功能绝对强大。Eigen包含了绝大部分你能用到的矩阵算法,同时提供许多第三方的接口。Eigen一个重要特点是没有什么依赖的库,本身仅有许多头文件组成,因此非常轻量而易于跨平台。你要做的就是把用到的头文件和你的代码放在一起就可以了。Eigen的一些特性:
    •支持整数、浮点数、复数,使用模板编程,可以为特殊的数据结构提供矩阵操作。比如在用ceres-solver进行做优化问题(比如bundle adjustment)的时候,有时候需要用模板编程写一个目标函数,ceres可以将模板自动替换为内部的一个可以自动求微分的特殊的double类型。而如果要在这个模板函数中进行矩阵计算,使用Eigen就会非常方便。
    •支持逐元素、分块、和整体的矩阵操作。
    •内含大量矩阵分解算法包括LU,LDLt,QR、SVD等等。
    •支持使用Intel MKL加速
    •部分功能支持多线程
    •稀疏矩阵支持良好,到今年新出的Eigen3.2,已经自带了SparseLU、SparseQR、共轭梯度(ConjugateGradient solver)、bi conjugate gradient stabilized solver等解稀疏矩阵的功能。同时提供SPQR、UmfPack等外部稀疏矩阵库的接口。
    •支持常用几何运算,包括旋转矩阵、四元数、矩阵变换、AngleAxis(欧拉角与Rodrigues变换)等等。
    •更新活跃,用户众多(Google、WilliowGarage也在用),使用Eigen的比较著名的开源项目有ROS(机器人操作系统)、PCL(点云处理库)、Google Ceres(优化算法)。OpenCV自带到Eigen的接口。

    Armadillo

    矩阵运算速度跟matlab一个量级
    目前使用比较广的C++矩阵运算库之一,是在C++下使用Matlab方式操作矩阵很好的选择,许多Matlab的矩阵操作函数都可以找到对应,这对习惯了Matlab的人来说实在是非常方便,另外如果要将Matlab下做研究的代码改写成C++,使用Armadillo也会很方便,这里有一个简易的Matlab到Armadillo的语法转换。下面列了一些Armadillo的特性:
    •支持整数,浮点数,和复数矩阵。
    •支持矩阵逐元素操作,包括abs · conj · conv_to · eps · imag/real · misc functions (exp, log, pow, sqrt, round, sign, …) · trigonometric functions (cos, sin, …)等等。
    •支持矩阵分块操作。
    •支持对整体矩阵的操作diagvec · min/max · prod · sum · statistics (mean, stddev, …) · accu · as_scalar · det · dot/cdot/norm_dot · log_det · norm · rank · trace等等。
    •Matlab用户,你甚至可以找到你熟悉的hist · histc · unique · cumsum · sort_index · find · repmat · linspace等函数。
    •除了自带的矩阵基本运算之外,可自动检测是否安装有BLAS,或更快的 OpenBLAS, Intel MKL, AMD ACML,并使用他们替代自带基本运算实现。
    •提供接口使用LAPACK进行矩阵分解运算,svd · qr · lu · fft等等。
    •提供了稀疏矩阵类,支持常用操作,但暂时没有矩阵分解的实现。
    •更新比较活跃,有一些计算机视觉、机器学习、物理方面的开源项目在使用,比如MLPACK (Machine Learning Library)。
    总体来讲很好用的矩阵库,速度上因为可以使用OpenBLAS等库进行加速,因此还是不错的

    Lapack (Linear Algebra PACKage)

    http://www.netlib.org/lapack/
    LAPACK,其名为Linear Algebra PACKage的缩写,是一以Fortran编程语言写就,用于数值计算的函式集。 LAPACK提供了丰富的工具函式,可用于诸如解多元线性方程式、线性系统方程组的最小平方解、计算特征向量、用于计算矩阵QR分解的Householder转换、以及奇异值分解等问题

    Pardiso

    http://www.pardiso-project.org/
    The package PARDISO is a thread-safe, high-performance, robust, memory efficient and easy to use software for solving large sparse symmetric and unsymmetric linear systems of equations on shared-memory and distributed-memory multiprocessors.

    Openblas

    是一个高性能多核 BLAS 库,是 GotoBLAS2 1.13 BSD 版本的衍生版。OpenBLAS 的编译依赖系统环境,并且没有原生单线程版本,在实验这哦那个,通过设置 OMP_NUM_THREADS=1 来模拟单线程版本,可能会带来一点点的性能下降。

    Intel mkl

    https://software.intel.com/en-us/intel-mkl
    Fastest and most used math library for Intel and compatible processors**
    • Features highly optimized, threaded and vectorized functions to maximize performance on each processor family
    • Utilizes de facto standard C and Fortran APIs for compatibility with BLAS, LAPACK and FFTW functions from other math libraries
    • Available with both free community-supported and paid support licenses

    来自 https://software.intel.com/en-us/intel-mkl

    Trilinos

    和PETSc同是美国能源部ODE2000支持开发的20多个ACTS工具箱之一,用于大规模计算。

    Blas(Basic Linear Algebra Subprograms)
    http://www.netlib.org/blas/

    参考意见:
    ○ 在windows上Eigen的使用或许比Armadillo方便些
    ○ Intel mkl 强大,速度快
    http://nghiaho.com/?p=954 (对比测试)
    http://www.leexiang.com/the-performance-of-matrix-multiplication-among-openblas-intel-mkl-and-eigen (对比测试)
    http://www.oschina.net/project/tag/239/Mathematics-computin?sort=view&lang=21&os=0 (其他软件合集)

    展开全文
  • 超好用的纯C语言矩阵运算

    千次阅读 2019-06-16 18:57:51
    超好用的纯C语言矩阵运算库 最近开发基于异构传感器的异步定位数据的卡尔曼滤波,因为最终要用在一个老爷DSP上,已有代码都是C写的,不想研究这种老爷设备下的C++调用,Eigen这超好用的矩阵运算库算是没法用了。随便...

    超好用的纯C语言矩阵运算库 easyMatrix


    最近开发基于异构传感器的异步定位数据的卡尔曼滤波,因为最终要用在一个老爷DSP上,已有代码都是C写的,不想研究这种老爷设备下的C++调用,Eigen这超好用的矩阵运算库算是没法用了。随便搜了几个C矩阵库发现不是超难用就是有bug,比如csdn某小哥写的matrix库,矩阵维数一大,逆矩阵计算就出错。项目结束后,想着可能很多嵌入式开发人员也许有类似的需求,决定自己还是写一个操作友好的matrix库,工具虽小,好用就行。
    支持矩阵的加、减、乘法运算,转置、逆矩阵、伴随矩阵、LU分解、行列式计算、线性方程求解等功能
    废话不说,上代码:
    github地址:https://github.com/fellylanma/easyMatrix

    //使用例子
    int main() {
    	//create matrix on stack
     	float val1[] = {1,2,3,4,5,6,7,8,9};
    	float val2[] = {1,1,1,1,1,1,1,1,1};
    	CREATE_MATRIX_ONSTACK(3,3,A,val1);
    	CREATE_MATRIX_ONSTACK(3,3,B,val2);
    	CREATE_MATRIX_ONSTACK(3,3,C,NULL);
    	addMatrix(&A,&B,&C);
    	dumpMatrix(&C);
        printf("det is:%f\n",detMatrix(&A));
    	int x = A.rows;
    	int y = A.cols;
    	printf("element[1,1] is:%f\n",A.element[1*x+1]);
    	//create matrix on heap and you need to delete the pointer
    	CREATE_DYNAMIC_MATRIX_ONHEAP(x,y,D,val1);
    	CREATE_DYNAMIC_MATRIX_ONHEAP(x,y,E,val2);
    	CREATE_DYNAMIC_MATRIX_ONHEAP(x,y,F,NULL);
    	multiMatrix(D,E,&C);
    	multiMatrix(D,E,F);
    
    	DELETE_DYNAMIC_MATRIX(D);
    	DELETE_DYNAMIC_MATRIX(E);
    	DELETE_DYNAMIC_MATRIX(F);
    }
    

    easyMatrix.h

    #ifndef _MAGRIDE_PLANNING_EASYMATRIX_H_
    #define _MAGRIDE_PLANNING_EASYMATRIX_H_
    
    #include <stdlib.h>
    typedef unsigned char uint8;
    typedef float DATA_TYPE;
    
    
    #define CREATE_MATRIX_ONSTACK(x,y,matrix,initval) \
    struct easyMatrix matrix;\
    DATA_TYPE val##x##N##y##N##matrix[x*y];\
        matrix.rows = x;\
        matrix.cols = y;\
        matrix.element = val##x##N##y##N##matrix;\
        if(initval!=NULL) setMatrix(initval, &(matrix))
    
    #define CREATE_DYNAMIC_MATRIX_ONHEAP(x,y,matrix,initval) \
    struct easyMatrix *matrix = (struct easyMatrix*)malloc(sizeof(struct easyMatrix));\
    matrix->rows = x;\
    matrix->cols = y;\
    matrix->element = (DATA_TYPE*) malloc(sizeof(DATA_TYPE)*(x)*(y));\
    if(initval!=NULL) setMatrix(initval, (matrix))
    
    #define DELETE_DYNAMIC_MATRIX(matrix) \
        free((matrix)->element);\
        free(matrix)
    
    struct easyMatrix {\
        uint8 rows,cols;\
        DATA_TYPE* element;
    };\
    
    struct easyMatrix* setMatrix(DATA_TYPE* const a,struct easyMatrix* c);
    struct easyMatrix* copyMatrix(struct easyMatrix* const a,struct easyMatrix* c);
    struct easyMatrix* transMatrix(struct easyMatrix* const a,struct easyMatrix* c);
    
    DATA_TYPE detMatrix(struct easyMatrix* const a);
    
    DATA_TYPE invMatrix(struct easyMatrix* const a, struct easyMatrix*b);
    
    struct easyMatrix* scaleMatrix(DATA_TYPE, struct easyMatrix* const a, struct easyMatrix*);
    
    struct easyMatrix* addMatrix(const struct easyMatrix* const a, const struct easyMatrix *const  b, struct easyMatrix * c);
    
    struct easyMatrix* leftMatrix(uint8, uint8, struct easyMatrix* const a, struct easyMatrix* b);
    
    struct easyMatrix* subMatrix(struct easyMatrix* const a, struct easyMatrix* const  b, struct easyMatrix* c);
    
    struct easyMatrix* multiMatrix(struct easyMatrix* const a, struct easyMatrix* const b, struct easyMatrix* c);
    
    struct easyMatrix* zerosMatrix(struct easyMatrix* e);
    
    struct easyMatrix* eyesMatrix(struct easyMatrix* e);
    
    void dumpMatrix(struct easyMatrix* const e);
    
    struct easyMatrix* adjMatrix(struct easyMatrix* const a,struct easyMatrix* c);
    
    struct easyMatrix* getLUMatrix(struct easyMatrix* const A, struct easyMatrix* L,struct easyMatrix* U) ;
    
    struct easyMatrix* invLMatrix(struct easyMatrix* const L, struct easyMatrix* L_inv) ;
    struct easyMatrix* invUMatrix(struct easyMatrix* const U, struct easyMatrix* U_inv) ;
    
    struct easyMatrix* solveEquationMatrix(const struct easyMatrix* const A,const struct easyMatrix* const Y, struct easyMatrix* X) ;
    
    
    DATA_TYPE fastDetMatrix(struct easyMatrix* const in) ;
    #endif//_MAGRIDE_PLANNING_EASYMATRIX_H_
    

    #easyMatrix.c

    #include <stdlib.h>
    #include <float.h>
    #include "easyMatrix.h"
    
    int isFiniteNumber(double d) {
        return (d<=DBL_MAX&&d>=-DBL_MAX);
    }
    struct easyMatrix* setMatrix(DATA_TYPE * const a,struct easyMatrix* c) {
        uint8 x = c->rows;
        uint8 y = c->cols;
        int t = x*y;
        for(int i=0;i<t;++i) {
            c->element[i] = a[i];
        }
        return c;
    }
    
    struct easyMatrix* copyMatrix(struct easyMatrix* const a,struct easyMatrix* c) {
        if(a->rows != c->rows) return NULL;
        if(a->cols != c->cols) return NULL;
        int t = a->rows*a->cols;
        for(int i=0;i<t;++i) {
            c->element[i] = a->element[i];
        }
        return c;
    }
    
    struct easyMatrix* transMatrix(struct easyMatrix* const a,struct easyMatrix* c) {
        if(a->rows != c->cols) return NULL;
        if(a->cols != c->rows) return NULL;
        int index = 0;
        int index_src = 0;
        for(uint8 ii=0;ii<a->cols;++ii) {
            index_src=ii;
            for(uint8 jj=0;jj<a->rows;++jj) {
                //c->element[index] = a->element[jj*a->cols+ii];
                c->element[index] = a->element[index_src];
                index++;
                index_src+=a->cols;
            }
        }
        return c;
    }
    
    struct easyMatrix* leftMatrix(uint8 x_i,uint8 y_i, struct easyMatrix* const in, struct easyMatrix* out) {
        if(in->rows != in->cols) return NULL;
        if(out->rows != out->cols) return NULL;
        if(in->rows != (out->rows+1)) return NULL;
        int index = 0;
        int index_src = 0;
        uint8 x =in->rows;
        uint8 y =in->cols;
        for(uint8 kk=0;kk<x;++kk) {
            for(uint8 ww=0;ww<y;++ww) {
                if(!(kk==x_i||ww==y_i)) {
                    //out->element[index] = in->element[kk*y+ww];
                    out->element[index] = in->element[index_src];
                    index++;
                }
                index_src++;
            }
        }
        return out;
    }
    struct easyMatrix* adjMatrix(struct easyMatrix* const in, struct easyMatrix* out) {
        if(in->rows != out->cols) return NULL;
        if(in->cols != out->rows) return NULL;
        int index = 0;
        uint8 x = in->rows;
        uint8 y = in->cols;
        CREATE_DYNAMIC_MATRIX_ONHEAP(x-1,y-1,ret,NULL);
        signed char sign1 = 1;
        signed char sign2 = 1;
        for(uint8 ii=0;ii<x;++ii) {
            sign2 = sign1;
            index = ii;
            for(uint8 jj=0;jj<y;++jj) {
                leftMatrix(ii,jj,in,ret);
                //out->element[jj*y+ii] = sign2*detMatrix(ret);
                out->element[index] = sign2*detMatrix(ret);
                sign2 = - sign2;    
                index+=y;
            }
            
            sign1 = - sign1;
        }
        DELETE_DYNAMIC_MATRIX(ret);
        return out;
    }
    
    DATA_TYPE invMatrix(struct easyMatrix *const in , struct easyMatrix * out) {
        if(in->cols!=in->rows) return 0;
        if(in->rows != out->cols) return 0;
        if(in->cols != out->rows) return 0;
        uint8 N = in->cols;
        CREATE_DYNAMIC_MATRIX_ONHEAP(N,N,L,NULL);
        CREATE_DYNAMIC_MATRIX_ONHEAP(N,N,LINV,NULL);
        CREATE_DYNAMIC_MATRIX_ONHEAP(N,N,U,NULL);
        CREATE_DYNAMIC_MATRIX_ONHEAP(N,N,UINV,NULL);
        getLUMatrix(in,L,U);
        invLMatrix(L,LINV);
        invUMatrix(U,UINV);
        multiMatrix(UINV,LINV,out);
        double s = 1;
        for(int i = 0;i<N;i++) 
            s *= U->element[i*N+i];   
        /*
        adjMatrix(in,out);
        DATA_TYPE scale = detMatrix(in);
        if(scale<1e-5&&scale>-1e-5) return 0.0;
        scale = 1/scale;
        scaleMatrix(scale,out,out);
    */
        DELETE_DYNAMIC_MATRIX(L);
        DELETE_DYNAMIC_MATRIX(U);
        DELETE_DYNAMIC_MATRIX(LINV);
        DELETE_DYNAMIC_MATRIX(UINV);
        return isFiniteNumber(s);
    }
    
    struct easyMatrix* getLUMatrix(struct easyMatrix* const A, struct easyMatrix* L,struct easyMatrix* U) {
        int row=0;
        DATA_TYPE s = 0;
        uint8 N = A->cols;
        int t = N*N;
        for(int i =0;i<t;i++) {
            L->element[i] = 1e-20;
            U->element[i] = 1e-20;
        }
        for(int i=0;i<N;i++) {
            L->element[i*N+i] = 1.0;
        }
        for(int i=0;i<N;i++) {
            for(int j=i;j<N;j++) {
                s = 0.0;
                for(int k=0;k<i;++k) {
                    s+=L->element[i*N+k]*U->element[k*N+j];
                }
                U->element[i*N+j]= A->element[i*N+j] - s; 
            }
            for (int j = i + 1;j < N;j++) {
                s = 0.0;
                for (int k = 0; k < i; k++)
                {
                    s += L->element[j*N+k] * U->element[k*N+i];
                }
                L->element[j*N+i] = (A->element[j*N+i] - s) / U->element[i*N+i];      //按列计算l值
            }
        }
        return L;
    
    }
    
    struct easyMatrix* invLMatrix(struct easyMatrix* const L, struct easyMatrix* L_inv) { 
        uint8 N = L->cols;
        DATA_TYPE s;
        int t = N*N;
        for(int i =0;i<t;i++) {
            L_inv->element[i] = 1e-13;
        }
        for (uint8 i = 0;i < N;i++)  {
            L_inv->element[i*N+i] = 1;
        }
        for (uint8 i= 1;i < N;i++) {
            for (uint8 j = 0;j < i;j++) {
                s = 0;
                for (uint8 k = 0;k < i;k++) {
                    s += L->element[i*N+k] * L_inv->element[k*N+j];
                }
                L_inv->element[i*N+j] = -s;
            }
        }
        return L_inv;
    }
    struct easyMatrix* invUMatrix(struct easyMatrix* const U, struct easyMatrix* U_inv) { 
        uint8 N = U->cols;
        DATA_TYPE s;
        int t = N*N;
        for(int i =0;i<t;i++) {
            U_inv->element[i] = 1e-13;
        }
     for (uint8 i = 0;i < N;i++)                    //按列序,列内按照从下到上,计算u的逆矩阵
        {
            U_inv->element[i*N+i] = 1 / U->element[i*N+i];
        }
        for (uint8 i = 1;i < N;i++) {
            for (int j = i - 1;j >=0;j--) {
                s = 0;
                for (uint8 k = j + 1;k <= i;k++) {
                    s += U->element[j*N+k] * U_inv->element[k*N+i];
                }
                U_inv->element[j*N+i] = -s / U->element[j*N+j];
            }
        }
        return U_inv;
    }
    DATA_TYPE fastDetMatrix(struct easyMatrix* const in) {
        uint8 x = in->rows;
        uint8 y = in->cols;
        if(x!=y) return 0;
        if(x==0 ) return 0;
        if(x==1 ) return in->element[0];
        DATA_TYPE *a =in->element;
        if(x==2) return(a[0]*a[3]-a[1]*a[2]);
        int N = x;
        CREATE_DYNAMIC_MATRIX_ONHEAP(N,N,L,NULL);
        CREATE_DYNAMIC_MATRIX_ONHEAP(N,N,U,NULL);
        getLUMatrix(in,L,U);
        double s = 1;
        for(int i = 0;i<N;i++) 
            s *= U->element[i*N+i];
        DELETE_DYNAMIC_MATRIX(L);
        DELETE_DYNAMIC_MATRIX(U);
        return s;
    }
    
    DATA_TYPE detMatrix(struct easyMatrix* const in) {
        uint8 x = in->rows;
        uint8 y = in->cols;
        if(x!=y) return 0;
        if(x==0 ) return 0;
        if(x==1 ) return in->element[0];
        DATA_TYPE *a =in->element;
        if(x==2) return(a[0]*a[3]-a[1]*a[2]);
    
        DATA_TYPE result = 0;
        signed char sign = 1;
        CREATE_DYNAMIC_MATRIX_ONHEAP(x-1,y-1,ret,NULL);
        for(uint8 i=0;i<x;++i) {
            leftMatrix(0,i,in,ret);
            result += sign*a[i]*detMatrix(ret);
            sign = - sign;
        }
        DELETE_DYNAMIC_MATRIX(ret);
        return result;
    }
    
    struct easyMatrix* addMatrix(const struct easyMatrix* const a,const struct easyMatrix* const b, struct easyMatrix* c) {
        if(a->cols != b->cols) return NULL;
        if(a->rows != b->rows) return NULL;
        struct easyMatrix* obj = (struct easyMatrix*)a;
        int t = obj->rows*obj->cols;
        for(int i=0;i<t;++i) {
            c->element[i] = obj->element[i]+b->element[i];
        }
        return c;
    }
    
    struct easyMatrix* subMatrix(struct easyMatrix* const a, struct easyMatrix* const b, struct easyMatrix* c) {
        if(a->cols != b->cols) return NULL;
        if(a->rows != b->rows) return NULL;
        struct easyMatrix* obj = (struct easyMatrix*)a;
        int t = obj->rows*obj->cols;
        for(int i=0;i<t;++i) {
            c->element[i] = obj->element[i]-b->element[i];
        }
        return c;
    }
    
    struct easyMatrix* scaleMatrix(DATA_TYPE scale, struct easyMatrix* const a, struct easyMatrix* b) {
        int t = a->cols*a->rows;
        for (int i = 0;i<t;++i) {
            b->element[i] = a->element[i]*scale;
        }
        return b;
    }
    
    struct easyMatrix* multiMatrix(struct easyMatrix* const a,struct easyMatrix* const b, struct easyMatrix* c) {
        if(NULL==c) return NULL;
        if(c == a || c == b) return NULL;
        if(a->cols != b->rows) return NULL;
        int count = 0;
        int t_cnt = 0;
        int z_cnt = 0;
        uint8 x = a->rows;
        uint8 y = a->cols;
        uint8 z = b->cols;
        for(uint8 i = 0;i<x;++i) {
            for(uint8 k = 0;k<z;++k) {
                c->element[count] = 0;
                z_cnt = 0;
                for(uint8 j = 0;j<y;++j) {
                    c->element[count] += a->element[t_cnt+j]*b->element[z_cnt+k];
                    z_cnt += z;
                }
                count++;
            }
            t_cnt+=y;
        }
        return c;
    }
    
    struct easyMatrix* zerosMatrix(struct easyMatrix* e) {
        int t = e->cols*e->rows;
        for(int i=0;i<t;++i) {
            e->element[i] = 0;
        }
        return e;
    }
    
    struct easyMatrix* eyesMatrix(struct easyMatrix* e) {
        if(e->rows != e->cols) return NULL;
        zerosMatrix(e);
        int index = 0;
        for(uint8 i=0;i<e->rows;++i) {
            e->element[index] = 1.0;
            index+=(e->cols);
            ++index;
        }
        return e;
    }
    
    struct easyMatrix* solveEquationMatrix(const struct easyMatrix* const A,const struct easyMatrix* const Y, struct easyMatrix* X) { 
        CREATE_DYNAMIC_MATRIX_ONHEAP(A->rows,A->cols,AINV,NULL);
        invMatrix(A,AINV);
        multiMatrix(AINV,Y,X);
        DELETE_DYNAMIC_MATRIX(AINV);
        return X;
    }
    void dumpMatrix(struct easyMatrix* const e) {
        int count = 0;
        int x = e->rows;
        int y = e->cols;
        printf("cols is:%d, rows is:%d\n",x,y);
        for(uint8 i = 0;i<x;++i) {
            for(uint8 j = 0;j<y;++j) {
                printf("%8f,",e->element[count]);
                ++count;
            }
            printf("\n");
        }
        return;
    }
    
    展开全文
  • 特别说明 此资料来自豆丁网) 您现在所看到的文档是使用下载器所生成的文档 此文档的原件位于 感谢您的支持 抱米花 /p-60444020.html
  • 工作学习上都需要用到C语言裸机下运算矩阵,找了一些库不是很理想; 浏览了上述博客中博主的程序,发现他 @shuoyueqishilove 写得很清晰, 简明易懂;不过原程序无法直接达到我想在stm32上运行的需求,安全性 也有很...
  • 异想家纯C语言矩阵运算

    千次阅读 2016-11-30 23:55:36
    所以只好用C语言写一些在高级语言里一个函数就解决的算法了,由于算法需要运用矩阵运算,自己就先用纯C语言写了个简单的矩阵运算库。  代码只实现了矩阵最基本的运算,包括矩阵的加、减、乘、数乘、转置、行列式、...

      Sandeepin最近做的项目中需要在嵌入式芯片里跑一些算法,而这些单片机性能不上不下,它能跑些简单的程序,但又还没到上Linux系统的地步。所以只好用C语言写一些在高级语言里一个函数就解决的算法了,由于算法需要运用矩阵运算,自己就先用纯C语言写了个简单的矩阵运算库。

      代码只实现了矩阵最基本的运算,包括矩阵的加、减、乘、数乘、转置、行列式、逆矩阵、代数余子式、伴随矩阵等运算。此外增加了一些实用函数,如显示矩阵、从csv文件读取保存矩阵等函数。具体的例子在主函数中体现,其中还用自己这个矩阵运算库做了一个简单的应用,利用公式β=(X'X)^(-1)X'Y来进行多元线性回归系数计算,大家可以参考参考,欢迎批评。

      JfzMatLib.c文件代码:

    #include "JfzMatLib.h"
    
    int main(int argc, char** argv)
    {
    	//矩阵的基本运算:加、减、乘、数乘、转置、行列式、逆矩阵、代数余子式、伴随矩阵
    	//初始实验矩阵
    	double A[] = { -3, 2, -5, -1, 0, -2, 3, -4, 1 };
    	double B[] = { 1, 4, 7, 3, 0, 5, -1, 9, 11 };
    	double C[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };
    	//计算结果矩阵
    	double *Add = (double*)malloc(sizeof(double) * 9);
    	double *Sub = (double*)malloc(sizeof(double) * 9);
    	double *Mul = (double*)malloc(sizeof(double) * 9);
    	double *kMat = (double*)malloc(sizeof(double) * 9);
    	double *CT = (double*)malloc(sizeof(double) * 9);
    	double *AInv = (double*)malloc(sizeof(double) * 9);
    	double *Adj = (double*)malloc(sizeof(double) * 9);
    	//显示矩阵A、B、C
    	printf("A:\n"); MatShow(A, 3, 3);
    	printf("B:\n"); MatShow(B, 3, 3);
    	printf("C:\n"); MatShow(C, 3, 3);
    	//矩阵相加
    	printf("A+B:\n");
    	Add = MatAdd(A, B, 3, 3); MatShow(Add, 3, 3);
    	//矩阵相减
    	printf("A-B:\n");
    	Sub = MatSub(A, B, 3, 3); MatShow(Sub, 3, 3);
    	//矩阵相乘
    	printf("A*B:\n");
    	Mul = MatMul(A, 3, 3, B, 3, 3); MatShow(Mul, 3, 3);
    	//矩阵数乘
    	printf("2*C:\n");
    	kMat = MatMulk(C, 3, 3, 2); MatShow(kMat, 3, 3);
    	//矩阵转置
    	printf("C的转置:\n");
    	CT = MatT(C, 3, 3); MatShow(CT, 3, 3);
    	//矩阵行列式值
    	printf("B的行列式值:\n");
    	printf("%16lf\n", MatDet(B, 3));
    	printf("C的行列式值:\n");
    	printf("%16lf\n", MatDet(C, 3));
    	//矩阵的逆
    	printf("A的逆:\n");
    	AInv = MatInv(A, 3, 3); MatShow(AInv, 3, 3);
    	printf("C的逆:\n");
    	MatInv(C, 3, 3);
    	//矩阵代数余子式
    	printf("C的(0,0)元素的代数余子式:\n");
    	printf("%16lf\n", MatACof(C, 3, 0, 0));
    	//矩阵伴随矩阵
    	printf("A的伴随矩阵:\n");
    	Adj = MatAdj(A, 3, 3); MatShow(Adj, 3, 3);
    
    	//蒋方正矩阵库应用:多元线性回归
    	//多元线性回归公式:β=(X'X)^(-1)X'Y
    	double X[15][5] = {
    		1, 316, 1536, 874, 981,//第一列要补1
    		1, 385, 1771, 777, 1386,
    		1, 299, 1565, 678, 1672,
    		1, 326, 1970, 785, 1864,
    		1, 441, 1890, 785, 2143,
    		1, 460, 2050, 709, 2176,
    		1, 470, 1873, 673, 1769,
    		1, 504, 1955, 793, 2207,
    		1, 348, 2016, 968, 2251,
    		1, 400, 2199, 944, 2390,
    		1, 496, 1328, 749, 2287,
    		1, 497, 1920, 952, 2388,
    		1, 533, 1400, 1452, 2093,
    		1, 506, 1612, 1587, 2083,
    		1, 458, 1613, 1485, 2390
    	};
    	double Y[15][1] = {
    		3894,
    		4628,
    		4569,
    		5340,
    		5449,
    		5599,
    		5010,
    		5694,
    		5792,
    		6126,
    		5025,
    		5924,
    		5657,
    		6019,
    		6141
    	};
    	//多元线性回归公式:β=(X'X)^(-1)X'Y
    	double *XT = (double*)malloc(sizeof(double) * 5 * 15);
    	double *XTX = (double*)malloc(sizeof(double) * 5 * 5);
    	double *InvXTX = (double*)malloc(sizeof(double) * 5 * 5);
    	double *InvXTXXT = (double*)malloc(sizeof(double) * 5 * 15);
    	double *InvXTXXTY = (double*)malloc(sizeof(double) * 5 * 1);
    	XT = MatT((double*)X, 15, 5);
    	XTX = MatMul(XT, 5, 15, (double*)X, 15, 5);
    	InvXTX = MatInv(XTX, 5, 5);
    	InvXTXXT = MatMul(InvXTX, 5, 5, XT, 5, 15);
    	InvXTXXTY = MatMul(InvXTXXT, 5, 15, (double*)Y, 15, 1);
    	printf("多元线性回归β系数:\n");
    	MatShow(InvXTXXTY, 5, 1);
    
    	//保存矩阵到csv
    	MatWrite("XTX.csv", XTX, 5, 5);
    	MatWrite("InvXTXXTY.csv", InvXTXXTY, 5, 1);
    	MatWrite("Fuck.csv", A, 3, 3);
    	MatWrite("Fuck2.csv", A, 9, 1);
    
    	//从csv读取矩阵
    	double *Fuck = (double*)malloc(sizeof(double) * 3 * 3);
    	Fuck = MatRead("Fuck.csv");
    	MatShow(Fuck, 3, 3);
    
    	double *Fuck2 = (double*)malloc(sizeof(double) * 9 * 1);
    	Fuck2 = MatRead("Fuck2.csv");
    	MatShow(Fuck2, 9, 1);
    
    	double *InvXTXXTYread = (double*)malloc(sizeof(double) * 5 * 1);
    	InvXTXXTYread = MatRead("InvXTXXTY.csv");
    	MatShow(InvXTXXTYread, 5, 1);
    
    	double *XTXread = (double*)malloc(sizeof(double) * 5 * 5);
    	XTXread = MatRead("XTX.csv");
    	MatShow(XTXread, 5, 5);
    	system("pause");
    }

      JfzMatLib.h头文件代码:

    #include <stdio.h>
    #include <stdlib.h>
    #include <string.h>
    
    void    MatShow(double* Mat, int row, int col);
    double* MatAdd(double* A, double* B, int row, int col);
    double* MatSub(double* A, double* B, int row, int col);
    double* MatMul(double* A, int Arow, int Acol, double* B, int Brow, int Bcol);
    double* MatMulk(double *A, int row, int col, double k);
    double* MatT(double* A, int row, int col);
    double  MatDet(double *A, int row);
    double* MatInv(double *A, int row, int col);
    double  MatACof(double *A, int row, int m, int n);
    double* MatAdj(double *A, int row, int col);
    double *MatRead(char *csvFileName, int row, int col);
    void MatWrite(char *A, int row, int col);
    
    // (det用)功能:求逆序对的个数
    int inver_order(int list[], int n)
    {
    	int ret = 0;
    	for (int i = 1; i < n; i++)
    		for (int j = 0; j < i; j++)
    			if (list[j] > list[i])
    				ret++;
    	return ret;
    }
    
    // (det用)功能:符号函数,正数返回1,负数返回-1
    int sgn(int order)
    {
    	return order % 2 ? -1 : 1;
    }
    
    // (det用)功能:交换两整数a、b的值
    void swap(int *a, int *b)
    {
    	int m;
    	m = *a;
    	*a = *b;
    	*b = m;
    }
    
    // 功能:求矩阵行列式的核心函数
    double det(double *p, int n, int k, int list[], double sum)
    {
    	if (k >= n)
    	{
    		int order = inver_order(list, n);
    		double item = (double)sgn(order);
    		for (int i = 0; i < n; i++)
    		{
    			//item *= p[i][list[i]];
    			item *= *(p + i * n + list[i]);
    		}
    		return sum + item;
    	}
    	else
    	{
    		for (int i = k; i < n; i++)
    		{
    			swap(&list[k], &list[i]);
    			sum = det(p, n, k + 1, list, sum);
    			swap(&list[k], &list[i]);
    		}
    	}
    	return sum;
    }
    
    // 功能:矩阵显示
    // 形参:(输入)矩阵首地址指针Mat,矩阵行数row和列数col。
    // 返回:无
    void MatShow(double* Mat, int row, int col)
    {
    	for (int i = 0; i < row*col; i++)
    	{
    		printf("%16lf ", Mat[i]);
    		if (0 == (i + 1) % col) printf("\n");
    	}
    }
    
    // 功能:矩阵相加
    // 形参:(输入)矩阵A首地址指针A,矩阵B首地址指针B,矩阵A(也是矩阵B)行数row和列数col
    // 返回:A+B
    double* MatAdd(double* A, double* B, int row, int col)
    {
    	double *Out = (double*)malloc(sizeof(double) * row * col);
    	for (int i = 0; i < row; i++)
    		for (int j = 0; j < col; j++)
    			Out[col*i + j] = A[col*i + j] + B[col*i + j];
    	return Out;
    }
    
    // 功能:矩阵相减
    // 形参:(输入)矩阵A,矩阵B,矩阵A(也是矩阵B)行数row和列数col
    // 返回:A-B
    double* MatSub(double* A, double* B, int row, int col)
    {
    	double *Out = (double*)malloc(sizeof(double) * row * col);
    	for (int i = 0; i < row; i++)
    		for (int j = 0; j < col; j++)
    			Out[col*i + j] = A[col*i + j] - B[col*i + j];
    	return Out;
    }
    
    // 功能:矩阵相乘
    // 形参:(输入)矩阵A,矩阵A行数row和列数col,矩阵B,矩阵B行数row和列数col
    // 返回:A*B
    double* MatMul(double* A, int Arow, int Acol, double* B, int Brow, int Bcol)
    {
    	double *Out = (double*)malloc(sizeof(double) * Arow * Bcol);
    	if (Acol != Brow)
    	{
    		printf("        Shit!矩阵不可乘!\n");
    		return NULL;
    	}
    	if (Acol == Brow)
    	{
    	int i, j, k;
    	for (i = 0; i < Arow; i++)
    		for (j = 0; j < Bcol; j++)
    		{
    			Out[Bcol*i + j] = 0;
    			for (k = 0; k < Acol; k++)
    				Out[Bcol*i + j] = Out[Bcol*i + j] + A[Acol*i + k] * B[Bcol*k + j];
    		}
    	return Out;
    	}
    }
    
    // 功能:矩阵数乘(实数k乘以矩阵A)
    // 形参:(输入)矩阵A首地址指针,矩阵行数row和列数col,实数k
    // 返回:kA
    double* MatMulk(double *A, int row, int col, double k)
    {
    	double *Out = (double*)malloc(sizeof(double) * row * col);
    	for (int i = 0; i < row * col; i++)
    	{
    		*Out = *A * k;
    		Out++;
    		A++;
    	}
    	Out = Out - row * col;
    	return Out;
    }
    
    // 功能:矩阵转置
    // 形参:(输入)矩阵A首地址指针A,行数row和列数col
    // 返回:A的转置
    double* MatT(double* A, int row, int col)
    {
    	double *Out = (double*)malloc(sizeof(double) * row * col);
    	for (int i = 0; i < row; i++)
    		for (int j = 0; j < col; j++)
    			Out[row*j + i] = A[col*i + j];
    	return Out;
    }
    
    // 功能:求行列式值
    // 形参:(输入)矩阵A首地址指针A,行数row
    // 返回:A的行列式值
    double MatDet(double *A, int row)
    {
    	int *list = (int*)malloc(sizeof(int) * row);
    	for (int i = 0; i < row; i++)
    		list[i] = i;
    	double Out = det(A, row, 0, list, 0.0);
    	free(list);
    	return Out;
    }
    
    // 功能:矩阵的逆
    // 形参:(输入)矩阵A首地址指针A,行数row和列数col
    // 返回:A的逆
    double *MatInv(double *A, int row, int col)
    {
    	double *Out = (double*)malloc(sizeof(double) * row * col);
    	double det = MatDet(A, row); //求行列式
    	if (det == 0)
    	{
    		printf("        Fuck!矩阵不可逆!\n");
    		return NULL;
    	}
    	if (det != 0)
    	{
    		Out = MatAdj(A, row, col); //求伴随矩阵
    		int len = row * row;
    		for (int i = 0; i < len; i++)
    			*(Out + i) /= det;
    		return Out;
    	}
    }
    
    // 功能:求代数余子式
    // 形参:(输入)矩阵A首地址指针A,矩阵行数row, 元素a的下标m,n(从0开始),
    // 返回:NxN 矩阵中元素A(mn)的代数余子式
    double MatACof(double *A, int row, int m, int n)
    {
    	int len = (row - 1) * (row - 1);
    	double *cofactor = (double*)malloc(sizeof(double) * len);
    
    	int count = 0;
    	int raw_len = row * row;
    	for (int i = 0; i < raw_len; i++)
    		if (i / row != m && i % row != n)
    			*(cofactor + count++) = *(A + i);
    	double ret = MatDet(cofactor, row - 1);
    	if ((m + n) % 2)
    		ret = -ret;
    	free(cofactor);
    	return ret;
    }
    
    // 功能:求伴随矩阵
    // 形参:(输入)矩阵A首地址指针A,行数row和列数col
    // 返回:A的伴随矩阵
    double *MatAdj(double *A, int row, int col)
    {
    	double *Out = (double*)malloc(sizeof(double) * row * col);
    	int len = row * row;
    	int count = 0;
    	for (int i = 0; i < len; i++)
    	{
    		*(Out + count++) = MatACof(A, row, i % row, i / row);
    	}
    	return Out;
    }
    
    // 读取文件行数
    int FileReadRow(const char *filename)
    {
    	FILE *f = fopen(filename, "r");
    	int i = 0;
    	char str[4096];
    	while (NULL != fgets(str, 4096, f))
    		++i;
    	printf("数组行数:%d\n", i);
    	return i;
    }
    
    // 读取文件每行数据数(逗号数+1)
    int FileReadCol(const char *filename)
    {
    	FILE *f = fopen(filename, "r");
    	int i = 0;
    	char str[4096];
    	fgets(str, 4096, f);
    	for (int j = 0; j < strlen(str); j++)
    	{
    		if (',' == str[j]) i++;
    	}
    	i++;// 数据数=逗号数+1
    	printf("数组列数:%d\n", i);
    	return i;
    }
    
    // 逗号间隔数据提取
    void GetCommaData(char str_In[4096], double double_Out[1024])
    {
    	int str_In_len = strlen(str_In); 
    	//printf("str_In_len:%d\n", str_In_len);
    	char str_Data_temp[128];
    	int j = 0;
    	int double_Out_num = 0;
    	for (int i = 0; i < str_In_len; i++)
    	{
    		//不是逗号,则是数据,存入临时数组中
    		if (',' != str_In[i]) str_Data_temp[j++] = str_In[i];
    		//是逗号或\n(最后一个数据),则数据转换为double,保存到输出数组
    		if (',' == str_In[i] || '\n' == str_In[i]) { str_Data_temp[j] = '\0'; j = 0; /*printf("str_Data_temp:%s\n", str_Data_temp); */double_Out[double_Out_num++] = atof(str_Data_temp); memset(str_Data_temp, 0, sizeof(str_Data_temp)); }
    	}
    }
    
    // 功能:从csv文件读矩阵,保存到指针中
    // 形参:(输入)csv文件名,预计行数row和列数col
    // 返回:矩阵指针A
    double *MatRead(char *csvFileName)
    {
    	int row = FileReadRow(csvFileName); 
    	int col = FileReadCol(csvFileName);
    	double *Out = (double*)malloc(sizeof(double) * row * col);
    	FILE *f = fopen(csvFileName, "r");
    	char buffer[4096];
    	while (fgets(buffer, sizeof(buffer), f))
    	{
    		//printf("buffer[%s]\n",buffer);
    		double double_Out[128] = { 0 };
    		GetCommaData(buffer, double_Out);
    		for (int i = 0; i < col; i++)
    		{
    			//printf("double_Out:%lf\n", double_Out[i]);
    			*Out = double_Out[i];
    			Out++;
    		}
    		
    	}
    	Out = Out - row * col;//指针移回数据开头
    	fclose(f);
    	return Out;
    }
    
    // 功能:将矩阵A存入csv文件中
    // 形参:(输入)保存的csv文件名,矩阵A首地址指针A,行数row和列数col
    // 返回:无
    void MatWrite(const char *csvFileName, double *A, int row, int col)
    {
    	FILE *DateFile;
    	double *Ap = A;
    	DateFile = fopen(csvFileName, "w");//追加的方式保存生成的时间戳
    	for (int i = 0; i < row*col; i++)
    	{
    		if ((i+1) % col == 0) fprintf(DateFile, "%lf\n", *Ap);//保存到文件,到列数换行
    		else fprintf(DateFile, "%lf,", *Ap);//加逗号
    		Ap++;
    	}
    	fclose(DateFile);
    }
    

      运行结果如图:


    展开全文
  • C语言矩阵运算库(Light Matrix)

    万次阅读 多人点赞 2017-08-09 04:14:37
    最近在做卡尔曼滤波和最小二乘的一些算法,都需要运用到矩阵运算,所以索性就写了个纯C的矩阵库(Light Matrix),只所以叫Light Matrix,因为目前只包含了矩阵的基本运算,尽量做到短小精悍。目前支持矩阵的加,减...

    最近在飞控上做卡尔曼滤波和最小二乘的一些算法,都需要运用到矩阵的运算,所以索性就写了个纯C的矩阵库(Light Matrix),之所以叫Light Matrix,因为目前只包含了矩阵的基本运算,尽量做到短小精悍。目前支持矩阵的加,减,乘,转置,行列式,伴随矩阵和逆矩阵,后续有时间会继续更新,可以关注我的github(https://github.com/zjc666/LightMatrix)
    如果发现问题欢迎随时交流,或者在git上提交request


    具体代码如下:

    main.c

    #include <stdio.h>
    #include "light_matrix.h"
    
    void main(void)
    {
        Mat a;
        float val[] = {
            1, 2, 3,
            4, 5, 6,
        };
        Mat b;
        float val2[] = {
            3, 6,
            8, 1,
            9, 2
        };
        Mat c;
        Mat d;
        float val3[] = {
            3, 2, -3,
            10, -3, 2,
            -3, 5, 9,
        };
        float det;
    
        MatCreate(&a, 2, 3);
        MatDump(MatSetVal(&a, val));
        MatCreate(&b, 3, 2);
        MatDump(MatSetVal(&b, val2));
        MatCreate(&c, 3, 3);
        MatCreate(&d, 3, 3);
        MatDump(MatSetVal(&d, val3));
    
        MatDump(MatMul(&b, &a, &c));
        MatDump(MatTrans(&a, &b));
        MatDump(MatInv(&d, &c));
        det = MatDet(&d);
        printf("det(d)=%.2f\n", det);
    
        MatDelete(&a);
        MatDelete(&b);
        MatDelete(&c);
        MatDelete(&d);
    }

    light_matrix.h

    /*
     * Light Matrix: C code implementation for basic matrix operation
     *
     * Copyright (C) 2017 Jiachi Zou
     *
     * This code is free software: you can redistribute it and/or modify
     * it under the terms of the GNU Lesser General Public License as
     * published by the Free Software Foundation, either version 3 of the
     * License, or (at your option) any later version.
     *
     * This code is distributed in the hope that it will be useful,
     * but WITHOUT ANY WARRANTY without even the implied warranty of
     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     * GNU General Public License for more details.
     *
     * You should have received a copy of the GNU Lesser General Public License
     * along with code.  If not, see <http:#www.gnu.org/licenses/>.
     */
    
    #ifndef __LIGHT_MATRIX__
    #define __LIGHT_MATRIX__
    
    typedef struct  {
        int row, col;
        float **element;
        unsigned char init;
    }Mat;
    
    Mat* MatCreate(Mat* mat, int row, int col);
    void MatDelete(Mat* mat);
    Mat* MatSetVal(Mat* mat, float* val);
    void MatDump(const Mat* mat);
    
    Mat* MatZeros(Mat* mat);
    Mat* MatEye(Mat* mat);
    
    Mat* MatAdd(Mat* src1, Mat* src2, Mat* dst);
    Mat* MatSub(Mat* src1, Mat* src2, Mat* dst);
    Mat* MatMul(Mat* src1, Mat* src2, Mat* dst);
    Mat* MatTrans(Mat* src, Mat* dst);
    float MatDet(Mat* mat);
    Mat* MatAdj(Mat* src, Mat* dst);
    Mat* MatInv(Mat* src, Mat* dst);
    
    #endif

    light_matrix.c

    /*
     * Light Matrix: C code implementation for basic matrix operation
     *
     * Copyright (C) 2017 Jiachi Zou
     *
     * This code is free software: you can redistribute it and/or modify
     * it under the terms of the GNU Lesser General Public License as
     * published by the Free Software Foundation, either version 3 of the
     * License, or (at your option) any later version.
     *
     * This code is distributed in the hope that it will be useful,
     * but WITHOUT ANY WARRANTY without even the implied warranty of
     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     * GNU General Public License for more details.
     *
     * You should have received a copy of the GNU Lesser General Public License
     * along with code.  If not, see <http:#www.gnu.org/licenses/>.
     */
    
    #include "light_matrix.h"
    #include <stdio.h>
    #include <stdlib.h>
    
    #define MAT_LEGAL_CHECKING
    #define MAT_INIT_FLAG   0x5C
    
    /************************************************************************/
    /*                          Private Function                            */
    /************************************************************************/
    
    void swap(int *a, int *b)
    {
      int m;
      m = *a;
      *a = *b;
      *b = m;
    }
    
    void perm(int list[], int k, int m, int* p, Mat* mat, float* det) 
    {
      int i;
    
      if(k > m)
      {
        float res = mat->element[0][list[0]];
    
        for(i = 1; i < mat->row ; i++){
            res *= mat->element[i][list[i]];
        }
        if(*p%2){
            //odd is negative
            *det -= res;
        }else{
            //even is positive
            *det += res;
        }
      }
      else
      {
        perm(list, k + 1, m, p, mat, det);
        for(i = k+1; i <= m; i++)
        {
          swap(&list[k], &list[i]);
          *p += 1;
          perm(list, k + 1, m, p, mat, det);
          swap(&list[k], &list[i]);
          *p -= 1; 
        }
      }
    }
    
    /************************************************************************/
    /*                           Public Function                            */
    /************************************************************************/
    
    Mat* MatCreate(Mat* mat, int row, int col)
    {
        int i;
    
    #ifdef MAT_LEGAL_CHECKING
        if(mat->init == MAT_INIT_FLAG){
            if(mat->row == row && mat->col == col)
                return mat;
            else
                MatDelete(mat);
        }
    #endif
    
        mat->element = (float**)malloc(row * sizeof(float*));
        for(i = 0 ; i < row ; i++){
            mat->element[i] = (float*)malloc(col * sizeof(float));  
        }
    
        if(mat->element == NULL){
            return NULL;
        }
        mat->row = row;
        mat->col = col;
        mat->init = MAT_INIT_FLAG;
    
        return mat;
    }
    
    void MatDelete(Mat* mat)
    {
        int i;
    
    #ifdef MAT_LEGAL_CHECKING
        if(mat->init != MAT_INIT_FLAG){
            return ;
        }
    #endif
    
        for(i = 0 ; i<mat->row ; i++)
            free(mat->element[i]);
        free(mat->element);
    
        mat->init = 0;
    }
    
    Mat* MatSetVal(Mat* mat, float* val)
    {
        int row,col;
    
        for(row = 0 ; row < mat->row ; row++){
            for(col = 0 ; col < mat->col ; col++){
                mat->element[row][col] = val[col + row * mat->col];
            }
        }
    
        return mat;
    }
    
    void MatDump(const Mat* mat)
    {
        int row,col;
    
    #ifdef MAT_LEGAL_CHECKING
        if(mat == NULL){
            return ;
        }
    #endif
    
        printf("Mat %dx%d:\n", mat->row, mat->col);
        for(row = 0 ; row < mat->row ; row++){
            for(col = 0 ; col < mat->col ; col++){
                printf("%.4f\t", mat->element[row][col]);
            }
            printf("\n");
        }
    }
    
    Mat* MatZeros(Mat* mat)
    {
        int row,col;
    
    #ifdef MAT_LEGAL_CHECKING
        if(mat->init != MAT_INIT_FLAG){
            printf("err check, none init matrix for MatZeros\n");
            return NULL;
        }
    #endif
    
        for(row = 0 ; row < mat->row ; row++){
            for(col = 0 ; col < mat->col ; col++){
                mat->element[row][col] = 0.0f;
            }
        }
    
        return mat;
    }
    
    Mat* MatEye(Mat* mat)
    {
        int i;
    
    #ifdef MAT_LEGAL_CHECKING
        if(mat->init != MAT_INIT_FLAG){
            printf("err check, none init matrix for MatEye\n");
            return NULL;
        }
    #endif
    
        for(i = 0 ; i < min(mat->row, mat->col) ; i++){
            mat->element[i][i] = 1.0f;
        }
    
        return mat;
    }
    
    /* dst = src1 + src2 */
    Mat* MatAdd(Mat* src1, Mat* src2, Mat* dst)
    {
        int row, col;
    
    #ifdef MAT_LEGAL_CHECKING
        if( !(src1->row == src2->row && src2->row == dst->row && src1->col == src2->col && src2->col == dst->col) ){
            printf("err check, unmatch matrix for MatAdd\n");
            MatDump(src1);
            MatDump(src2);
            MatDump(dst);
            return NULL;
        }
    #endif
    
        for(row = 0 ; row < src1->row ; row++){
            for(col = 0 ; col < src1->col ; col++){
                dst->element[row][col] = src1->element[row][col] + src2->element[row][col];
            }
        }
    
        return dst;
    }
    
    /* dst = src1 - src2 */
    Mat* MatSub(Mat* src1, Mat* src2, Mat* dst)
    {
        int row, col;
    
    #ifdef MAT_LEGAL_CHECKING
        if( !(src1->row == src2->row && src2->row == dst->row && src1->col == src2->col && src2->col == dst->col) ){
            printf("err check, unmatch matrix for MatSub\n");
            MatDump(src1);
            MatDump(src2);
            MatDump(dst);
            return NULL;
        }
    #endif
    
        for(row = 0 ; row < src1->row ; row++){
            for(col = 0 ; col < src1->col ; col++){
                dst->element[row][col] = src1->element[row][col] - src2->element[row][col];
            }
        }
    
        return dst;
    }
    
    /* dst = src1 * src2 */
    Mat* MatMul(Mat* src1, Mat* src2, Mat* dst)
    {
        int row, col;
        int i;
        float temp;
    
    #ifdef MAT_LEGAL_CHECKING
        if( src1->col != src2->row || src1->row != dst->row || src2->col != dst->col ){
            printf("err check, unmatch matrix for MatMul\n");
            MatDump(src1);
            MatDump(src2);
            MatDump(dst);
            return NULL;
        }
    #endif
    
        for(row = 0 ; row < dst->row ; row++){
            for(col = 0 ; col < dst->col ; col++){
                temp = 0.0f;
                for(i = 0 ; i < src1->col ; i++){
                    temp += src1->element[row][i] * src2->element[i][col];
                }
                dst->element[row][col] = temp;
            }
        }
    
        return dst;
    }
    
    /* dst = src' */
    Mat* MatTrans(Mat* src, Mat* dst)
    {
        int row, col;
    
    #ifdef MAT_LEGAL_CHECKING
        if( src->row != dst->col || src->col != dst->row ){
            printf("err check, unmatch matrix for MatTranspose\n");
            MatDump(src);
            MatDump(dst);
            return NULL;
        }
    #endif
    
        for(row = 0 ; row < dst->row ; row++){
            for(col = 0 ; col < dst->col ; col++){
                dst->element[row][col] = src->element[col][row];
            }
        }
    
        return dst;
    }
    
    // return det(mat)
    float MatDet(Mat* mat)
    {
        float det = 0.0f;
        int plarity = 0;
        int *list;
        int i;
    
    #ifdef MAT_LEGAL_CHECKING
        if( mat->row != mat->col){
            printf("err check, not a square matrix for MatDetermine\n");
            MatDump(mat);
            return 0.0f;
        }
    #endif
    
        list = (int*)malloc(sizeof(int)*mat->col);
        for(i = 0 ; i < mat->col ; i++)
            list[i] = i;
    
        perm(list, 0, mat->row-1, &plarity, mat, &det);
        free(list);
    
        return det;
    }
    
    // dst = adj(src)
    Mat* MatAdj(Mat* src, Mat* dst)
    {
        Mat smat;
        int row, col;
        int i,j,r,c;
        float det;
    
    #ifdef MAT_LEGAL_CHECKING
        if( src->row != src->col || src->row != dst->row || src->col != dst->col){
            printf("err check, not a square matrix for MatAdj\n");
            MatDump(src);
            MatDump(dst);
            return NULL;
        }
    #endif
    
        MatCreate(&smat, src->row-1, src->col-1);
    
        for(row = 0 ; row < src->row ; row++){
            for(col = 0 ; col < src->col ; col++){
                r = 0;
                for(i = 0 ; i < src->row ; i++){
                    if(i == row)
                        continue;
                    c = 0;
                    for(j = 0; j < src->col ; j++){
                        if(j == col)
                            continue;
                        smat.element[r][c] = src->element[i][j];
                        c++;
                    }
                    r++;
                }
                det = MatDet(&smat);
                if((row+col)%2)
                    det = -det;
                dst->element[col][row] = det;
            }
        }
    
        MatDelete(&smat);
    
        return dst;
    }
    
    // dst = src^(-1)
    Mat* MatInv(Mat* src, Mat* dst)
    {
        Mat adj_mat;
        float det;
        int row, col;
    
    #ifdef MAT_LEGAL_CHECKING
        if( src->row != src->col || src->row != dst->row || src->col != dst->col){
            printf("err check, not a square matrix for MatInv\n");
            MatDump(src);
            MatDump(dst);
            return NULL;
        }
    #endif
    
        MatCreate(&adj_mat, src->row, src->col);
        MatAdj(src, &adj_mat);
        det = MatDet(src);
    
        for(row = 0 ; row < src->row ; row++){
            for(col = 0 ; col < src->col ; col++)
                dst->element[row][col] = adj_mat.element[row][col]/det;
        }
    
        MatDelete(&adj_mat);
    
        return dst;
    }
    展开全文
  • C语言实现矩阵运算

    万次阅读 多人点赞 2019-09-04 21:28:06
    最近在学习机器人运动控制学,用到了矩阵运算,并用C语言实现之 一个矩阵最基本的有行数line,列数row和 行数乘以列数个数据(row*line), 所以用一个最基本的结构体变量来表示一个矩阵; 矩阵的结构体: typedef struct...
  • 自写的C语言矩阵简易运算

    千次阅读 多人点赞 2019-07-22 19:52:30
    ROS中的矩阵运算也是基于Eigen库的,但是我目前想自己做一做这个底层驱动,涉及正逆运动学、关节速度规划、空间姿态插补算法等,而我现有的单片机不支持这个Eigen库,所以就写了一个简单的基于C语言矩阵运算库,...
  • 基于C语言矩阵运算

    万次阅读 多人点赞 2018-05-23 22:46:43
    最近本着锻炼自己编程能力的目的,捣鼓了一下矩阵运算,看到网上很多的矩阵运算库都是基于C++的,很少有基于C语言的,于是自己想要写一个纯C语言矩阵运算库。大神看到了不要喷我,我只是个小白。 个人感觉矩阵...
  • c语言实现矩阵运算

    千次阅读 2020-06-12 20:08:49
    1.一个矩阵要有的信息包括它的行数列数以及每一个元素的信息,所以我们设置了一个结构体来包含这些信息,结构体的成员就是整型的行数和列数和一个表示元素内容的二级指针。其中二级指针比较难以理解,二级指针是指向...
  • #include,n? ?float?a[20][20],b[20][20],c[20][20]? ?int?i,j? ? printf"请输入矩阵行数? ?scanf%d&m? ?printf"请输入矩阵列数? ?scanf%d&n? ?printf"请输入第一个矩阵? ?for(i=0;i;i
  • //任务二矩阵的基本运算 #include<stdio.h> #include<stdlib.h> #define R1 4//矩阵 MA 行数可以按具体情况修改 #define C1 4//矩阵 MA 列数可以按具体情况修改 #define R2 4//矩阵 MB 行数可以按具体情况修改 #...
  • .z .z 任务二矩阵的基本运算 #in clude<stdio.h> #in clude<stdlib.h> #define R1 4矩阵MA行数可以按具体情况修改 #define C1 4//矩阵MA列数可以按具体情况修改 #define R2 4//矩阵MB行数可以按具体情况修改 #define...
  • //任务二矩阵的基本运算 #include<stdio.h> #include<stdlib.h> #define R1 4//矩阵MA行数可以按具体情况修改 #define C1 4//矩阵MA列数可以按具体情况修改 #define R2 4//矩阵MB行数可以按具体情况修改 #define C2 ...
  • C语言实现的矩阵运算

    热门讨论 2013-01-13 11:05:23
    这是我在VS2010 环境下用C语言写的几个有用的矩阵运算的算法。包括求矩阵的逆、转置、行列式和乘法运算。
  • c语言实现矩阵及其运算,gcc编译器

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 803
精华内容 321
关键字:

c语言矩阵运算

c语言 订阅