Just another update

Task/LU-decomposition/00DESCRIPTION

@@ -1,10 +1,16 @@
-Every square matrix <math>A</math> can be decomposed into a product of a lower triangular matrix <math>L</math> and a upper triangular matrix <math>U</math>, as described in [[wp:LU decomposition|LU decomposition]].
+Every square matrix <math>A</math> can be decomposed into a product of a lower triangular matrix <math>L</math> and a upper triangular matrix <math>U</math>,
+as described in [[wp:LU decomposition|LU decomposition]].
 
 :<math>A = LU</math>
 
-It is a modified form of Gaussian elimination. While the [[Cholesky decomposition]] only works for symmetric, positive definite matrices, the more general LU decomposition works for any square matrix.
+It is a modified form of Gaussian elimination.
+While the [[Cholesky decomposition]] only works for symmetric,
+positive definite matrices, the more general LU decomposition
+works for any square matrix.
 
-There are several algorithms for calculating L and U. To derive ''Crout's algorithm'' for a 3x3 example, we have to solve the following system:
+There are several algorithms for calculating L and U.
+To derive ''Crout's algorithm'' for a 3x3 example,
+we have to solve the following system:
 
 :<math>
 A =
@@ -88,7 +94,7 @@ and then for <math>L</math>
 
 :<math>l_{ij} = \frac{1}{u_{jj}} (a_{ij} - \sum_{k=1}^{j-1} u_{kj}l_{ik})</math>
 
-We see in the second formula that to get the <math>l_{ij}</math> below the diagonal, we have to divide by the diagonal element (pivot) <math>u_{ij}</math>, so we get problems when <math>u_{ij}</math> is either 0 or very small, which leads to numerical instability.
+We see in the second formula that to get the <math>l_{ij}</math> below the diagonal, we have to divide by the diagonal element (pivot) <math>u_{jj}</math>, so we get problems when <math>u_{jj}</math> is either 0 or very small, which leads to numerical instability.
 
 The solution to this problem is ''pivoting'' <math>A</math>, which means rearranging the rows of <math>A</math>, prior to the <math>LU</math> decomposition, in a way that the largest element of each column gets onto the diagonal of <math>A</math>. Rearranging the columns means to multiply <math>A</math> by a permutation matrix <math>P</math>:
 
@@ -119,7 +125,9 @@ The decomposition algorithm is then applied on the rearranged matrix so that
 
 '''Task description'''
 
-The task is to implement a routine which will take a square nxn matrix <math>A</math> and return a lower triangular matrix <math>L</math>, a upper triangular matrix <math>U</math> and a permutation matrix <math>P</math>, so that the above equation is fullfilled. You should then test it on the following two examples and include your output.
+The task is to implement a routine which will take a square nxn matrix <math>A</math> and return a lower triangular matrix <math>L</math>, a upper triangular matrix <math>U</math> and a permutation matrix <math>P</math>,
+so that the above equation is fullfilled.
+You should then test it on the following two examples and include your output.
 
 Example 1:
 <pre>

Task/LU-decomposition/D/lu-decomposition.d

@@ -1,11 +1,11 @@
 import std.stdio, std.algorithm, std.typecons, std.numeric,
        std.array, std.conv, std.string, std.range;
 
-bool isRectangular(T)(in T[][] m) pure nothrow {
+bool isRectangular(T)(in T[][] m) pure nothrow @nogc {
     return m.all!(r => r.length == m[0].length);
 }
 
-bool isSquare(T)(in T[][] m) pure nothrow {
+bool isSquare(T)(in T[][] m) pure nothrow @nogc {
     return m.isRectangular && m[0].length == m.length;
 }
 
@@ -34,7 +34,7 @@ in {
 } body {
     immutable n = m.length;
     auto id = iota(n)
-              .map!((in j) => n.iota.map!(i => cast(T)(i == j)).array)
+              .map!((in j) => n.iota.map!(i => T(i == j)).array)
               .array;
 
     foreach (immutable i; 0 .. n) {

Task/LU-decomposition/Fortran/lu-decomposition.f

@@ -1,75 +1,66 @@
 program lu1
-implicit none
-
-real*8     :: a1(3,3), a2(4,4)
-
-a1 = reshape((/real*8::1,2,1,3,4,1,5,7,0/),(/3,3/))
-call check(a1)
-
-a2 = reshape((/real*8::11,1,3,2,9,5,17,5,24,2,18,7,2,6,1,1/),(/4,4/))
-call check(a2)
+  implicit none
+  call check( reshape([real(8)::1,2,1,3,4,1,5,7,0                  ],[3,3]) )
+  call check( reshape([real(8)::11,1,3,2,9,5,17,5,24,2,18,7,2,6,1,1],[4,4]) )
 
 contains
 
-subroutine lu(a,p)
-! in situ decomposition, correspondes to LAPACK's dgebtrf
-implicit none
-real*8, intent(inout) :: a(:,:)
-integer, intent(out)  :: p(:)
-integer               :: n, i,j,k,ii
-n = ubound(a,1)
-p = (/ ( i, i=1,n ) /)
-do k = 1,n-1
-    ii = k-1+maxloc(abs(a(p(k:),k)),1)
-    if (ii /= k )  then
-        p((/k, ii/)) = p((/ii, k/))
-    end if
-    a(p(k+1:),k) = a(p(k+1:),k)/a(p(k),k)
-    forall (j = k+1:n)
-        a(p(k+1:),j) = a(p(k+1:),j)-a(p(k+1:),k)*a(p(k),j)
+  subroutine check(a)
+    real(8), intent(in)   :: a(:,:)
+    integer               :: i,j,n
+    real(8), allocatable  :: aa(:,:),l(:,:),u(:,:)
+    integer, allocatable  :: p(:,:)
+    integer, allocatable  :: ipiv(:)
+    n = size(a,1)
+    allocate(aa(n,n),l(n,n),u(n,n),p(n,n),ipiv(n))
+    forall (j=1:n,i=1:n)
+      aa(i,j) = a(i,j)
+      u (i,j) = 0d0
+      p (i,j) = merge(1  ,0  ,i.eq.j)
+      l (i,j) = merge(1d0,0d0,i.eq.j)
     end forall
-end do
-end subroutine
+    call lu(aa, ipiv)
+    do i = 1,n
+       l(i, :i-1) = aa(ipiv(i), :i-1)
+       u(i,i:   ) = aa(ipiv(i),i:   )
+    end do
+    p(ipiv,:) = p
+    call mat_print('a',a)
+    call mat_print('p',p)
+    call mat_print('l',l)
+    call mat_print('u',u)
+    print *, "residual"
+    print *, "|| P.A - L.U || =  ", maxval(abs(matmul(p,a)-matmul(l,u)))
+  end subroutine
 
-subroutine check(a)
-implicit none
-real*8, intent(in) :: a(:,:)
-real*8             :: aa(ubound(a,1), ubound(a,2))
-real*8             :: l(ubound(a,1), ubound(a,2))
-real*8             :: u(ubound(a,1), ubound(a,2))
-integer            :: p(ubound(a,1), ubound(a,2)), ipiv(ubound(a,1))
-integer            :: i,  n
-character(len=100) :: fmt
+  subroutine lu(a,p)
+!   in situ decomposition, corresponds to LAPACK's dgebtrf
+    real(8), intent(inout) :: a(:,:)
+    integer, intent(out  ) :: p(:)
+    integer                :: n, i,j,k,kmax
+    n = size(a,1)
+    p = [ ( i, i=1,n ) ]
+    do k = 1,n-1
+        kmax = maxloc(abs(a(p(k:),k)),1) + k-1
+        if (kmax /= k ) p([k, kmax]) = p([kmax, k])
+        a(p(k+1:),k) = a(p(k+1:),k) / a(p(k),k)
+        forall (j=k+1:n) a(p(k+1:),j) = a(p(k+1:),j) - a(p(k+1:),k) * a(p(k),j)
+    end do
+  end subroutine
 
-n = ubound(a,1)
-aa = a ! work with a copy
-p = 0; l=0; u = 0
-forall (i=1:n)
-    p(i,i) = 1d0; l(i,i) = 1d0 ! convert permutation vector a matrix
-end forall
+  subroutine mat_print(amsg,a)
+    character(*), intent(in) :: amsg
+    class    (*), intent(in) :: a(:,:)
+    integer                  :: i
+    print*,' '
+    print*,amsg
+    do i=1,size(a,1)
+      select type (a)
+        type is (real(8)) ; print'(100f8.2)',a(i,:)
+        type is (integer) ; print'(100i8  )',a(i,:)
+      end select
+    end do
+    print*,' '
+  end subroutine
 
-call lu(aa, ipiv)
-do i = 1,n
-   l(i,:i-1) = aa(ipiv(i),:i-1)
-   u(i,i:) = aa(ipiv(i),i:)
-end do
-p(ipiv,:) = p
-write (fmt,"(a,i1,a)") "(",n,"(f8.2,1x))"
-
-print *, "a"
-print fmt, transpose(a)
-
-print *, "p"
-print fmt, transpose(dble(p))
-
-print *, "l"
-print fmt, transpose(l)
-
-print *, "u"
-print fmt, transpose(u)
-
-print *, "residual"
-print *, "|| P.A - L.U || =  ", maxval(abs(matmul(p,a)-matmul(l,u)))
-
-end subroutine
-end program>
+end program

Task/LU-decomposition/Ruby/lu-decomposition-1.rb

@@ -10,10 +10,10 @@ class Matrix
       (0 ... row_size).each do |j|
         if j >= i
           # upper
-          u[i][j] = tmp[i,j] - (0 .. i-1).inject(0.0) {|sum, k| sum + u[k][j] * l[i][k]}
+          u[i][j] = tmp[i,j] - (0 ... i).inject(0.0) {|sum, k| sum + u[k][j] * l[i][k]}
         else
           # lower
-          l[i][j] = (tmp[i,j] - (0 .. j-1).inject(0.0) {|sum, k| sum + u[k][j] * l[i][k]}) / u[j][j]
+          l[i][j] = (tmp[i,j] - (0 ... j).inject(0.0) {|sum, k| sum + u[k][j] * l[i][k]}) / u[j][j]
         end
       end
     end
@@ -37,29 +37,29 @@ class Matrix
     Matrix[*id]
   end
 
-  def pretty_print(format)
-    each_with_index do |val, i, j|
-      print "#{format} " % val
-      puts "" if j==column_size-1
-    end
+  def pretty_print(format, head=nil)
+    puts head if head
+    puts each_slice(column_size).map{|row| format*row_size % row}
   end
 end
 
+puts "Example 1:"
 a = Matrix[[1,  3,  5],
            [2,  4,  7],
            [1,  1,  0]]
-puts "A"; a.pretty_print("%2d")
+a.pretty_print(" %2d", "A")
 l, u, p = a.lu_decomposition
-puts "U"; u.pretty_print("%8.5f")
-puts "L"; l.pretty_print("%8.5f")
-puts "P"; p.pretty_print("%d")
+l.pretty_print(" %8.5f", "L")
+u.pretty_print(" %8.5f", "U")
+p.pretty_print(" %d",    "P")
 
+puts "\nExample 2:"
 a = Matrix[[11, 9,24,2],
            [ 1, 5, 2,6],
            [ 3,17,18,1],
            [ 2, 5, 7,1]]
-puts "A"; a.pretty_print("%2d")
+a.pretty_print(" %2d", "A")
 l, u, p = a.lu_decomposition
-puts "U"; u.pretty_print("%8.5f")
-puts "L"; l.pretty_print("%8.5f")
-puts "P"; p.pretty_print("%d")
+l.pretty_print(" %8.5f", "L")
+u.pretty_print(" %8.5f", "U")
+p.pretty_print(" %d",    "P")

Task/LU-decomposition/Ruby/lu-decomposition-2.rb

@@ -0,0 +1,4 @@
+l, u, p = a.lup_decomposition
+l.pretty_print(" %8.5f", "L")
+u.pretty_print(" %8.5f", "U")
+p.pretty_print(" %d",    "P")