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.
9 // Dlaln2 solves a linear equation or a system of 2 linear equations of the form
10 // (ca A - w D) X = scale B, if trans == false,
11 // (ca A^T - w D) X = scale B, if trans == true,
12 // where A is a na×na real matrix, ca is a real scalar, D is a na×na diagonal
13 // real matrix, w is a scalar, real if nw == 1, complex if nw == 2, and X and B
14 // are na×1 matrices, real if w is real, complex if w is complex.
16 // If w is complex, X and B are represented as na×2 matrices, the first column
17 // of each being the real part and the second being the imaginary part.
19 // na and nw must be 1 or 2, otherwise Dlaln2 will panic.
21 // d1 and d2 are the diagonal elements of D. d2 is not used if na == 1.
23 // wr and wi represent the real and imaginary part, respectively, of the scalar
24 // w. wi is not used if nw == 1.
26 // smin is the desired lower bound on the singular values of A. This should be
27 // a safe distance away from underflow or overflow, say, between
28 // (underflow/machine precision) and (overflow*machine precision).
30 // If both singular values of (ca A - w D) are less than smin, smin*identity
31 // will be used instead of (ca A - w D). If only one singular value is less than
32 // smin, one element of (ca A - w D) will be perturbed enough to make the
33 // smallest singular value roughly smin. If both singular values are at least
34 // smin, (ca A - w D) will not be perturbed. In any case, the perturbation will
35 // be at most some small multiple of max(smin, ulp*norm(ca A - w D)). The
36 // singular values are computed by infinity-norm approximations, and thus will
37 // only be correct to a factor of 2 or so.
39 // All input quantities are assumed to be smaller than overflow by a reasonable
42 // scale is a scaling factor less than or equal to 1 which is chosen so that X
43 // can be computed without overflow. X is further scaled if necessary to assure
44 // that norm(ca A - w D)*norm(X) is less than overflow.
46 // xnorm contains the infinity-norm of X when X is regarded as a na×nw real
49 // ok will be false if (ca A - w D) had to be perturbed to make its smallest
50 // singular value greater than smin, otherwise ok will be true.
52 // Dlaln2 is an internal routine. It is exported for testing purposes.
53 func (impl Implementation) Dlaln2(trans bool, na, nw int, smin, ca float64, a []float64, lda int, d1, d2 float64, b []float64, ldb int, wr, wi float64, x []float64, ldx int) (scale, xnorm float64, ok bool) {
54 // TODO(vladimir-ch): Consider splitting this function into two, one
55 // handling the real case (nw == 1) and the other handling the complex
56 // case (nw == 2). Given that Go has complex types, their signatures
57 // would be simpler and more natural, and the implementation not as
60 if na != 1 && na != 2 {
61 panic("lapack: invalid value of na")
63 if nw != 1 && nw != 2 {
64 panic("lapack: invalid value of nw")
66 checkMatrix(na, na, a, lda)
67 checkMatrix(na, nw, b, ldb)
68 checkMatrix(na, nw, x, ldx)
72 smini := math.Max(smin, smlnum)
78 // 1×1 (i.e., scalar) system C X = B.
84 csr := ca*a[0] - wr*d1
85 cnorm := math.Abs(csr)
87 // If |C| < smini, use C = smini.
94 // Check scaling for X = B / C.
95 bnorm := math.Abs(b[0])
96 if cnorm < 1 && bnorm > math.Max(1, bignum*cnorm) {
101 x[0] = b[0] * scale / csr
102 xnorm = math.Abs(x[0])
104 return scale, xnorm, ok
107 // Complex 1×1 system (w is complex).
110 csr := ca*a[0] - wr*d1
112 cnorm := math.Abs(csr) + math.Abs(csi)
114 // If |C| < smini, use C = smini.
122 // Check scaling for X = B / C.
123 bnorm := math.Abs(b[0]) + math.Abs(b[1])
124 if cnorm < 1 && bnorm > math.Max(1, bignum*cnorm) {
129 cx := complex(scale*b[0], scale*b[1]) / complex(csr, csi)
130 x[0], x[1] = real(cx), imag(cx)
131 xnorm = math.Abs(x[0]) + math.Abs(x[1])
133 return scale, xnorm, ok
138 // Compute the real part of
161 // Real 2×2 system (w is real).
163 // Find the largest element in C.
166 for j, v := range crv {
174 // If norm(C) < smini, use smini*identity.
176 bnorm := math.Max(math.Abs(b[0]), math.Abs(b[ldb]))
177 if smini < 1 && bnorm > math.Max(1, bignum*smini) {
180 temp := scale / smini
182 x[ldx] = temp * b[ldb]
186 return scale, xnorm, ok
189 // Gaussian elimination with complete pivoting.
190 // Form upper triangular matrix
194 ur12 := crv[pivot[icmax][1]]
195 cr21 := crv[pivot[icmax][2]]
196 cr22 := crv[pivot[icmax][3]]
199 ur22 := cr22 - ur12*lr21
201 // If smaller pivot < smini, use smini.
202 if math.Abs(ur22) < smini {
209 // If the pivot lies in the second row, swap the rows.
216 br2 -= lr21 * br1 // Apply the Gaussian elimination step to the right-hand side.
218 bbnd := math.Max(math.Abs(ur22*ur11r*br1), math.Abs(br2))
219 if bbnd > 1 && math.Abs(ur22) < 1 && bbnd >= bignum*math.Abs(ur22) {
223 // Solve the linear system ur*xr=br.
224 xr2 := br2 * scale / ur22
225 xr1 := scale*br1*ur11r - ur11r*ur12*xr2
227 // If the pivot lies in the second column, swap the components of the solution.
234 xnorm = math.Max(math.Abs(xr1), math.Abs(xr2))
236 // Further scaling if norm(A)*norm(X) > overflow.
237 if xnorm > 1 && cmax > 1 && xnorm > bignum/cmax {
238 temp := cmax / bignum
245 return scale, xnorm, ok
248 // Complex 2×2 system (w is complex).
250 // Find the largest element in C.
259 for j, v := range crv {
261 if v+math.Abs(civ[j]) > cmax {
262 cmax = v + math.Abs(civ[j])
267 // If norm(C) < smini, use smini*identity.
269 br1 := math.Abs(b[0]) + math.Abs(b[1])
270 br2 := math.Abs(b[ldb]) + math.Abs(b[ldb+1])
271 bnorm := math.Max(br1, br2)
272 if smini < 1 && bnorm > 1 && bnorm > bignum*smini {
275 temp := scale / smini
278 x[ldb] = temp * b[ldb]
279 x[ldb+1] = temp * b[ldb+1]
283 return scale, xnorm, ok
286 // Gaussian elimination with complete pivoting.
289 ur12 := crv[pivot[icmax][1]]
290 ui12 := civ[pivot[icmax][1]]
291 cr21 := crv[pivot[icmax][2]]
292 ci21 := civ[pivot[icmax][2]]
293 cr22 := crv[pivot[icmax][3]]
294 ci22 := civ[pivot[icmax][3]]
301 if icmax == 0 || icmax == 3 {
302 // Off-diagonals of pivoted C are real.
303 if math.Abs(ur11) > math.Abs(ui11) {
305 ur11r = 1 / (ur11 * (1 + temp*temp))
306 ui11r = -temp * ur11r
309 ui11r = -1 / (ui11 * (1 + temp*temp))
310 ur11r = -temp * ui11r
316 ur22 = cr22 - ur12*lr21
317 ui22 = ci22 - ur12*li21
319 // Diagonals of pivoted C are real.
321 // ui11r is already 0.
326 ur22 = cr22 - ur12*lr21 + ui12*li21
327 ui22 = -ur12*li21 - ui12*lr21
329 u22abs := math.Abs(ur22) + math.Abs(ui22)
331 // If smaller pivot < smini, use smini.
341 // If the pivot lies in the second row, swap the rows.
352 br2 += -lr21*br1 + li21*bi1
353 bi2 += -li21*br1 - lr21*bi1
355 bbnd1 := u22abs * (math.Abs(ur11r) + math.Abs(ui11r)) * (math.Abs(br1) + math.Abs(bi1))
356 bbnd2 := math.Abs(br2) + math.Abs(bi2)
357 bbnd := math.Max(bbnd1, bbnd2)
358 if bbnd > 1 && u22abs < 1 && bbnd >= bignum*u22abs {
366 cx2 := complex(br2, bi2) / complex(ur22, ui22)
367 xr2, xi2 := real(cx2), imag(cx2)
368 xr1 := ur11r*br1 - ui11r*bi1 - ur12s*xr2 + ui12s*xi2
369 xi1 := ui11r*br1 + ur11r*bi1 - ui12s*xr2 - ur12s*xi2
371 // If the pivot lies in the second column, swap the components of the solution.
382 xnorm = math.Max(math.Abs(xr1)+math.Abs(xi1), math.Abs(xr2)+math.Abs(xi2))
384 // Further scaling if norm(A)*norm(X) > overflow.
385 if xnorm > 1 && cmax > 1 && xnorm > bignum/cmax {
386 temp := cmax / bignum
395 return scale, xnorm, ok