X-Git-Url: http://git.osdn.net/view?p=bytom%2Fvapor.git;a=blobdiff_plain;f=vendor%2Fgonum.org%2Fv1%2Fgonum%2Flapack%2Fgonum%2Fdlatrs.go;fp=vendor%2Fgonum.org%2Fv1%2Fgonum%2Flapack%2Fgonum%2Fdlatrs.go;h=0000000000000000000000000000000000000000;hp=f0c94764e6a9db41739f3e7737bba6077cfc1baf;hb=2cf5801b2e693a45de9b51ec9aa9c1f787d57105;hpb=0dff3fcf4fbd306176d561d721c1c31e58d90742 diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlatrs.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlatrs.go deleted file mode 100644 index f0c94764..00000000 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dlatrs.go +++ /dev/null @@ -1,350 +0,0 @@ -// Copyright ©2015 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" - - "gonum.org/v1/gonum/blas" - "gonum.org/v1/gonum/blas/blas64" -) - -// Dlatrs solves a triangular system of equations scaled to prevent overflow. It -// solves -// A * x = scale * b if trans == blas.NoTrans -// A^T * x = scale * b if trans == blas.Trans -// where the scale s is set for numeric stability. -// -// A is an n×n triangular matrix. On entry, the slice x contains the values of -// of b, and on exit it contains the solution vector x. -// -// If normin == true, cnorm is an input and cnorm[j] contains the norm of the off-diagonal -// part of the j^th column of A. If trans == blas.NoTrans, cnorm[j] must be greater -// than or equal to the infinity norm, and greater than or equal to the one-norm -// otherwise. If normin == false, then cnorm is treated as an output, and is set -// to contain the 1-norm of the off-diagonal part of the j^th column of A. -// -// Dlatrs is an internal routine. It is exported for testing purposes. -func (impl Implementation) Dlatrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, normin bool, n int, a []float64, lda int, x []float64, cnorm []float64) (scale float64) { - if uplo != blas.Upper && uplo != blas.Lower { - panic(badUplo) - } - if trans != blas.Trans && trans != blas.NoTrans { - panic(badTrans) - } - if diag != blas.Unit && diag != blas.NonUnit { - panic(badDiag) - } - upper := uplo == blas.Upper - noTrans := trans == blas.NoTrans - nonUnit := diag == blas.NonUnit - - if n < 0 { - panic(nLT0) - } - checkMatrix(n, n, a, lda) - checkVector(n, x, 1) - checkVector(n, cnorm, 1) - - if n == 0 { - return 0 - } - smlnum := dlamchS / dlamchP - bignum := 1 / smlnum - scale = 1 - bi := blas64.Implementation() - if !normin { - if upper { - cnorm[0] = 0 - for j := 1; j < n; j++ { - cnorm[j] = bi.Dasum(j, a[j:], lda) - } - } else { - for j := 0; j < n-1; j++ { - cnorm[j] = bi.Dasum(n-j-1, a[(j+1)*lda+j:], lda) - } - cnorm[n-1] = 0 - } - } - // Scale the column norms by tscal if the maximum element in cnorm is greater than bignum. - imax := bi.Idamax(n, cnorm, 1) - tmax := cnorm[imax] - var tscal float64 - if tmax <= bignum { - tscal = 1 - } else { - tscal = 1 / (smlnum * tmax) - bi.Dscal(n, tscal, cnorm, 1) - } - - // Compute a bound on the computed solution vector to see if bi.Dtrsv can be used. - j := bi.Idamax(n, x, 1) - xmax := math.Abs(x[j]) - xbnd := xmax - var grow float64 - var jfirst, jlast, jinc int - if noTrans { - if upper { - jfirst = n - 1 - jlast = -1 - jinc = -1 - } else { - jfirst = 0 - jlast = n - jinc = 1 - } - // Compute the growth in A * x = b. - if tscal != 1 { - grow = 0 - goto Solve - } - if nonUnit { - grow = 1 / math.Max(xbnd, smlnum) - xbnd = grow - for j := jfirst; j != jlast; j += jinc { - if grow <= smlnum { - goto Solve - } - tjj := math.Abs(a[j*lda+j]) - xbnd = math.Min(xbnd, math.Min(1, tjj)*grow) - if tjj+cnorm[j] >= smlnum { - grow *= tjj / (tjj + cnorm[j]) - } else { - grow = 0 - } - } - grow = xbnd - } else { - grow = math.Min(1, 1/math.Max(xbnd, smlnum)) - for j := jfirst; j != jlast; j += jinc { - if grow <= smlnum { - goto Solve - } - grow *= 1 / (1 + cnorm[j]) - } - } - } else { - if upper { - jfirst = 0 - jlast = n - jinc = 1 - } else { - jfirst = n - 1 - jlast = -1 - jinc = -1 - } - if tscal != 1 { - grow = 0 - goto Solve - } - if nonUnit { - grow = 1 / (math.Max(xbnd, smlnum)) - xbnd = grow - for j := jfirst; j != jlast; j += jinc { - if grow <= smlnum { - goto Solve - } - xj := 1 + cnorm[j] - grow = math.Min(grow, xbnd/xj) - tjj := math.Abs(a[j*lda+j]) - if xj > tjj { - xbnd *= tjj / xj - } - } - grow = math.Min(grow, xbnd) - } else { - grow = math.Min(1, 1/math.Max(xbnd, smlnum)) - for j := jfirst; j != jlast; j += jinc { - if grow <= smlnum { - goto Solve - } - xj := 1 + cnorm[j] - grow /= xj - } - } - } - -Solve: - if grow*tscal > smlnum { - // Use the Level 2 BLAS solve if the reciprocal of the bound on - // elements of X is not too small. - bi.Dtrsv(uplo, trans, diag, n, a, lda, x, 1) - if tscal != 1 { - bi.Dscal(n, 1/tscal, cnorm, 1) - } - return scale - } - - // Use a Level 1 BLAS solve, scaling intermediate results. - if xmax > bignum { - scale = bignum / xmax - bi.Dscal(n, scale, x, 1) - xmax = bignum - } - if noTrans { - for j := jfirst; j != jlast; j += jinc { - xj := math.Abs(x[j]) - var tjj, tjjs float64 - if nonUnit { - tjjs = a[j*lda+j] * tscal - } else { - tjjs = tscal - if tscal == 1 { - goto Skip1 - } - } - tjj = math.Abs(tjjs) - if tjj > smlnum { - if tjj < 1 { - if xj > tjj*bignum { - rec := 1 / xj - bi.Dscal(n, rec, x, 1) - scale *= rec - xmax *= rec - } - } - x[j] /= tjjs - xj = math.Abs(x[j]) - } else if tjj > 0 { - if xj > tjj*bignum { - rec := (tjj * bignum) / xj - if cnorm[j] > 1 { - rec /= cnorm[j] - } - bi.Dscal(n, rec, x, 1) - scale *= rec - xmax *= rec - } - x[j] /= tjjs - xj = math.Abs(x[j]) - } else { - for i := 0; i < n; i++ { - x[i] = 0 - } - x[j] = 1 - xj = 1 - scale = 0 - xmax = 0 - } - Skip1: - if xj > 1 { - rec := 1 / xj - if cnorm[j] > (bignum-xmax)*rec { - rec *= 0.5 - bi.Dscal(n, rec, x, 1) - scale *= rec - } - } else if xj*cnorm[j] > bignum-xmax { - bi.Dscal(n, 0.5, x, 1) - scale *= 0.5 - } - if upper { - if j > 0 { - bi.Daxpy(j, -x[j]*tscal, a[j:], lda, x, 1) - i := bi.Idamax(j, x, 1) - xmax = math.Abs(x[i]) - } - } else { - if j < n-1 { - bi.Daxpy(n-j-1, -x[j]*tscal, a[(j+1)*lda+j:], lda, x[j+1:], 1) - i := j + bi.Idamax(n-j-1, x[j+1:], 1) - xmax = math.Abs(x[i]) - } - } - } - } else { - for j := jfirst; j != jlast; j += jinc { - xj := math.Abs(x[j]) - uscal := tscal - rec := 1 / math.Max(xmax, 1) - var tjjs float64 - if cnorm[j] > (bignum-xj)*rec { - rec *= 0.5 - if nonUnit { - tjjs = a[j*lda+j] * tscal - } else { - tjjs = tscal - } - tjj := math.Abs(tjjs) - if tjj > 1 { - rec = math.Min(1, rec*tjj) - uscal /= tjjs - } - if rec < 1 { - bi.Dscal(n, rec, x, 1) - scale *= rec - xmax *= rec - } - } - var sumj float64 - if uscal == 1 { - if upper { - sumj = bi.Ddot(j, a[j:], lda, x, 1) - } else if j < n-1 { - sumj = bi.Ddot(n-j-1, a[(j+1)*lda+j:], lda, x[j+1:], 1) - } - } else { - if upper { - for i := 0; i < j; i++ { - sumj += (a[i*lda+j] * uscal) * x[i] - } - } else if j < n { - for i := j + 1; i < n; i++ { - sumj += (a[i*lda+j] * uscal) * x[i] - } - } - } - if uscal == tscal { - x[j] -= sumj - xj := math.Abs(x[j]) - var tjjs float64 - if nonUnit { - tjjs = a[j*lda+j] * tscal - } else { - tjjs = tscal - if tscal == 1 { - goto Skip2 - } - } - tjj := math.Abs(tjjs) - if tjj > smlnum { - if tjj < 1 { - if xj > tjj*bignum { - rec = 1 / xj - bi.Dscal(n, rec, x, 1) - scale *= rec - xmax *= rec - } - } - x[j] /= tjjs - } else if tjj > 0 { - if xj > tjj*bignum { - rec = (tjj * bignum) / xj - bi.Dscal(n, rec, x, 1) - scale *= rec - xmax *= rec - } - x[j] /= tjjs - } else { - for i := 0; i < n; i++ { - x[i] = 0 - } - x[j] = 1 - scale = 0 - xmax = 0 - } - } else { - x[j] = x[j]/tjjs - sumj - } - Skip2: - xmax = math.Max(xmax, math.Abs(x[j])) - } - } - scale /= tscal - if tscal != 1 { - bi.Dscal(n, 1/tscal, cnorm, 1) - } - return scale -}