func (m *Matrix) L() *Matrix {
	rows, cols := m.Size()
	m2 := New(rows, cols)
	for row := 0; row < rows; row++ {
		for col := row; col < cols; col++ {
			m2.Set(row, col, m.Get(row, col))
		}
	}
	return m2
}
L returns the matrix with zeros below the diagonal.
package matrix
import (
	"bytes"
	"errors"
	"fmt"
	"math"
)
const (
	DefaultEpsilon = 0.000001
)
var (
	ErrDimensionMismatch = errors.New("dimension mismatch")
	ErrSingularValue = errors.New("singular value")
)
func New(rows, cols int, values ...float64) *Matrix {
	if len(values) == 0 {
		return &Matrix{
			stride:   cols,
			epsilon:  DefaultEpsilon,
			elements: make([]float64, rows*cols),
		}
	}
	elems := make([]float64, rows*cols)
	copy(elems, values)
	return &Matrix{
		stride:   cols,
		epsilon:  DefaultEpsilon,
		elements: elems,
	}
}
func Identity(order int) *Matrix {
	m := New(order, order)
	for i := 0; i < order; i++ {
		m.Set(i, i, 1)
	}
	return m
}
func Zero(rows, cols int) *Matrix {
	return New(rows, cols)
}
func Ones(rows, cols int) *Matrix {
	ones := make([]float64, rows*cols)
	for i := 0; i < (rows * cols); i++ {
		ones[i] = 1
	}
	return &Matrix{
		stride:   cols,
		epsilon:  DefaultEpsilon,
		elements: ones,
	}
}
func Eye(n int) *Matrix {
	m := Zero(n, n)
	for i := 0; i < len(m.elements); i += n + 1 {
		m.elements[i] = 1
	}
	return m
}
func NewFromArrays(a [][]float64) *Matrix {
	rows := len(a)
	if rows == 0 {
		return nil
	}
	cols := len(a[0])
	m := New(rows, cols)
	for row := 0; row < rows; row++ {
		for col := 0; col < cols; col++ {
			m.Set(row, col, a[row][col])
		}
	}
	return m
}
type Matrix struct {
	epsilon  float64
	elements []float64
	stride   int
}
func (m *Matrix) String() string {
	buffer := bytes.NewBuffer(nil)
	rows, cols := m.Size()
	for row := 0; row < rows; row++ {
		for col := 0; col < cols; col++ {
			buffer.WriteString(f64s(m.Get(row, col)))
			buffer.WriteRune(' ')
		}
		buffer.WriteRune('\n')
	}
	return buffer.String()
}
func (m *Matrix) Epsilon() float64 {
	return m.epsilon
}
func (m *Matrix) WithEpsilon(epsilon float64) *Matrix {
	m.epsilon = epsilon
	return m
}
func (m *Matrix) Each(action func(row, col int, value float64)) {
	rows, cols := m.Size()
	for row := 0; row < rows; row++ {
		for col := 0; col < cols; col++ {
			action(row, col, m.Get(row, col))
		}
	}
}
func (m *Matrix) Round() *Matrix {
	rows, cols := m.Size()
	for row := 0; row < rows; row++ {
		for col := 0; col < cols; col++ {
			m.Set(row, col, roundToEpsilon(m.Get(row, col), m.epsilon))
		}
	}
	return m
}
func (m *Matrix) Arrays() [][]float64 {
	rows, cols := m.Size()
	a := make([][]float64, rows)
	for row := 0; row < rows; row++ {
		a[row] = make([]float64, cols)
		for col := 0; col < cols; col++ {
			a[row][col] = m.Get(row, col)
		}
	}
	return a
}
func (m *Matrix) Size() (rows, cols int) {
	rows = len(m.elements) / m.stride
	cols = m.stride
	return
}
func (m *Matrix) IsSquare() bool {
	return m.stride == (len(m.elements) / m.stride)
}
func (m *Matrix) IsSymmetric() bool {
	rows, cols := m.Size()
	if rows != cols {
		return false
	}
	for i := 0; i < rows; i++ {
		for j := 0; j < i; j++ {
			if m.Get(i, j) != m.Get(j, i) {
				return false
			}
		}
	}
	return true
}
func (m *Matrix) Get(row, col int) float64 {
	index := (m.stride * row) + col
	return m.elements[index]
}
func (m *Matrix) Set(row, col int, val float64) {
	index := (m.stride * row) + col
	m.elements[index] = val
}
func (m *Matrix) Col(col int) Vector {
	rows, _ := m.Size()
	values := make([]float64, rows)
	for row := 0; row < rows; row++ {
		values[row] = m.Get(row, col)
	}
	return Vector(values)
}
func (m *Matrix) Row(row int) Vector {
	_, cols := m.Size()
	values := make([]float64, cols)
	for col := 0; col < cols; col++ {
		values[col] = m.Get(row, col)
	}
	return Vector(values)
}
func (m *Matrix) SubMatrix(i, j, rows, cols int) *Matrix {
	return &Matrix{
		elements: m.elements[i*m.stride+j : i*m.stride+j+(rows-1)*m.stride+cols],
		stride:   m.stride,
		epsilon:  m.epsilon,
	}
}
func (m *Matrix) ScaleRow(row int, scale float64) {
	startIndex := row * m.stride
	for i := startIndex; i < m.stride; i++ {
		m.elements[i] = m.elements[i] * scale
	}
}
func (m *Matrix) scaleAddRow(rd int, rs int, f float64) {
	indexd := rd * m.stride
	indexs := rs * m.stride
	for col := 0; col < m.stride; col++ {
		m.elements[indexd] += f * m.elements[indexs]
		indexd++
		indexs++
	}
}
func (m *Matrix) SwapRows(i, j int) {
	var vi, vj float64
	for col := 0; col < m.stride; col++ {
		vi = m.Get(i, col)
		vj = m.Get(j, col)
		m.Set(i, col, vj)
		m.Set(j, col, vi)
	}
}
func (m *Matrix) Augment(m2 *Matrix) (*Matrix, error) {
	mr, mc := m.Size()
	m2r, m2c := m2.Size()
	if mr != m2r {
		return nil, ErrDimensionMismatch
	}
	m3 := Zero(mr, mc+m2c)
	for row := 0; row < mr; row++ {
		for col := 0; col < mc; col++ {
			m3.Set(row, col, m.Get(row, col))
		}
		for col := 0; col < m2c; col++ {
			m3.Set(row, mc+col, m2.Get(row, col))
		}
	}
	return m3, nil
}
func (m *Matrix) Copy() *Matrix {
	m2 := &Matrix{stride: m.stride, epsilon: m.epsilon, elements: make([]float64, len(m.elements))}
	copy(m2.elements, m.elements)
	return m2
}
func (m *Matrix) DiagonalVector() Vector {
	rows, cols := m.Size()
	rank := minInt(rows, cols)
	values := make([]float64, rank)
	for index := 0; index < rank; index++ {
		values[index] = m.Get(index, index)
	}
	return Vector(values)
}
func (m *Matrix) Diagonal() *Matrix {
	rows, cols := m.Size()
	rank := minInt(rows, cols)
	m2 := New(rank, rank)
	for index := 0; index < rank; index++ {
		m2.Set(index, index, m.Get(index, index))
	}
	return m2
}
func (m *Matrix) Equals(other *Matrix) bool {
	if other == nil && m != nil {
		return false
	} else if other == nil {
		return true
	}
	if m.stride != other.stride {
		return false
	}
	msize := len(m.elements)
	m2size := len(other.elements)
	if msize != m2size {
		return false
	}
	for i := 0; i < msize; i++ {
		if m.elements[i] != other.elements[i] {
			return false
		}
	}
	return true
}