|
| 1 | +package vecxt |
| 2 | + |
| 3 | +import dev.ludovic.netlib.lapack.JavaLAPACK |
| 4 | +import org.netlib.util.intW |
| 5 | +import vecxt.matrix.Matrix |
| 6 | +import vecxt.MatrixInstance.* |
| 7 | +import vecxt.MatrixHelper.zeros |
| 8 | +import vecxt.BoundsCheck.BoundsCheck |
| 9 | + |
| 10 | +/** LU decomposition with partial pivoting using LAPACK. |
| 11 | + * |
| 12 | + * Decomposes a matrix A into the product P*A = L*U where: |
| 13 | + * - P is a permutation matrix |
| 14 | + * - L is a lower triangular matrix with unit diagonal |
| 15 | + * - U is an upper triangular matrix |
| 16 | + */ |
| 17 | +object LU: |
| 18 | + private lazy final val lapack = JavaLAPACK.getInstance() |
| 19 | + |
| 20 | + /** Computes the LU decomposition with partial pivoting of a matrix. |
| 21 | + * |
| 22 | + * Uses LAPACK's dgetrf routine to compute the factorization P*A = L*U where: |
| 23 | + * - P is a permutation matrix (represented as a pivot array) |
| 24 | + * - L is a lower triangular matrix with ones on the diagonal |
| 25 | + * - U is an upper triangular matrix |
| 26 | + * |
| 27 | + * The decomposition allows efficient solving of linear systems and computation of determinants. |
| 28 | + * |
| 29 | + * Note: The factorization is computed even for singular matrices. In such cases, at least one diagonal element of U |
| 30 | + * will be zero (or near zero), indicating the matrix is rank-deficient. |
| 31 | + * |
| 32 | + * @param m |
| 33 | + * The input matrix to decompose. Can be square or rectangular (m x n). |
| 34 | + * @return |
| 35 | + * A named tuple containing: |
| 36 | + * - L: Lower triangular matrix with unit diagonal (rows × min(rows,cols)) |
| 37 | + * - U: Upper triangular matrix (min(rows,cols) × cols) |
| 38 | + * - P: Permutation array representing successive row interchanges. At step i, row i was swapped with row P(i). The |
| 39 | + * permutation matrix can be reconstructed by applying these swaps in order. |
| 40 | + * @throws IllegalArgumentException |
| 41 | + * if matrix is empty or if an argument to LAPACK is invalid |
| 42 | + */ |
| 43 | + inline def lu(m: Matrix[Double])(using |
| 44 | + inline bc: BoundsCheck |
| 45 | + ): (L: Matrix[Double], U: Matrix[Double], P: Array[Int]) = |
| 46 | + nonEmptyMatCheck(m) |
| 47 | + |
| 48 | + val rows = m.rows |
| 49 | + val cols = m.cols |
| 50 | + val minDim = math.min(rows, cols) |
| 51 | + |
| 52 | + // Copy input matrix as LAPACK will overwrite it |
| 53 | + val A = m.deepCopy |
| 54 | + |
| 55 | + // Allocate pivot array - LAPACK uses 1-based indexing |
| 56 | + val ipiv = Array.ofDim[Int](minDim) |
| 57 | + |
| 58 | + val info = new intW(0) |
| 59 | + |
| 60 | + // Compute LU factorization |
| 61 | + lapack.dgetrf( |
| 62 | + rows, // M: number of rows |
| 63 | + cols, // N: number of columns |
| 64 | + A.raw, // A: matrix to factor (overwritten with L and U) |
| 65 | + math.max(1, rows), // LDA: leading dimension |
| 66 | + ipiv, // IPIV: pivot indices |
| 67 | + info // INFO: status |
| 68 | + ) |
| 69 | + |
| 70 | + // Check for errors |
| 71 | + if info.`val` < 0 then |
| 72 | + throw new IllegalArgumentException(s"LU decomposition failed: argument ${-info.`val`} had an illegal value") |
| 73 | + end if |
| 74 | + |
| 75 | + // Note: info.val > 0 means U(info.val, info.val) is exactly zero |
| 76 | + // The factorization has been completed, but U is singular |
| 77 | + // We allow this and let the caller handle it |
| 78 | + |
| 79 | + // Extract L and U from the factored matrix |
| 80 | + // LAPACK stores L and U in the same matrix: |
| 81 | + // - U is in the upper triangle (including diagonal) |
| 82 | + // - L is in the lower triangle (diagonal is implicitly 1) |
| 83 | + |
| 84 | + val L = Matrix.zeros[Double](rows, minDim) |
| 85 | + val U = Matrix.zeros[Double](minDim, cols) |
| 86 | + |
| 87 | + // Extract L (lower triangle with unit diagonal) |
| 88 | + var i = 0 |
| 89 | + while i < rows do |
| 90 | + var j = 0 |
| 91 | + while j < minDim do |
| 92 | + if i > j then |
| 93 | + // Below diagonal: copy from A |
| 94 | + L(i, j) = A(i, j) |
| 95 | + else if i == j then |
| 96 | + // On diagonal: set to 1 |
| 97 | + L(i, j) = 1.0 |
| 98 | + // else: above diagonal is 0 (already initialized) |
| 99 | + end if |
| 100 | + j += 1 |
| 101 | + end while |
| 102 | + i += 1 |
| 103 | + end while |
| 104 | + |
| 105 | + // Extract U (upper triangle) |
| 106 | + i = 0 |
| 107 | + while i < minDim do |
| 108 | + var j = 0 |
| 109 | + while j < cols do |
| 110 | + if i <= j then |
| 111 | + // On or above diagonal: copy from A |
| 112 | + U(i, j) = A(i, j) |
| 113 | + // else: below diagonal is 0 (already initialized) |
| 114 | + end if |
| 115 | + j += 1 |
| 116 | + end while |
| 117 | + i += 1 |
| 118 | + end while |
| 119 | + |
| 120 | + // Convert LAPACK's 1-based pivot indices to 0-based |
| 121 | + val pivotArray = Array.ofDim[Int](minDim) |
| 122 | + i = 0 |
| 123 | + while i < minDim do |
| 124 | + pivotArray(i) = ipiv(i) - 1 // Convert from 1-based to 0-based |
| 125 | + i += 1 |
| 126 | + end while |
| 127 | + |
| 128 | + (L = L, U = U, P = pivotArray) |
| 129 | + end lu |
| 130 | + |
| 131 | +end LU |
0 commit comments