//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /**::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: * @author Robert Davis, John Miller * @version 1.0 * @date Sun Sep 16 14:09:25 EDT 2012 * @see LICENSE (MIT style license file). */ package scalation.linalgebra import math.{max, min} import scalation.math.Complex import scalation.math.Complex.{abs, _0} import scalation.util.Error //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SymTriMatrixC` class stores and operates on symmetric tridiagonal matrices. * The elements are of type of Complex. A matrix is stored as two vectors: * the diagonal vector and the sub-diagonal vector. * @param d1 the first/row dimension (symmetric => d2 = d1) */ class SymTriMatrixC (val d1: Int) extends Matric with Error with Serializable { // Note: implementations for the following methods are from the Matric trait: // foreach, mag, rank, sameDimensions, leDimensions, sameCrossDimensions, // isSquare, isSymmetric /** Dimension 1 */ lazy val dim1 = d1 /** Dimension 2 */ lazy val dim2 = d1 /** Size of the sub-diagonal */ private val d1_1 = d1 - 1 /** Range for the diagonal */ private val range_d = 0 until d1 /** Range for the sub-diagonal */ private val range_s = 0 until d1_1 /** Diagonal of the matrix */ private var _dg: VectorC = new VectorC (d1) /** Sub-diagonal (also same for sup-diagonal) of the matrix */ private var _sd: VectorC = new VectorC (d1_1) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Construct a symmetric tridiagonal matrix with the given diagonal and sub-diagonal. * @param v1 the diagonal vector * @param v2 the sub-diagonal vector */ def this (v1: VectorC, v2: VectorC) { this (v1.dim) for (i <- range_d) _dg(i) = v1(i) for (i <- range_s) _sd(i) = v2(i) } // constructor //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Construct a symmetric tridiagonal matrix from the given matrix. * @param a the matrix of values to assign */ def this (a: Matric) { this (a.dim1) for (i <- range_d) _dg(i) = a(i, i) for (i <- range_s) _sd(i) = a(i, i+1) } // constructor //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get the diagonal of the matrix. */ def dg: VectorC = _dg //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set the diagonal of the matrix. * @param v the vector to assign to the diagonal */ def dg_ (v: VectorC) { _dg = v } //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get the sub-diagonal of the matrix. */ def sd: VectorC = _sd //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set the sub-diagonal of the matrix. * @param v the vector to assign to the sub-diagonal */ def sd_ (v: VectorC) { _sd = v } //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = { if (i == j) _dg(i) // on diagonal else if (i == j + 1) _sd(j) // on sub-diagonal (below diagonal) else if (i + 1 == j) _sd(i) // on sup-diagonal (above diagonal) else throw new Exception ("apply: element not on tridiagonal") } // apply //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get this matrix's element at the i,j-th index position, returning 0. if * off tridiagonal. * @param i the row index * @param j the column index */ def at (i: Int, j: Int): Complex = { if (i < 0 || j < 0 || i >= d1 || j >= d1) _0 else if (i == j) _dg(i) // on diagonal else if (i == j + 1) _sd(j) // on sub-diagonal (below diagonal) else if (i + 1 == j) _sd(i) // on sup-diagonal (above diagonal) else _0 } // at //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Get this matrix's vector at the i-th index position (i-th row). * @param i the row index */ def apply (i: Int): VectorC = { val u = new VectorC (d1) u(i) = _dg(i) if (i > 0) u(i-1) = _sd(i-1) if (i < d1_1) u(i+1) = _sd(i) u } // apply //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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): SymTriMatrixC = { if (ir != jr) flaw ("apply", "requires same ranges to maintain symmetry") slice (ir.start, ir.end) } // apply //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = col(j)(ir) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = this(i)(jr) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set this matrix's element at the i,j-th index position to the scalar x. * Only store x if it is non-zero. * @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) { if (i == j) _dg(i) = x else if (i == j + 1) _sd(j) = x else if (i + 1 == j) _sd(i) = x else flaw ("update", "element not on tridiagonal") } // update //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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) { _dg(i) = u(i) if (i > 0) _sd(i-1) = u(i-1) if (i < d1_1) _sd(i) = u(i+1) } // update //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set a slice this matrix row-wise on range ir and column-wise on range jr. * Ex: a(2..4, 3..5) = b * @param ir the row range * @param jr the column range * @param b the matrix to assign */ def update (ir: Range, jr: Range, b: SymTriMatrixC) { if (ir != jr) flaw ("update", "requires same ranges to maintain symmetry") for (i <- ir) { _dg(i) = b.dg(i - ir.start) if (i > ir.start) _sd(i-1) = b.sd(i - ir.start - 1) } // for } // update //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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) { col(j)(ir) = u } //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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) { this(i)(jr) = u } //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set all the elements in this matrix to the scalar x. * @param x the scalar value to assign */ def set (x: Complex) { for (i <- range1) { _dg(i) = x if (i > 0) _sd(i) = x } // for } // set //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set all 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]]) { throw new NoSuchMethodException ("values for SymTriMatrixC should be diagonal") } // set //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set this matrix's ith row starting at 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) { throw new NoSuchMethodException ("values for SymTriMatrixC should be diagonal") } // set //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Slice this matrix row-wise from to end. * @param from the start row of the slice (inclusive) * @param end the end row of the slice (exclusive) */ def slice (from: Int, end: Int): SymTriMatrixC = { val c = new SymTriMatrixC (end - from) for (i <- c.range1) { c._dg(i) = _dg(i + from) if (i > 0) c._sd(i - 1) = _sd(i + from - 1) } // for c } // slice //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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): SymTriMatrixC = { throw new NoSuchMethodException ("SymTriMatrixC must be symmetric") } // slice //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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): SymTriMatrixC = { throw new NoSuchMethodException ("SymTriMatrixC does not support sliceExclude") } // sliceExclude //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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]): SymTriMatrixC = { throw new NoSuchMethodException ("SymTriMatrixC does not support selectRows") } // selectRows //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = { val u = new VectorC (dim1 - from) for (i <- (from max col-1) until (dim1 min col+2)) u(i-from) = this(i, col) u } // col //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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) { _dg(col) = u(col) if (col > 0) _sd(col-1) = u(col-1) } // setCol //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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]): SymTriMatrixC = { throw new NoSuchMethodException ("SymTriMatrixC does not support selectCol") } // selectCols //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Transpose this matrix (rows => columns). Note, since the matrix is * symmetric, it returns itself. */ def t: SymTriMatrixC = this //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Concatenate this matrix and vector u. * @param u the vector to be concatenated as the new last row in matrix */ def ++ (u: VectorC): SymTriMatrixC = { throw new NoSuchMethodException ("SymTriMatrixC does not support ++") } // ++ //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Add this matrix and matrix b. * @param b the matrix to add (requires leDimensions) */ def + (b: Matric): SymTriMatrixC = { val trid = b.asInstanceOf [SymTriMatrixC] if (d1 == trid.d1) { new SymTriMatrixC (_dg + trid.dg, _sd + trid.sd) } else { flaw ("+", "matrix b has the wrong dimensions") null } // if } // + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Add this matrix and scalar x. * @param x the scalar to add */ def + (x: Complex): SymTriMatrixC = { new SymTriMatrixC (_dg + x, _sd + x) } // + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Add in-place this matrix and matrix b. * @param b the matrix to add (requires leDimensions) */ def += (b: Matric): SymTriMatrixC = { val trid = b.asInstanceOf [SymTriMatrixC] if (d1 == trid.d1) { _dg += trid.dg _sd += trid.sd } else { flaw ("+=", "matrix b has the wrong dimensions") } // if this } // += //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Add in-place this matrix and scalar x. * @param x the scalar to add */ def += (x: Complex): SymTriMatrixC = { _dg += x; _sd += x; this } // += //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** From this matrix subtract matrix b. * @param b the matrix to subtract (requires leDimensions) */ def - (b: Matric): SymTriMatrixC = { val trid = b.asInstanceOf [SymTriMatrixC] if (d1 == trid.d1) { new SymTriMatrixC (_dg - trid.dg, _sd - trid.sd) } else { flaw ("-", "matrix b has the wrong dimensions") null } // if this } // - //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** From this matrix subtract scalar x. * @param x the scalar to subtract */ def - (x: Complex): SymTriMatrixC = { new SymTriMatrixC (_dg - x, _sd - x) } // - //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** From this matrix subtract in-place matrix b. * @param b the matrix to subtract (requires leDimensions) */ def -= (b: Matric): SymTriMatrixC = { val trid = b.asInstanceOf [SymTriMatrixC] if (d1 == trid.d1) { _dg -= trid.dg _sd -= trid.sd } else { flaw ("-=", "matrix b has the wrong dimensions") } // if this } // -= //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** From this matrix subtract in-place scalar x. * @param x the scalar to subtract */ def -= (x: Complex): SymTriMatrixC = { _dg -= x; _sd -= x; this } // -= //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by matrix b. * @param b the matrix to multiply by */ def * (b: Matric): SymTriMatrixC = { throw new NoSuchMethodException ("SymTriMatrixC does not support * with general matrices") } // * //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by matrix b. Requires b to have type SymTriMatrixD, * but returns a more general type of matrix. * @param b the matrix to multiply by */ def * (b: SymTriMatrixC): MatrixC = { val c = new MatrixC (d1) for (i <- 0 until d1; j <- (i-2 max 0) to (i+2 min d1_1)) { var sum = _0 val k1 = ((i min j) - 1) max 0 val k2 = ((i max j) + 1) min d1_1 for (k <- k1 to k2) sum += at(i, k) * b.at(k, j) c(i, j) = sum } // for c } // * //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by vector u. * @param u the vector to multiply by */ def * (u: VectorC): VectorC = { val c = new VectorC (d1) c(0) = _dg(0) * u(0) + _sd(0) * u(1) for (i <- 1 until d1_1) { c(i) = _sd(i-1) * u(i-1) + _dg(i) * u(i) + _sd(i) * u(i+1) } // for c(d1-1) = _sd(d1-2) * u(d1-2) + _dg(d1-1) * u(d1-1) c } // * //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by scalar x. * @param x the scalar to multiply by */ def * (x: Complex): SymTriMatrixC = { new SymTriMatrixC (_dg * x, _sd * x) } // * //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply in-place this matrix by matrix b * @param b the matrix to multiply by */ def *= (b: Matric): SymTriMatrixC = { throw new NoSuchMethodException ("inplace matrix multiplication not yet implemented") } // *= //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply in-place this matrix by scalar x. * @param x the scalar to multiply by */ def *= (x: Complex): SymTriMatrixC = { _dg *= x; _sd *= x; this } // *= //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply this matrix by vector u to produce another matrix (a_ij * u_j) * @param u the vector to multiply by */ def ** (u: VectorC): SymTriMatrixC = { throw new NoSuchMethodException ("matrix * vector -> matrix not implemented yet") } // ** //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Multiply in-place this matrix by vector u to produce another matrix (a_ij * u_j) * @param u the vector to multiply by */ def **= (u: VectorC): SymTriMatrixC = { throw new NoSuchMethodException ("inplace matrix * vector -> matrix not implemented yet") } // **= //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Divide this matrix by scalar x. * @param x the scalar to divide by */ def / (x: Complex): SymTriMatrixC = { new SymTriMatrixC (_dg / x, _sd / x) } // / //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Divide in-place this matrix by scalar x. * @param x the scalar to divide by */ def /= (x: Complex): SymTriMatrixC = { _dg /= x; _sd /= x; this } // /= //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Raise this matrix to the pth power (for some integer p >= 2). * @param p the power to raise this matrix to */ def ~^ (p: Int): MatrixC = { throw new NoSuchMethodException ("matrix power function (~^) not implemented yet") } // ~^ //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Find the maximum element in this matrix. * @param e the ending row index (exclusive) for the search */ def max (e: Int = dim1): Complex = _dg(0 until e).max() max _sd(0 until e).max() //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Find the minimum element in this matrix. * @param e the ending row index (exclusive) for the search */ def min (e: Int = dim1): Complex = _dg(0 until e).min() min _sd(0 until e).min() //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = { if (k == 0) _dg else if (math.abs (k) == 1) _sd else { flaw ("getDiag", "nothing stored for diagonal " + k); null } } // getDiag //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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) { if (k == 0) _dg = u else if (math.abs (k) == 1) _sd = u else flaw ("setDiag", "nothing stored for diagonal " + k) } // setDiag //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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) { _dg.set (x) } //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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): SymTriMatrixC = { val s = if (relative) mag else Complex._1 // use matrix magnitude or 1 for (i <- range_d) if (abs (_dg(i)) <= s * thres) _dg(i) = _0 for (i <- range_s) if (abs (_sd(i)) <= s * thres) _sd(i) = _0 this } // clean //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the determinant of this matrix. */ def det: Complex = detHelper (d1 - 1) //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Solve for x in the equation a*x = d where a is this matrix * @param d the constant vector. */ def solve (d: VectorC): VectorC = { val x = new VectorC (d1) val c = new VectorC (d1) val a = new VectorC (d1) val b = new VectorC (d1) for (i <- range_s) { c(i) = _sd(i); a(i) = _sd(i) } for (i <- range_d) b(i) = _dg(i) c(0) = c(0) / b(0) d(0) = d(0) / b(0) for (i <- 1 until d1) { val id = Complex._1 / (b(i) - c(i-1) * a(i)) c(i) = c(i) * id d(i) = (d(i) - d(i-1) * a(i)) * id } // for x(d1-1) = d(d1-1) for (i <- d1 - 2 to 0 by -1) x(i) = d(i) - c(i) * x(i+1) x } // solve //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Helper method for computing the determinant of this matrix. * @param n the current dimension */ private def detHelper (n: Int): Complex = { if (n == 0) _dg(0) else if (n == 1) _dg(0) * _dg(1) - _sd(0) * _sd(0) else _dg(n) * detHelper (n - 1) - _sd(n-1) * _sd(n-1) * detHelper (n - 2) } // detHelper //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = { if (dim2 != dim1 + 1) flaw ("nullspace", "requires n (columns) = m (rows) + 1") reduce.col(dim2 - 1) * Complex._1n ++ Complex._1 } // nullspace //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = { if (dim2 != dim1 + 1) flaw ("nullspace", "requires n (columns) = m (rows) + 1") reduce_ip col(dim2 - 1) * Complex._1n ++ Complex._1 } // nullspace_ip //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = _dg.sum //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the sum of this matrix, i.e., the sum of its elements. */ def sum: Complex = _dg.sum + _sd.sum + _sd.sum //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the sum of the lower triangular region of this matrix. */ def sumLower: Complex = _sd.sum //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = _dg.norm1 + _sd.norm1 + _sd.norm1 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 = { val c = new VectorC (dim2) for (j <- range2) c(j) = col(j).norm1 c.max () } // norm1 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix is rectangular (all rows have the same number * of columns). */ def isRectangular: Boolean = true //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Check whether this matrix is nonnegative (has no negative elements). */ def isNonnegative: Boolean = _dg.isNonnegative && _sd.isNonnegative //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Convert this symmetric tridiagonal matrix to a string showing the diagonal * vector followed by the sub-diagonal vector. */ override def toString: String = "\nSymTriMatrixC(\t" + _dg + ", \n\t\t" + _sd + ")" //-------------------------------------------------------------------------- // The following methods are not useful for Symmetric Tridiagonal matrices: //-------------------------------------------------------------------------- def lud: Tuple2 [Matric, Matric] = { throw new NoSuchMethodException ("not implemented") } // lud def lud_ip: Tuple2 [Matric, Matric] = { throw new NoSuchMethodException ("not implemented") } // lud_ip def solve (l: Matric, u: Matric, b: VectorC): VectorC = { throw new NoSuchMethodException ("not implemented") } // solve def solve (lu: Tuple2 [Matric, Matric], b: VectorC): VectorC = { throw new NoSuchMethodException ("not implemented") } // solve def diag (b: Matric): SymTriMatrixC = { throw new NoSuchMethodException ("not implemented") } // diag def diag (p: Int, q: Int): SymTriMatrixC = { throw new NoSuchMethodException("not implemented") } // diag def inverse: SymTriMatrixC = { throw new NoSuchMethodException("not implemented") } // inverse def inverse_ip: SymTriMatrixC = { throw new NoSuchMethodException("not implemented") } // inverse_ip def reduce: SymTriMatrixC = { throw new NoSuchMethodException("not implemented") } // reduce def reduce_ip { throw new NoSuchMethodException("not implemented") } // reduce_ip } // SymTriMatrixC class //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SymTriMatrixCTest` object is used to test the `SymTriMatrixC` class. */ object SymTriMatrixCTest extends App { val a = new SymTriMatrixC (VectorC (1.0, 2.0, 3.0), VectorC (4.0, 5.0)) val b = new SymTriMatrixC (VectorC (2.0, 3.0, 4.0), VectorC (5.0, 6.0)) val u = VectorC (5.0, 3.0, 6.0) val c = new MatrixC ((3, 3), 1.0, 4.0, 0.0, 4.0, 2.0, 5.0, 0.0, 5.0, 3.0) val d = new MatrixC ((3, 3), 2.0, 5.0, 0.0, 5.0, 3.0, 6.0, 0.0, 6.0, 4.0) println ("a = " + a) println ("b = " + b) println ("a.det = " + a.det) println ("a + b = " + (a + b)) println ("a - b = " + (a - b)) println ("a * b = " + (a * b)) println ("a * u = " + (a * u)) println ("c * d = " + (c * d)) } // SymTriMatrixCTest object