| // Copyright ©2016 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" |
| |
| // Dlaln2 solves a linear equation or a system of 2 linear equations of the form |
| // (ca A - w D) X = scale B, if trans == false, |
| // (ca A^T - w D) X = scale B, if trans == true, |
| // where A is a na×na real matrix, ca is a real scalar, D is a na×na diagonal |
| // real matrix, w is a scalar, real if nw == 1, complex if nw == 2, and X and B |
| // are na×1 matrices, real if w is real, complex if w is complex. |
| // |
| // If w is complex, X and B are represented as na×2 matrices, the first column |
| // of each being the real part and the second being the imaginary part. |
| // |
| // na and nw must be 1 or 2, otherwise Dlaln2 will panic. |
| // |
| // d1 and d2 are the diagonal elements of D. d2 is not used if na == 1. |
| // |
| // wr and wi represent the real and imaginary part, respectively, of the scalar |
| // w. wi is not used if nw == 1. |
| // |
| // smin is the desired lower bound on the singular values of A. This should be |
| // a safe distance away from underflow or overflow, say, between |
| // (underflow/machine precision) and (overflow*machine precision). |
| // |
| // If both singular values of (ca A - w D) are less than smin, smin*identity |
| // will be used instead of (ca A - w D). If only one singular value is less than |
| // smin, one element of (ca A - w D) will be perturbed enough to make the |
| // smallest singular value roughly smin. If both singular values are at least |
| // smin, (ca A - w D) will not be perturbed. In any case, the perturbation will |
| // be at most some small multiple of max(smin, ulp*norm(ca A - w D)). The |
| // singular values are computed by infinity-norm approximations, and thus will |
| // only be correct to a factor of 2 or so. |
| // |
| // All input quantities are assumed to be smaller than overflow by a reasonable |
| // factor. |
| // |
| // scale is a scaling factor less than or equal to 1 which is chosen so that X |
| // can be computed without overflow. X is further scaled if necessary to assure |
| // that norm(ca A - w D)*norm(X) is less than overflow. |
| // |
| // xnorm contains the infinity-norm of X when X is regarded as a na×nw real |
| // matrix. |
| // |
| // ok will be false if (ca A - w D) had to be perturbed to make its smallest |
| // singular value greater than smin, otherwise ok will be true. |
| // |
| // Dlaln2 is an internal routine. It is exported for testing purposes. |
| 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) { |
| // TODO(vladimir-ch): Consider splitting this function into two, one |
| // handling the real case (nw == 1) and the other handling the complex |
| // case (nw == 2). Given that Go has complex types, their signatures |
| // would be simpler and more natural, and the implementation not as |
| // convoluted. |
| |
| if na != 1 && na != 2 { |
| panic("lapack: invalid value of na") |
| } |
| if nw != 1 && nw != 2 { |
| panic("lapack: invalid value of nw") |
| } |
| checkMatrix(na, na, a, lda) |
| checkMatrix(na, nw, b, ldb) |
| checkMatrix(na, nw, x, ldx) |
| |
| smlnum := 2 * dlamchS |
| bignum := 1 / smlnum |
| smini := math.Max(smin, smlnum) |
| |
| ok = true |
| scale = 1 |
| |
| if na == 1 { |
| // 1×1 (i.e., scalar) system C X = B. |
| |
| if nw == 1 { |
| // Real 1×1 system. |
| |
| // C = ca A - w D. |
| csr := ca*a[0] - wr*d1 |
| cnorm := math.Abs(csr) |
| |
| // If |C| < smini, use C = smini. |
| if cnorm < smini { |
| csr = smini |
| cnorm = smini |
| ok = false |
| } |
| |
| // Check scaling for X = B / C. |
| bnorm := math.Abs(b[0]) |
| if cnorm < 1 && bnorm > math.Max(1, bignum*cnorm) { |
| scale = 1 / bnorm |
| } |
| |
| // Compute X. |
| x[0] = b[0] * scale / csr |
| xnorm = math.Abs(x[0]) |
| |
| return scale, xnorm, ok |
| } |
| |
| // Complex 1×1 system (w is complex). |
| |
| // C = ca A - w D. |
| csr := ca*a[0] - wr*d1 |
| csi := -wi * d1 |
| cnorm := math.Abs(csr) + math.Abs(csi) |
| |
| // If |C| < smini, use C = smini. |
| if cnorm < smini { |
| csr = smini |
| csi = 0 |
| cnorm = smini |
| ok = false |
| } |
| |
| // Check scaling for X = B / C. |
| bnorm := math.Abs(b[0]) + math.Abs(b[1]) |
| if cnorm < 1 && bnorm > math.Max(1, bignum*cnorm) { |
| scale = 1 / bnorm |
| } |
| |
| // Compute X. |
| cx := complex(scale*b[0], scale*b[1]) / complex(csr, csi) |
| x[0], x[1] = real(cx), imag(cx) |
| xnorm = math.Abs(x[0]) + math.Abs(x[1]) |
| |
| return scale, xnorm, ok |
| } |
| |
| // 2×2 system. |
| |
| // Compute the real part of |
| // C = ca A - w D |
| // or |
| // C = ca A^T - w D. |
| crv := [4]float64{ |
| ca*a[0] - wr*d1, |
| ca * a[1], |
| ca * a[lda], |
| ca*a[lda+1] - wr*d2, |
| } |
| if trans { |
| crv[1] = ca * a[lda] |
| crv[2] = ca * a[1] |
| } |
| |
| pivot := [4][4]int{ |
| {0, 1, 2, 3}, |
| {1, 0, 3, 2}, |
| {2, 3, 0, 1}, |
| {3, 2, 1, 0}, |
| } |
| |
| if nw == 1 { |
| // Real 2×2 system (w is real). |
| |
| // Find the largest element in C. |
| var cmax float64 |
| var icmax int |
| for j, v := range crv { |
| v = math.Abs(v) |
| if v > cmax { |
| cmax = v |
| icmax = j |
| } |
| } |
| |
| // If norm(C) < smini, use smini*identity. |
| if cmax < smini { |
| bnorm := math.Max(math.Abs(b[0]), math.Abs(b[ldb])) |
| if smini < 1 && bnorm > math.Max(1, bignum*smini) { |
| scale = 1 / bnorm |
| } |
| temp := scale / smini |
| x[0] = temp * b[0] |
| x[ldx] = temp * b[ldb] |
| xnorm = temp * bnorm |
| ok = false |
| |
| return scale, xnorm, ok |
| } |
| |
| // Gaussian elimination with complete pivoting. |
| // Form upper triangular matrix |
| // [ur11 ur12] |
| // [ 0 ur22] |
| ur11 := crv[icmax] |
| ur12 := crv[pivot[icmax][1]] |
| cr21 := crv[pivot[icmax][2]] |
| cr22 := crv[pivot[icmax][3]] |
| ur11r := 1 / ur11 |
| lr21 := ur11r * cr21 |
| ur22 := cr22 - ur12*lr21 |
| |
| // If smaller pivot < smini, use smini. |
| if math.Abs(ur22) < smini { |
| ur22 = smini |
| ok = false |
| } |
| |
| var br1, br2 float64 |
| if icmax > 1 { |
| // If the pivot lies in the second row, swap the rows. |
| br1 = b[ldb] |
| br2 = b[0] |
| } else { |
| br1 = b[0] |
| br2 = b[ldb] |
| } |
| br2 -= lr21 * br1 // Apply the Gaussian elimination step to the right-hand side. |
| |
| bbnd := math.Max(math.Abs(ur22*ur11r*br1), math.Abs(br2)) |
| if bbnd > 1 && math.Abs(ur22) < 1 && bbnd >= bignum*math.Abs(ur22) { |
| scale = 1 / bbnd |
| } |
| |
| // Solve the linear system ur*xr=br. |
| xr2 := br2 * scale / ur22 |
| xr1 := scale*br1*ur11r - ur11r*ur12*xr2 |
| if icmax&0x1 != 0 { |
| // If the pivot lies in the second column, swap the components of the solution. |
| x[0] = xr2 |
| x[ldx] = xr1 |
| } else { |
| x[0] = xr1 |
| x[ldx] = xr2 |
| } |
| xnorm = math.Max(math.Abs(xr1), math.Abs(xr2)) |
| |
| // Further scaling if norm(A)*norm(X) > overflow. |
| if xnorm > 1 && cmax > 1 && xnorm > bignum/cmax { |
| temp := cmax / bignum |
| x[0] *= temp |
| x[ldx] *= temp |
| xnorm *= temp |
| scale *= temp |
| } |
| |
| return scale, xnorm, ok |
| } |
| |
| // Complex 2×2 system (w is complex). |
| |
| // Find the largest element in C. |
| civ := [4]float64{ |
| -wi * d1, |
| 0, |
| 0, |
| -wi * d2, |
| } |
| var cmax float64 |
| var icmax int |
| for j, v := range crv { |
| v := math.Abs(v) |
| if v+math.Abs(civ[j]) > cmax { |
| cmax = v + math.Abs(civ[j]) |
| icmax = j |
| } |
| } |
| |
| // If norm(C) < smini, use smini*identity. |
| if cmax < smini { |
| br1 := math.Abs(b[0]) + math.Abs(b[1]) |
| br2 := math.Abs(b[ldb]) + math.Abs(b[ldb+1]) |
| bnorm := math.Max(br1, br2) |
| if smini < 1 && bnorm > 1 && bnorm > bignum*smini { |
| scale = 1 / bnorm |
| } |
| temp := scale / smini |
| x[0] = temp * b[0] |
| x[1] = temp * b[1] |
| x[ldb] = temp * b[ldb] |
| x[ldb+1] = temp * b[ldb+1] |
| xnorm = temp * bnorm |
| ok = false |
| |
| return scale, xnorm, ok |
| } |
| |
| // Gaussian elimination with complete pivoting. |
| ur11 := crv[icmax] |
| ui11 := civ[icmax] |
| ur12 := crv[pivot[icmax][1]] |
| ui12 := civ[pivot[icmax][1]] |
| cr21 := crv[pivot[icmax][2]] |
| ci21 := civ[pivot[icmax][2]] |
| cr22 := crv[pivot[icmax][3]] |
| ci22 := civ[pivot[icmax][3]] |
| var ( |
| ur11r, ui11r float64 |
| lr21, li21 float64 |
| ur12s, ui12s float64 |
| ur22, ui22 float64 |
| ) |
| if icmax == 0 || icmax == 3 { |
| // Off-diagonals of pivoted C are real. |
| if math.Abs(ur11) > math.Abs(ui11) { |
| temp := ui11 / ur11 |
| ur11r = 1 / (ur11 * (1 + temp*temp)) |
| ui11r = -temp * ur11r |
| } else { |
| temp := ur11 / ui11 |
| ui11r = -1 / (ui11 * (1 + temp*temp)) |
| ur11r = -temp * ui11r |
| } |
| lr21 = cr21 * ur11r |
| li21 = cr21 * ui11r |
| ur12s = ur12 * ur11r |
| ui12s = ur12 * ui11r |
| ur22 = cr22 - ur12*lr21 |
| ui22 = ci22 - ur12*li21 |
| } else { |
| // Diagonals of pivoted C are real. |
| ur11r = 1 / ur11 |
| // ui11r is already 0. |
| lr21 = cr21 * ur11r |
| li21 = ci21 * ur11r |
| ur12s = ur12 * ur11r |
| ui12s = ui12 * ur11r |
| ur22 = cr22 - ur12*lr21 + ui12*li21 |
| ui22 = -ur12*li21 - ui12*lr21 |
| } |
| u22abs := math.Abs(ur22) + math.Abs(ui22) |
| |
| // If smaller pivot < smini, use smini. |
| if u22abs < smini { |
| ur22 = smini |
| ui22 = 0 |
| ok = false |
| } |
| |
| var br1, bi1 float64 |
| var br2, bi2 float64 |
| if icmax > 1 { |
| // If the pivot lies in the second row, swap the rows. |
| br1 = b[ldb] |
| bi1 = b[ldb+1] |
| br2 = b[0] |
| bi2 = b[1] |
| } else { |
| br1 = b[0] |
| bi1 = b[1] |
| br2 = b[ldb] |
| bi2 = b[ldb+1] |
| } |
| br2 += -lr21*br1 + li21*bi1 |
| bi2 += -li21*br1 - lr21*bi1 |
| |
| bbnd1 := u22abs * (math.Abs(ur11r) + math.Abs(ui11r)) * (math.Abs(br1) + math.Abs(bi1)) |
| bbnd2 := math.Abs(br2) + math.Abs(bi2) |
| bbnd := math.Max(bbnd1, bbnd2) |
| if bbnd > 1 && u22abs < 1 && bbnd >= bignum*u22abs { |
| scale = 1 / bbnd |
| br1 *= scale |
| bi1 *= scale |
| br2 *= scale |
| bi2 *= scale |
| } |
| |
| cx2 := complex(br2, bi2) / complex(ur22, ui22) |
| xr2, xi2 := real(cx2), imag(cx2) |
| xr1 := ur11r*br1 - ui11r*bi1 - ur12s*xr2 + ui12s*xi2 |
| xi1 := ui11r*br1 + ur11r*bi1 - ui12s*xr2 - ur12s*xi2 |
| if icmax&0x1 != 0 { |
| // If the pivot lies in the second column, swap the components of the solution. |
| x[0] = xr2 |
| x[1] = xi2 |
| x[ldx] = xr1 |
| x[ldx+1] = xi1 |
| } else { |
| x[0] = xr1 |
| x[1] = xi1 |
| x[ldx] = xr2 |
| x[ldx+1] = xi2 |
| } |
| xnorm = math.Max(math.Abs(xr1)+math.Abs(xi1), math.Abs(xr2)+math.Abs(xi2)) |
| |
| // Further scaling if norm(A)*norm(X) > overflow. |
| if xnorm > 1 && cmax > 1 && xnorm > bignum/cmax { |
| temp := cmax / bignum |
| x[0] *= temp |
| x[1] *= temp |
| x[ldx] *= temp |
| x[ldx+1] *= temp |
| xnorm *= temp |
| scale *= temp |
| } |
| |
| return scale, xnorm, ok |
| } |