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

   1 // Copyright ©2017 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         "gonum.org/v1/gonum/lapack"
  12 )
  13
  14 // Dggsvd3 computes the generalized singular value decomposition (GSVD)
  15 // of an m×n matrix A and p×n matrix B:
  16 //  U^T*A*Q = D1*[ 0 R ]
  17 //
  18 //  V^T*B*Q = D2*[ 0 R ]
  19 // where U, V and Q are orthogonal matrices.
  20 //
  21 // Dggsvd3 returns k and l, the dimensions of the sub-blocks. k+l
  22 // is the effective numerical rank of the (m+p)×n matrix [ A^T B^T ]^T.
  23 // R is a (k+l)×(k+l) nonsingular upper triangular matrix, D1 and
  24 // D2 are m×(k+l) and p×(k+l) diagonal matrices and of the following
  25 // structures, respectively:
  26 //
  27 // If m-k-l >= 0,
  28 //
  29 //                    k  l
  30 //       D1 =     k [ I  0 ]
  31 //                l [ 0  C ]
  32 //            m-k-l [ 0  0 ]
  33 //
  34 //                  k  l
  35 //       D2 = l   [ 0  S ]
  36 //            p-l [ 0  0 ]
  37 //
  38 //               n-k-l  k    l
  39 //  [ 0 R ] = k [  0   R11  R12 ] k
  40 //            l [  0    0   R22 ] l
  41 //
  42 // where
  43 //
  44 //  C = diag( alpha_k, ... , alpha_{k+l} ),
  45 //  S = diag( beta_k,  ... , beta_{k+l} ),
  46 //  C^2 + S^2 = I.
  47 //
  48 // R is stored in
  49 //  A[0:k+l, n-k-l:n]
  50 // on exit.
  51 //
  52 // If m-k-l < 0,
  53 //
  54 //                 k m-k k+l-m
  55 //      D1 =   k [ I  0    0  ]
  56 //           m-k [ 0  C    0  ]
  57 //
  58 //                   k m-k k+l-m
  59 //      D2 =   m-k [ 0  S    0  ]
  60 //           k+l-m [ 0  0    I  ]
  61 //             p-l [ 0  0    0  ]
  62 //
  63 //                 n-k-l  k   m-k  k+l-m
  64 //  [ 0 R ] =    k [ 0    R11  R12  R13 ]
  65 //             m-k [ 0     0   R22  R23 ]
  66 //           k+l-m [ 0     0    0   R33 ]
  67 //
  68 // where
  69 //  C = diag( alpha_k, ... , alpha_m ),
  70 //  S = diag( beta_k,  ... , beta_m ),
  71 //  C^2 + S^2 = I.
  72 //
  73 //  R = [ R11 R12 R13 ] is stored in A[1:m, n-k-l+1:n]
  74 //      [  0  R22 R23 ]
  75 // and R33 is stored in
  76 //  B[m-k:l, n+m-k-l:n] on exit.
  77 //
  78 // Dggsvd3 computes C, S, R, and optionally the orthogonal transformation
  79 // matrices U, V and Q.
  80 //
  81 // jobU, jobV and jobQ are options for computing the orthogonal matrices. The behavior
  82 // is as follows
  83 //  jobU == lapack.GSVDU        Compute orthogonal matrix U
  84 //  jobU == lapack.GSVDNone     Do not compute orthogonal matrix.
  85 // The behavior is the same for jobV and jobQ with the exception that instead of
  86 // lapack.GSVDU these accept lapack.GSVDV and lapack.GSVDQ respectively.
  87 // The matrices U, V and Q must be m×m, p×p and n×n respectively unless the
  88 // relevant job parameter is lapack.GSVDNone.
  89 //
  90 // alpha and beta must have length n or Dggsvd3 will panic. On exit, alpha and
  91 // beta contain the generalized singular value pairs of A and B
  92 //   alpha[0:k] = 1,
  93 //   beta[0:k]  = 0,
  94 // if m-k-l >= 0,
  95 //   alpha[k:k+l] = diag(C),
  96 //   beta[k:k+l]  = diag(S),
  97 // if m-k-l < 0,
  98 //   alpha[k:m]= C, alpha[m:k+l]= 0
  99 //   beta[k:m] = S, beta[m:k+l] = 1.
 100 // if k+l < n,
 101 //   alpha[k+l:n] = 0 and
 102 //   beta[k+l:n]  = 0.
 103 //
 104 // On exit, iwork contains the permutation required to sort alpha descending.
 105 //
 106 // iwork must have length n, work must have length at least max(1, lwork), and
 107 // lwork must be -1 or greater than n, otherwise Dggsvd3 will panic. If
 108 // lwork is -1, work[0] holds the optimal lwork on return, but Dggsvd3 does
 109 // not perform the GSVD.
 110 func (impl Implementation) Dggsvd3(jobU, jobV, jobQ lapack.GSVDJob, m, n, p int, a []float64, lda int, b []float64, ldb int, alpha, beta, u []float64, ldu int, v []float64, ldv int, q []float64, ldq int, work []float64, lwork int, iwork []int) (k, l int, ok bool) {
 111         checkMatrix(m, n, a, lda)
 112         checkMatrix(p, n, b, ldb)
 113
 114         wantu := jobU == lapack.GSVDU
 115         if wantu {
 116                 checkMatrix(m, m, u, ldu)
 117         } else if jobU != lapack.GSVDNone {
 118                 panic(badGSVDJob + "U")
 119         }
 120         wantv := jobV == lapack.GSVDV
 121         if wantv {
 122                 checkMatrix(p, p, v, ldv)
 123         } else if jobV != lapack.GSVDNone {
 124                 panic(badGSVDJob + "V")
 125         }
 126         wantq := jobQ == lapack.GSVDQ
 127         if wantq {
 128                 checkMatrix(n, n, q, ldq)
 129         } else if jobQ != lapack.GSVDNone {
 130                 panic(badGSVDJob + "Q")
 131         }
 132
 133         if len(alpha) != n {
 134                 panic(badAlpha)
 135         }
 136         if len(beta) != n {
 137                 panic(badBeta)
 138         }
 139
 140         if lwork != -1 && lwork <= n {
 141                 panic(badWork)
 142         }
 143         if len(work) < max(1, lwork) {
 144                 panic(shortWork)
 145         }
 146         if len(iwork) < n {
 147                 panic(badWork)
 148         }
 149
 150         // Determine optimal work length.
 151         impl.Dggsvp3(jobU, jobV, jobQ,
 152                 m, p, n,
 153                 a, lda,
 154                 b, ldb,
 155                 0, 0,
 156                 u, ldu,
 157                 v, ldv,
 158                 q, ldq,
 159                 iwork,
 160                 work, work, -1)
 161         lwkopt := n + int(work[0])
 162         lwkopt = max(lwkopt, 2*n)
 163         lwkopt = max(lwkopt, 1)
 164         work[0] = float64(lwkopt)
 165         if lwork == -1 {
 166                 return 0, 0, true
 167         }
 168
 169         // Compute the Frobenius norm of matrices A and B.
 170         anorm := impl.Dlange(lapack.NormFrob, m, n, a, lda, nil)
 171         bnorm := impl.Dlange(lapack.NormFrob, p, n, b, ldb, nil)
 172
 173         // Get machine precision and set up threshold for determining
 174         // the effective numerical rank of the matrices A and B.
 175         tola := float64(max(m, n)) * math.Max(anorm, dlamchS) * dlamchP
 176         tolb := float64(max(p, n)) * math.Max(bnorm, dlamchS) * dlamchP
 177
 178         // Preprocessing.
 179         k, l = impl.Dggsvp3(jobU, jobV, jobQ,
 180                 m, p, n,
 181                 a, lda,
 182                 b, ldb,
 183                 tola, tolb,
 184                 u, ldu,
 185                 v, ldv,
 186                 q, ldq,
 187                 iwork,
 188                 work[:n], work[n:], lwork-n)
 189
 190         // Compute the GSVD of two upper "triangular" matrices.
 191         _, ok = impl.Dtgsja(jobU, jobV, jobQ,
 192                 m, p, n,
 193                 k, l,
 194                 a, lda,
 195                 b, ldb,
 196                 tola, tolb,
 197                 alpha, beta,
 198                 u, ldu,
 199                 v, ldv,
 200                 q, ldq,
 201                 work)
 202
 203         // Sort the singular values and store the pivot indices in iwork
 204         // Copy alpha to work, then sort alpha in work.
 205         bi := blas64.Implementation()
 206         bi.Dcopy(n, alpha, 1, work[:n], 1)
 207         ibnd := min(l, m-k)
 208         for i := 0; i < ibnd; i++ {
 209                 // Scan for largest alpha_{k+i}.
 210                 isub := i
 211                 smax := work[k+i]
 212                 for j := i + 1; j < ibnd; j++ {
 213                         if v := work[k+j]; v > smax {
 214                                 isub = j
 215                                 smax = v
 216                         }
 217                 }
 218                 if isub != i {
 219                         work[k+isub] = work[k+i]
 220                         work[k+i] = smax
 221                         iwork[k+i] = k + isub
 222                 } else {
 223                         iwork[k+i] = k + i
 224                 }
 225         }
 226
 227         work[0] = float64(lwkopt)
 228
 229         return k, l, ok
 230 }