vendor/gonum.org/v1/gonum/lapack/gonum/dggsvp3.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"
  11         "gonum.org/v1/gonum/lapack"
  12 )
  13
  14 // Dggsvp3 computes orthogonal matrices U, V and Q such that
  15 //
  16 //                  n-k-l  k    l
  17 //  U^T*A*Q =    k [ 0    A12  A13 ] if m-k-l >= 0;
  18 //               l [ 0     0   A23 ]
  19 //           m-k-l [ 0     0    0  ]
  20 //
  21 //                  n-k-l  k    l
  22 //  U^T*A*Q =    k [ 0    A12  A13 ] if m-k-l < 0;
  23 //             m-k [ 0     0   A23 ]
  24 //
  25 //                  n-k-l  k    l
  26 //  V^T*B*Q =    l [ 0     0   B13 ]
  27 //             p-l [ 0     0    0  ]
  28 //
  29 // where the k×k matrix A12 and l×l matrix B13 are non-singular
  30 // upper triangular. A23 is l×l upper triangular if m-k-l >= 0,
  31 // otherwise A23 is (m-k)×l upper trapezoidal.
  32 //
  33 // Dggsvp3 returns k and l, the dimensions of the sub-blocks. k+l
  34 // is the effective numerical rank of the (m+p)×n matrix [ A^T B^T ]^T.
  35 //
  36 // jobU, jobV and jobQ are options for computing the orthogonal matrices. The behavior
  37 // is as follows
  38 //  jobU == lapack.GSVDU        Compute orthogonal matrix U
  39 //  jobU == lapack.GSVDNone     Do not compute orthogonal matrix.
  40 // The behavior is the same for jobV and jobQ with the exception that instead of
  41 // lapack.GSVDU these accept lapack.GSVDV and lapack.GSVDQ respectively.
  42 // The matrices U, V and Q must be m×m, p×p and n×n respectively unless the
  43 // relevant job parameter is lapack.GSVDNone.
  44 //
  45 // tola and tolb are the convergence criteria for the Jacobi-Kogbetliantz
  46 // iteration procedure. Generally, they are the same as used in the preprocessing
  47 // step, for example,
  48 //  tola = max(m, n)*norm(A)*eps,
  49 //  tolb = max(p, n)*norm(B)*eps.
  50 // Where eps is the machine epsilon.
  51 //
  52 // iwork must have length n, work must have length at least max(1, lwork), and
  53 // lwork must be -1 or greater than zero, otherwise Dggsvp3 will panic.
  54 //
  55 // Dggsvp3 is an internal routine. It is exported for testing purposes.
  56 func (impl Implementation) Dggsvp3(jobU, jobV, jobQ lapack.GSVDJob, m, p, n int, a []float64, lda int, b []float64, ldb int, tola, tolb float64, u []float64, ldu int, v []float64, ldv int, q []float64, ldq int, iwork []int, tau, work []float64, lwork int) (k, l int) {
  57         const forward = true
  58
  59         checkMatrix(m, n, a, lda)
  60         checkMatrix(p, n, b, ldb)
  61
  62         wantu := jobU == lapack.GSVDU
  63         if !wantu && jobU != lapack.GSVDNone {
  64                 panic(badGSVDJob + "U")
  65         }
  66         if jobU != lapack.GSVDNone {
  67                 checkMatrix(m, m, u, ldu)
  68         }
  69
  70         wantv := jobV == lapack.GSVDV
  71         if !wantv && jobV != lapack.GSVDNone {
  72                 panic(badGSVDJob + "V")
  73         }
  74         if jobV != lapack.GSVDNone {
  75                 checkMatrix(p, p, v, ldv)
  76         }
  77
  78         wantq := jobQ == lapack.GSVDQ
  79         if !wantq && jobQ != lapack.GSVDNone {
  80                 panic(badGSVDJob + "Q")
  81         }
  82         if jobQ != lapack.GSVDNone {
  83                 checkMatrix(n, n, q, ldq)
  84         }
  85
  86         if len(iwork) != n {
  87                 panic(badWork)
  88         }
  89         if lwork != -1 && lwork < 1 {
  90                 panic(badWork)
  91         }
  92         if len(work) < max(1, lwork) {
  93                 panic(badWork)
  94         }
  95
  96         var lwkopt int
  97         impl.Dgeqp3(p, n, b, ldb, iwork, tau, work, -1)
  98         lwkopt = int(work[0])
  99         if wantv {
 100                 lwkopt = max(lwkopt, p)
 101         }
 102         lwkopt = max(lwkopt, min(n, p))
 103         lwkopt = max(lwkopt, m)
 104         if wantq {
 105                 lwkopt = max(lwkopt, n)
 106         }
 107         impl.Dgeqp3(m, n, a, lda, iwork, tau, work, -1)
 108         lwkopt = max(lwkopt, int(work[0]))
 109         lwkopt = max(1, lwkopt)
 110         if lwork == -1 {
 111                 work[0] = float64(lwkopt)
 112                 return 0, 0
 113         }
 114
 115         // tau check must come after lwkopt query since
 116         // the Dggsvd3 call for lwkopt query may have
 117         // lwork == -1, and tau is provided by work.
 118         if len(tau) < n {
 119                 panic(badTau)
 120         }
 121
 122         // QR with column pivoting of B: B*P = V*[ S11 S12 ].
 123         //                                       [  0   0  ]
 124         for i := range iwork[:n] {
 125                 iwork[i] = 0
 126         }
 127         impl.Dgeqp3(p, n, b, ldb, iwork, tau, work, lwork)
 128
 129         // Update A := A*P.
 130         impl.Dlapmt(forward, m, n, a, lda, iwork)
 131
 132         // Determine the effective rank of matrix B.
 133         for i := 0; i < min(p, n); i++ {
 134                 if math.Abs(b[i*ldb+i]) > tolb {
 135                         l++
 136                 }
 137         }
 138
 139         if wantv {
 140                 // Copy the details of V, and form V.
 141                 impl.Dlaset(blas.All, p, p, 0, 0, v, ldv)
 142                 if p > 1 {
 143                         impl.Dlacpy(blas.Lower, p-1, min(p, n), b[ldb:], ldb, v[ldv:], ldv)
 144                 }
 145                 impl.Dorg2r(p, p, min(p, n), v, ldv, tau, work)
 146         }
 147
 148         // Clean up B.
 149         for i := 1; i < l; i++ {
 150                 r := b[i*ldb : i*ldb+i]
 151                 for j := range r {
 152                         r[j] = 0
 153                 }
 154         }
 155         if p > l {
 156                 impl.Dlaset(blas.All, p-l, n, 0, 0, b[l*ldb:], ldb)
 157         }
 158
 159         if wantq {
 160                 // Set Q = I and update Q := Q*P.
 161                 impl.Dlaset(blas.All, n, n, 0, 1, q, ldq)
 162                 impl.Dlapmt(forward, n, n, q, ldq, iwork)
 163         }
 164
 165         if p >= l && n != l {
 166                 // RQ factorization of [ S11 S12 ]: [ S11 S12 ] = [ 0 S12 ]*Z.
 167                 impl.Dgerq2(l, n, b, ldb, tau, work)
 168
 169                 // Update A := A*Z^T.
 170                 impl.Dormr2(blas.Right, blas.Trans, m, n, l, b, ldb, tau, a, lda, work)
 171
 172                 if wantq {
 173                         // Update Q := Q*Z^T.
 174                         impl.Dormr2(blas.Right, blas.Trans, n, n, l, b, ldb, tau, q, ldq, work)
 175                 }
 176
 177                 // Clean up B.
 178                 impl.Dlaset(blas.All, l, n-l, 0, 0, b, ldb)
 179                 for i := 1; i < l; i++ {
 180                         r := b[i*ldb+n-l : i*ldb+i+n-l]
 181                         for j := range r {
 182                                 r[j] = 0
 183                         }
 184                 }
 185         }
 186
 187         // Let              N-L     L
 188         //            A = [ A11    A12 ] M,
 189         //
 190         // then the following does the complete QR decomposition of A11:
 191         //
 192         //          A11 = U*[  0  T12 ]*P1^T.
 193         //                  [  0   0  ]
 194         for i := range iwork[:n-l] {
 195                 iwork[i] = 0
 196         }
 197         impl.Dgeqp3(m, n-l, a, lda, iwork[:n-l], tau, work, lwork)
 198
 199         // Determine the effective rank of A11.
 200         for i := 0; i < min(m, n-l); i++ {
 201                 if math.Abs(a[i*lda+i]) > tola {
 202                         k++
 203                 }
 204         }
 205
 206         // Update A12 := U^T*A12, where A12 = A[0:m, n-l:n].
 207         impl.Dorm2r(blas.Left, blas.Trans, m, l, min(m, n-l), a, lda, tau, a[n-l:], lda, work)
 208
 209         if wantu {
 210                 // Copy the details of U, and form U.
 211                 impl.Dlaset(blas.All, m, m, 0, 0, u, ldu)
 212                 if m > 1 {
 213                         impl.Dlacpy(blas.Lower, m-1, min(m, n-l), a[lda:], lda, u[ldu:], ldu)
 214                 }
 215                 impl.Dorg2r(m, m, min(m, n-l), u, ldu, tau, work)
 216         }
 217
 218         if wantq {
 219                 // Update Q[0:n, 0:n-l] := Q[0:n, 0:n-l]*P1.
 220                 impl.Dlapmt(forward, n, n-l, q, ldq, iwork[:n-l])
 221         }
 222
 223         // Clean up A: set the strictly lower triangular part of
 224         // A[0:k, 0:k] = 0, and A[k:m, 0:n-l] = 0.
 225         for i := 1; i < k; i++ {
 226                 r := a[i*lda : i*lda+i]
 227                 for j := range r {
 228                         r[j] = 0
 229                 }
 230         }
 231         if m > k {
 232                 impl.Dlaset(blas.All, m-k, n-l, 0, 0, a[k*lda:], lda)
 233         }
 234
 235         if n-l > k {
 236                 // RQ factorization of [ T11 T12 ] = [ 0 T12 ]*Z1.
 237                 impl.Dgerq2(k, n-l, a, lda, tau, work)
 238
 239                 if wantq {
 240                         // Update Q[0:n, 0:n-l] := Q[0:n, 0:n-l]*Z1^T.
 241                         impl.Dorm2r(blas.Right, blas.Trans, n, n-l, k, a, lda, tau, q, ldq, work)
 242                 }
 243
 244                 // Clean up A.
 245                 impl.Dlaset(blas.All, k, n-l-k, 0, 0, a, lda)
 246                 for i := 1; i < k; i++ {
 247                         r := a[i*lda+n-k-l : i*lda+i+n-k-l]
 248                         for j := range r {
 249                                 a[j] = 0
 250                         }
 251                 }
 252         }
 253
 254         if m > k {
 255                 // QR factorization of A[k:m, n-l:n].
 256                 impl.Dgeqr2(m-k, l, a[k*lda+n-l:], lda, tau, work)
 257                 if wantu {
 258                         // Update U[:, k:m) := U[:, k:m]*U1.
 259                         impl.Dorm2r(blas.Right, blas.NoTrans, m, m-k, min(m-k, l), a[k*lda+n-l:], lda, tau, u[k:], ldu, work)
 260                 }
 261
 262                 // Clean up A.
 263                 for i := k + 1; i < m; i++ {
 264                         r := a[i*lda+n-l : i*lda+min(n-l+i-k, n)]
 265                         for j := range r {
 266                                 r[j] = 0
 267                         }
 268                 }
 269         }
 270
 271         work[0] = float64(lwkopt)
 272         return k, l
 273 }