C++矩阵库：newmat

矩阵运算基础（线性代数）：http://www2.edu-edu.com.cn/lesson_crs78/self/j_0022/soft/ch0605.html
逆矩阵：https://baike.baidu.com/item/逆矩阵/10481136?fr=aladdin
三角矩阵：https://baike.baidu.com/item/三角矩阵/10232350?fr=aladdin
乔里斯基（cholesky）分解是将一个正定的厄米特（hermite）矩阵分解为一个下三角矩阵和它的共轭转置矩阵。
厄米特矩阵（共轭矩阵），矩阵对称的两个数字互相共轭。实数的共轭是自己本身，复数的共轭-实部不变虚部相反。

一下代码为C语言版本，源代码转自：http://blog.csdn.net/jdh99/article/details/7402710 作者：jdh
C++版本的源代码：http://blog.csdn.net/jdh99/article/details/7360091 作者：jdh

还有一篇关于C实现矩阵运算的文章：http://blog.csdn.net/linaijunix/article/details/50358617
该篇文章中实现了转置、加减法、乘法、求逆和乔里斯基分解。

源代码:
matrix.h:
#ifndef _MATRIX_H   
#define _MATRIX_H 
 
//头文件   
#include "stdio.h"
#include "stdlib.h"
   
//矩阵数据结构   
//二维矩阵   
struct _Matrix   
{    
    int m;   
    int n;   
    float *arr;   
};  
 
//矩阵方法 
//设置m   
void matrix_set_m(struct _Matrix *m,int mm);   
//设置n   
void matrix_set_n(struct _Matrix *m,int nn);   
//初始化   
void matrix_init(struct _Matrix *m);   
//释放   
void matrix_free(struct _Matrix *m);   
//读取i,j坐标的数据   
//失败返回-31415,成功返回值   
float matrix_read(struct _Matrix *m,int i,int j);   
//写入i,j坐标的数据   
//失败返回-1,成功返回1   
int matrix_write(struct _Matrix *m,int i,int j,float val);  
 
//矩阵运算 
//成功返回1,失败返回-1   
int matrix_add(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);   
//C = A - B   
//成功返回1,失败返回-1   
int matrix_subtract(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);   
//C = A * B   
//成功返回1,失败返回-1   
int matrix_multiply(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);   
//行列式的值,只能计算2 * 2,3 * 3   
//失败返回-31415,成功返回值   
float matrix_det(struct _Matrix *A);   
//求转置矩阵,B = AT   
//成功返回1,失败返回-1   
int matrix_transpos(struct _Matrix *A,struct _Matrix *B);   
//求逆矩阵,B = A^(-1)   
//成功返回1,失败返回-1   
int matrix_inverse(struct _Matrix *A,struct _Matrix *B);   
   
#endif   

matrix.c:
#include "matrix.h" 
 
//矩阵方法 
//设置m   
void matrix_set_m(struct _Matrix *m,int mm) 
{ 
    m->m = mm;  
} 
 
//设置n   
void matrix_set_n(struct _Matrix *m,int nn) 
{ 
    m->n = nn; 
} 
 
//初始化   
void matrix_init(struct _Matrix *m) 
{ 
    m->arr = (float *)malloc(m->m * m->n * sizeof(float)); 
} 
 
//释放   
void matrix_free(struct _Matrix *m) 
{ 
    free(m->arr); 
} 
 
//读取i,j坐标的数据   
//失败返回-31415,成功返回值   
float matrix_read(struct _Matrix *m,int i,int j) 
{ 
    if (i >= m->m || j >= m->n)   
    {   
        return -31415;   
    }   
       
    return *(m->arr + i * m->n + j);   
} 
 
//写入i,j坐标的数据   
//失败返回-1,成功返回1   
int matrix_write(struct _Matrix *m,int i,int j,float val) 
{ 
    if (i >= m->m || j >= m->n)   
    {   
        return -1;   
    }   
       
    *(m->arr + i * m->n + j) = val;   
    return 1;   
} 
 
//矩阵运算 
//成功返回1,失败返回-1   
int matrix_add(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C) 
{ 
    int i = 0;   
    int j = 0;   
       
    //判断是否可以运算   
    if (A->m != B->m || A->n != B->n || \ 
        A->m != C->m || A->n != C->n)   
    {   
        return -1;   
    }   
    //运算   
    for (i = 0;i < C->m;i++)   
    {   
        for (j = 0;j < C->n;j++)   
        {   
            matrix_write(C,i,j,matrix_read(A,i,j) + matrix_read(B,i,j));   
        }   
    }   
       
    return 1;  
} 
 
//C = A - B   
//成功返回1,失败返回-1   
int matrix_subtract(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C) 
{ 
    int i = 0;   
    int j = 0;   
       
    //判断是否可以运算   
    if (A->m != B->m || A->n != B->n || \ 
        A->m != C->m || A->n != C->n)   
    {   
        return -1;   
    }   
    //运算   
    for (i = 0;i < C->m;i++)   
    {   
        for (j = 0;j < C->n;j++)   
        {   
            matrix_write(C,i,j,matrix_read(A,i,j) - matrix_read(B,i,j));   
        }   
    }   
       
    return 1;  
} 
 
//C = A * B   
//成功返回1,失败返回-1   
int matrix_multiply(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C) 
{ 
    int i = 0;   
    int j = 0;   
    int k = 0;   
    float temp = 0;   
       
    //判断是否可以运算   
    if (A->m != C->m || B->n != C->n || \ 
        A->n != B->m)   
    {   
        return -1;   
    }   
    //运算   
    for (i = 0;i < C->m;i++)   
    {   
        for (j = 0;j < C->n;j++)   
        {   
            temp = 0;   
            for (k = 0;k < A->n;k++)   
            {   
                temp += matrix_read(A,i,k) * matrix_read(B,k,j);   
            }   
            matrix_write(C,i,j,temp);   
        }   
    }   
       
    return 1;  
} 
 
//行列式的值,只能计算2 * 2,3 * 3   
//失败返回-31415,成功返回值   
float matrix_det(struct _Matrix *A) 
{ 
    float value = 0;   
       
    //判断是否可以运算   
    if (A->m != A->n || (A->m != 2 && A->m != 3))   
    {   
        return -31415;   
    }   
    //运算   
    if (A->m == 2)   
    {   
        value = matrix_read(A,0,0) * matrix_read(A,1,1) - matrix_read(A,0,1) * matrix_read(A,1,0);   
    }   
    else   
    {   
        value = matrix_read(A,0,0) * matrix_read(A,1,1) * matrix_read(A,2,2) + \ 
                matrix_read(A,0,1) * matrix_read(A,1,2) * matrix_read(A,2,0) + \ 
                matrix_read(A,0,2) * matrix_read(A,1,0) * matrix_read(A,2,1) - \ 
                matrix_read(A,0,0) * matrix_read(A,1,2) * matrix_read(A,2,1) - \ 
                matrix_read(A,0,1) * matrix_read(A,1,0) * matrix_read(A,2,2) - \ 
                matrix_read(A,0,2) * matrix_read(A,1,1) * matrix_read(A,2,0);   
    }   
       
    return value;  
} 
 
//求转置矩阵,B = AT   
//成功返回1,失败返回-1   
int matrix_transpos(struct _Matrix *A,struct _Matrix *B) 
{ 
    int i = 0;   
    int j = 0;   
       
    //判断是否可以运算   
    if (A->m != B->n || A->n != B->m)   
    {   
        return -1;   
    }   
    //运算   
    for (i = 0;i < B->m;i++)   
    {   
        for (j = 0;j < B->n;j++)   
        {   
            matrix_write(B,i,j,matrix_read(A,j,i));   
        }   
    }   
       
    return 1;   
} 
 
//求逆矩阵,B = A^(-1)   
//成功返回1,失败返回-1   
int matrix_inverse(struct _Matrix *A,struct _Matrix *B) 
{ 
    int i = 0;   
    int j = 0;   
    int k = 0;   
    struct _Matrix m;   
    float temp = 0;   
    float b = 0;   
       
    //判断是否可以运算   
    if (A->m != A->n || B->m != B->n || A->m != B->m)   
    {   
        return -1;   
    }   
       
       
       
    //增广矩阵m = A | B初始化    
    matrix_set_m(&m,A->m); 
    matrix_set_n(&m,2 * A->m); 
    matrix_init(&m); 
    for (i = 0;i < m.m;i++)   
    {   
        for (j = 0;j < m.n;j++)   
        {   
            if (j <= A->n - 1)   
            {   
                matrix_write(&m,i,j,matrix_read(A,i,j));   
            }   
            else   
            {   
                if (i == j - A->n)   
                {   
                    matrix_write(&m,i,j,1);   
                }   
                else   
                {   
                    matrix_write(&m,i,j,0);   
                }   
            }   
        }   
    }   
       
    //高斯消元   
    //变换下三角   
    for (k = 0;k < m.m - 1;k++)   
    {   
        //如果坐标为k,k的数为0,则行变换   
        if (matrix_read(&m,k,k) == 0)   
        {   
            for (i = k + 1;i < m.m;i++)   
            {   
                if (matrix_read(&m,i,k) != 0)   
                {   
                    break;   
                }   
            }   
            if (i >= m.m)   
            {   
                return -1;   
            }   
            else   
            {   
                //交换行   
                for (j = 0;j < m.n;j++)   
                {   
                    temp = matrix_read(&m,k,j);   
                    matrix_write(&m,k,j,matrix_read(&m,k + 1,j));   
                    matrix_write(&m,k + 1,j,temp);   
                }   
            }   
        }   
           
        //消元   
        for (i = k + 1;i < m.m;i++)   
        {   
            //获得倍数   
            b = matrix_read(&m,i,k) / matrix_read(&m,k,k);   
            //行变换   
            for (j = 0;j < m.n;j++)   
            {   
                temp = matrix_read(&m,i,j) - b * matrix_read(&m,k,j);   
                matrix_write(&m,i,j,temp);   
            }   
        }   
    }   
    //变换上三角   
    for (k = m.m - 1;k > 0;k--)   
    {   
        //如果坐标为k,k的数为0,则行变换   
        if (matrix_read(&m,k,k) == 0)   
        {   
            for (i = k + 1;i < m.m;i++)   
            {   
                if (matrix_read(&m,i,k) != 0)   
                {   
                    break;   
                }   
            }   
            if (i >= m.m)   
            {   
                return -1;   
            }   
            else   
            {   
                //交换行   
                for (j = 0;j < m.n;j++)   
                {   
                    temp = matrix_read(&m,k,j);   
                    matrix_write(&m,k,j,matrix_read(&m,k + 1,j));   
                    matrix_write(&m,k + 1,j,temp);   
                }   
            }   
        }   
           
        //消元   
        for (i = k - 1;i >= 0;i--)   
        {   
            //获得倍数   
            b = matrix_read(&m,i,k) / matrix_read(&m,k,k);   
            //行变换   
            for (j = 0;j < m.n;j++)   
            {   
                temp = matrix_read(&m,i,j) - b * matrix_read(&m,k,j);   
                matrix_write(&m,i,j,temp);   
            }   
        }   
    }   
    //将左边方阵化为单位矩阵   
    for (i = 0;i < m.m;i++)   
    {   
        if (matrix_read(&m,i,i) != 1)   
        {   
            //获得倍数   
            b = 1 / matrix_read(&m,i,i);   
            //行变换   
            for (j = 0;j < m.n;j++)   
            {   
                temp = matrix_read(&m,i,j) * b;   
                matrix_write(&m,i,j,temp);   
            }   
        }   
    }   
    //求得逆矩阵   
    for (i = 0;i < B->m;i++)   
    {   
        for (j = 0;j < B->m;j++)   
        {   
            matrix_write(B,i,j,matrix_read(&m,i,j + m.m));   
        }   
    }   
    //释放增广矩阵   
    matrix_free(&m);   
       
    return 1;  
} 

main.c:(测试代码)
#include  
#include "matrix.h" 
 
//打印2维矩阵 
void printf_matrix(struct _Matrix *A) 
{ 
    int i = 0; 
    int j = 0; 
    int m = 0; 
    int n = 0; 
     
    m = A->m; 
    n = A->n; 
    for (i = 0;i < m;i++) 
    { 
        for (j = 0;j < n;j++) 
        { 
            printf("%f\t",matrix_read(A,i,j)); 
        } 
        printf("\n"); 
    } 
} 
 
int main() 
{ 
    int i = 0; 
    int j = 0; 
    int k = 0; 
    struct _Matrix m1; 
    struct _Matrix m2; 
    struct _Matrix m3; 
     
    //初始化内存 
    matrix_set_m(&m1,3); 
    matrix_set_n(&m1,3); 
    matrix_init(&m1); 
    matrix_set_m(&m2,3); 
    matrix_set_n(&m2,3); 
    matrix_init(&m2); 
    matrix_set_m(&m3,3); 
    matrix_set_n(&m3,3); 
    matrix_init(&m3); 
     
    //初始化数据 
    k = 1; 
    for (i = 0;i < m1.m;i++) 
    { 
        for (j = 0;j < m1.n;j++) 
        { 
            matrix_write(&m1,i,j,k++); 
        } 
    } 
     
    for (i = 0;i < m2.m;i++) 
    { 
        for (j = 0;j < m2.n;j++) 
        { 
            matrix_write(&m2,i,j,k++); 
        } 
    } 
     
    //原数据 
    printf("A:\n"); 
    printf_matrix(&m1); 
    printf("B:\n"); 
    printf_matrix(&m2); 
     
    printf("A:行列式的值%f\n",matrix_det(&m1)); 
     
    //C = A + B 
    if (matrix_add(&m1,&m2,&m3) > 0) 
    { 
        printf("C = A + B:\n"); 
        printf_matrix(&m3); 
    } 
     
    //C = A - B 
    if (matrix_subtract(&m1,&m2,&m3) > 0) 
    { 
        printf("C = A - B:\n"); 
        printf_matrix(&m3); 
    } 
     
    //C = A * B 
    if (matrix_multiply(&m1,&m2,&m3) > 0) 
    { 
        printf("C = A * B:\n"); 
        printf_matrix(&m3); 
    } 
     
    //C = AT 
    if (matrix_transpos(&m1,&m3) > 0) 
    { 
        printf("C = AT:\n"); 
        printf_matrix(&m3); 
    } 
     
    if (matrix_inverse(&m1,&m3) > 0) 
    { 
        printf("C = A^(-1):\n"); 
        printf_matrix(&m3); 
    } 
     
    getchar(); 
    return 0; 
}