Re: [問題] 矩陣!!限定只能用兩個矩陣...

[ADF (............￾ ￾NN￾)]

※ 引述《LPH66 ((short)(-15074))》之銘言：
: 剛剛想到了XD
: 把 B 做 LUP decomposition
: 也就是求得一個下三角矩陣 L 一個上三角矩陣 U 及一個排列矩陣 P
: 使得 B = LUP
: (講起來好像很複雜  其實做一次高斯消去法把它「紀錄」下來就是了 XD
:  看維基: http://zh-tw.wikipedia.org/wiki/LU%E5%88%86%E8%A7%A3
:  做法修改一下還可以在 B 裡面 in-place 完成!)
: 所以 AB 就成了 ALUP
: 而乘一個三角矩陣是可以 in-place 來做的
: (也就是把 乘B 分成了三步 乘L 乘U 和 乘P
:  其中 乘L 和 乘U 都是乘三角矩陣  可以 in-place
:  乘P 只代表我們要把 A 的 column 間換來換去 XD)

無聊 幫忙寫實作碼 XD

#include <iostream>

using namespace std;
#define SIZE 3

void mulUpMat( float M[SIZE][SIZE] , float U[SIZE][SIZE] )
{

        for ( int i = 0 ; i < SIZE ; ++ i)
        for ( int j = SIZE - 1 ; j >= 0 ; --j )
        {
                float result = 0.0f;
                for ( int k = 0 ; k <= j ; ++k)
                        result += M[i][k] * U[k][j];
                M[i][j] = result;
        }
}


void mulLowMat( float M[SIZE][SIZE] , float L[SIZE][SIZE] )
{
        for ( int i = 0 ; i < SIZE ; ++i )
        for ( int j = 0 ; j < SIZE ; ++j )
        {
                float result = 0.0f;
                for ( int k = j ; k < SIZE ; ++k)
                {
                        if ( k != j )
                                result += M[i][k] * L[k][j];
                        else
                                result += M[i][k];
                }
                M[i][j] = result;
        }
}

void Mul( float S[SIZE][SIZE] , float M[SIZE][SIZE] , float U[SIZE][SIZE]  )
{
        for ( int i = 0 ; i < SIZE ; ++i )
        for ( int j = 0 ; j < SIZE ; ++j )
        {
                float result = 0.0f;
                for ( int k = 0 ; k < SIZE ; ++k)
                        result += M[i][k] * U[k][j];
                S[i][j] = result;
        }
}

void print( float M[SIZE][SIZE] )
{
        for( int i = 0 ; i < SIZE ; ++ i )
        {
                for( int j = 0 ; j < SIZE ; ++ j )
                {
                        cout << M[i][j] << " ";
                }
                cout << endl;
        }
}

void exchangeCol( float M[SIZE][SIZE] , int i1 , int i2 )
{
        for ( int j = 0 ; j < SIZE ; ++j )
        {
                std::swap( M[i1][j] , M[i2][j] );
        }
}

void LUDec( float A[SIZE][SIZE] , float L[SIZE][SIZE] ,  float P[SIZE][SIZE] )
{
        for ( int n = 0 ; n < SIZE ; ++n )
        {
                int ch = n;
                while ( A[ch][n] == 0 )
                {
                        ++ch;
                        if ( ch == SIZE)
                                break;
                }

                if ( ch == SIZE )
                        continue;

                if ( ch != n )
                {
                        exchangeCol( A , n , ch );
                        exchangeCol( P , n , ch );
                }

                if ( &A[0][0] != &L[0][0] )
                        L[n][n] = 1.0f;

                for ( int i = n + 1 ;i < SIZE ; ++i  )
                {
                        float s  =  A[i][n] / A[n][n];
                        A[i][n] = 0;
                        for( int j = n + 1 ; j < SIZE ; ++ j )
                        {
                                A[i][j] -= s * A[n][j];
                        }
                        L[i][n] = s;
                }
        }
}



int main()
{
        float A[ SIZE ][ SIZE ]=
        {
                0, 3, 1,
                5, 3, 5,
                7, 8 ,9,
        };

        float B[ SIZE ][ SIZE ] =
        {
                1, 3, 2,
                3, 2, 9,
                4, 7, 8,
        };

        float S[ SIZE ][ SIZE ];

        Mul( S , A , B );
        print( S );

        LUDec( B , B , A );
        mulLowMat( A , B );
        mulUpMat( A , B );
        print( A );

        return 0;
}

--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 140.112.233.54