+++ /dev/null
-// 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/blas64"
-)
-
-// Dlacn2 estimates the 1-norm of an n×n matrix A using sequential updates with
-// matrix-vector products provided externally.
-//
-// Dlacn2 is called sequentially and it returns the value of est and kase to be
-// used on the next call.
-// On the initial call, kase must be 0.
-// In between calls, x must be overwritten by
-// A * X if kase was returned as 1,
-// A^T * X if kase was returned as 2,
-// and all other parameters must not be changed.
-// On the final return, kase is returned as 0, v contains A*W where W is a
-// vector, and est = norm(V)/norm(W) is a lower bound for 1-norm of A.
-//
-// v, x, and isgn must all have length n and n must be at least 1, otherwise
-// Dlacn2 will panic. isave is used for temporary storage.
-//
-// Dlacn2 is an internal routine. It is exported for testing purposes.
-func (impl Implementation) Dlacn2(n int, v, x []float64, isgn []int, est float64, kase int, isave *[3]int) (float64, int) {
- if n < 1 {
- panic("lapack: non-positive n")
- }
- checkVector(n, x, 1)
- checkVector(n, v, 1)
- if len(isgn) < n {
- panic("lapack: insufficient isgn length")
- }
- if isave[0] < 0 || isave[0] > 5 {
- panic("lapack: bad isave value")
- }
- if isave[0] == 0 && kase != 0 {
- panic("lapack: bad isave value")
- }
- itmax := 5
- bi := blas64.Implementation()
- if kase == 0 {
- for i := 0; i < n; i++ {
- x[i] = 1 / float64(n)
- }
- kase = 1
- isave[0] = 1
- return est, kase
- }
- switch isave[0] {
- default:
- panic("unreachable")
- case 1:
- if n == 1 {
- v[0] = x[0]
- est = math.Abs(v[0])
- kase = 0
- return est, kase
- }
- est = bi.Dasum(n, x, 1)
- for i := 0; i < n; i++ {
- x[i] = math.Copysign(1, x[i])
- isgn[i] = int(x[i])
- }
- kase = 2
- isave[0] = 2
- return est, kase
- case 2:
- isave[1] = bi.Idamax(n, x, 1)
- isave[2] = 2
- for i := 0; i < n; i++ {
- x[i] = 0
- }
- x[isave[1]] = 1
- kase = 1
- isave[0] = 3
- return est, kase
- case 3:
- bi.Dcopy(n, x, 1, v, 1)
- estold := est
- est = bi.Dasum(n, v, 1)
- sameSigns := true
- for i := 0; i < n; i++ {
- if int(math.Copysign(1, x[i])) != isgn[i] {
- sameSigns = false
- break
- }
- }
- if !sameSigns && est > estold {
- for i := 0; i < n; i++ {
- x[i] = math.Copysign(1, x[i])
- isgn[i] = int(x[i])
- }
- kase = 2
- isave[0] = 4
- return est, kase
- }
- case 4:
- jlast := isave[1]
- isave[1] = bi.Idamax(n, x, 1)
- if x[jlast] != math.Abs(x[isave[1]]) && isave[2] < itmax {
- isave[2] += 1
- for i := 0; i < n; i++ {
- x[i] = 0
- }
- x[isave[1]] = 1
- kase = 1
- isave[0] = 3
- return est, kase
- }
- case 5:
- tmp := 2 * (bi.Dasum(n, x, 1)) / float64(3*n)
- if tmp > est {
- bi.Dcopy(n, x, 1, v, 1)
- est = tmp
- }
- kase = 0
- return est, kase
- }
- // Iteration complete. Final stage
- altsgn := 1.0
- for i := 0; i < n; i++ {
- x[i] = altsgn * (1 + float64(i)/float64(n-1))
- altsgn *= -1
- }
- kase = 1
- isave[0] = 5
- return est, kase
-}
OSDN Git Service
