//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** @author John Miller * @version 1.0 * @date Sun Sep 16 14:09:25 EDT 2012 * @see LICENSE (MIT style license file). */ package scalation.linalgebra import scalation.math.Complex import scalation.math.Complex.abs import scalation.util.Error //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `Matric` trait specifies the operations to be defined by three concrete * implemeting classes: MatrixC, SparseMatrixC and SymTriMatrixC. */ trait Matric extends Error { /** Matrix dimension 1 (# rows) */ val dim1: Int /** Matrix dimension 2 (# columns) */ val dim2: Int /** Range for the storage array on dimension 1 (rows) */ protected val range1 = 0 until dim1 /** Range for the storage array on dimension 2 (columns) */ protected val range2 = 0 until dim2 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get this matrix's element at the i,j-th index position. * @param i the row index * @param j the column index */ def apply (i: Int, j: Int): Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get this matrix's vector at the i-th index position (i-th row). * @param i the row index */ def apply (i: Int): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get a slice this matrix row-wise on range ir and column-wise on range jr. * Ex: b = a(2..4, 3..5) * @param ir the row range * @param jr the column range */ def apply (ir: Range, jr: Range): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get a slice this matrix row-wise on range ir and column-wise at index j. * Ex: u = a(2..4, 3) * @param ir the row range * @param j the column index */ def apply (ir: Range, j: Int): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get a slice this matrix row-wise at index i and column-wise on range jr. * Ex: u = a(2, 3..5) * @param i the row index * @param jr the column range */ def apply (i: Int, jr: Range): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set this matrix's element at the i,j-th index position to the scalar x. * @param i the row index * @param j the column index * @param x the scalar value to assign */ def update (i: Int, j: Int, x: Complex) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set this matrix's row at the i-th index position to the vector u. * @param i the row index * @param u the vector value to assign */ def update (i: Int, u: VectorC) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set a slice this matrix row-wise on range ir and column-wise at index j. * Ex: a(2..4, 3) = u * @param ir the row range * @param j the column index * @param u the vector to assign */ def update (ir: Range, j: Int, u: VectorC) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set a slice this matrix row-wise at index i and column-wise on range jr. * Ex: a(2, 3..5) = u * @param i the row index * @param jr the column range * @param u the vector to assign */ def update (i: Int, jr: Range, u: VectorC) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set all the elements in this matrix to the scalar x. * @param x the scalar value to assign */ def set (x: Complex) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set all the elements in this matrix to the scalar x. * @param x the scalar value to assign */ def set (x: Double) { set (Complex (x)) } //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set the values in this matrix as copies of the values in 2D array u. * @param u the 2D array of values to assign */ def set (u: Array [Array [Complex]]) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set this matrix's ith row starting a column j to the vector u. * @param i the row index * @param u the vector value to assign * @param j the starting column index */ def set (i: Int, u: VectorC, j: Int = 0) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set all the lower triangular elements in this matrix to the scalar x. * @param x the scalar value to assign */ def setLower (x: Complex) { for (i <- range1; j <- 0 until i) this(i, j) = x } //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Iterate over the matrix row by row. * @param f the function to apply */ def foreach [U] (f: Array [Complex] => U) { var i = 0 while (i < dim1) { f (this(i)()); i += 1 } } // foreach //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Slice this matrix row-wise from to end. * @param from the start of the slice * @param end the end of the slice */ def slice (from: Int, end: Int): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Slice this matrix row-wise r_from to r_end and column-wise c_from to c_end. * @param r_from the start of the row slice * @param r_end the end of the row slice * @param c_from the start of the column slice * @param c_end the end of the column slice */ def slice (r_from: Int, r_end: Int, c_from: Int, c_end: Int): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Slice this matrix excluding the given row and column. * @param row the row to exclude * @param col the column to exclude */ def sliceExclude (row: Int, col: Int): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Select rows from this matrix according to the given index/basis. * @param rowIndex the row index positions (e.g., (0, 2, 5)) */ def selectRows (rowIndex: Array [Int]): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get column 'col' from the matrix, returning it as a vector. * @param col the column to extract from the matrix * @param from the position to start extracting from */ def col (col: Int, from: Int = 0): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set column 'col' of the matrix to a vector. * @param col the column to set * @param u the vector to assign to the column */ def setCol (col: Int, u: VectorC) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Select columns from this matrix according to the given index/basis. * Ex: Can be used to divide a matrix into a basis and a non-basis. * @param colIndex the column index positions (e.g., (0, 2, 5)) */ def selectCols (colIndex: Array [Int]): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Return 1 if the condition is true else 0 * @param cond the condition to evaluate */ def oneIf (cond: Boolean): Int = if (cond) 1 else 0 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Transpose this matrix (rows => columns). */ def t: Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Concatenate this matrix and vector u. * @param u the vector to be concatenated as the new last row in matrix */ def ++ (u: VectorC): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Add this matrix and scalar x. * @param x the scalar to add */ def + (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Add in-place this matrix and scalar x. * @param x the scalar to add */ def += (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** From this matrix subtract scalar x. * @param x the scalar to subtract */ def - (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** From this matrix subtract in-place scalar x. * @param x the scalar to subtract */ def -= (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by vector u. * @param u the vector to multiply by */ def * (u: VectorC): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by scalar x. * @param x the scalar to multiply by */ def * (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply in-place this matrix by scalar x. * @param x the scalar to multiply by */ def *= (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by vector u to produce another matrix (a_ij * b_j) * @param b the vector to multiply by */ def ** (u: VectorC): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply in-place this matrix by vector u to produce another matrix (a_ij * b_j) * @param u the vector to multiply by */ def **= (u: VectorC): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Divide this matrix by scalar x. * @param x the scalar to divide by */ def / (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Divide in-place this matrix by scalar x. * @param x the scalar to divide by */ def /= (x: Complex): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Raise this matrix to the pth power (for some integer p >= 2). * @param p the power to raise this matrix to */ def ~^ (p: Int): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Find the maximum element in this matrix. * @param e the ending row index (exclusive) for the search */ def max (e: Int = dim1): Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Find the minimum element in this matrix. * @param e the ending row index (exclusive) for the search */ def min (e: Int = dim1): Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Find the magnitude of this matrix, the element value farthest from zero. */ def mag: Complex = abs (max ()) max abs (min ()) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Decompose this matrix into the product of lower and upper triangular * matrices (l, u) using the LU Decomposition algorithm. This version uses * partial pivoting. */ def lud: Tuple2 [Matric, Matric] //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Decompose in-place this matrix into the product of lower and upper triangular * matrices (l, u) using the LU Decomposition algorithm. This version uses * partial pivoting. */ def lud_ip: Tuple2 [Matric, Matric] //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Solve for x in the equation l*u*x = b (see lud above). * @param l the lower triangular matrix * @param u the upper triangular matrix * @param b the constant vector */ def solve (l: Matric, u: Matric, b: VectorC): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Solve for x in the equation l*u*x = b (see lud above). * @param lu the lower and upper triangular matrices * @param b the constant vector */ def solve (lu: Tuple2 [Matric, Matric], b: VectorC): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Solve for x in the equation a*x = b where a is this matrix (see lud above). * @param b the constant vector. */ def solve (b: VectorC): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Determine the rank of this m by n matrix by taking the upper triangular * matrix from the LU Decomposition and counting the number of non-zero * diagonal elements. */ def rank: Int = { val max = if (dim1 < dim2) dim1 else dim2 // rank <= min (m, n) val u = lud._2 // upper triangular matrix var count = 0 for (i <- 0 until max if this(i, i) != 0.0) count += 1 count } // rank //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Form a matrix [Ip, this, Iq] where Ir is a r by r identity matrix, by * positioning the three matrices Ip, this and Iq along the diagonal. * @param p the size of identity matrix Ip * @param q the size of identity matrix Iq */ def diag (p: Int, q: Int): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get the kth diagonal of this matrix. Assumes dim2 >= dim1. * @param k how far above the main diagonal, e.g., (-1, 0, 1) for (sub, main, super) */ def getDiag (k: Int = 0): VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set the kth diagonal of this matrix to the vector u. Assumes dim2 >= dim1. * @param u the vector to set the diagonal to * @param k how far above the main diagonal, e.g., (-1, 0, 1) for (sub, main, super) */ def setDiag (u: VectorC, k: Int = 0) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set the main diagonal of this matrix to the scalar x. Assumes dim2 >= dim1. * @param x the scalar to set the diagonal to */ def setDiag (x: Complex) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Invert this matrix (requires a squareMatrix) and use partial pivoting. */ def inverse: Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Invert in-place this matrix (requires a squareMatrix) and use partial pivoting. */ def inverse_ip: Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Use Gauss-Jordan reduction on this matrix to make the left part embed an * identity matrix. A constraint on this m by n matrix is that n >= m. */ def reduce: Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Use Gauss-Jordan reduction in-place on this matrix to make the left part * embed an identity matrix. A constraint on this m by n matrix is that n >= m. */ def reduce_ip //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Clean values in matrix at or below the threshold by setting them to zero. * Iterative algorithms give approximate values and if very close to zero, * may throw off other calculations, e.g., in computing eigenvectors. * @param thres the cutoff threshold (a small value) * @param relative whether to use relative or absolute cutoff */ def clean (thres: Double, relative: Boolean = true): Matric //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the (right) nullspace of this m by n matrix (requires n = m + 1) * by performing Gauss-Jordan reduction and extracting the negation of the * last column augmented by 1. The nullspace of matrix a is "this vector v * times any scalar s", i.e., a*(v*s) = 0. The left nullspace of matrix a is * the same as the right nullspace of a.t (a transpose). */ def nullspace: VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the (right) nullspace in-place of this m by n matrix (requires n = m + 1) * by performing Gauss-Jordan reduction and extracting the negation of the * last column augmented by 1. The nullspace of matrix a is "this vector v * times any scalar s", i.e., a*(v*s) = 0. The left nullspace of matrix a is * the same as the right nullspace of a.t (a transpose). */ def nullspace_ip: VectorC //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the trace of this matrix, i.e., the sum of the elements on the * main diagonal. Should also equal the sum of the eigenvalues. * @see Eigen.scala */ def trace: Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the sum of this matrix, i.e., the sum of its elements. */ def sum: Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the sum of the lower triangular region of this matrix. */ def sumLower: Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the abs sum of this matrix, i.e., the sum of the absolute value * of its elements. This is useful for comparing matrices (a - b).sumAbs */ def sumAbs: Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the 1-norm of this matrix, i.e., the maximum 1-norm of the * column vectors. This is useful for comparing matrices (a - b).norm1 */ def norm1: Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the determinant of this matrix. */ def det: Complex //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix and the other Matrix have the same dimensions. * @param b the other matrix */ def sameDimensions (b: Matric): Boolean = dim1 == b.dim1 && dim2 == b.dim2 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix dimensions are less than or equal to (le) those * of the other Matrix. * @param b the other matrix */ def leDimensions (b: Matric): Boolean = dim1 <= b.dim1 && dim2 <= b.dim2 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix and the other matrix have the same cross dimensions. * @param b the other matrix */ def sameCrossDimensions (b: Matric): Boolean = dim2 == b.dim1 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix is rectangular (all rows have the same number * of columns). */ def isRectangular: Boolean //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix is square (same row and column dimensions). */ def isSquare: Boolean = dim1 == dim2 && isRectangular //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix is nonnegative (has no negative elements). */ def isNonnegative: Boolean //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix is symmetric. */ def isSymmetric: Boolean = { for (i <- 0 to dim1 - 2; j <- i + 1 until dim2 if this(i, j) != this(j, i)) return false true } // isSymmetric } // Matric class