github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/iparmq.go (about)

     1  // Copyright ©2016 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 "math"
     8  
     9  // Iparmq returns problem and machine dependent parameters useful for Dhseqr and
    10  // related subroutines for eigenvalue problems.
    11  //
    12  // ispec specifies the parameter to return:
    13  //  12: Crossover point between Dlahqr and Dlaqr0. Will be at least 11.
    14  //  13: Deflation window size.
    15  //  14: Nibble crossover point. Determines when to skip a multi-shift QR sweep.
    16  //  15: Number of simultaneous shifts in a multishift QR iteration.
    17  //  16: Select structured matrix multiply.
    18  // For other values of ispec Iparmq will panic.
    19  //
    20  // name is the name of the calling function. name must be in uppercase but this
    21  // is not checked.
    22  //
    23  // opts is not used and exists for future use.
    24  //
    25  // n is the order of the Hessenberg matrix H.
    26  //
    27  // ilo and ihi specify the block [ilo:ihi+1,ilo:ihi+1] that is being processed.
    28  //
    29  // lwork is the amount of workspace available.
    30  //
    31  // Except for ispec input parameters are not checked.
    32  //
    33  // Iparmq is an internal routine. It is exported for testing purposes.
    34  func (Implementation) Iparmq(ispec int, name, opts string, n, ilo, ihi, lwork int) int {
    35  	nh := ihi - ilo + 1
    36  	ns := 2
    37  	switch {
    38  	case nh >= 30:
    39  		ns = 4
    40  	case nh >= 60:
    41  		ns = 10
    42  	case nh >= 150:
    43  		ns = max(10, nh/int(math.Log(float64(nh))/math.Ln2))
    44  	case nh >= 590:
    45  		ns = 64
    46  	case nh >= 3000:
    47  		ns = 128
    48  	case nh >= 6000:
    49  		ns = 256
    50  	}
    51  	ns = max(2, ns-(ns%2))
    52  
    53  	switch ispec {
    54  	default:
    55  		panic(badIspec)
    56  
    57  	case 12:
    58  		// Matrices of order smaller than nmin get sent to Dlahqr, the
    59  		// classic double shift algorithm. This must be at least 11.
    60  		const nmin = 75
    61  		return nmin
    62  
    63  	case 13:
    64  		const knwswp = 500
    65  		if nh <= knwswp {
    66  			return ns
    67  		}
    68  		return 3 * ns / 2
    69  
    70  	case 14:
    71  		// Skip a computationally expensive multi-shift QR sweep with
    72  		// Dlaqr5 whenever aggressive early deflation finds at least
    73  		// nibble*(window size)/100 deflations. The default, small,
    74  		// value reflects the expectation that the cost of looking
    75  		// through the deflation window with Dlaqr3 will be
    76  		// substantially smaller.
    77  		const nibble = 14
    78  		return nibble
    79  
    80  	case 15:
    81  		return ns
    82  
    83  	case 16:
    84  		if len(name) != 6 {
    85  			panic(badName)
    86  		}
    87  		const (
    88  			k22min = 14
    89  			kacmin = 14
    90  		)
    91  		var acc22 int
    92  		switch {
    93  		case name[1:] == "GGHRD" || name[1:] == "GGHD3":
    94  			acc22 = 1
    95  			if nh >= k22min {
    96  				acc22 = 2
    97  			}
    98  		case name[3:] == "EXC":
    99  			if nh >= kacmin {
   100  				acc22 = 1
   101  			}
   102  			if nh >= k22min {
   103  				acc22 = 2
   104  			}
   105  		case name[1:] == "HSEQR" || name[1:5] == "LAQR":
   106  			if ns >= kacmin {
   107  				acc22 = 1
   108  			}
   109  			if ns >= k22min {
   110  				acc22 = 2
   111  			}
   112  		}
   113  		return acc22
   114  	}
   115  }