gonum.org/v1/gonum@v0.14.0/lapack/gonum/dlacn2.go (about)

     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  //
    21  //	A * X    if kase was returned as 1,
    22  //	Aᵀ * X   if kase was returned as 2,
    23  //
    24  // and all other parameters must not be changed.
    25  // On the final return, kase is returned as 0, v contains A*W where W is a
    26  // vector, and est = norm(V)/norm(W) is a lower bound for 1-norm of A.
    27  //
    28  // v, x, and isgn must all have length n and n must be at least 1, otherwise
    29  // Dlacn2 will panic. isave is used for temporary storage.
    30  //
    31  // Dlacn2 is an internal routine. It is exported for testing purposes.
    32  func (impl Implementation) Dlacn2(n int, v, x []float64, isgn []int, est float64, kase int, isave *[3]int) (float64, int) {
    33  	switch {
    34  	case n < 1:
    35  		panic(nLT1)
    36  	case len(v) < n:
    37  		panic(shortV)
    38  	case len(x) < n:
    39  		panic(shortX)
    40  	case len(isgn) < n:
    41  		panic(shortIsgn)
    42  	case isave[0] < 0 || 5 < isave[0]:
    43  		panic(badIsave)
    44  	case isave[0] == 0 && kase != 0:
    45  		panic(badIsave)
    46  	}
    47  
    48  	const itmax = 5
    49  	bi := blas64.Implementation()
    50  
    51  	if kase == 0 {
    52  		for i := 0; i < n; i++ {
    53  			x[i] = 1 / float64(n)
    54  		}
    55  		kase = 1
    56  		isave[0] = 1
    57  		return est, kase
    58  	}
    59  	switch isave[0] {
    60  	case 1:
    61  		if n == 1 {
    62  			v[0] = x[0]
    63  			est = math.Abs(v[0])
    64  			kase = 0
    65  			return est, kase
    66  		}
    67  		est = bi.Dasum(n, x, 1)
    68  		for i := 0; i < n; i++ {
    69  			x[i] = math.Copysign(1, x[i])
    70  			isgn[i] = int(x[i])
    71  		}
    72  		kase = 2
    73  		isave[0] = 2
    74  		return est, kase
    75  	case 2:
    76  		isave[1] = bi.Idamax(n, x, 1)
    77  		isave[2] = 2
    78  		for i := 0; i < n; i++ {
    79  			x[i] = 0
    80  		}
    81  		x[isave[1]] = 1
    82  		kase = 1
    83  		isave[0] = 3
    84  		return est, kase
    85  	case 3:
    86  		bi.Dcopy(n, x, 1, v, 1)
    87  		estold := est
    88  		est = bi.Dasum(n, v, 1)
    89  		sameSigns := true
    90  		for i := 0; i < n; i++ {
    91  			if int(math.Copysign(1, x[i])) != isgn[i] {
    92  				sameSigns = false
    93  				break
    94  			}
    95  		}
    96  		if !sameSigns && est > estold {
    97  			for i := 0; i < n; i++ {
    98  				x[i] = math.Copysign(1, x[i])
    99  				isgn[i] = int(x[i])
   100  			}
   101  			kase = 2
   102  			isave[0] = 4
   103  			return est, kase
   104  		}
   105  	case 4:
   106  		jlast := isave[1]
   107  		isave[1] = bi.Idamax(n, x, 1)
   108  		if x[jlast] != math.Abs(x[isave[1]]) && isave[2] < itmax {
   109  			isave[2] += 1
   110  			for i := 0; i < n; i++ {
   111  				x[i] = 0
   112  			}
   113  			x[isave[1]] = 1
   114  			kase = 1
   115  			isave[0] = 3
   116  			return est, kase
   117  		}
   118  	case 5:
   119  		tmp := 2 * (bi.Dasum(n, x, 1)) / float64(3*n)
   120  		if tmp > est {
   121  			bi.Dcopy(n, x, 1, v, 1)
   122  			est = tmp
   123  		}
   124  		kase = 0
   125  		return est, kase
   126  	}
   127  	// Iteration complete. Final stage
   128  	altsgn := 1.0
   129  	for i := 0; i < n; i++ {
   130  		x[i] = altsgn * (1 + float64(i)/float64(n-1))
   131  		altsgn *= -1
   132  	}
   133  	kase = 1
   134  	isave[0] = 5
   135  	return est, kase
   136  }