aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.hgignore2
-rw-r--r--Eigen/src/Core/util/Memory.h2
-rw-r--r--Eigen/src/SVD/JacobiSVD.h13
-rw-r--r--Eigen/src/SparseCore/CompressedStorage.h2
-rw-r--r--Eigen/src/SparseLU/SparseLU.h5
-rw-r--r--Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h4
-rw-r--r--lapack/complex_double.cpp3
-rw-r--r--lapack/complex_single.cpp3
-rw-r--r--lapack/double.cpp3
-rw-r--r--lapack/eigenvalues.cpp23
-rw-r--r--lapack/lapack_common.h7
-rw-r--r--lapack/single.cpp3
-rw-r--r--lapack/svd.cpp138
-rw-r--r--test/sparse_solver.h33
-rw-r--r--test/sparselu.cpp3
-rw-r--r--unsupported/Eigen/src/BDCSVD/BDCSVD.h30
16 files changed, 214 insertions, 60 deletions
diff --git a/.hgignore b/.hgignore
index e33ba2e9d..769a47f1f 100644
--- a/.hgignore
+++ b/.hgignore
@@ -30,3 +30,5 @@ log
patch
a
a.*
+lapack/testing
+lapack/reference
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index 7b39285ea..b49a18cc7 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2008-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009 Kenneth Riddile <kfriddile@yahoo.com>
// Copyright (C) 2010 Hauke Heibel <hauke.heibel@gmail.com>
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index 59f88a15a..7a8aa8d3f 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -692,21 +692,20 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
// limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
+ // Scaling factor to reduce over/under-flows
+ RealScalar scale = matrix.cwiseAbs().maxCoeff();
+ if(scale==RealScalar(0)) scale = RealScalar(1);
+
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
- if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
+ if(!m_qr_precond_morecols.run(*this, matrix/scale) && !m_qr_precond_morerows.run(*this, matrix/scale))
{
- m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
+ m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
}
-
- // Scaling factor to reduce over/under-flows
- RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff();
- if(scale==RealScalar(0)) scale = RealScalar(1);
- m_workMatrix /= scale;
/*** step 2. The main Jacobi SVD iteration. ***/
diff --git a/Eigen/src/SparseCore/CompressedStorage.h b/Eigen/src/SparseCore/CompressedStorage.h
index 2845e44e7..e09ad97f0 100644
--- a/Eigen/src/SparseCore/CompressedStorage.h
+++ b/Eigen/src/SparseCore/CompressedStorage.h
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index b02796293..d72d7f150 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -249,14 +249,13 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
// Initialize with the determinant of the row matrix
Scalar det = Scalar(1.);
- //Note that the diagonal blocks of U are stored in supernodes,
+ // Note that the diagonal blocks of U are stored in supernodes,
// which are available in the L part :)
for (Index j = 0; j < this->cols(); ++j)
{
for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
{
- if(it.row() < j) continue;
- if(it.row() == j)
+ if(it.index() == j)
{
det *= abs(it.value());
break;
diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
index 6b0b83e0b..e8ee35a94 100644
--- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
+++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
@@ -189,8 +189,8 @@ class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator
m_idval(mat.colIndexPtr()[outer]),
m_startidval(m_idval),
m_endidval(mat.colIndexPtr()[outer+1]),
- m_idrow(mat.rowIndexPtr()[outer]),
- m_endidrow(mat.rowIndexPtr()[outer+1])
+ m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
+ m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
{}
inline InnerIterator& operator++()
{
diff --git a/lapack/complex_double.cpp b/lapack/complex_double.cpp
index 424d2b8ca..c9c575273 100644
--- a/lapack/complex_double.cpp
+++ b/lapack/complex_double.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
+#include "svd.cpp"
diff --git a/lapack/complex_single.cpp b/lapack/complex_single.cpp
index c0b2d29ae..6d11b26cd 100644
--- a/lapack/complex_single.cpp
+++ b/lapack/complex_single.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
+#include "svd.cpp"
diff --git a/lapack/double.cpp b/lapack/double.cpp
index d86549e19..ea78bb662 100644
--- a/lapack/double.cpp
+++ b/lapack/double.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
#include "eigenvalues.cpp"
+#include "svd.cpp"
diff --git a/lapack/eigenvalues.cpp b/lapack/eigenvalues.cpp
index 6141032ab..921c51569 100644
--- a/lapack/eigenvalues.cpp
+++ b/lapack/eigenvalues.cpp
@@ -7,10 +7,10 @@
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-#include "common.h"
+#include "lapack_common.h"
#include <Eigen/Eigenvalues>
-// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges
+// computes eigen values and vectors of a general N-by-N matrix A
EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info))
{
// TODO exploit the work buffer
@@ -22,24 +22,7 @@ EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Sca
else if(*n<0) *info = -3;
else if(*lda<std::max(1,*n)) *info = -5;
else if((!query_size) && *lwork<std::max(1,3**n-1)) *info = -8;
-
-// if(*info==0)
-// {
-// int nb = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 )
-// LWKOPT = MAX( 1, ( NB+2 )*N )
-// WORK( 1 ) = LWKOPT
-// *
-// IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
-// $ INFO = -8
-// END IF
-// *
-// IF( INFO.NE.0 ) THEN
-// CALL XERBLA( 'SSYEV ', -INFO )
-// RETURN
-// ELSE IF( LQUERY ) THEN
-// RETURN
-// END IF
-
+
if(*info!=0)
{
int e = -*info;
diff --git a/lapack/lapack_common.h b/lapack/lapack_common.h
index e558c1409..a93598784 100644
--- a/lapack/lapack_common.h
+++ b/lapack/lapack_common.h
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2010-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -18,6 +18,11 @@
typedef Eigen::Map<Eigen::Transpositions<Eigen::Dynamic,Eigen::Dynamic,int> > PivotsType;
+#if ISCOMPLEX
+#define EIGEN_LAPACK_ARG_IF_COMPLEX(X) X,
+#else
+#define EIGEN_LAPACK_ARG_IF_COMPLEX(X)
+#endif
#endif // EIGEN_LAPACK_COMMON_H
diff --git a/lapack/single.cpp b/lapack/single.cpp
index a64ed44e1..c7da3effa 100644
--- a/lapack/single.cpp
+++ b/lapack/single.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
#include "eigenvalues.cpp"
+#include "svd.cpp"
diff --git a/lapack/svd.cpp b/lapack/svd.cpp
new file mode 100644
index 000000000..ecac3bab1
--- /dev/null
+++ b/lapack/svd.cpp
@@ -0,0 +1,138 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "lapack_common.h"
+#include <Eigen/SVD>
+#include <unsupported/Eigen/BDCSVD>
+
+// computes the singular values/vectors a general M-by-N matrix A using divide-and-conquer
+EIGEN_LAPACK_FUNC(gesdd,(char *jobz, int *m, int* n, Scalar* a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt, Scalar* /*work*/, int* lwork,
+ EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int * /*iwork*/, int *info))
+{
+ // TODO exploit the work buffer
+ bool query_size = *lwork==-1;
+ int diag_size = (std::min)(*m,*n);
+
+ *info = 0;
+ if(*jobz!='A' && *jobz!='S' && *jobz!='O' && *jobz!='N') *info = -1;
+ else if(*m<0) *info = -2;
+ else if(*n<0) *info = -3;
+ else if(*lda<std::max(1,*m)) *info = -5;
+ else if(*lda<std::max(1,*m)) *info = -8;
+ else if(*ldu <1 || (*jobz=='A' && *ldu <*m)
+ || (*jobz=='O' && *m<*n && *ldu<*m)) *info = -8;
+ else if(*ldvt<1 || (*jobz=='A' && *ldvt<*n)
+ || (*jobz=='S' && *ldvt<diag_size)
+ || (*jobz=='O' && *m>=*n && *ldvt<*n)) *info = -10;
+
+ if(*info!=0)
+ {
+ int e = -*info;
+ return xerbla_(SCALAR_SUFFIX_UP"GESDD ", &e, 6);
+ }
+
+ if(query_size)
+ {
+ *lwork = 0;
+ return 0;
+ }
+
+ if(*n==0 || *m==0)
+ return 0;
+
+ PlainMatrixType mat(*m,*n);
+ mat = matrix(a,*m,*n,*lda);
+
+ int option = *jobz=='A' ? ComputeFullU|ComputeFullV
+ : *jobz=='S' ? ComputeThinU|ComputeThinV
+ : *jobz=='O' ? ComputeThinU|ComputeThinV
+ : 0;
+
+ BDCSVD<PlainMatrixType> svd(mat,option);
+
+ make_vector(s,diag_size) = svd.singularValues().head(diag_size);
+
+ if(*jobz=='A')
+ {
+ matrix(u,*m,*m,*ldu) = svd.matrixU();
+ matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
+ }
+ else if(*jobz=='S')
+ {
+ matrix(u,*m,diag_size,*ldu) = svd.matrixU();
+ matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
+ }
+ else if(*jobz=='O' && *m>=*n)
+ {
+ matrix(a,*m,*n,*lda) = svd.matrixU();
+ matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
+ }
+ else if(*jobz=='O')
+ {
+ matrix(u,*m,*m,*ldu) = svd.matrixU();
+ matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
+ }
+
+ return 0;
+}
+
+// computes the singular values/vectors a general M-by-N matrix A using two sided jacobi algorithm
+EIGEN_LAPACK_FUNC(gesvd,(char *jobu, char *jobv, int *m, int* n, Scalar* a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt, Scalar* /*work*/, int* lwork,
+ EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int *info))
+{
+ // TODO exploit the work buffer
+ bool query_size = *lwork==-1;
+ int diag_size = (std::min)(*m,*n);
+
+ *info = 0;
+ if( *jobu!='A' && *jobu!='S' && *jobu!='O' && *jobu!='N') *info = -1;
+ else if((*jobv!='A' && *jobv!='S' && *jobv!='O' && *jobv!='N')
+ || (*jobu=='O' && *jobv=='O')) *info = -2;
+ else if(*m<0) *info = -3;
+ else if(*n<0) *info = -4;
+ else if(*lda<std::max(1,*m)) *info = -6;
+ else if(*ldu <1 || ((*jobu=='A' || *jobu=='S') && *ldu<*m)) *info = -9;
+ else if(*ldvt<1 || (*jobv=='A' && *ldvt<*n)
+ || (*jobv=='S' && *ldvt<diag_size)) *info = -11;
+
+ if(*info!=0)
+ {
+ int e = -*info;
+ return xerbla_(SCALAR_SUFFIX_UP"GESVD ", &e, 6);
+ }
+
+ if(query_size)
+ {
+ *lwork = 0;
+ return 0;
+ }
+
+ if(*n==0 || *m==0)
+ return 0;
+
+ PlainMatrixType mat(*m,*n);
+ mat = matrix(a,*m,*n,*lda);
+
+ int option = (*jobu=='A' ? ComputeFullU : *jobu=='S' || *jobu=='O' ? ComputeThinU : 0)
+ | (*jobv=='A' ? ComputeFullV : *jobv=='S' || *jobv=='O' ? ComputeThinV : 0);
+
+ JacobiSVD<PlainMatrixType> svd(mat,option);
+
+ make_vector(s,diag_size) = svd.singularValues().head(diag_size);
+
+ if(*jobu=='A') matrix(u,*m,*m,*ldu) = svd.matrixU();
+ else if(*jobu=='S') matrix(u,*m,diag_size,*ldu) = svd.matrixU();
+ else if(*jobu=='O') matrix(a,*m,diag_size,*lda) = svd.matrixU();
+
+ if(*jobv=='A') matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
+ else if(*jobv=='S') matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
+ else if(*jobv=='O') matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
+
+ return 0;
+}
diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index 8c8d7f939..b10eaebcc 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -133,7 +133,23 @@ void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType&
Scalar refDet = dA.determinant();
VERIFY_IS_APPROX(refDet,solver.determinant());
}
+template<typename Solver, typename DenseMat>
+void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
+{
+ using std::abs;
+ typedef typename Solver::MatrixType Mat;
+ typedef typename Mat::Scalar Scalar;
+
+ solver.compute(A);
+ if (solver.info() != Success)
+ {
+ std::cerr << "sparse solver testing: factorization failed (check_sparse_abs_determinant)\n";
+ return;
+ }
+ Scalar refDet = abs(dA.determinant());
+ VERIFY_IS_APPROX(refDet,solver.absDeterminant());
+}
template<typename Solver, typename DenseMat>
int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300)
@@ -333,3 +349,20 @@ template<typename Solver> void check_sparse_square_determinant(Solver& solver)
check_sparse_determinant(solver, A, dA);
}
}
+
+template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver)
+{
+ typedef typename Solver::MatrixType Mat;
+ typedef typename Mat::Scalar Scalar;
+ typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+
+ // generate the problem
+ Mat A;
+ DenseMatrix dA;
+ generate_sparse_square_problem(solver, A, dA, 30);
+ A.makeCompressed();
+ for (int i = 0; i < g_repeat; i++) {
+ check_sparse_abs_determinant(solver, A, dA);
+ }
+}
+
diff --git a/test/sparselu.cpp b/test/sparselu.cpp
index 37980defc..52371cb12 100644
--- a/test/sparselu.cpp
+++ b/test/sparselu.cpp
@@ -44,6 +44,9 @@ template<typename T> void test_sparselu_T()
check_sparse_square_solving(sparselu_colamd);
check_sparse_square_solving(sparselu_amd);
check_sparse_square_solving(sparselu_natural);
+
+ check_sparse_square_abs_determinant(sparselu_colamd);
+ check_sparse_square_abs_determinant(sparselu_amd);
}
void test_sparselu()
diff --git a/unsupported/Eigen/src/BDCSVD/BDCSVD.h b/unsupported/Eigen/src/BDCSVD/BDCSVD.h
index e8e795589..e2551cf88 100644
--- a/unsupported/Eigen/src/BDCSVD/BDCSVD.h
+++ b/unsupported/Eigen/src/BDCSVD/BDCSVD.h
@@ -194,7 +194,7 @@ public:
}; //end class BDCSVD
-// Methode to allocate ans initialize matrix and attributes
+// Method to allocate and initialize matrix and attributes
template<typename MatrixType>
void BDCSVD<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
{
@@ -215,24 +215,6 @@ void BDCSVD<MatrixType>::allocate(Index rows, Index cols, unsigned int computati
if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize);
}// end allocate
-// Methode which compute the BDCSVD for the int
-template<>
-BDCSVD<Matrix<int, Dynamic, Dynamic> >& BDCSVD<Matrix<int, Dynamic, Dynamic> >::compute(const MatrixType& matrix, unsigned int computationOptions)
-{
- allocate(matrix.rows(), matrix.cols(), computationOptions);
- m_nonzeroSingularValues = 0;
- m_computed = Matrix<int, Dynamic, Dynamic>::Zero(rows(), cols());
-
- m_singularValues.head(m_diagSize).setZero();
-
- if (m_computeFullU) m_matrixU.setZero(rows(), rows());
- if (m_computeFullV) m_matrixV.setZero(cols(), cols());
- m_isInitialized = true;
- return *this;
-}
-
-
-// Methode which compute the BDCSVD
template<typename MatrixType>
BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions)
{
@@ -612,6 +594,8 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
#endif
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
+ assert(U.allFinite());
+ assert(V.allFinite());
assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n);
assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n);
assert(m_naiveU.allFinite());
@@ -836,8 +820,12 @@ void BDCSVD<MatrixType>::perturbCol0
using std::sqrt;
Index n = col0.size();
Index m = perm.size();
- Index last = perm(m-1);
-
+ if(m==0)
+ {
+ zhat.setZero();
+ return;
+ }
+ Index last = last = perm(m-1);
// The offset permits to skip deflated entries while computing zhat
for (Index k = 0; k < n; ++k)
{
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-27 20:47:38 +0000