LU decomposition: Difference between revisions

=={{header|Kotlin}}==

<lang scala>// version 1.1.4-3

typealias Vector = DoubleArray

typealias Matrix = Array<Vector>

operator fun Matrix.times(other: Matrix): Matrix {

val rows1 = this.size

val cols1 = this[0].size

val rows2 = other.size

val cols2 = other[0].size

require(cols1 == rows2)

val result = Matrix(rows1) { Vector(cols2) }

for (i in 0 until rows1) {

for (j in 0 until cols2) {

for (k in 0 until rows2) {

result[i][j] += this[i][k] * other[k][j]

}

}

}

return result

}

fun pivotize(m: Matrix): Matrix {

val n = m.size

val im = Array(n) { Vector(n) }

for (i in 0 until n) im[i][i] = 1.0

for (i in 0 until n) {

var max = m[i][i]

var row = i

for (j in i until n) {

if (m[j][i] > max) {

max = m[j][i]

row = j

}

}

if (i != row) {

val t = im[i]

im[i] = im[row]

im[row] = t

}

}

return im

}

fun lu(a: Matrix): Array<Matrix> {

val n = a.size

val l = Array(n) { Vector(n) }

val u = Array(n) { Vector(n) }

val p = pivotize(a)

val a2 = p * a

for (j in 0 until n) {

l[j][j] = 1.0

for (i in 0 until j + 1) {

var sum = 0.0

for (k in 0 until i) sum += u[k][j] * l[i][k]

u[i][j] = a2[i][j] - sum

}

for (i in j until n) {

var sum2 = 0.0

for(k in 0 until j) sum2 += u[k][j] * l[i][k]

l[i][j] = (a2[i][j] - sum2) / u[j][j]

}

}

return arrayOf(l, u, p)

}

fun printMatrix(title: String, m: Matrix, f: String) {

val n = m.size

println("\n$title\n")

for (i in 0 until n) {

for (j in 0 until n) print("${f.format(m[i][j])} ")

println()

}

}

fun main(args: Array<String>) {

val a1 = arrayOf(

doubleArrayOf( 1.0, 3.0, 5.0),

doubleArrayOf( 2.0, 4.0, 7.0),

doubleArrayOf( 1.0, 1.0, 0.0)

)

val (l1, u1, p1) = lu(a1)

println("EXAMPLE 1:-")

printMatrix("A:", a1, "%1.0f")

printMatrix("L:", l1, "% 7.5f")

printMatrix("U:", u1, "% 8.5f")

printMatrix("P:", p1, "%1.0f")

val a2 = arrayOf(

doubleArrayOf(11.0, 9.0, 24.0, 2.0),

doubleArrayOf( 1.0, 5.0, 2.0, 6.0),

doubleArrayOf( 3.0, 17.0, 18.0, 1.0),

doubleArrayOf( 2.0, 5.0, 7.0, 1.0)

)

val (l2, u2, p2) = lu(a2)

println("\nEXAMPLE 2:-")

printMatrix("A:", a2, "%2.0f")

printMatrix("L:", l2, "%7.5f")

printMatrix("U:", u2, "%8.5f")

printMatrix("P:", p2, "%1.0f")

}</lang>

{{out}}

<pre>

EXAMPLE 1:-

A:

1 3 5

2 4 7

1 1 0

L:

1.00000 0.00000 0.00000

0.50000 1.00000 0.00000

0.50000 -1.00000 1.00000

U:

2.00000 4.00000 7.00000

0.00000 1.00000 1.50000

0.00000 0.00000 -2.00000

P:

0 1 0

1 0 0

0 0 1

EXAMPLE 2:-

A:

11 9 24 2

1 5 2 6

3 17 18 1

2 5 7 1

L:

1.00000 0.00000 0.00000 0.00000

0.27273 1.00000 0.00000 0.00000

0.09091 0.28750 1.00000 0.00000

0.18182 0.23125 0.00360 1.00000

U:

11.00000 9.00000 24.00000 2.00000

0.00000 14.54545 11.45455 0.45455

0.00000 0.00000 -3.47500 5.68750

0.00000 0.00000 0.00000 0.51079

P:

1 0 0 0

0 0 1 0

0 1 0 0

0 0 0 1

</pre>