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