github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/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  	"github.com/jingcheng-WU/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ᵀ * 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  	switch {
    32  	case n < 1:
    33  		panic(nLT1)
    34  	case len(v) < n:
    35  		panic(shortV)
    36  	case len(x) < n:
    37  		panic(shortX)
    38  	case len(isgn) < n:
    39  		panic(shortIsgn)
    40  	case isave[0] < 0 || 5 < isave[0]:
    41  		panic(badIsave)
    42  	case isave[0] == 0 && kase != 0:
    43  		panic(badIsave)
    44  	}
    45  
    46  	const itmax = 5
    47  	bi := blas64.Implementation()
    48  
    49  	if kase == 0 {
    50  		for i := 0; i < n; i++ {
    51  			x[i] = 1 / float64(n)
    52  		}
    53  		kase = 1
    54  		isave[0] = 1
    55  		return est, kase
    56  	}
    57  	switch isave[0] {
    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  }