Chinaunix首页 | 论坛 | 博客
  • 博客访问: 553436
  • 博文数量: 92
  • 博客积分: 2511
  • 博客等级: 少校
  • 技术积分: 932
  • 用 户 组: 普通用户
  • 注册时间: 2008-10-19 10:10
文章分类
文章存档

2011年(6)

2010年(27)

2009年(37)

2008年(22)

我的朋友

分类: WINDOWS

2009-12-22 10:41:07

A是一个方块矩阵ALU分解是将它分解成如下形式:

 A = LU, \,

其中LU分别是下三角矩阵和上三角矩阵。

例如对于一个3 \times 3的矩阵,就有

 
        \begin{bmatrix}
           a_{11} & a_{12} & a_{13} \\
           a_{21} & a_{22} & a_{23} \\
           a_{31} & a_{32} & a_{33} \\
        \end{bmatrix} =
      \begin{bmatrix}
           l_{11} & 0 & 0 \\
           l_{21} & l_{22} & 0 \\
           l_{31} & l_{32} & l_{33} \\
        \end{bmatrix}
        \begin{bmatrix}
           u_{11} & u_{12} & u_{13} \\
           0 & u_{22} & u_{23} \\
           0 & 0 & u_{33} \\
        \end{bmatrix}

一个LDU分解是一个如下形式的分解:

A = LDU

其中D对角矩阵LU单位三角矩阵(对角线上全是1的三角矩阵)。

LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵,其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法(Doolittle algorithm):从下至上地对矩阵A做初等行变换,将对角线左下方的元素变成零,然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵,这一系列单位下三角矩阵的乘积的逆就是L矩阵,它也是一个单位下三角矩阵。

这类算法的复杂度一般在\frac{2n^3}{3}左右,对充分消元的分解则不然。

[编辑]杜尔里特算法

对给定的N × N矩阵

A = (an,n)

A(0): = A

然后定义对于n = 1,...,N-1的情况如下:

在第n步,消去矩阵A(n-1)的第n列主对角线下的元素:将A(n-1)的第n行乘以l_{i,n} := -\frac{a_{i,n}^{(n-1)}}{a_{n,n}^{(n-1)}}之后加到第i行上去。其中i = n+1,\ldots,N

这相当于在A(n-1)的左边乘上一个单位下三角矩阵:

 
L_n =
\begin{pmatrix}
     1 &        &           &         &         & 0 \\
       & \ddots &           &         &         &   \\
       &        &         1 &         &         &   \\
       &        & l_{n+1,n} &  \ddots &         &   \\
       &        &    \vdots &         &  \ddots &   \\
     0 &        &   l_{N,n} &         &         & 1 \\
\end{pmatrix}.

于是,定义为:设

A(n): = LnA(n − 1).

经过N-1轮操作后,所有在主对角线下的系数都为0了,于是我们得到了一个上三角矩阵:A(N-1),这时就有:

 
A = L_{1}^{-1} L_{1} A^{(0)}
= L_{1}^{-1} A^{(1)} = L_{1}^{-1} L_{2}^{-1} L_{2} A^{(1)} = 
L_{1}^{-1}L_{2}^{-1} A^{(2)} =\ldots = L_{1}^{-1} \ldots L_{N-1}^{-1} A^{(N-1)}.

这时,矩阵A(N-1) 就是UL=L_{1}^{-1} \ldots L_{N-1}^{-1}。于是我们得到分解:A = LU

显然,要是算法成立,在每步操作时必须有a_{n,n}^{(n-1)} \neq 0。如果这一条件不成立,就要将第n行和另一行交换,由此就会出现一个置换矩阵P。这就是为什么一般来说LU分解里会带有一个置换矩阵的原因。


!Author:RainMan
!2009-12-20
!LU分解法解线性方程组AX=B
!该程序从文件输入,格式为
!矩阵的维数是 3
!A矩阵
!4 2 1
!8 7 4
!16 17 15
!B矩阵
!11 34 95

PROGRAM solve_equation
    REAL,DIMENSION(:,:),ALLOCATABLE::A
    REAL,DIMENSION(:),ALLOCATABLE::B
    CHARACTER*(80):: inputfile = 'dengjin.txt'
    INTEGER DIM
    
CALL GET_DIM(inputfile,DIM)
    ALLOCATE(A(DIM,DIM),B(DIM))
   
 CALL INPUT_MATRIX(A,B,DIM,inputfile)
    
CALL LU(A,DIM)
    
CALL solve(A,B,B,DIM)
    PRINT*,B
END

SUBROUTINE GET_DIM(inputfile,n)
    CHARACTER*(*) inputfile
    CHARACTER*(80) str
    INTEGER n
    OPEN(10,FILE = inputfile)
    READ(10,*)str,n
    CLOSE(10)
END

SUBROUTINE INPUT_MATRIX(A,B,DIM,inputfile)
    CHARACTER*(*) inputfile
    INTEGER DIM
    INTEGER A(DIM,DIM),B(DIM)
    
    OPEN(20,FILE=inputfile)
    READ(20,*)
    READ(20,*)
    READ(20,*) ((A(I,J),J=1,DIM),I=1,DIM)
    READ(20,*)
    READ(20,*) (B(I),I=1,DIM)
    CLOSE(20)
END

SUBROUTINE LU(A,DIM)
    INTEGER DIM
    REAL A(DIM,DIM)
    
    DO I=2,DIM
        DO J=I,DIM
            A(J,I-1) = A(J,I-1)*1.0/A(I-1,I-1)
            m = A(J,I-1)
            DO K=I,DIM
                A(J,K) = A(J,K)-m*A(I-1,K)
            ENDDO
        ENDDO
    ENDDO

END

SUBROUTINE solve(A,B,X,DIM)
    INTEGER:: DIM
    REAL A(DIM,DIM),B(DIM),X(DIM)
!----------LY = B
    DO I=2,DIM
        DO J=1,I-1
            B(I) = B(I) - A(I,J)*B(J)
        ENDDO
    ENDDO

    B(DIM) = B(DIM)/A(DIM,DIM)
!---------UX=Y Y stores in B,same as X
    DO I=DIM-1,1,-1
        DO J=I+1,DIM
            B(I) = B(I) - A(I,J)*B(J)
        ENDDO
        B(I) = B(I)/A(I,I)
    ENDDO
END


阅读(9556) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~