矩阵相乘算法

  这学期学校的课程多了很多的计算数学的课，像《地球物理中的有限元法》《地球物理中的边界元法》和《数值分析》之类的，搞得我需要写很多数学的算法，刚开始写个简单的吧，就是矩阵的乘法，大家都知道

设  矩阵 ,  矩阵 , 则由元素

 

构成的  矩阵

称为矩阵  与  的乘积, 记为 .

可能是我逻辑不是很好，照这个定义写的算法嵌套太多了，我写不出来呵呵。所以就像找另外的方法，结果我想出来一个比这个简单得多的办法，我的算法如下，注释也在其中：

//Author:dengjin_CN

//Time:12:06 2009-10-14



#include<stdio.h>
#include<stdlib.h>

/*输入的两个矩阵第一个存储在m1中，第二个存储在m2中，结
果保存在result中，第一个矩阵为x行y列，第二个矩阵为y行z列*/
void matrix_multiply(int *m1,int *m2,int *result,int x,int y,int z);

int main(void){
    
    int result[6] = {0};//此数组的维数等于x*z，每次运行需要修改

    int m1[] = {1,3,2,4,5,2,3,4,3,2,3,4};
    int m2[] = {2,3,1,6,0,7,8,0};
    
    //For testing:

    matrix_multiply(m1,m2,result,3,4,2);

    return EXIT_SUCCESS;    
}

void
matrix_multiply(int *m1,int *m2,int *result,int x,int y,int z){
    
    int l_m1[x][y],l_m2[y][z],l_result[y][x*z];//把原来的一维数组转换到相应行列的二维数组，方便运算

    int i=0,j=0,k=0;
    int loop = 0;

    int part_of_result[y][x*z];//多次运算值

    int (*p_part_of_result)[x*z] = &part_of_result[0];//用指针运算，以下凡是p_part_of_result均可以直接换成part_of_result

    int *p_result = result;
    
    //converting

    for(i=0;i<x;i++)
        for(j=0;j<y;j++)
            l_m1[i][j] = *m1++;
    
    //converting

    for(i=0;i<y;i++)
        for(j=0;j<z;j++)
            l_m2[i][j] = *m2++;
    
    /*算法部分：
    假定有如下两个矩阵相乘
    [a b c d] [A B]
    [e f g h] X [C D]
    [i j k l] [E F]
         [G H]
    把第一个矩阵的每一列的每一个值均与第二个矩阵的对应行相乘，比如
    a*A+a*B+e*A+e*B+i*A+i*B存在第一个part_of_result中，相应的
    b*C+b*D+f*C+f*D+j*C+j*D存在第二个part_of_result中，以此类推(一共y次)
    .
    .
    .
    d*G+d*H+h*G+h*G+l*G+l*H存在最后一个part_of_result中，然后把对应的列相加存到result中，即
    a*A+b*C+...+d*G存在result[0]
    a*B+b*D+...+d*H存在result[1],以此类推
    存在result中的值即为最后结果，格式化输出
    */
    while(loop < y){
        k = 0;
        for(i=0;i<x;i++)
            for(j=0;j<z;j++){
                p_part_of_result[loop][k] = l_m1[i][loop]*l_m2[loop][j];
                k++;
            }
        loop++;
    }

    
    for(i=0;i<x*z;i++){
        for(j=0;j<y;j++)
            *p_result += p_part_of_result[j][i];
        p_result++;
    }
        
    
    //output format

    for(i=0;i<x;i++){
        putchar('[');
        
        for(j=0;j<z;j++)
            printf("%4d ",*result++);
            
        putchar(']');
        putchar('\n');
        
    }
    
}


结果呢，我今天发这个博文的时候Google了一下，发现在wikipedia上一些解释，可以看看：

一般矩阵乘积

矩阵相乘最重要的方法当然是一般矩阵乘积了，它只有在第一个矩阵的列数(column)和第二个矩阵的行数(row)相同时才有定义。一般单指矩阵乘积时，指的便是一般矩阵乘积。若A为m×n矩阵，B为n×p矩阵，则他们的乘积AB(有时记做A · B)会是一个m×p矩阵。其乘积矩阵的元素如下面式子得出：

以上是用矩阵单元的代数系统来说明这类乘法的抽象性质。

[编辑]由定义直接计算

左边的图表示出要如何计算AB的(1,2)和(3,3)元素，当A是个4×2矩阵和B是个2×3矩阵时。分别来自两个矩阵的元素都依箭头方向而两两配对，把每一对中的两个元素相乘，再把这些乘积加总起来，最后得到的值即为箭头相交位置的值。

[编辑]系数－向量方法

这种矩阵乘积亦可由稍微不同的观点来思考：把向量和各系数相乘后相加起来。设A和B是两个给定如下的矩阵：

 and 

则

举个例子来说: