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

   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 "gonum.org/v1/gonum/lapack"
   8
   9 // Dtrexc reorders the real Schur factorization of a n×n real matrix
  10 //  A = Q*T*Q^T
  11 // so that the diagonal block of T with row index ifst is moved to row ilst.
  12 //
  13 // On entry, T must be in Schur canonical form, that is, block upper triangular
  14 // with 1×1 and 2×2 diagonal blocks; each 2×2 diagonal block has its diagonal
  15 // elements equal and its off-diagonal elements of opposite sign.
  16 //
  17 // On return, T will be reordered by an orthogonal similarity transformation Z
  18 // as Z^T*T*Z, and will be again in Schur canonical form.
  19 //
  20 // If compq is lapack.UpdateSchur, on return the matrix Q of Schur vectors will be
  21 // updated by postmultiplying it with Z.
  22 // If compq is lapack.None, the matrix Q is not referenced and will not be
  23 // updated.
  24 // For other values of compq Dtrexc will panic.
  25 //
  26 // ifst and ilst specify the reordering of the diagonal blocks of T. The block
  27 // with row index ifst is moved to row ilst, by a sequence of transpositions
  28 // between adjacent blocks.
  29 //
  30 // If ifst points to the second row of a 2×2 block, ifstOut will point to the
  31 // first row, otherwise it will be equal to ifst.
  32 //
  33 // ilstOut will point to the first row of the block in its final position. If ok
  34 // is true, ilstOut may differ from ilst by +1 or -1.
  35 //
  36 // It must hold that
  37 //  0 <= ifst < n, and  0 <= ilst < n,
  38 // otherwise Dtrexc will panic.
  39 //
  40 // If ok is false, two adjacent blocks were too close to swap because the
  41 // problem is very ill-conditioned. T may have been partially reordered, and
  42 // ilstOut will point to the first row of the block at the position to which it
  43 // has been moved.
  44 //
  45 // work must have length at least n, otherwise Dtrexc will panic.
  46 //
  47 // Dtrexc is an internal routine. It is exported for testing purposes.
  48 func (impl Implementation) Dtrexc(compq lapack.EVComp, n int, t []float64, ldt int, q []float64, ldq int, ifst, ilst int, work []float64) (ifstOut, ilstOut int, ok bool) {
  49         checkMatrix(n, n, t, ldt)
  50         var wantq bool
  51         switch compq {
  52         default:
  53                 panic("lapack: bad value of compq")
  54         case lapack.None:
  55                 // Nothing to do because wantq is already false.
  56         case lapack.UpdateSchur:
  57                 wantq = true
  58                 checkMatrix(n, n, q, ldq)
  59         }
  60         if (ifst < 0 || n <= ifst) && n > 0 {
  61                 panic("lapack: ifst out of range")
  62         }
  63         if (ilst < 0 || n <= ilst) && n > 0 {
  64                 panic("lapack: ilst out of range")
  65         }
  66         if len(work) < n {
  67                 panic(badWork)
  68         }
  69
  70         ok = true
  71
  72         // Quick return if possible.
  73         if n <= 1 {
  74                 return ifst, ilst, true
  75         }
  76
  77         // Determine the first row of specified block
  78         // and find out it is 1×1 or 2×2.
  79         if ifst > 0 && t[ifst*ldt+ifst-1] != 0 {
  80                 ifst--
  81         }
  82         nbf := 1 // Size of the first block.
  83         if ifst+1 < n && t[(ifst+1)*ldt+ifst] != 0 {
  84                 nbf = 2
  85         }
  86         // Determine the first row of the final block
  87         // and find out it is 1×1 or 2×2.
  88         if ilst > 0 && t[ilst*ldt+ilst-1] != 0 {
  89                 ilst--
  90         }
  91         nbl := 1 // Size of the last block.
  92         if ilst+1 < n && t[(ilst+1)*ldt+ilst] != 0 {
  93                 nbl = 2
  94         }
  95
  96         switch {
  97         case ifst == ilst:
  98                 return ifst, ilst, true
  99
 100         case ifst < ilst:
 101                 // Update ilst.
 102                 switch {
 103                 case nbf == 2 && nbl == 1:
 104                         ilst--
 105                 case nbf == 1 && nbl == 2:
 106                         ilst++
 107                 }
 108                 here := ifst
 109                 for here < ilst {
 110                         // Swap block with next one below.
 111                         if nbf == 1 || nbf == 2 {
 112                                 // Current block either 1×1 or 2×2.
 113                                 nbnext := 1 // Size of the next block.
 114                                 if here+nbf+1 < n && t[(here+nbf+1)*ldt+here+nbf] != 0 {
 115                                         nbnext = 2
 116                                 }
 117                                 ok = impl.Dlaexc(wantq, n, t, ldt, q, ldq, here, nbf, nbnext, work)
 118                                 if !ok {
 119                                         return ifst, here, false
 120                                 }
 121                                 here += nbnext
 122                                 // Test if 2×2 block breaks into two 1×1 blocks.
 123                                 if nbf == 2 && t[(here+1)*ldt+here] == 0 {
 124                                         nbf = 3
 125                                 }
 126                                 continue
 127                         }
 128
 129                         // Current block consists of two 1×1 blocks each of
 130                         // which must be swapped individually.
 131                         nbnext := 1 // Size of the next block.
 132                         if here+3 < n && t[(here+3)*ldt+here+2] != 0 {
 133                                 nbnext = 2
 134                         }
 135                         ok = impl.Dlaexc(wantq, n, t, ldt, q, ldq, here+1, 1, nbnext, work)
 136                         if !ok {
 137                                 return ifst, here, false
 138                         }
 139                         if nbnext == 1 {
 140                                 // Swap two 1×1 blocks, no problems possible.
 141                                 impl.Dlaexc(wantq, n, t, ldt, q, ldq, here, 1, nbnext, work)
 142                                 here++
 143                                 continue
 144                         }
 145                         // Recompute nbnext in case 2×2 split.
 146                         if t[(here+2)*ldt+here+1] == 0 {
 147                                 nbnext = 1
 148                         }
 149                         if nbnext == 2 {
 150                                 // 2×2 block did not split.
 151                                 ok = impl.Dlaexc(wantq, n, t, ldt, q, ldq, here, 1, nbnext, work)
 152                                 if !ok {
 153                                         return ifst, here, false
 154                                 }
 155                         } else {
 156                                 // 2×2 block did split.
 157                                 impl.Dlaexc(wantq, n, t, ldt, q, ldq, here, 1, 1, work)
 158                                 impl.Dlaexc(wantq, n, t, ldt, q, ldq, here+1, 1, 1, work)
 159                         }
 160                         here += 2
 161                 }
 162                 return ifst, here, true
 163
 164         default: // ifst > ilst
 165                 here := ifst
 166                 for here > ilst {
 167                         // Swap block with next one above.
 168                         if nbf == 1 || nbf == 2 {
 169                                 // Current block either 1×1 or 2×2.
 170                                 nbnext := 1
 171                                 if here-2 >= 0 && t[(here-1)*ldt+here-2] != 0 {
 172                                         nbnext = 2
 173                                 }
 174                                 ok = impl.Dlaexc(wantq, n, t, ldt, q, ldq, here-nbnext, nbnext, nbf, work)
 175                                 if !ok {
 176                                         return ifst, here, false
 177                                 }
 178                                 here -= nbnext
 179                                 // Test if 2×2 block breaks into two 1×1 blocks.
 180                                 if nbf == 2 && t[(here+1)*ldt+here] == 0 {
 181                                         nbf = 3
 182                                 }
 183                                 continue
 184                         }
 185
 186                         // Current block consists of two 1×1 blocks each of
 187                         // which must be swapped individually.
 188                         nbnext := 1
 189                         if here-2 >= 0 && t[(here-1)*ldt+here-2] != 0 {
 190                                 nbnext = 2
 191                         }
 192                         ok = impl.Dlaexc(wantq, n, t, ldt, q, ldq, here-nbnext, nbnext, 1, work)
 193                         if !ok {
 194                                 return ifst, here, false
 195                         }
 196                         if nbnext == 1 {
 197                                 // Swap two 1×1 blocks, no problems possible.
 198                                 impl.Dlaexc(wantq, n, t, ldt, q, ldq, here, nbnext, 1, work)
 199                                 here--
 200                                 continue
 201                         }
 202                         // Recompute nbnext in case 2×2 split.
 203                         if t[here*ldt+here-1] == 0 {
 204                                 nbnext = 1
 205                         }
 206                         if nbnext == 2 {
 207                                 // 2×2 block did not split.
 208                                 ok = impl.Dlaexc(wantq, n, t, ldt, q, ldq, here-1, 2, 1, work)
 209                                 if !ok {
 210                                         return ifst, here, false
 211                                 }
 212                         } else {
 213                                 // 2×2 block did split.
 214                                 impl.Dlaexc(wantq, n, t, ldt, q, ldq, here, 1, 1, work)
 215                                 impl.Dlaexc(wantq, n, t, ldt, q, ldq, here-1, 1, 1, work)
 216                         }
 217                         here -= 2
 218                 }
 219                 return ifst, here, true
 220         }
 221 }