1// Copyright ©2017 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
5package gonum
6
7import "math"
8
9// Dlags2 computes 2-by-2 orthogonal matrices U, V and Q with the
10// triangles of A and B specified by upper.
11//
12// If upper is true
13//
14//  U^T*A*Q = U^T*[ a1 a2 ]*Q = [ x  0 ]
15//                [ 0  a3 ]     [ x  x ]
16// and
17//  V^T*B*Q = V^T*[ b1 b2 ]*Q = [ x  0 ]
18//                [ 0  b3 ]     [ x  x ]
19//
20// otherwise
21//
22//  U^T*A*Q = U^T*[ a1 0  ]*Q = [ x  x ]
23//                [ a2 a3 ]     [ 0  x ]
24// and
25//  V^T*B*Q = V^T*[ b1 0  ]*Q = [ x  x ]
26//                [ b2 b3 ]     [ 0  x ].
27//
28// The rows of the transformed A and B are parallel, where
29//
30//  U = [  csu  snu ], V = [  csv snv ], Q = [  csq   snq ]
31//      [ -snu  csu ]      [ -snv csv ]      [ -snq   csq ]
32//
33// Dlags2 is an internal routine. It is exported for testing purposes.
34func (impl Implementation) Dlags2(upper bool, a1, a2, a3, b1, b2, b3 float64) (csu, snu, csv, snv, csq, snq float64) {
35	if upper {
36		// Input matrices A and B are upper triangular matrices.
37		//
38		// Form matrix C = A*adj(B) = [ a b ]
39		//                            [ 0 d ]
40		a := a1 * b3
41		d := a3 * b1
42		b := a2*b1 - a1*b2
43
44		// The SVD of real 2-by-2 triangular C.
45		//
46		//  [ csl -snl ]*[ a b ]*[  csr  snr ] = [ r 0 ]
47		//  [ snl  csl ] [ 0 d ] [ -snr  csr ]   [ 0 t ]
48		_, _, snr, csr, snl, csl := impl.Dlasv2(a, b, d)
49
50		if math.Abs(csl) >= math.Abs(snl) || math.Abs(csr) >= math.Abs(snr) {
51			// Compute the [0, 0] and [0, 1] elements of U^T*A and V^T*B,
52			// and [0, 1] element of |U|^T*|A| and |V|^T*|B|.
53
54			ua11r := csl * a1
55			ua12 := csl*a2 + snl*a3
56
57			vb11r := csr * b1
58			vb12 := csr*b2 + snr*b3
59
60			aua12 := math.Abs(csl)*math.Abs(a2) + math.Abs(snl)*math.Abs(a3)
61			avb12 := math.Abs(csr)*math.Abs(b2) + math.Abs(snr)*math.Abs(b3)
62
63			// Zero [0, 1] elements of U^T*A and V^T*B.
64			if math.Abs(ua11r)+math.Abs(ua12) != 0 {
65				if aua12/(math.Abs(ua11r)+math.Abs(ua12)) <= avb12/(math.Abs(vb11r)+math.Abs(vb12)) {
66					csq, snq, _ = impl.Dlartg(-ua11r, ua12)
67				} else {
68					csq, snq, _ = impl.Dlartg(-vb11r, vb12)
69				}
70			} else {
71				csq, snq, _ = impl.Dlartg(-vb11r, vb12)
72			}
73
74			csu = csl
75			snu = -snl
76			csv = csr
77			snv = -snr
78		} else {
79			// Compute the [1, 0] and [1, 1] elements of U^T*A and V^T*B,
80			// and [1, 1] element of |U|^T*|A| and |V|^T*|B|.
81
82			ua21 := -snl * a1
83			ua22 := -snl*a2 + csl*a3
84
85			vb21 := -snr * b1
86			vb22 := -snr*b2 + csr*b3
87
88			aua22 := math.Abs(snl)*math.Abs(a2) + math.Abs(csl)*math.Abs(a3)
89			avb22 := math.Abs(snr)*math.Abs(b2) + math.Abs(csr)*math.Abs(b3)
90
91			// Zero [1, 1] elements of U^T*A and V^T*B, and then swap.
92			if math.Abs(ua21)+math.Abs(ua22) != 0 {
93				if aua22/(math.Abs(ua21)+math.Abs(ua22)) <= avb22/(math.Abs(vb21)+math.Abs(vb22)) {
94					csq, snq, _ = impl.Dlartg(-ua21, ua22)
95				} else {
96					csq, snq, _ = impl.Dlartg(-vb21, vb22)
97				}
98			} else {
99				csq, snq, _ = impl.Dlartg(-vb21, vb22)
100			}
101
102			csu = snl
103			snu = csl
104			csv = snr
105			snv = csr
106		}
107	} else {
108		// Input matrices A and B are lower triangular matrices
109		//
110		// Form matrix C = A*adj(B) = [ a 0 ]
111		//                            [ c d ]
112		a := a1 * b3
113		d := a3 * b1
114		c := a2*b3 - a3*b2
115
116		// The SVD of real 2-by-2 triangular C
117		//
118		// [ csl -snl ]*[ a 0 ]*[  csr  snr ] = [ r 0 ]
119		// [ snl  csl ] [ c d ] [ -snr  csr ]   [ 0 t ]
120		_, _, snr, csr, snl, csl := impl.Dlasv2(a, c, d)
121
122		if math.Abs(csr) >= math.Abs(snr) || math.Abs(csl) >= math.Abs(snl) {
123			// Compute the [1, 0] and [1, 1] elements of U^T*A and V^T*B,
124			// and [1, 0] element of |U|^T*|A| and |V|^T*|B|.
125
126			ua21 := -snr*a1 + csr*a2
127			ua22r := csr * a3
128
129			vb21 := -snl*b1 + csl*b2
130			vb22r := csl * b3
131
132			aua21 := math.Abs(snr)*math.Abs(a1) + math.Abs(csr)*math.Abs(a2)
133			avb21 := math.Abs(snl)*math.Abs(b1) + math.Abs(csl)*math.Abs(b2)
134
135			// Zero [1, 0] elements of U^T*A and V^T*B.
136			if (math.Abs(ua21) + math.Abs(ua22r)) != 0 {
137				if aua21/(math.Abs(ua21)+math.Abs(ua22r)) <= avb21/(math.Abs(vb21)+math.Abs(vb22r)) {
138					csq, snq, _ = impl.Dlartg(ua22r, ua21)
139				} else {
140					csq, snq, _ = impl.Dlartg(vb22r, vb21)
141				}
142			} else {
143				csq, snq, _ = impl.Dlartg(vb22r, vb21)
144			}
145
146			csu = csr
147			snu = -snr
148			csv = csl
149			snv = -snl
150		} else {
151			// Compute the [0, 0] and [0, 1] elements of U^T *A and V^T *B,
152			// and [0, 0] element of |U|^T*|A| and |V|^T*|B|.
153
154			ua11 := csr*a1 + snr*a2
155			ua12 := snr * a3
156
157			vb11 := csl*b1 + snl*b2
158			vb12 := snl * b3
159
160			aua11 := math.Abs(csr)*math.Abs(a1) + math.Abs(snr)*math.Abs(a2)
161			avb11 := math.Abs(csl)*math.Abs(b1) + math.Abs(snl)*math.Abs(b2)
162
163			// Zero [0, 0] elements of U^T*A and V^T*B, and then swap.
164			if (math.Abs(ua11) + math.Abs(ua12)) != 0 {
165				if aua11/(math.Abs(ua11)+math.Abs(ua12)) <= avb11/(math.Abs(vb11)+math.Abs(vb12)) {
166					csq, snq, _ = impl.Dlartg(ua12, ua11)
167				} else {
168					csq, snq, _ = impl.Dlartg(vb12, vb11)
169				}
170			} else {
171				csq, snq, _ = impl.Dlartg(vb12, vb11)
172			}
173
174			csu = snr
175			snu = csr
176			csv = snl
177			snv = csl
178		}
179	}
180
181	return csu, snu, csv, snv, csq, snq
182}
183