C語言科學(xué)計算入門之矩陣乘法的相關(guān)計算
1.矩陣相乘
矩陣相乘應(yīng)滿足的條件:
(1) 矩陣A的列數(shù)必須等于矩陣B的行數(shù),矩陣A與矩陣B才能相乘;
(2) 矩陣C的行數(shù)等于矩陣A的行數(shù),矩陣C的列數(shù)等于矩陣B的列數(shù);
(3) 矩陣C中第i行第j列的元素等于矩陣A的第i行元素與矩陣B的第j列元素對應(yīng)乘積之和,即
如:
則:
2. 常用矩陣相乘算法
用A的第i行分別和B的第j列的各個元素相乘求和,求得C的第i行j列的元素,這種算法中,B的訪問是按列進(jìn)行訪問的,代碼如下:
void arymul(int a[4][5], int b[5][3], int c[4][3]) { int i, j, k; int temp; for(i = 0; i < 4; i++){ for(j = 0; j < 3; j++){ temp = 0; for(k = 0; k < 5; k++){ temp += a[i][k] * b[k][j]; } c[i][j] = temp; printf("%d/t", c[i][j]); } printf("%d/n"); } }
3. 改進(jìn)的算法
矩陣A、B、C都按行(數(shù)據(jù)的存儲順序)訪問,以提高存儲器訪問效率,對于A的第i行中,第j列的元素分別和B的第j行的元素相乘,對于B中相同的列k在上述計算過程中求和,從而得到C第i行k列的數(shù)據(jù),代碼如下:
void arymul1(int a[4][5], int b[5][3], int c[4][3]) { int i, j, k; int temp[3] = {0}; for(i = 0; i < 4; i++){ for(k = 0; k < 3; k ++) temp[k] = 0; for(j = 0; j < 5; j++){//當(dāng)前行的每個元素 for(k = 0; k < 3; k++){ temp[k] += a[i][j] * b[j][k]; } } for(k = 0; k < 3; k++){ c[i][k] = temp[k]; printf("%d/t", c[i][k]); } printf("%d/n"); } }
這種算法很容易轉(zhuǎn)到稀疏矩陣的相乘算法。
PS:斯特拉森算法的實現(xiàn)
斯特拉森方法,是由v.斯特拉森在1969年提出的一個方法。
我們先討論二階矩陣的計算方法。
對于二階矩陣
a11 a12 b11 b12 A = a21 a22 B = b21 b22
先計算下面7個量(1)
x1 = (a11 + a22) * (b11 + b22); x2 = (a21 + a22) * b11; x3 = a11 * (b12 - b22); x4 = a22 * (b21 - b11); x5 = (a11 + a12) * b22; x6 = (a21 - a11) * (b11 + b12); x7 = (a12 - a22) * (b21 + b22);
再設(shè)C = AB。根據(jù)矩陣相乘的規(guī)則,C的各元素為(2)
c11 = a11 * b11 + a12 * b21 c12 = a11 * b12 + a12 * b22 c21 = a21 * b11 + a22 * b21 c22 = a21 * b12 + a22 * b22
比較(1)(2),C的各元素可以表示為(3)
c11 = x1 + x4 - x5 + x7 c12 = x3 + x5 c21 = x2 + x4 c22 = x1 + x3 - x2 + x6
根據(jù)以上的方法,我們就可以計算4階矩陣了,先將4階矩陣A和B劃分成四塊2階矩陣,分別利用公式計算它們的乘積,再使用(1)(3)來計算出最后結(jié)果。
ma11 ma12 mb11 mb12 A4 = ma21 ma22 B4 = mb21 mb22
其中
a11 a12 a13 a14 b11 b12 b13 b14 ma11 = a21 a22 ma12 = a23 a24 mb11 = b21 b22 mb12 = b23 b24 a31 a32 a33 a34 b31 b32 b33 b34 ma21 = a41 a42 ma22 = a43 a44 mb21 = b41 b42 mb22 = b43 b44
實現(xiàn)
// 計算2X2矩陣 void Multiply2X2(float& fOut_11, float& fOut_12, float& fOut_21, float& fOut_22, float f1_11, float f1_12, float f1_21, float f1_22, float f2_11, float f2_12, float f2_21, float f2_22) { const float x1((f1_11 + f1_22) * (f2_11 + f2_22)); const float x2((f1_21 + f1_22) * f2_11); const float x3(f1_11 * (f2_12 - f2_22)); const float x4(f1_22 * (f2_21 - f2_11)); const float x5((f1_11 + f1_12) * f2_22); const float x6((f1_21 - f1_11) * (f2_11 + f2_12)); const float x7((f1_12 - f1_22) * (f2_21 + f2_22)); fOut_11 = x1 + x4 - x5 + x7; fOut_12 = x3 + x5; fOut_21 = x2 + x4; fOut_22 = x1 - x2 + x3 + x6; } // 計算4X4矩陣 void Multiply(CLAYMATRIX& mOut, const CLAYMATRIX& m1, const CLAYMATRIX& m2) { float fTmp[7][4]; // (ma11 + ma22) * (mb11 + mb22) Multiply2X2(fTmp[0][0], fTmp[0][1], fTmp[0][2], fTmp[0][3], m1._11 + m1._33, m1._12 + m1._34, m1._21 + m1._43, m1._22 + m1._44, m2._11 + m2._33, m2._12 + m2._34, m2._21 + m2._43, m2._22 + m2._44); // (ma21 + ma22) * mb11 Multiply2X2(fTmp[1][0], fTmp[1][1], fTmp[1][2], fTmp[1][3], m1._31 + m1._33, m1._32 + m1._34, m1._41 + m1._43, m1._42 + m1._44, m2._11, m2._12, m2._21, m2._22); // ma11 * (mb12 - mb22) Multiply2X2(fTmp[2][0], fTmp[2][1], fTmp[2][2], fTmp[2][3], m1._11, m1._12, m1._21, m1._22, m2._13 - m2._33, m2._14 - m2._34, m2._23 - m2._43, m2._24 - m2._44); // ma22 * (mb21 - mb11) Multiply2X2(fTmp[3][0], fTmp[3][1], fTmp[3][2], fTmp[3][3], m1._33, m1._34, m1._43, m1._44, m2._31 - m2._11, m2._32 - m2._12, m2._41 - m2._21, m2._42 - m2._22); // (ma11 + ma12) * mb22 Multiply2X2(fTmp[4][0], fTmp[4][1], fTmp[4][2], fTmp[4][3], m1._11 + m1._13, m1._12 + m1._14, m1._21 + m1._23, m1._22 + m1._24, m2._33, m2._34, m2._43, m2._44); // (ma21 - ma11) * (mb11 + mb12) Multiply2X2(fTmp[5][0], fTmp[5][1], fTmp[5][2], fTmp[5][3], m1._31 - m1._11, m1._32 - m1._12, m1._41 - m1._21, m1._42 - m1._22, m2._11 + m2._13, m2._12 + m2._14, m2._21 + m2._23, m2._22 + m2._24); // (ma12 - ma22) * (mb21 + mb22) Multiply2X2(fTmp[6][0], fTmp[6][1], fTmp[6][2], fTmp[6][3], m1._13 - m1._33, m1._14 - m1._34, m1._23 - m1._43, m1._24 - m1._44, m2._31 + m2._33, m2._32 + m2._34, m2._41 + m2._43, m2._42 + m2._44); // 第一塊 mOut._11 = fTmp[0][0] + fTmp[3][0] - fTmp[4][0] + fTmp[6][0]; mOut._12 = fTmp[0][1] + fTmp[3][1] - fTmp[4][1] + fTmp[6][1]; mOut._21 = fTmp[0][2] + fTmp[3][2] - fTmp[4][2] + fTmp[6][2]; mOut._22 = fTmp[0][3] + fTmp[3][3] - fTmp[4][3] + fTmp[6][3]; // 第二塊 mOut._13 = fTmp[2][0] + fTmp[4][0]; mOut._14 = fTmp[2][1] + fTmp[4][1]; mOut._23 = fTmp[2][2] + fTmp[4][2]; mOut._24 = fTmp[2][3] + fTmp[4][3]; // 第三塊 mOut._31 = fTmp[1][0] + fTmp[3][0]; mOut._32 = fTmp[1][1] + fTmp[3][1]; mOut._41 = fTmp[1][2] + fTmp[3][2]; mOut._42 = fTmp[1][3] + fTmp[3][3]; // 第四塊 mOut._33 = fTmp[0][0] - fTmp[1][0] + fTmp[2][0] + fTmp[5][0]; mOut._34 = fTmp[0][1] - fTmp[1][1] + fTmp[2][1] + fTmp[5][1]; mOut._43 = fTmp[0][2] - fTmp[1][2] + fTmp[2][2] + fTmp[5][2]; mOut._44 = fTmp[0][3] - fTmp[1][3] + fTmp[2][3] + fTmp[5][3]; }
比較
在標(biāo)準(zhǔn)的定義算法中我們需要進(jìn)行n * n * n次乘法運算,新算法中我們需要進(jìn)行7log2n次乘法,對于最常用的4階矩陣: 原算法 新算法
加法次數(shù) 48 72(48次加法,24次減法)
乘法次數(shù) 64 49
需要額外空間 16 * sizeof(float) 28 * sizeof(float)
新算法要比原算法多了24次減法運算,少了15次乘法。但因為浮點乘法的運算速度要遠(yuǎn)遠(yuǎn)慢于加/減法運算,所以新算法的整體速度有所提高。
欄 目:C語言
下一篇:C++編程中的const關(guān)鍵字常見用法總結(jié)
本文標(biāo)題:C語言科學(xué)計算入門之矩陣乘法的相關(guān)計算
本文地址:http://mengdiqiu.com.cn/a1/Cyuyan/2631.html
您可能感興趣的文章
- 04-02c語言函數(shù)調(diào)用后清空內(nèi)存 c語言調(diào)用函數(shù)刪除字符
- 04-02c語言的正則匹配函數(shù) c語言正則表達(dá)式函數(shù)庫
- 04-02func函數(shù)+在C語言 func函數(shù)在c語言中
- 04-02c語言中對數(shù)函數(shù)的表達(dá)式 c語言中對數(shù)怎么表達(dá)
- 04-02c語言用函數(shù)寫分段 用c語言表示分段函數(shù)
- 04-02c語言編寫函數(shù)冒泡排序 c語言冒泡排序法函數(shù)
- 04-02c語言沒有round函數(shù) round c語言
- 04-02c語言分段函數(shù)怎么求 用c語言求分段函數(shù)
- 04-02C語言中怎么打出三角函數(shù) c語言中怎么打出三角函數(shù)的值
- 04-02c語言調(diào)用函數(shù)求fibo C語言調(diào)用函數(shù)求階乘


閱讀排行
本欄相關(guān)
- 04-02c語言函數(shù)調(diào)用后清空內(nèi)存 c語言調(diào)用
- 04-02func函數(shù)+在C語言 func函數(shù)在c語言中
- 04-02c語言的正則匹配函數(shù) c語言正則表達(dá)
- 04-02c語言用函數(shù)寫分段 用c語言表示分段
- 04-02c語言中對數(shù)函數(shù)的表達(dá)式 c語言中對
- 04-02c語言編寫函數(shù)冒泡排序 c語言冒泡排
- 04-02c語言沒有round函數(shù) round c語言
- 04-02c語言分段函數(shù)怎么求 用c語言求分段
- 04-02C語言中怎么打出三角函數(shù) c語言中怎
- 04-02c語言調(diào)用函數(shù)求fibo C語言調(diào)用函數(shù)求
隨機(jī)閱讀
- 01-11Mac OSX 打開原生自帶讀寫NTFS功能(圖文
- 01-10C#中split用法實例總結(jié)
- 01-11ajax實現(xiàn)頁面的局部加載
- 01-10delphi制作wav文件的方法
- 01-10使用C語言求解撲克牌的順子及n個骰子
- 01-10SublimeText編譯C開發(fā)環(huán)境設(shè)置
- 08-05dedecms(織夢)副欄目數(shù)量限制代碼修改
- 04-02jquery與jsp,用jquery
- 08-05織夢dedecms什么時候用欄目交叉功能?
- 08-05DEDE織夢data目錄下的sessions文件夾有什