vendor/gonum.org/v1/gonum/lapack/gonum/dlacn2.go

   1 // Copyright ©2015 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.
   4
   5 package gonum
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas/blas64"
  11 )
  12
  13 // Dlacn2 estimates the 1-norm of an n×n matrix A using sequential updates with
  14 // matrix-vector products provided externally.
  15 //
  16 // Dlacn2 is called sequentially and it returns the value of est and kase to be
  17 // used on the next call.
  18 // On the initial call, kase must be 0.
  19 // In between calls, x must be overwritten by
  20 //  A * X    if kase was returned as 1,
  21 //  A^T * X  if kase was returned as 2,
  22 // and all other parameters must not be changed.
  23 // On the final return, kase is returned as 0, v contains A*W where W is a
  24 // vector, and est = norm(V)/norm(W) is a lower bound for 1-norm of A.
  25 //
  26 // v, x, and isgn must all have length n and n must be at least 1, otherwise
  27 // Dlacn2 will panic. isave is used for temporary storage.
  28 //
  29 // Dlacn2 is an internal routine. It is exported for testing purposes.
  30 func (impl Implementation) Dlacn2(n int, v, x []float64, isgn []int, est float64, kase int, isave *[3]int) (float64, int) {
  31         if n < 1 {
  32                 panic("lapack: non-positive n")
  33         }
  34         checkVector(n, x, 1)
  35         checkVector(n, v, 1)
  36         if len(isgn) < n {
  37                 panic("lapack: insufficient isgn length")
  38         }
  39         if isave[0] < 0 || isave[0] > 5 {
  40                 panic("lapack: bad isave value")
  41         }
  42         if isave[0] == 0 && kase != 0 {
  43                 panic("lapack: bad isave value")
  44         }
  45         itmax := 5
  46         bi := blas64.Implementation()
  47         if kase == 0 {
  48                 for i := 0; i < n; i++ {
  49                         x[i] = 1 / float64(n)
  50                 }
  51                 kase = 1
  52                 isave[0] = 1
  53                 return est, kase
  54         }
  55         switch isave[0] {
  56         default:
  57                 panic("unreachable")
  58         case 1:
  59                 if n == 1 {
  60                         v[0] = x[0]
  61                         est = math.Abs(v[0])
  62                         kase = 0
  63                         return est, kase
  64                 }
  65                 est = bi.Dasum(n, x, 1)
  66                 for i := 0; i < n; i++ {
  67                         x[i] = math.Copysign(1, x[i])
  68                         isgn[i] = int(x[i])
  69                 }
  70                 kase = 2
  71                 isave[0] = 2
  72                 return est, kase
  73         case 2:
  74                 isave[1] = bi.Idamax(n, x, 1)
  75                 isave[2] = 2
  76                 for i := 0; i < n; i++ {
  77                         x[i] = 0
  78                 }
  79                 x[isave[1]] = 1
  80                 kase = 1
  81                 isave[0] = 3
  82                 return est, kase
  83         case 3:
  84                 bi.Dcopy(n, x, 1, v, 1)
  85                 estold := est
  86                 est = bi.Dasum(n, v, 1)
  87                 sameSigns := true
  88                 for i := 0; i < n; i++ {
  89                         if int(math.Copysign(1, x[i])) != isgn[i] {
  90                                 sameSigns = false
  91                                 break
  92                         }
  93                 }
  94                 if !sameSigns && est > estold {
  95                         for i := 0; i < n; i++ {
  96                                 x[i] = math.Copysign(1, x[i])
  97                                 isgn[i] = int(x[i])
  98                         }
  99                         kase = 2
 100                         isave[0] = 4
 101                         return est, kase
 102                 }
 103         case 4:
 104                 jlast := isave[1]
 105                 isave[1] = bi.Idamax(n, x, 1)
 106                 if x[jlast] != math.Abs(x[isave[1]]) && isave[2] < itmax {
 107                         isave[2] += 1
 108                         for i := 0; i < n; i++ {
 109                                 x[i] = 0
 110                         }
 111                         x[isave[1]] = 1
 112                         kase = 1
 113                         isave[0] = 3
 114                         return est, kase
 115                 }
 116         case 5:
 117                 tmp := 2 * (bi.Dasum(n, x, 1)) / float64(3*n)
 118                 if tmp > est {
 119                         bi.Dcopy(n, x, 1, v, 1)
 120                         est = tmp
 121                 }
 122                 kase = 0
 123                 return est, kase
 124         }
 125         // Iteration complete. Final stage
 126         altsgn := 1.0
 127         for i := 0; i < n; i++ {
 128                 x[i] = altsgn * (1 + float64(i)/float64(n-1))
 129                 altsgn *= -1
 130         }
 131         kase = 1
 132         isave[0] = 5
 133         return est, kase
 134 }