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.
10 "gonum.org/v1/gonum/blas"
11 "gonum.org/v1/gonum/blas/blas64"
12 "gonum.org/v1/gonum/lapack"
15 // Dtgsja computes the generalized singular value decomposition (GSVD)
16 // of two real upper triangular or trapezoidal matrices A and B.
18 // A and B have the following forms, which may be obtained by the
19 // preprocessing subroutine Dggsvp from a general m×n matrix A and p×n
23 // A = k [ 0 A12 A13 ] if m-k-l >= 0;
28 // A = k [ 0 A12 A13 ] if m-k-l < 0;
35 // where the k×k matrix A12 and l×l matrix B13 are non-singular
36 // upper triangular. A23 is l×l upper triangular if m-k-l >= 0,
37 // otherwise A23 is (m-k)×l upper trapezoidal.
41 // U^T*A*Q = D1*[ 0 R ], V^T*B*Q = D2*[ 0 R ],
43 // where U, V and Q are orthogonal matrices.
44 // R is a non-singular upper triangular matrix, and D1 and D2 are
45 // diagonal matrices, which are of the following structures:
59 // [ 0 R ] = k [ 0 R11 R12 ] k
64 // C = diag( alpha_k, ... , alpha_{k+l} ),
65 // S = diag( beta_k, ... , beta_{k+l} ),
84 // [ 0 R ] = k [ 0 R11 R12 R13 ]
85 // m-k [ 0 0 R22 R23 ]
86 // k+l-m [ 0 0 0 R33 ]
89 // C = diag( alpha_k, ... , alpha_m ),
90 // S = diag( beta_k, ... , beta_m ),
93 // R = [ R11 R12 R13 ] is stored in A[0:m, n-k-l:n]
95 // and R33 is stored in
96 // B[m-k:l, n+m-k-l:n] on exit.
98 // The computation of the orthogonal transformation matrices U, V or Q
99 // is optional. These matrices may either be formed explicitly, or they
100 // may be post-multiplied into input matrices U1, V1, or Q1.
102 // Dtgsja essentially uses a variant of Kogbetliantz algorithm to reduce
103 // min(l,m-k)×l triangular or trapezoidal matrix A23 and l×l
104 // matrix B13 to the form:
106 // U1^T*A13*Q1 = C1*R1; V1^T*B13*Q1 = S1*R1,
108 // where U1, V1 and Q1 are orthogonal matrices. C1 and S1 are diagonal
109 // matrices satisfying
113 // and R1 is an l×l non-singular upper triangular matrix.
115 // jobU, jobV and jobQ are options for computing the orthogonal matrices. The behavior
117 // jobU == lapack.GSVDU Compute orthogonal matrix U
118 // jobU == lapack.GSVDUnit Use unit-initialized matrix
119 // jobU == lapack.GSVDNone Do not compute orthogonal matrix.
120 // The behavior is the same for jobV and jobQ with the exception that instead of
121 // lapack.GSVDU these accept lapack.GSVDV and lapack.GSVDQ respectively.
122 // The matrices U, V and Q must be m×m, p×p and n×n respectively unless the
123 // relevant job parameter is lapack.GSVDNone.
125 // k and l specify the sub-blocks in the input matrices A and B:
126 // A23 = A[k:min(k+l,m), n-l:n) and B13 = B[0:l, n-l:n]
127 // of A and B, whose GSVD is going to be computed by Dtgsja.
129 // tola and tolb are the convergence criteria for the Jacobi-Kogbetliantz
130 // iteration procedure. Generally, they are the same as used in the preprocessing
131 // step, for example,
132 // tola = max(m, n)*norm(A)*eps,
133 // tolb = max(p, n)*norm(B)*eps,
134 // where eps is the machine epsilon.
136 // work must have length at least 2*n, otherwise Dtgsja will panic.
138 // alpha and beta must have length n or Dtgsja will panic. On exit, alpha and
139 // beta contain the generalized singular value pairs of A and B
143 // alpha[k:k+l] = diag(C),
144 // beta[k:k+l] = diag(S),
146 // alpha[k:m]= C, alpha[m:k+l]= 0
147 // beta[k:m] = S, beta[m:k+l] = 1.
149 // alpha[k+l:n] = 0 and
152 // On exit, A[n-k:n, 0:min(k+l,m)] contains the triangular matrix R or part of R
153 // and if necessary, B[m-k:l, n+m-k-l:n] contains a part of R.
155 // Dtgsja returns whether the routine converged and the number of iteration cycles
158 // Dtgsja is an internal routine. It is exported for testing purposes.
159 func (impl Implementation) Dtgsja(jobU, jobV, jobQ lapack.GSVDJob, m, p, n, k, l int, a []float64, lda int, b []float64, ldb int, tola, tolb float64, alpha, beta, u []float64, ldu int, v []float64, ldv int, q []float64, ldq int, work []float64) (cycles int, ok bool) {
162 checkMatrix(m, n, a, lda)
163 checkMatrix(p, n, b, ldb)
172 initu := jobU == lapack.GSVDUnit
173 wantu := initu || jobU == lapack.GSVDU
174 if !initu && !wantu && jobU != lapack.GSVDNone {
175 panic(badGSVDJob + "U")
177 if jobU != lapack.GSVDNone {
178 checkMatrix(m, m, u, ldu)
181 initv := jobV == lapack.GSVDUnit
182 wantv := initv || jobV == lapack.GSVDV
183 if !initv && !wantv && jobV != lapack.GSVDNone {
184 panic(badGSVDJob + "V")
186 if jobV != lapack.GSVDNone {
187 checkMatrix(p, p, v, ldv)
190 initq := jobQ == lapack.GSVDUnit
191 wantq := initq || jobQ == lapack.GSVDQ
192 if !initq && !wantq && jobQ != lapack.GSVDNone {
193 panic(badGSVDJob + "Q")
195 if jobQ != lapack.GSVDNone {
196 checkMatrix(n, n, q, ldq)
203 // Initialize U, V and Q, if necessary
205 impl.Dlaset(blas.All, m, m, 0, 1, u, ldu)
208 impl.Dlaset(blas.All, p, p, 0, 1, v, ldv)
211 impl.Dlaset(blas.All, n, n, 0, 1, q, ldq)
214 bi := blas64.Implementation()
215 minTol := math.Min(tola, tolb)
217 // Loop until convergence.
219 for cycles = 1; cycles <= maxit; cycles++ {
222 for i := 0; i < l-1; i++ {
223 for j := i + 1; j < l; j++ {
224 var a1, a2, a3 float64
226 a1 = a[(k+i)*lda+n-l+i]
229 a3 = a[(k+j)*lda+n-l+j]
238 a2 = a[(k+i)*lda+n-l+j]
243 a2 = a[(k+j)*lda+n-l+i]
248 csu, snu, csv, snv, csq, snq := impl.Dlags2(upper, a1, a2, a3, b1, b2, b3)
250 // Update (k+i)-th and (k+j)-th rows of matrix A: U^T*A.
252 bi.Drot(l, a[(k+j)*lda+n-l:], 1, a[(k+i)*lda+n-l:], 1, csu, snu)
255 // Update i-th and j-th rows of matrix B: V^T*B.
256 bi.Drot(l, b[j*ldb+n-l:], 1, b[i*ldb+n-l:], 1, csv, snv)
258 // Update (n-l+i)-th and (n-l+j)-th columns of matrices
259 // A and B: A*Q and B*Q.
260 bi.Drot(min(k+l, m), a[n-l+j:], lda, a[n-l+i:], lda, csq, snq)
261 bi.Drot(l, b[n-l+j:], ldb, b[n-l+i:], ldb, csq, snq)
265 a[(k+i)*lda+n-l+j] = 0
270 a[(k+j)*lda+n-l+i] = 0
275 // Update orthogonal matrices U, V, Q, if desired.
276 if wantu && k+j < m {
277 bi.Drot(m, u[k+j:], ldu, u[k+i:], ldu, csu, snu)
280 bi.Drot(p, v[j:], ldv, v[i:], ldv, csv, snv)
283 bi.Drot(n, q[n-l+j:], ldq, q[n-l+i:], ldq, csq, snq)
289 // The matrices A13 and B13 were lower triangular at the start
290 // of the cycle, and are now upper triangular.
292 // Convergence test: test the parallelism of the corresponding
295 for i := 0; i < min(l, m-k); i++ {
296 bi.Dcopy(l-i, a[(k+i)*lda+n-l+i:], 1, work, 1)
297 bi.Dcopy(l-i, b[i*ldb+n-l+i:], 1, work[l:], 1)
298 ssmin := impl.Dlapll(l-i, work, 1, work[l:], 1)
299 error = math.Max(error, ssmin)
301 if math.Abs(error) <= minTol {
302 // The algorithm has converged.
303 // Compute the generalized singular value pairs (alpha, beta)
304 // and set the triangular matrix R to array A.
305 for i := 0; i < k; i++ {
310 for i := 0; i < min(l, m-k); i++ {
311 a1 := a[(k+i)*lda+n-l+i]
317 // Change sign if necessary.
319 bi.Dscal(l-i, -1, b[i*ldb+n-l+i:], 1)
321 bi.Dscal(p, -1, v[i:], ldv)
324 beta[k+i], alpha[k+i], _ = impl.Dlartg(math.Abs(gamma), 1)
326 if alpha[k+i] >= beta[k+i] {
327 bi.Dscal(l-i, 1/alpha[k+i], a[(k+i)*lda+n-l+i:], 1)
329 bi.Dscal(l-i, 1/beta[k+i], b[i*ldb+n-l+i:], 1)
330 bi.Dcopy(l-i, b[i*ldb+n-l+i:], 1, a[(k+i)*lda+n-l+i:], 1)
335 bi.Dcopy(l-i, b[i*ldb+n-l+i:], 1, a[(k+i)*lda+n-l+i:], 1)
339 for i := m; i < k+l; i++ {
344 for i := k + l; i < n; i++ {
355 // The algorithm has not converged after maxit cycles.