1 // Copyright ©2015 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.
8 "gonum.org/v1/gonum/blas"
9 "gonum.org/v1/gonum/blas/blas64"
10 "gonum.org/v1/gonum/lapack"
11 "gonum.org/v1/gonum/lapack/gonum"
14 var lapack64 lapack.Float64 = gonum.Implementation{}
16 // Use sets the LAPACK float64 implementation to be used by subsequent BLAS calls.
17 // The default implementation is native.Implementation.
18 func Use(l lapack.Float64) {
22 // Potrf computes the Cholesky factorization of a.
23 // The factorization has the form
24 // A = U^T * U if a.Uplo == blas.Upper, or
25 // A = L * L^T if a.Uplo == blas.Lower,
26 // where U is an upper triangular matrix and L is lower triangular.
27 // The triangular matrix is returned in t, and the underlying data between
28 // a and t is shared. The returned bool indicates whether a is positive
29 // definite and the factorization could be finished.
30 func Potrf(a blas64.Symmetric) (t blas64.Triangular, ok bool) {
31 ok = lapack64.Dpotrf(a.Uplo, a.N, a.Data, a.Stride)
40 // Gecon estimates the reciprocal of the condition number of the n×n matrix A
41 // given the LU decomposition of the matrix. The condition number computed may
42 // be based on the 1-norm or the ∞-norm.
44 // a contains the result of the LU decomposition of A as computed by Getrf.
46 // anorm is the corresponding 1-norm or ∞-norm of the original matrix A.
48 // work is a temporary data slice of length at least 4*n and Gecon will panic otherwise.
50 // iwork is a temporary data slice of length at least n and Gecon will panic otherwise.
51 func Gecon(norm lapack.MatrixNorm, a blas64.General, anorm float64, work []float64, iwork []int) float64 {
52 return lapack64.Dgecon(norm, a.Cols, a.Data, a.Stride, anorm, work, iwork)
55 // Gels finds a minimum-norm solution based on the matrices A and B using the
56 // QR or LQ factorization. Gels returns false if the matrix
57 // A is singular, and true if this solution was successfully found.
59 // The minimization problem solved depends on the input parameters.
61 // 1. If m >= n and trans == blas.NoTrans, Gels finds X such that || A*X - B||_2
63 // 2. If m < n and trans == blas.NoTrans, Gels finds the minimum norm solution of
65 // 3. If m >= n and trans == blas.Trans, Gels finds the minimum norm solution of
67 // 4. If m < n and trans == blas.Trans, Gels finds X such that || A*X - B||_2
69 // Note that the least-squares solutions (cases 1 and 3) perform the minimization
70 // per column of B. This is not the same as finding the minimum-norm matrix.
72 // The matrix A is a general matrix of size m×n and is modified during this call.
73 // The input matrix B is of size max(m,n)×nrhs, and serves two purposes. On entry,
74 // the elements of b specify the input matrix B. B has size m×nrhs if
75 // trans == blas.NoTrans, and n×nrhs if trans == blas.Trans. On exit, the
76 // leading submatrix of b contains the solution vectors X. If trans == blas.NoTrans,
77 // this submatrix is of size n×nrhs, and of size m×nrhs otherwise.
79 // Work is temporary storage, and lwork specifies the usable memory length.
80 // At minimum, lwork >= max(m,n) + max(m,n,nrhs), and this function will panic
81 // otherwise. A longer work will enable blocked algorithms to be called.
82 // In the special case that lwork == -1, work[0] will be set to the optimal working
84 func Gels(trans blas.Transpose, a blas64.General, b blas64.General, work []float64, lwork int) bool {
85 return lapack64.Dgels(trans, a.Rows, a.Cols, b.Cols, a.Data, a.Stride, b.Data, b.Stride, work, lwork)
88 // Geqrf computes the QR factorization of the m×n matrix A using a blocked
89 // algorithm. A is modified to contain the information to construct Q and R.
90 // The upper triangle of a contains the matrix R. The lower triangular elements
91 // (not including the diagonal) contain the elementary reflectors. tau is modified
92 // to contain the reflector scales. tau must have length at least min(m,n), and
93 // this function will panic otherwise.
95 // The ith elementary reflector can be explicitly constructed by first extracting
99 // v[j] = a[j*lda+i] j > i
100 // and computing H_i = I - tau[i] * v * v^T.
102 // The orthonormal matrix Q can be constucted from a product of these elementary
103 // reflectors, Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n).
105 // Work is temporary storage, and lwork specifies the usable memory length.
106 // At minimum, lwork >= m and this function will panic otherwise.
107 // Geqrf is a blocked QR factorization, but the block size is limited
108 // by the temporary space available. If lwork == -1, instead of performing Geqrf,
109 // the optimal work length will be stored into work[0].
110 func Geqrf(a blas64.General, tau, work []float64, lwork int) {
111 lapack64.Dgeqrf(a.Rows, a.Cols, a.Data, a.Stride, tau, work, lwork)
114 // Gelqf computes the LQ factorization of the m×n matrix A using a blocked
115 // algorithm. A is modified to contain the information to construct L and Q. The
116 // lower triangle of a contains the matrix L. The elements above the diagonal
117 // and the slice tau represent the matrix Q. tau is modified to contain the
118 // reflector scales. tau must have length at least min(m,n), and this function
119 // will panic otherwise.
121 // See Geqrf for a description of the elementary reflectors and orthonormal
122 // matrix Q. Q is constructed as a product of these elementary reflectors,
123 // Q = H_{k-1} * ... * H_1 * H_0.
125 // Work is temporary storage, and lwork specifies the usable memory length.
126 // At minimum, lwork >= m and this function will panic otherwise.
127 // Gelqf is a blocked LQ factorization, but the block size is limited
128 // by the temporary space available. If lwork == -1, instead of performing Gelqf,
129 // the optimal work length will be stored into work[0].
130 func Gelqf(a blas64.General, tau, work []float64, lwork int) {
131 lapack64.Dgelqf(a.Rows, a.Cols, a.Data, a.Stride, tau, work, lwork)
134 // Gesvd computes the singular value decomposition of the input matrix A.
136 // The singular value decomposition is
137 // A = U * Sigma * V^T
138 // where Sigma is an m×n diagonal matrix containing the singular values of A,
139 // U is an m×m orthogonal matrix and V is an n×n orthogonal matrix. The first
140 // min(m,n) columns of U and V are the left and right singular vectors of A
143 // jobU and jobVT are options for computing the singular vectors. The behavior
145 // jobU == lapack.SVDAll All m columns of U are returned in u
146 // jobU == lapack.SVDInPlace The first min(m,n) columns are returned in u
147 // jobU == lapack.SVDOverwrite The first min(m,n) columns of U are written into a
148 // jobU == lapack.SVDNone The columns of U are not computed.
149 // The behavior is the same for jobVT and the rows of V^T. At most one of jobU
150 // and jobVT can equal lapack.SVDOverwrite, and Gesvd will panic otherwise.
152 // On entry, a contains the data for the m×n matrix A. During the call to Gesvd
153 // the data is overwritten. On exit, A contains the appropriate singular vectors
154 // if either job is lapack.SVDOverwrite.
156 // s is a slice of length at least min(m,n) and on exit contains the singular
157 // values in decreasing order.
159 // u contains the left singular vectors on exit, stored columnwise. If
160 // jobU == lapack.SVDAll, u is of size m×m. If jobU == lapack.SVDInPlace u is
161 // of size m×min(m,n). If jobU == lapack.SVDOverwrite or lapack.SVDNone, u is
164 // vt contains the left singular vectors on exit, stored rowwise. If
165 // jobV == lapack.SVDAll, vt is of size n×m. If jobVT == lapack.SVDInPlace vt is
166 // of size min(m,n)×n. If jobVT == lapack.SVDOverwrite or lapack.SVDNone, vt is
169 // work is a slice for storing temporary memory, and lwork is the usable size of
170 // the slice. lwork must be at least max(5*min(m,n), 3*min(m,n)+max(m,n)).
171 // If lwork == -1, instead of performing Gesvd, the optimal work length will be
172 // stored into work[0]. Gesvd will panic if the working memory has insufficient
175 // Gesvd returns whether the decomposition successfully completed.
176 func Gesvd(jobU, jobVT lapack.SVDJob, a, u, vt blas64.General, s, work []float64, lwork int) (ok bool) {
177 return lapack64.Dgesvd(jobU, jobVT, a.Rows, a.Cols, a.Data, a.Stride, s, u.Data, u.Stride, vt.Data, vt.Stride, work, lwork)
180 // Getrf computes the LU decomposition of the m×n matrix A.
181 // The LU decomposition is a factorization of A into
183 // where P is a permutation matrix, L is a unit lower triangular matrix, and
184 // U is a (usually) non-unit upper triangular matrix. On exit, L and U are stored
187 // ipiv is a permutation vector. It indicates that row i of the matrix was
188 // changed with ipiv[i]. ipiv must have length at least min(m,n), and will panic
189 // otherwise. ipiv is zero-indexed.
191 // Getrf is the blocked version of the algorithm.
193 // Getrf returns whether the matrix A is singular. The LU decomposition will
194 // be computed regardless of the singularity of A, but division by zero
195 // will occur if the false is returned and the result is used to solve a
196 // system of equations.
197 func Getrf(a blas64.General, ipiv []int) bool {
198 return lapack64.Dgetrf(a.Rows, a.Cols, a.Data, a.Stride, ipiv)
201 // Getri computes the inverse of the matrix A using the LU factorization computed
202 // by Getrf. On entry, a contains the PLU decomposition of A as computed by
203 // Getrf and on exit contains the reciprocal of the original matrix.
205 // Getri will not perform the inversion if the matrix is singular, and returns
206 // a boolean indicating whether the inversion was successful.
208 // Work is temporary storage, and lwork specifies the usable memory length.
209 // At minimum, lwork >= n and this function will panic otherwise.
210 // Getri is a blocked inversion, but the block size is limited
211 // by the temporary space available. If lwork == -1, instead of performing Getri,
212 // the optimal work length will be stored into work[0].
213 func Getri(a blas64.General, ipiv []int, work []float64, lwork int) (ok bool) {
214 return lapack64.Dgetri(a.Cols, a.Data, a.Stride, ipiv, work, lwork)
217 // Getrs solves a system of equations using an LU factorization.
218 // The system of equations solved is
219 // A * X = B if trans == blas.Trans
220 // A^T * X = B if trans == blas.NoTrans
221 // A is a general n×n matrix with stride lda. B is a general matrix of size n×nrhs.
223 // On entry b contains the elements of the matrix B. On exit, b contains the
224 // elements of X, the solution to the system of equations.
226 // a and ipiv contain the LU factorization of A and the permutation indices as
227 // computed by Getrf. ipiv is zero-indexed.
228 func Getrs(trans blas.Transpose, a blas64.General, b blas64.General, ipiv []int) {
229 lapack64.Dgetrs(trans, a.Cols, b.Cols, a.Data, a.Stride, ipiv, b.Data, b.Stride)
232 // Ggsvd3 computes the generalized singular value decomposition (GSVD)
233 // of an m×n matrix A and p×n matrix B:
234 // U^T*A*Q = D1*[ 0 R ]
236 // V^T*B*Q = D2*[ 0 R ]
237 // where U, V and Q are orthogonal matrices.
239 // Ggsvd3 returns k and l, the dimensions of the sub-blocks. k+l
240 // is the effective numerical rank of the (m+p)×n matrix [ A^T B^T ]^T.
241 // R is a (k+l)×(k+l) nonsingular upper triangular matrix, D1 and
242 // D2 are m×(k+l) and p×(k+l) diagonal matrices and of the following
243 // structures, respectively:
257 // [ 0 R ] = k [ 0 R11 R12 ] k
262 // C = diag( alpha_k, ... , alpha_{k+l} ),
263 // S = diag( beta_k, ... , beta_{k+l} ),
277 // D2 = m-k [ 0 S 0 ]
282 // [ 0 R ] = k [ 0 R11 R12 R13 ]
283 // m-k [ 0 0 R22 R23 ]
284 // k+l-m [ 0 0 0 R33 ]
287 // C = diag( alpha_k, ... , alpha_m ),
288 // S = diag( beta_k, ... , beta_m ),
291 // R = [ R11 R12 R13 ] is stored in A[1:m, n-k-l+1:n]
293 // and R33 is stored in
294 // B[m-k:l, n+m-k-l:n] on exit.
296 // Ggsvd3 computes C, S, R, and optionally the orthogonal transformation
297 // matrices U, V and Q.
299 // jobU, jobV and jobQ are options for computing the orthogonal matrices. The behavior
301 // jobU == lapack.GSVDU Compute orthogonal matrix U
302 // jobU == lapack.GSVDNone Do not compute orthogonal matrix.
303 // The behavior is the same for jobV and jobQ with the exception that instead of
304 // lapack.GSVDU these accept lapack.GSVDV and lapack.GSVDQ respectively.
305 // The matrices U, V and Q must be m×m, p×p and n×n respectively unless the
306 // relevant job parameter is lapack.GSVDNone.
308 // alpha and beta must have length n or Ggsvd3 will panic. On exit, alpha and
309 // beta contain the generalized singular value pairs of A and B
313 // alpha[k:k+l] = diag(C),
314 // beta[k:k+l] = diag(S),
316 // alpha[k:m]= C, alpha[m:k+l]= 0
317 // beta[k:m] = S, beta[m:k+l] = 1.
319 // alpha[k+l:n] = 0 and
322 // On exit, iwork contains the permutation required to sort alpha descending.
324 // iwork must have length n, work must have length at least max(1, lwork), and
325 // lwork must be -1 or greater than n, otherwise Ggsvd3 will panic. If
326 // lwork is -1, work[0] holds the optimal lwork on return, but Ggsvd3 does
327 // not perform the GSVD.
328 func Ggsvd3(jobU, jobV, jobQ lapack.GSVDJob, a, b blas64.General, alpha, beta []float64, u, v, q blas64.General, work []float64, lwork int, iwork []int) (k, l int, ok bool) {
329 return lapack64.Dggsvd3(jobU, jobV, jobQ, a.Rows, a.Cols, b.Rows, a.Data, a.Stride, b.Data, b.Stride, alpha, beta, u.Data, u.Stride, v.Data, v.Stride, q.Data, q.Stride, work, lwork, iwork)
332 // Lange computes the matrix norm of the general m×n matrix A. The input norm
333 // specifies the norm computed.
334 // lapack.MaxAbs: the maximum absolute value of an element.
335 // lapack.MaxColumnSum: the maximum column sum of the absolute values of the entries.
336 // lapack.MaxRowSum: the maximum row sum of the absolute values of the entries.
337 // lapack.Frobenius: the square root of the sum of the squares of the entries.
338 // If norm == lapack.MaxColumnSum, work must be of length n, and this function will panic otherwise.
339 // There are no restrictions on work for the other matrix norms.
340 func Lange(norm lapack.MatrixNorm, a blas64.General, work []float64) float64 {
341 return lapack64.Dlange(norm, a.Rows, a.Cols, a.Data, a.Stride, work)
344 // Lansy computes the specified norm of an n×n symmetric matrix. If
345 // norm == lapack.MaxColumnSum or norm == lapackMaxRowSum work must have length
346 // at least n and this function will panic otherwise.
347 // There are no restrictions on work for the other matrix norms.
348 func Lansy(norm lapack.MatrixNorm, a blas64.Symmetric, work []float64) float64 {
349 return lapack64.Dlansy(norm, a.Uplo, a.N, a.Data, a.Stride, work)
352 // Lantr computes the specified norm of an m×n trapezoidal matrix A. If
353 // norm == lapack.MaxColumnSum work must have length at least n and this function
354 // will panic otherwise. There are no restrictions on work for the other matrix norms.
355 func Lantr(norm lapack.MatrixNorm, a blas64.Triangular, work []float64) float64 {
356 return lapack64.Dlantr(norm, a.Uplo, a.Diag, a.N, a.N, a.Data, a.Stride, work)
359 // Lapmt rearranges the columns of the m×n matrix X as specified by the
360 // permutation k_0, k_1, ..., k_{n-1} of the integers 0, ..., n-1.
362 // If forward is true a forward permutation is performed:
364 // X[0:m, k[j]] is moved to X[0:m, j] for j = 0, 1, ..., n-1.
366 // otherwise a backward permutation is performed:
368 // X[0:m, j] is moved to X[0:m, k[j]] for j = 0, 1, ..., n-1.
370 // k must have length n, otherwise Lapmt will panic. k is zero-indexed.
371 func Lapmt(forward bool, x blas64.General, k []int) {
372 lapack64.Dlapmt(forward, x.Rows, x.Cols, x.Data, x.Stride, k)
375 // Ormlq multiplies the matrix C by the othogonal matrix Q defined by
376 // A and tau. A and tau are as returned from Gelqf.
377 // C = Q * C if side == blas.Left and trans == blas.NoTrans
378 // C = Q^T * C if side == blas.Left and trans == blas.Trans
379 // C = C * Q if side == blas.Right and trans == blas.NoTrans
380 // C = C * Q^T if side == blas.Right and trans == blas.Trans
381 // If side == blas.Left, A is a matrix of side k×m, and if side == blas.Right
382 // A is of size k×n. This uses a blocked algorithm.
384 // Work is temporary storage, and lwork specifies the usable memory length.
385 // At minimum, lwork >= m if side == blas.Left and lwork >= n if side == blas.Right,
386 // and this function will panic otherwise.
387 // Ormlq uses a block algorithm, but the block size is limited
388 // by the temporary space available. If lwork == -1, instead of performing Ormlq,
389 // the optimal work length will be stored into work[0].
391 // Tau contains the Householder scales and must have length at least k, and
392 // this function will panic otherwise.
393 func Ormlq(side blas.Side, trans blas.Transpose, a blas64.General, tau []float64, c blas64.General, work []float64, lwork int) {
394 lapack64.Dormlq(side, trans, c.Rows, c.Cols, a.Rows, a.Data, a.Stride, tau, c.Data, c.Stride, work, lwork)
397 // Ormqr multiplies an m×n matrix C by an orthogonal matrix Q as
398 // C = Q * C, if side == blas.Left and trans == blas.NoTrans,
399 // C = Q^T * C, if side == blas.Left and trans == blas.Trans,
400 // C = C * Q, if side == blas.Right and trans == blas.NoTrans,
401 // C = C * Q^T, if side == blas.Right and trans == blas.Trans,
402 // where Q is defined as the product of k elementary reflectors
403 // Q = H_0 * H_1 * ... * H_{k-1}.
405 // If side == blas.Left, A is an m×k matrix and 0 <= k <= m.
406 // If side == blas.Right, A is an n×k matrix and 0 <= k <= n.
407 // The ith column of A contains the vector which defines the elementary
408 // reflector H_i and tau[i] contains its scalar factor. tau must have length k
409 // and Ormqr will panic otherwise. Geqrf returns A and tau in the required
412 // work must have length at least max(1,lwork), and lwork must be at least n if
413 // side == blas.Left and at least m if side == blas.Right, otherwise Ormqr will
416 // work is temporary storage, and lwork specifies the usable memory length. At
417 // minimum, lwork >= m if side == blas.Left and lwork >= n if side ==
418 // blas.Right, and this function will panic otherwise. Larger values of lwork
419 // will generally give better performance. On return, work[0] will contain the
420 // optimal value of lwork.
422 // If lwork is -1, instead of performing Ormqr, the optimal workspace size will
423 // be stored into work[0].
424 func Ormqr(side blas.Side, trans blas.Transpose, a blas64.General, tau []float64, c blas64.General, work []float64, lwork int) {
425 lapack64.Dormqr(side, trans, c.Rows, c.Cols, a.Cols, a.Data, a.Stride, tau, c.Data, c.Stride, work, lwork)
428 // Pocon estimates the reciprocal of the condition number of a positive-definite
429 // matrix A given the Cholesky decmposition of A. The condition number computed
430 // is based on the 1-norm and the ∞-norm.
432 // anorm is the 1-norm and the ∞-norm of the original matrix A.
434 // work is a temporary data slice of length at least 3*n and Pocon will panic otherwise.
436 // iwork is a temporary data slice of length at least n and Pocon will panic otherwise.
437 func Pocon(a blas64.Symmetric, anorm float64, work []float64, iwork []int) float64 {
438 return lapack64.Dpocon(a.Uplo, a.N, a.Data, a.Stride, anorm, work, iwork)
441 // Syev computes all eigenvalues and, optionally, the eigenvectors of a real
442 // symmetric matrix A.
444 // w contains the eigenvalues in ascending order upon return. w must have length
445 // at least n, and Syev will panic otherwise.
447 // On entry, a contains the elements of the symmetric matrix A in the triangular
448 // portion specified by uplo. If jobz == lapack.ComputeEV a contains the
449 // orthonormal eigenvectors of A on exit, otherwise on exit the specified
450 // triangular region is overwritten.
452 // Work is temporary storage, and lwork specifies the usable memory length. At minimum,
453 // lwork >= 3*n-1, and Syev will panic otherwise. The amount of blocking is
454 // limited by the usable length. If lwork == -1, instead of computing Syev the
455 // optimal work length is stored into work[0].
456 func Syev(jobz lapack.EVJob, a blas64.Symmetric, w, work []float64, lwork int) (ok bool) {
457 return lapack64.Dsyev(jobz, a.Uplo, a.N, a.Data, a.Stride, w, work, lwork)
460 // Trcon estimates the reciprocal of the condition number of a triangular matrix A.
461 // The condition number computed may be based on the 1-norm or the ∞-norm.
463 // work is a temporary data slice of length at least 3*n and Trcon will panic otherwise.
465 // iwork is a temporary data slice of length at least n and Trcon will panic otherwise.
466 func Trcon(norm lapack.MatrixNorm, a blas64.Triangular, work []float64, iwork []int) float64 {
467 return lapack64.Dtrcon(norm, a.Uplo, a.Diag, a.N, a.Data, a.Stride, work, iwork)
470 // Trtri computes the inverse of a triangular matrix, storing the result in place
473 // Trtri will not perform the inversion if the matrix is singular, and returns
474 // a boolean indicating whether the inversion was successful.
475 func Trtri(a blas64.Triangular) (ok bool) {
476 return lapack64.Dtrtri(a.Uplo, a.Diag, a.N, a.Data, a.Stride)
479 // Trtrs solves a triangular system of the form A * X = B or A^T * X = B. Trtrs
480 // returns whether the solve completed successfully. If A is singular, no solve is performed.
481 func Trtrs(trans blas.Transpose, a blas64.Triangular, b blas64.General) (ok bool) {
482 return lapack64.Dtrtrs(a.Uplo, trans, a.Diag, a.N, b.Cols, a.Data, a.Stride, b.Data, b.Stride)
485 // Geev computes the eigenvalues and, optionally, the left and/or right
486 // eigenvectors for an n×n real nonsymmetric matrix A.
488 // The right eigenvector v_j of A corresponding to an eigenvalue λ_j
491 // and the left eigenvector u_j corresponding to an eigenvalue λ_j is defined by
492 // u_j^H A = λ_j u_j^H,
493 // where u_j^H is the conjugate transpose of u_j.
495 // On return, A will be overwritten and the left and right eigenvectors will be
496 // stored, respectively, in the columns of the n×n matrices VL and VR in the
497 // same order as their eigenvalues. If the j-th eigenvalue is real, then
500 // and if it is not real, then j and j+1 form a complex conjugate pair and the
501 // eigenvectors can be recovered as
502 // u_j = VL[:,j] + i*VL[:,j+1],
503 // u_{j+1} = VL[:,j] - i*VL[:,j+1],
504 // v_j = VR[:,j] + i*VR[:,j+1],
505 // v_{j+1} = VR[:,j] - i*VR[:,j+1],
506 // where i is the imaginary unit. The computed eigenvectors are normalized to
507 // have Euclidean norm equal to 1 and largest component real.
509 // Left eigenvectors will be computed only if jobvl == lapack.ComputeLeftEV,
510 // otherwise jobvl must be lapack.None.
511 // Right eigenvectors will be computed only if jobvr == lapack.ComputeRightEV,
512 // otherwise jobvr must be lapack.None.
513 // For other values of jobvl and jobvr Geev will panic.
515 // On return, wr and wi will contain the real and imaginary parts, respectively,
516 // of the computed eigenvalues. Complex conjugate pairs of eigenvalues appear
517 // consecutively with the eigenvalue having the positive imaginary part first.
518 // wr and wi must have length n, and Geev will panic otherwise.
520 // work must have length at least lwork and lwork must be at least max(1,4*n) if
521 // the left or right eigenvectors are computed, and at least max(1,3*n) if no
522 // eigenvectors are computed. For good performance, lwork must generally be
523 // larger. On return, optimal value of lwork will be stored in work[0].
525 // If lwork == -1, instead of performing Geev, the function only calculates the
526 // optimal vaule of lwork and stores it into work[0].
528 // On return, first will be the index of the first valid eigenvalue.
529 // If first == 0, all eigenvalues and eigenvectors have been computed.
530 // If first is positive, Geev failed to compute all the eigenvalues, no
531 // eigenvectors have been computed and wr[first:] and wi[first:] contain those
532 // eigenvalues which have converged.
533 func Geev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob, a blas64.General, wr, wi []float64, vl, vr blas64.General, work []float64, lwork int) (first int) {
536 panic("lapack64: matrix not square")
538 if jobvl == lapack.ComputeLeftEV && (vl.Rows != n || vl.Cols != n) {
539 panic("lapack64: bad size of VL")
541 if jobvr == lapack.ComputeRightEV && (vr.Rows != n || vr.Cols != n) {
542 panic("lapack64: bad size of VR")
544 return lapack64.Dgeev(jobvl, jobvr, n, a.Data, a.Stride, wr, wi, vl.Data, vl.Stride, vr.Data, vr.Stride, work, lwork)