blob: 1622275d6bf6002e0ac3bb50a5a34c6cc57b6dbe [file] [log] [blame]
// Copyright ©2017 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package gonum
import "math"
// Dlags2 computes 2-by-2 orthogonal matrices U, V and Q with the
// triangles of A and B specified by upper.
//
// If upper is true
//
// Uᵀ*A*Q = Uᵀ*[ a1 a2 ]*Q = [ x 0 ]
// [ 0 a3 ] [ x x ]
// and
// Vᵀ*B*Q = Vᵀ*[ b1 b2 ]*Q = [ x 0 ]
// [ 0 b3 ] [ x x ]
//
// otherwise
//
// Uᵀ*A*Q = Uᵀ*[ a1 0 ]*Q = [ x x ]
// [ a2 a3 ] [ 0 x ]
// and
// Vᵀ*B*Q = Vᵀ*[ b1 0 ]*Q = [ x x ]
// [ b2 b3 ] [ 0 x ].
//
// The rows of the transformed A and B are parallel, where
//
// U = [ csu snu ], V = [ csv snv ], Q = [ csq snq ]
// [ -snu csu ] [ -snv csv ] [ -snq csq ]
//
// Dlags2 is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlags2(upper bool, a1, a2, a3, b1, b2, b3 float64) (csu, snu, csv, snv, csq, snq float64) {
if upper {
// Input matrices A and B are upper triangular matrices.
//
// Form matrix C = A*adj(B) = [ a b ]
// [ 0 d ]
a := a1 * b3
d := a3 * b1
b := a2*b1 - a1*b2
// The SVD of real 2-by-2 triangular C.
//
// [ csl -snl ]*[ a b ]*[ csr snr ] = [ r 0 ]
// [ snl csl ] [ 0 d ] [ -snr csr ] [ 0 t ]
_, _, snr, csr, snl, csl := impl.Dlasv2(a, b, d)
if math.Abs(csl) >= math.Abs(snl) || math.Abs(csr) >= math.Abs(snr) {
// Compute the [0, 0] and [0, 1] elements of Uᵀ*A and Vᵀ*B,
// and [0, 1] element of |U|ᵀ*|A| and |V|ᵀ*|B|.
ua11r := csl * a1
ua12 := csl*a2 + snl*a3
vb11r := csr * b1
vb12 := csr*b2 + snr*b3
aua12 := math.Abs(csl)*math.Abs(a2) + math.Abs(snl)*math.Abs(a3)
avb12 := math.Abs(csr)*math.Abs(b2) + math.Abs(snr)*math.Abs(b3)
// Zero [0, 1] elements of Uᵀ*A and Vᵀ*B.
if math.Abs(ua11r)+math.Abs(ua12) != 0 {
if aua12/(math.Abs(ua11r)+math.Abs(ua12)) <= avb12/(math.Abs(vb11r)+math.Abs(vb12)) {
csq, snq, _ = impl.Dlartg(-ua11r, ua12)
} else {
csq, snq, _ = impl.Dlartg(-vb11r, vb12)
}
} else {
csq, snq, _ = impl.Dlartg(-vb11r, vb12)
}
csu = csl
snu = -snl
csv = csr
snv = -snr
} else {
// Compute the [1, 0] and [1, 1] elements of Uᵀ*A and Vᵀ*B,
// and [1, 1] element of |U|ᵀ*|A| and |V|ᵀ*|B|.
ua21 := -snl * a1
ua22 := -snl*a2 + csl*a3
vb21 := -snr * b1
vb22 := -snr*b2 + csr*b3
aua22 := math.Abs(snl)*math.Abs(a2) + math.Abs(csl)*math.Abs(a3)
avb22 := math.Abs(snr)*math.Abs(b2) + math.Abs(csr)*math.Abs(b3)
// Zero [1, 1] elements of Uᵀ*A and Vᵀ*B, and then swap.
if math.Abs(ua21)+math.Abs(ua22) != 0 {
if aua22/(math.Abs(ua21)+math.Abs(ua22)) <= avb22/(math.Abs(vb21)+math.Abs(vb22)) {
csq, snq, _ = impl.Dlartg(-ua21, ua22)
} else {
csq, snq, _ = impl.Dlartg(-vb21, vb22)
}
} else {
csq, snq, _ = impl.Dlartg(-vb21, vb22)
}
csu = snl
snu = csl
csv = snr
snv = csr
}
} else {
// Input matrices A and B are lower triangular matrices
//
// Form matrix C = A*adj(B) = [ a 0 ]
// [ c d ]
a := a1 * b3
d := a3 * b1
c := a2*b3 - a3*b2
// The SVD of real 2-by-2 triangular C
//
// [ csl -snl ]*[ a 0 ]*[ csr snr ] = [ r 0 ]
// [ snl csl ] [ c d ] [ -snr csr ] [ 0 t ]
_, _, snr, csr, snl, csl := impl.Dlasv2(a, c, d)
if math.Abs(csr) >= math.Abs(snr) || math.Abs(csl) >= math.Abs(snl) {
// Compute the [1, 0] and [1, 1] elements of Uᵀ*A and Vᵀ*B,
// and [1, 0] element of |U|ᵀ*|A| and |V|ᵀ*|B|.
ua21 := -snr*a1 + csr*a2
ua22r := csr * a3
vb21 := -snl*b1 + csl*b2
vb22r := csl * b3
aua21 := math.Abs(snr)*math.Abs(a1) + math.Abs(csr)*math.Abs(a2)
avb21 := math.Abs(snl)*math.Abs(b1) + math.Abs(csl)*math.Abs(b2)
// Zero [1, 0] elements of Uᵀ*A and Vᵀ*B.
if (math.Abs(ua21) + math.Abs(ua22r)) != 0 {
if aua21/(math.Abs(ua21)+math.Abs(ua22r)) <= avb21/(math.Abs(vb21)+math.Abs(vb22r)) {
csq, snq, _ = impl.Dlartg(ua22r, ua21)
} else {
csq, snq, _ = impl.Dlartg(vb22r, vb21)
}
} else {
csq, snq, _ = impl.Dlartg(vb22r, vb21)
}
csu = csr
snu = -snr
csv = csl
snv = -snl
} else {
// Compute the [0, 0] and [0, 1] elements of Uᵀ *A and Vᵀ *B,
// and [0, 0] element of |U|ᵀ*|A| and |V|ᵀ*|B|.
ua11 := csr*a1 + snr*a2
ua12 := snr * a3
vb11 := csl*b1 + snl*b2
vb12 := snl * b3
aua11 := math.Abs(csr)*math.Abs(a1) + math.Abs(snr)*math.Abs(a2)
avb11 := math.Abs(csl)*math.Abs(b1) + math.Abs(snl)*math.Abs(b2)
// Zero [0, 0] elements of Uᵀ*A and Vᵀ*B, and then swap.
if (math.Abs(ua11) + math.Abs(ua12)) != 0 {
if aua11/(math.Abs(ua11)+math.Abs(ua12)) <= avb11/(math.Abs(vb11)+math.Abs(vb12)) {
csq, snq, _ = impl.Dlartg(ua12, ua11)
} else {
csq, snq, _ = impl.Dlartg(vb12, vb11)
}
} else {
csq, snq, _ = impl.Dlartg(vb12, vb11)
}
csu = snr
snu = csr
csv = snl
snv = csl
}
}
return csu, snu, csv, snv, csq, snq
}