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.
8 "gonum.org/v1/gonum/blas"
9 "gonum.org/v1/gonum/blas/blas64"
10 "gonum.org/v1/gonum/lapack/lapack64"
13 // Solve finds a minimum-norm solution to a system of linear equations defined
14 // by the matrices a and b. If A is singular or near-singular, a Condition error
15 // is returned. See the documentation for Condition for more information.
17 // The minimization problem solved depends on the input parameters:
18 // - if m >= n, find X such that ||A*X - B||_2 is minimized,
19 // - if m < n, find the minimum norm solution of A * X = B.
20 // The solution matrix, X, is stored in-place into the receiver.
21 func (m *Dense) Solve(a, b Matrix) error {
29 // TODO(btracey): Add special cases for SymDense, etc.
30 aU, aTrans := untranspose(a)
31 bU, bTrans := untranspose(b)
32 switch rma := aU.(type) {
40 switch rm := bU.(type) {
42 if m != bU || bTrans {
43 if m == bU || m.checkOverlap(rm.RawMatrix()) {
44 tmp := getWorkspace(br, bc, false)
56 // m and b share data so Copy cannot be used directly.
57 tmp := getWorkspace(br, bc, false)
64 rm := rma.RawTriangular()
65 blas64.Trsm(side, tA, 1, rm, m.mat)
66 work := getFloats(3*rm.N, false)
67 iwork := getInts(rm.N, false)
68 cond := lapack64.Trcon(CondNorm, rm, work, iwork)
71 if cond > ConditionTolerance {
72 return Condition(cond)
85 for i := 0; i < ar; i++ {
86 v := m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+ac]
94 return lu.Solve(m, false, b)
98 return qr.Solve(m, false, b)
102 return lq.Solve(m, false, b)
106 // SolveVec finds a minimum-norm solution to a system of linear equations defined
107 // by the matrix a and the right-hand side column vector b. If A is singular or
108 // near-singular, a Condition error is returned. See the documentation for
109 // Dense.Solve for more information.
110 func (v *VecDense) SolveVec(a Matrix, b Vector) error {
111 if _, bc := b.Dims(); bc != 1 {
116 // The Solve implementation is non-trivial, so rather than duplicate the code,
117 // instead recast the VecDenses as Dense and call the matrix code.
119 if rv, ok := b.(RawVectorer); ok {
120 bmat := rv.RawVector()
126 // We conditionally create bm as m when b and v are identical
127 // to prevent the overlap detection code from identifying m
128 // and bm as overlapping but not identical.
131 b := VecDense{mat: bmat, n: b.Len()}
134 return m.Solve(a, bm)