github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dlaexc.go (about) 1 // Copyright ©2016 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 gonum 6 7 import ( 8 "math" 9 10 "github.com/jingcheng-WU/gonum/blas" 11 "github.com/jingcheng-WU/gonum/blas/blas64" 12 "github.com/jingcheng-WU/gonum/lapack" 13 ) 14 15 // Dlaexc swaps two adjacent diagonal blocks of order 1 or 2 in an n×n upper 16 // quasi-triangular matrix T by an orthogonal similarity transformation. 17 // 18 // T must be in Schur canonical form, that is, block upper triangular with 1×1 19 // and 2×2 diagonal blocks; each 2×2 diagonal block has its diagonal elements 20 // equal and its off-diagonal elements of opposite sign. On return, T will 21 // contain the updated matrix again in Schur canonical form. 22 // 23 // If wantq is true, the transformation is accumulated in the n×n matrix Q, 24 // otherwise Q is not referenced. 25 // 26 // j1 is the index of the first row of the first block. n1 and n2 are the order 27 // of the first and second block, respectively. 28 // 29 // work must have length at least n, otherwise Dlaexc will panic. 30 // 31 // If ok is false, the transformed matrix T would be too far from Schur form. 32 // The blocks are not swapped, and T and Q are not modified. 33 // 34 // If n1 and n2 are both equal to 1, Dlaexc will always return true. 35 // 36 // Dlaexc is an internal routine. It is exported for testing purposes. 37 func (impl Implementation) Dlaexc(wantq bool, n int, t []float64, ldt int, q []float64, ldq int, j1, n1, n2 int, work []float64) (ok bool) { 38 switch { 39 case n < 0: 40 panic(nLT0) 41 case ldt < max(1, n): 42 panic(badLdT) 43 case wantq && ldt < max(1, n): 44 panic(badLdQ) 45 case j1 < 0 || n <= j1: 46 panic(badJ1) 47 case len(work) < n: 48 panic(shortWork) 49 case n1 < 0 || 2 < n1: 50 panic(badN1) 51 case n2 < 0 || 2 < n2: 52 panic(badN2) 53 } 54 55 if n == 0 || n1 == 0 || n2 == 0 { 56 return true 57 } 58 59 switch { 60 case len(t) < (n-1)*ldt+n: 61 panic(shortT) 62 case wantq && len(q) < (n-1)*ldq+n: 63 panic(shortQ) 64 } 65 66 if j1+n1 >= n { 67 // TODO(vladimir-ch): Reference LAPACK does this check whether 68 // the start of the second block is in the matrix T. It returns 69 // true if it is not and moreover it does not check whether the 70 // whole second block fits into T. This does not feel 71 // satisfactory. The only caller of Dlaexc is Dtrexc, so if the 72 // caller makes sure that this does not happen, we could be 73 // stricter here. 74 return true 75 } 76 77 j2 := j1 + 1 78 j3 := j1 + 2 79 80 bi := blas64.Implementation() 81 82 if n1 == 1 && n2 == 1 { 83 // Swap two 1×1 blocks. 84 t11 := t[j1*ldt+j1] 85 t22 := t[j2*ldt+j2] 86 87 // Determine the transformation to perform the interchange. 88 cs, sn, _ := impl.Dlartg(t[j1*ldt+j2], t22-t11) 89 90 // Apply transformation to the matrix T. 91 if n-j3 > 0 { 92 bi.Drot(n-j3, t[j1*ldt+j3:], 1, t[j2*ldt+j3:], 1, cs, sn) 93 } 94 if j1 > 0 { 95 bi.Drot(j1, t[j1:], ldt, t[j2:], ldt, cs, sn) 96 } 97 98 t[j1*ldt+j1] = t22 99 t[j2*ldt+j2] = t11 100 101 if wantq { 102 // Accumulate transformation in the matrix Q. 103 bi.Drot(n, q[j1:], ldq, q[j2:], ldq, cs, sn) 104 } 105 106 return true 107 } 108 109 // Swapping involves at least one 2×2 block. 110 // 111 // Copy the diagonal block of order n1+n2 to the local array d and 112 // compute its norm. 113 nd := n1 + n2 114 var d [16]float64 115 const ldd = 4 116 impl.Dlacpy(blas.All, nd, nd, t[j1*ldt+j1:], ldt, d[:], ldd) 117 dnorm := impl.Dlange(lapack.MaxAbs, nd, nd, d[:], ldd, work) 118 119 // Compute machine-dependent threshold for test for accepting swap. 120 eps := dlamchP 121 thresh := math.Max(10*eps*dnorm, dlamchS/eps) 122 123 // Solve T11*X - X*T22 = scale*T12 for X. 124 var x [4]float64 125 const ldx = 2 126 scale, _, _ := impl.Dlasy2(false, false, -1, n1, n2, d[:], ldd, d[n1*ldd+n1:], ldd, d[n1:], ldd, x[:], ldx) 127 128 // Swap the adjacent diagonal blocks. 129 switch { 130 case n1 == 1 && n2 == 2: 131 // Generate elementary reflector H so that 132 // ( scale, X11, X12 ) H = ( 0, 0, * ) 133 u := [3]float64{scale, x[0], 1} 134 _, tau := impl.Dlarfg(3, x[1], u[:2], 1) 135 t11 := t[j1*ldt+j1] 136 137 // Perform swap provisionally on diagonal block in d. 138 impl.Dlarfx(blas.Left, 3, 3, u[:], tau, d[:], ldd, work) 139 impl.Dlarfx(blas.Right, 3, 3, u[:], tau, d[:], ldd, work) 140 141 // Test whether to reject swap. 142 if math.Max(math.Abs(d[2*ldd]), math.Max(math.Abs(d[2*ldd+1]), math.Abs(d[2*ldd+2]-t11))) > thresh { 143 return false 144 } 145 146 // Accept swap: apply transformation to the entire matrix T. 147 impl.Dlarfx(blas.Left, 3, n-j1, u[:], tau, t[j1*ldt+j1:], ldt, work) 148 impl.Dlarfx(blas.Right, j2+1, 3, u[:], tau, t[j1:], ldt, work) 149 150 t[j3*ldt+j1] = 0 151 t[j3*ldt+j2] = 0 152 t[j3*ldt+j3] = t11 153 154 if wantq { 155 // Accumulate transformation in the matrix Q. 156 impl.Dlarfx(blas.Right, n, 3, u[:], tau, q[j1:], ldq, work) 157 } 158 159 case n1 == 2 && n2 == 1: 160 // Generate elementary reflector H so that: 161 // H ( -X11 ) = ( * ) 162 // ( -X21 ) = ( 0 ) 163 // ( scale ) = ( 0 ) 164 u := [3]float64{1, -x[ldx], scale} 165 _, tau := impl.Dlarfg(3, -x[0], u[1:], 1) 166 t33 := t[j3*ldt+j3] 167 168 // Perform swap provisionally on diagonal block in D. 169 impl.Dlarfx(blas.Left, 3, 3, u[:], tau, d[:], ldd, work) 170 impl.Dlarfx(blas.Right, 3, 3, u[:], tau, d[:], ldd, work) 171 172 // Test whether to reject swap. 173 if math.Max(math.Abs(d[ldd]), math.Max(math.Abs(d[2*ldd]), math.Abs(d[0]-t33))) > thresh { 174 return false 175 } 176 177 // Accept swap: apply transformation to the entire matrix T. 178 impl.Dlarfx(blas.Right, j3+1, 3, u[:], tau, t[j1:], ldt, work) 179 impl.Dlarfx(blas.Left, 3, n-j1-1, u[:], tau, t[j1*ldt+j2:], ldt, work) 180 181 t[j1*ldt+j1] = t33 182 t[j2*ldt+j1] = 0 183 t[j3*ldt+j1] = 0 184 185 if wantq { 186 // Accumulate transformation in the matrix Q. 187 impl.Dlarfx(blas.Right, n, 3, u[:], tau, q[j1:], ldq, work) 188 } 189 190 default: // n1 == 2 && n2 == 2 191 // Generate elementary reflectors H_1 and H_2 so that: 192 // H_2 H_1 ( -X11 -X12 ) = ( * * ) 193 // ( -X21 -X22 ) ( 0 * ) 194 // ( scale 0 ) ( 0 0 ) 195 // ( 0 scale ) ( 0 0 ) 196 u1 := [3]float64{1, -x[ldx], scale} 197 _, tau1 := impl.Dlarfg(3, -x[0], u1[1:], 1) 198 199 temp := -tau1 * (x[1] + u1[1]*x[ldx+1]) 200 u2 := [3]float64{1, -temp * u1[2], scale} 201 _, tau2 := impl.Dlarfg(3, -temp*u1[1]-x[ldx+1], u2[1:], 1) 202 203 // Perform swap provisionally on diagonal block in D. 204 impl.Dlarfx(blas.Left, 3, 4, u1[:], tau1, d[:], ldd, work) 205 impl.Dlarfx(blas.Right, 4, 3, u1[:], tau1, d[:], ldd, work) 206 impl.Dlarfx(blas.Left, 3, 4, u2[:], tau2, d[ldd:], ldd, work) 207 impl.Dlarfx(blas.Right, 4, 3, u2[:], tau2, d[1:], ldd, work) 208 209 // Test whether to reject swap. 210 m1 := math.Max(math.Abs(d[2*ldd]), math.Abs(d[2*ldd+1])) 211 m2 := math.Max(math.Abs(d[3*ldd]), math.Abs(d[3*ldd+1])) 212 if math.Max(m1, m2) > thresh { 213 return false 214 } 215 216 // Accept swap: apply transformation to the entire matrix T. 217 j4 := j1 + 3 218 impl.Dlarfx(blas.Left, 3, n-j1, u1[:], tau1, t[j1*ldt+j1:], ldt, work) 219 impl.Dlarfx(blas.Right, j4+1, 3, u1[:], tau1, t[j1:], ldt, work) 220 impl.Dlarfx(blas.Left, 3, n-j1, u2[:], tau2, t[j2*ldt+j1:], ldt, work) 221 impl.Dlarfx(blas.Right, j4+1, 3, u2[:], tau2, t[j2:], ldt, work) 222 223 t[j3*ldt+j1] = 0 224 t[j3*ldt+j2] = 0 225 t[j4*ldt+j1] = 0 226 t[j4*ldt+j2] = 0 227 228 if wantq { 229 // Accumulate transformation in the matrix Q. 230 impl.Dlarfx(blas.Right, n, 3, u1[:], tau1, q[j1:], ldq, work) 231 impl.Dlarfx(blas.Right, n, 3, u2[:], tau2, q[j2:], ldq, work) 232 } 233 } 234 235 if n2 == 2 { 236 // Standardize new 2×2 block T11. 237 a, b := t[j1*ldt+j1], t[j1*ldt+j2] 238 c, d := t[j2*ldt+j1], t[j2*ldt+j2] 239 var cs, sn float64 240 t[j1*ldt+j1], t[j1*ldt+j2], t[j2*ldt+j1], t[j2*ldt+j2], _, _, _, _, cs, sn = impl.Dlanv2(a, b, c, d) 241 if n-j1-2 > 0 { 242 bi.Drot(n-j1-2, t[j1*ldt+j1+2:], 1, t[j2*ldt+j1+2:], 1, cs, sn) 243 } 244 if j1 > 0 { 245 bi.Drot(j1, t[j1:], ldt, t[j2:], ldt, cs, sn) 246 } 247 if wantq { 248 bi.Drot(n, q[j1:], ldq, q[j2:], ldq, cs, sn) 249 } 250 } 251 if n1 == 2 { 252 // Standardize new 2×2 block T22. 253 j3 := j1 + n2 254 j4 := j3 + 1 255 a, b := t[j3*ldt+j3], t[j3*ldt+j4] 256 c, d := t[j4*ldt+j3], t[j4*ldt+j4] 257 var cs, sn float64 258 t[j3*ldt+j3], t[j3*ldt+j4], t[j4*ldt+j3], t[j4*ldt+j4], _, _, _, _, cs, sn = impl.Dlanv2(a, b, c, d) 259 if n-j3-2 > 0 { 260 bi.Drot(n-j3-2, t[j3*ldt+j3+2:], 1, t[j4*ldt+j3+2:], 1, cs, sn) 261 } 262 bi.Drot(j3, t[j3:], ldt, t[j4:], ldt, cs, sn) 263 if wantq { 264 bi.Drot(n, q[j3:], ldq, q[j4:], ldq, cs, sn) 265 } 266 } 267 268 return true 269 }