blob: d97e5edf865730f43d9073840ab0701a2c5fe411 [file] [log] [blame]
// 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ᵀ * 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) {
switch {
case n < 1:
panic(nLT1)
case len(v) < n:
panic(shortV)
case len(x) < n:
panic(shortX)
case len(isgn) < n:
panic(shortIsgn)
case isave[0] < 0 || 5 < isave[0]:
panic(badIsave)
case isave[0] == 0 && kase != 0:
panic(badIsave)
}
const 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] {
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
}