next | previous | forward | backward | up | top | index | toc | Macaulay2 web site

LUdecomposition -- LU decomposition

Synopsis

Description

The output matrices are mutable exactly when the input matrix is, but the matrix A is not modified

If Q is the m by m permutation matrix such that Q_(P_i,i) = 1, and all other entries are zero, then A = QLU.

There are several restrictions. The first is that there are only a limited number of rings for which this function is implemented. Second, if A is a mutable matrix defined over RR or CC, then A must be densely encoded. This restriction is not present if A is a matrix.

i1 : kk = ZZ/101

o1 = kk

o1 : QuotientRing
i2 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk

o2 = | 1  2 3  4  |
     | 1  3 6  10 |
     | 19 7 11 13 |

              3        4
o2 : Matrix kk  <--- kk
i3 : (P,L,U) = LUdecomposition A

o3 = ({0, 1, 2}, | 1  0   0 |, | 1 2 3  4  |)
                 | 1  1   0 |  | 0 1 3  6  |
                 | 19 -31 1 |  | 0 0 47 22 |

o3 : Sequence
i4 : Q = id_(kk^3) _ P

o4 = | 1 0 0 |
     | 0 1 0 |
     | 0 0 1 |

              3        3
o4 : Matrix kk  <--- kk
i5 : Q * L * U == matrix A

o5 = true
For matrices over RR or CC, this function calls the lapack routines, which are restricted to 53 bits of precision.
i6 : A = matrix"1,2,3,4,5,6;1,3,6,12,13,16;19,7,11,47,48,21" ** RR

o6 = | 1  2 3  4  5  6  |
     | 1  3 6  12 13 16 |
     | 19 7 11 47 48 21 |

                3          6
o6 : Matrix RR    <--- RR
              53         53
i7 : (P,L,U) = LUdecomposition A

o7 = ({2, 1, 0}, | 1        0   0 |, | 19 7       11      47      48     
                 | .0526316 1   0 |  | 0  2.63158 5.42105 9.52632 10.4737
                 | .0526316 .62 1 |  | 0  0       -.94    -4.38   -4.02  
     ------------------------------------------------------------------------
     21      |)
     14.8947 |
     -4.34   |

o7 : Sequence
i8 : Q = id_ (RR^3) _ P

o8 = | 0 0 1 |
     | 0 1 0 |
     | 1 0 0 |

                3          3
o8 : Matrix RR    <--- RR
              53         53
i9 : Q * L * U - A

o9 = | 0 -2.22045e-16 0 0 0 0 |
     | 0 0            0 0 0 0 |
     | 0 0            0 0 0 0 |

                3          6
o9 : Matrix RR    <--- RR
              53         53
i10 : clean(1e-15,oo)

o10 = 0

                 3          6
o10 : Matrix RR    <--- RR
               53         53
Mutable matrices can sometimes be useful for speed, and/or space. If A is a mutable matrix, it must be densely encoded (which is the default).
i11 : A = mutableMatrix(CC,5,10, Dense=>true)

o11 = 0

o11 : MutableMatrix
i12 : printingPrecision = 2

o12 = 2
i13 : setRandomSeed 0

o13 = 0
i14 : fillMatrix A

o14 = | .83+.27i  .64+.14i  .014+.33i .12+.86i .76+.02i .41+.088i .12+.78i .98+.79i .4+.12i    .06+.57i |
      | .44+.48i  .45+.077i .12+.72i  .8+.77i  .19+.73i .014+.18i .97+.42i .05+.55i .18+.53i   1+.32i   |
      | .084+.19i .4+.52i   .13+.28i  .88+.03i .3+.49i  .86+.04i  .7+.5i   .7+.78i  .37+.41i   .56+.16i |
      | .83+.04i  .34+.26i  .39+.97i  .4+.6i   .77+.28i .9+.84i   .73+.88i .23+.99i .56+.89i   .3+i     |
      | .88+.44i  .24+.77i  .82+.62i  .32+.7i  .39+.49i .58+.45i  .94+.87i .19+.85i .093+.017i .34+.83i |

o14 : MutableMatrix
i15 : (P,L,U) = LUdecomposition A;
i16 : Q = id_(CC^5) _ P

o16 = | 0 0 0 0 1 |
      | 0 1 0 0 0 |
      | 0 0 0 1 0 |
      | 0 0 1 0 0 |
      | 1 0 0 0 0 |

                 5          5
o16 : Matrix CC    <--- CC
               53         53
i17 : matrix Q * matrix L * matrix U - matrix A

o17 = | -1.1e-16          0        -2.2e-16-5.6e-17i -1.1e-16+1.1e-16i
      | 0                 0        5.6e-17-1.1e-16i  0                
      | -2.8e-17i         1.1e-16i 5.6e-17           1.1e-16          
      | -1.1e-16-5.6e-17i 0        0                 -5.6e-17         
      | 0                 0        0                 0                
      -----------------------------------------------------------------------
      -1.1e-16i        1.1e-16i -2.2e-16+1.1e-16i 0         5.6e-17i         
      0                5.6e-17i -1.1e-16i         0         0                
      5.6e-17-1.1e-16i 2.8e-17i 0                 -1.1e-16  -5.6e-17+5.6e-17i
      0                1.1e-16  -2.2e-16-1.1e-16i -1.1e-16i 0                
      0                0        0                 0         0                
      -----------------------------------------------------------------------
      -2.2e-16+1.1e-16i |
      1.1e-16-1.1e-16i  |
      -5.6e-17i         |
      -1.1e-16i         |
      0                 |

                 5          10
o17 : Matrix CC    <--- CC
               53         53
i18 : clean(1e-15,oo)

o18 = 0

                 5          10
o18 : Matrix CC    <--- CC
               53         53

Caveat

This function is limited in scope, but is sometimes useful for very large matrices

See also

Ways to use LUdecomposition :