1// Copyright ©2015 The Gonum Authors. All rights reserved. 2// Use of this source code is governed by a BSD-style 3// license that can be found in the LICENSE file. 4 5// Package mat provides implementations of float64 and complex128 matrix 6// structures and linear algebra operations on them. 7// 8// Overview 9// 10// This section provides a quick overview of the mat package. The following 11// sections provide more in depth commentary. 12// 13// mat provides: 14// - Interfaces for Matrix classes (Matrix, Symmetric, Triangular) 15// - Concrete implementations (Dense, SymDense, TriDense) 16// - Methods and functions for using matrix data (Add, Trace, SymRankOne) 17// - Types for constructing and using matrix factorizations (QR, LU) 18// - The complementary types for complex matrices, CMatrix, CSymDense, etc. 19// 20// A matrix may be constructed through the corresponding New function. If no 21// backing array is provided the matrix will be initialized to all zeros. 22// // Allocate a zeroed real matrix of size 3×5 23// zero := mat.NewDense(3, 5, nil) 24// If a backing data slice is provided, the matrix will have those elements. 25// Matrices are all stored in row-major format. 26// // Generate a 6×6 matrix of random values. 27// data := make([]float64, 36) 28// for i := range data { 29// data[i] = rand.NormFloat64() 30// } 31// a := mat.NewDense(6, 6, data) 32// Operations involving matrix data are implemented as functions when the values 33// of the matrix remain unchanged 34// tr := mat.Trace(a) 35// and are implemented as methods when the operation modifies the receiver. 36// zero.Copy(a) 37// 38// Receivers must be the correct size for the matrix operations, otherwise the 39// operation will panic. As a special case for convenience, a zero-value matrix 40// will be modified to have the correct size, allocating data if necessary. 41// var c mat.Dense // construct a new zero-sized matrix 42// c.Mul(a, a) // c is automatically adjusted to be 6×6 43// 44// Zero-value of a matrix 45// 46// A zero-value matrix is either the Go language definition of a zero-value or 47// is a zero-sized matrix with zero-length stride. Matrix implementations may have 48// a Reset method to revert the receiver into a zero-valued matrix and an IsZero 49// method that returns whether the matrix is zero-valued. 50// So the following will all result in a zero-value matrix. 51// - var a mat.Dense 52// - a := NewDense(0, 0, make([]float64, 0, 100)) 53// - a.Reset() 54// A zero-value matrix can not be sliced even if it does have an adequately sized 55// backing data slice, but can be expanded using its Grow method if it exists. 56// 57// The Matrix Interfaces 58// 59// The Matrix interface is the common link between the concrete types of real 60// matrices, The Matrix interface is defined by three functions: Dims, which 61// returns the dimensions of the Matrix, At, which returns the element in the 62// specified location, and T for returning a Transpose (discussed later). All of 63// the concrete types can perform these behaviors and so implement the interface. 64// Methods and functions are designed to use this interface, so in particular the method 65// func (m *Dense) Mul(a, b Matrix) 66// constructs a *Dense from the result of a multiplication with any Matrix types, 67// not just *Dense. Where more restrictive requirements must be met, there are also the 68// Symmetric and Triangular interfaces. For example, in 69// func (s *SymDense) AddSym(a, b Symmetric) 70// the Symmetric interface guarantees a symmetric result. 71// 72// The CMatrix interface plays the same role for complex matrices. The difference 73// is that the CMatrix type has the H method instead T, for returning the conjugate 74// transpose. 75// 76// (Conjugate) Transposes 77// 78// The T method is used for transposition on real matrices, and H is used for 79// conjugate transposition on complex matrices. For example, c.Mul(a.T(), b) computes 80// c = a^T * b. The mat types implement this method implicitly — 81// see the Transpose and Conjugate types for more details. Note that some 82// operations have a transpose as part of their definition, as in *SymDense.SymOuterK. 83// 84// Matrix Factorization 85// 86// Matrix factorizations, such as the LU decomposition, typically have their own 87// specific data storage, and so are each implemented as a specific type. The 88// factorization can be computed through a call to Factorize 89// var lu mat.LU 90// lu.Factorize(a) 91// The elements of the factorization can be extracted through methods on the 92// factorized type, i.e. *LU.UTo. The factorization types can also be used directly, 93// as in *Dense.SolveCholesky. Some factorizations can be updated directly, 94// without needing to update the original matrix and refactorize, 95// as in *LU.RankOne. 96// 97// BLAS and LAPACK 98// 99// BLAS and LAPACK are the standard APIs for linear algebra routines. Many 100// operations in mat are implemented using calls to the wrapper functions 101// in gonum/blas/blas64 and gonum/lapack/lapack64 and their complex equivalents. 102// By default, blas64 and lapack64 call the native Go implementations of the 103// routines. Alternatively, it is possible to use C-based implementations of the 104// APIs through the respective cgo packages and "Use" functions. The Go 105// implementation of LAPACK (used by default) makes calls 106// through blas64, so if a cgo BLAS implementation is registered, the lapack64 107// calls will be partially executed in Go and partially executed in C. 108// 109// Type Switching 110// 111// The Matrix abstraction enables efficiency as well as interoperability. Go's 112// type reflection capabilities are used to choose the most efficient routine 113// given the specific concrete types. For example, in 114// c.Mul(a, b) 115// if a and b both implement RawMatrixer, that is, they can be represented as a 116// blas64.General, blas64.Gemm (general matrix multiplication) is called, while 117// instead if b is a RawSymmetricer blas64.Symm is used (general-symmetric 118// multiplication), and if b is a *VecDense blas64.Gemv is used. 119// 120// There are many possible type combinations and special cases. No specific guarantees 121// are made about the performance of any method, and in particular, note that an 122// abstract matrix type may be copied into a concrete type of the corresponding 123// value. If there are specific special cases that are needed, please submit a 124// pull-request or file an issue. 125// 126// Invariants 127// 128// Matrix input arguments to functions are never directly modified. If an operation 129// changes Matrix data, the mutated matrix will be the receiver of a function. 130// 131// For convenience, a matrix may be used as both a receiver and as an input, e.g. 132// a.Pow(a, 6) 133// v.SolveVec(a.T(), v) 134// though in many cases this will cause an allocation (see Element Aliasing). 135// An exception to this rule is Copy, which does not allow a.Copy(a.T()). 136// 137// Element Aliasing 138// 139// Most methods in mat modify receiver data. It is forbidden for the modified 140// data region of the receiver to overlap the used data area of the input 141// arguments. The exception to this rule is when the method receiver is equal to one 142// of the input arguments, as in the a.Pow(a, 6) call above, or its implicit transpose. 143// 144// This prohibition is to help avoid subtle mistakes when the method needs to read 145// from and write to the same data region. There are ways to make mistakes using the 146// mat API, and mat functions will detect and complain about those. 147// There are many ways to make mistakes by excursion from the mat API via 148// interaction with raw matrix values. 149// 150// If you need to read the rest of this section to understand the behavior of 151// your program, you are being clever. Don't be clever. If you must be clever, 152// blas64 and lapack64 may be used to call the behavior directly. 153// 154// mat will use the following rules to detect overlap between the receiver and one 155// of the inputs: 156// - the input implements one of the Raw methods, and 157// - the address ranges of the backing data slices overlap, and 158// - the strides differ or there is an overlap in the used data elements. 159// If such an overlap is detected, the method will panic. 160// 161// The following cases will not panic: 162// - the data slices do not overlap, 163// - there is pointer identity between the receiver and input values after 164// the value has been untransposed if necessary. 165// 166// mat will not attempt to detect element overlap if the input does not implement a 167// Raw method. Method behavior is undefined if there is undetected overlap. 168// 169package mat // import "gonum.org/v1/gonum/mat" 170