This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
 Lines 87-102 void MatrixBase::makeHouseholde (-)a/Eigen/src/Householder/Householder.h (+110 lines) 87 ` beta = sqrt(numext::abs2(c0) + tailSqNorm);` 87 ` beta = sqrt(numext::abs2(c0) + tailSqNorm);` 88 ` if (numext::real(c0)>=RealScalar(0))` 88 ` if (numext::real(c0)>=RealScalar(0))` 89 ` beta = -beta;` 89 ` beta = -beta;` 90 ` essential = tail / (c0 - beta);` 90 ` essential = tail / (c0 - beta);` 91 ` tau = conj((beta - c0) / beta);` 91 ` tau = conj((beta - c0) / beta);` 92 ` }` 92 ` }` 93 `}` 93 `}` 94 94 95 `/** Computes the stable elementary reflector H in a stable way such that:` 96 ` * \f\$ H *this = [ beta 0 ... 0]^T \f\$` 97 ` * where the transformation H is:` 98 ` * \f\$ H = I - tau v v^*\f\$` 99 ` * and the vector v is:` 100 ` * \f\$ v^T = [1 essential^T] \f\$` 101 ` *` 102 ` * \note This routine provides, not only a stable elementary reflector,` 103 ` * but a compatible reflector with Lapack routines.` 104 ` * For instance here, beta is ensured to be greater than zero.` 105 ` *` 106 ` * On output:` 107 ` * \param essential the essential part of the vector \c v` 108 ` * \param tau the scaling factor of the Householder transformation` 109 ` * \param beta ( > 0) the result of H * \c *this` 110 ` *` 111 ` * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),` 112 ` * MatrixBase::applyHouseholderOnTheLeft(), MatrixBase::applyHouseholderOnTheRight()` 113 ` */` 114 `template` 115 `void MatrixBase::makeStableHouseholderInPlace(Scalar& tau, RealScalar& beta)` 116 `{` 117 ` VectorBlock::ret> essentialPart(derived(), 1, size()-1);` 118 ` makeStableHouseholder(essentialPart, tau, beta);` 119 `}` 120 121 122 `/** Computes the stable elementary reflector H in a stable way such that:` 123 ` * \f\$ H *this = [ beta 0 ... 0]^T \f\$` 124 ` * where the transformation H is:` 125 ` * \f\$ H = I - tau v v^*\f\$` 126 ` * and the vector v is:` 127 ` * \f\$ v^T = [1 essential^T] \f\$` 128 ` *` 129 ` * \note This routine provides, not only a stable elementary reflector,` 130 ` * but a compatible reflector with Lapack routines.` 131 ` * For instance here, beta is ensured to be greater than zero.` 132 ` *` 133 ` * On output:` 134 ` * \param essential the essential part of the vector \c v` 135 ` * \param tau the scaling factor of the Householder transformation` 136 ` * \param beta ( > 0) the result of H * \c *this` 137 ` *` 138 ` * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),` 139 ` * MatrixBase::applyHouseholderOnTheLeft(), MatrixBase::applyHouseholderOnTheRight()` 140 ` */` 141 `template` 142 `template` 143 `void MatrixBase::makeStableHouseholder(` 144 ` EssentialPart& essential,` 145 ` Scalar& tau,` 146 ` RealScalar& beta) const` 147 `{` 148 ` using std::sqrt;` 149 ` using numext::conj;` 150 151 ` EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)` 152 ` VectorBlock tail(derived(), 1, size()-1);` 153 ` RealScalar tailNorm = size()==1 ? RealScalar(0) : tail.blueNorm();` 154 ` Scalar c0 = coeff(0);` 155 ` if(tailNorm == RealScalar(0) && numext::imag(c0)==RealScalar(0))` 156 ` {` 157 ` beta = numext::real(c0);` 158 ` if(beta >= 0)` 159 ` tau = RealScalar(0);` 160 ` else` 161 ` {` 162 ` tau = 2;` 163 ` essential.setZero();` 164 ` beta = -beta;` 165 ` }` 166 ` }` 167 ` else` 168 ` {` 169 ` // Compute beta = sqrt(c0*c0 + tailNorm * tailNorm) in a stable way` 170 ` RealScalar c0_abs = (std::abs)(c0);` 171 ` RealScalar max_abs = (std::max)(c0_abs, tailNorm);` 172 ` RealScalar min_abs = (std::min)(c0_abs, tailNorm);` 173 ` if(min_abs == 0) beta = max_abs;` 174 ` else beta = max_abs * sqrt(1. + (min_abs / max_abs)*(min_abs / max_abs));` 175 ` if (numext::real(c0) < RealScalar(0)) beta = -beta;` 176 177 ` if(beta < 0)` 178 ` {` 179 ` beta = -beta;` 180 ` tau = numext::conj(beta-c0)/beta;` 181 ` essential = tail / (c0-beta);` 182 ` }` 183 ` else` 184 ` {` 185 ` RealScalar gamma = numext::real(c0) + beta;` 186 ` RealScalar delta = (numext::imag(c0)*(numext::imag(c0)/gamma) + tailNorm * (tailNorm/gamma));` 187 188 ` if(NumTraits::IsComplex)` 189 ` {` 190 `// std::complex delta_out(delta, numext::imag(c0));` 191 ` Scalar delta_out = c0;` 192 ` numext::real_ref(delta_out) = delta;` 193 ` tau = - delta_out/beta;` 194 ` essential = tail/(delta_out);` 195 ` }` 196 ` else` 197 ` {` 198 ` tau = delta / beta;` 199 ` essential = tail / (-delta);` 200 ` }` 201 ` }` 202 203 ` }` 204 `}` 95 `/** Apply the elementary reflector H given by` 205 `/** Apply the elementary reflector H given by` 96 ` * \f\$ H = I - tau v v^*\f\$` 206 ` * \f\$ H = I - tau v v^*\f\$` 97 ` * with` 207 ` * with` 98 ` * \f\$ v^T = [1 essential^T] \f\$` 208 ` * \f\$ v^T = [1 essential^T] \f\$` 99 ` * from the left to a vector or matrix.` 209 ` * from the left to a vector or matrix.` 100 ` *` 210 ` *` 101 ` * On input:` 211 ` * On input:` 102 ` * \param essential the essential part of the vector \c v` 212 ` * \param essential the essential part of the vector \c v`
 Lines 407-518 typename MatrixType::RealScalar ColPivHo (-)a/Eigen/src/QR/ColPivHouseholderQR.h (-86 / +101 lines) 407 `template` 407 `template` 408 `typename MatrixType::RealScalar ColPivHouseholderQR::logAbsDeterminant() const` 408 `typename MatrixType::RealScalar ColPivHouseholderQR::logAbsDeterminant() const` 409 `{` 409 `{` 410 ` eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");` 410 ` eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");` 411 ` eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");` 411 ` eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");` 412 ` return m_qr.diagonal().cwiseAbs().array().log().sum();` 412 ` return m_qr.diagonal().cwiseAbs().array().log().sum();` 413 `}` 413 `}` 414 414 415 `template` 416 `void colpivhouseholder_qr_inplace(MatrixType& qr, HCoeffsType& hCoeffs, RealRowVectorType& colSqNorms, ` 417 ` IntRowVectorType& colsTranspositions, PermutationType& colsPermutation,` 418 ` RowVectorType& temp, typename MatrixType::RealScalar& maxpivot, ` 419 ` typename MatrixType::Index& nonzero_pivots, typename MatrixType::Index& det_pq)` 420 `{` 421 ` typedef typename MatrixType::Index Index;` 422 ` typedef typename MatrixType::Scalar Scalar;` 423 ` typedef typename MatrixType::RealScalar RealScalar;` 424 ` typedef typename PermutationType::Index PermIndexType;` 425 ` using std::abs;` 426 ` Index rows = qr.rows();` 427 ` Index cols = qr.cols();` 428 ` Index size = qr.diagonalSize();` 429 ` ` 430 ` // the column permutation is stored as int indices, so just to be sure:` 431 ` eigen_assert(cols<=NumTraits::highest());` 432 433 ` hCoeffs.resize(size);` 434 435 ` temp.resize(cols); //FIXME Why is this not a local variable` 436 437 ` colsTranspositions.resize(qr.cols());` 438 ` Index number_of_transpositions = 0;` 439 440 ` colSqNorms.resize(cols);` 441 ` for(Index k = 0; k < cols; ++k)` 442 ` colSqNorms.coeffRef(k) = qr.col(k).squaredNorm();` 443 444 ` RealScalar threshold_helper = colSqNorms.maxCoeff() * numext::abs2(NumTraits::epsilon()) / RealScalar(rows);` 445 446 ` nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)` 447 ` maxpivot = RealScalar(0);` 448 449 ` for(Index k = 0; k < size; ++k)` 450 ` {` 451 ` // first, we look up in our table colSqNorms which column has the biggest squared norm` 452 ` Index biggest_col_index;` 453 ` RealScalar biggest_col_sq_norm = colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);` 454 ` biggest_col_index += k;` 455 456 ` // since our table colSqNorms accumulates imprecision at every step, we must now recompute` 457 ` // the actual squared norm of the selected column.` 458 ` // Note that not doing so does result in solve() sometimes returning inf/nan values` 459 ` // when running the unit test with 1000 repetitions.` 460 ` biggest_col_sq_norm = qr.col(biggest_col_index).tail(rows-k).squaredNorm();` 461 462 ` // we store that back into our table: it can't hurt to correct our table.` 463 ` colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;` 464 465 ` // if the current biggest column is smaller than epsilon times the initial biggest column,` 466 ` // terminate to avoid generating nan/inf values.` 467 ` // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)` 468 ` // repetitions of the unit test, with the result of solve() filled with large values of the order` 469 ` // of 1/(size*epsilon).` 470 ` if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))` 471 ` {` 472 ` nonzero_pivots = k;` 473 ` hCoeffs.tail(size-k).setZero();` 474 ` qr.bottomRightCorner(rows-k,cols-k)` 475 ` .template triangularView()` 476 ` .setZero();` 477 ` break;` 478 ` }` 479 480 ` // apply the transposition to the columns` 481 ` colsTranspositions.coeffRef(k) = biggest_col_index;` 482 ` if(k != biggest_col_index) {` 483 ` qr.col(k).swap(qr.col(biggest_col_index));` 484 ` std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_index));` 485 ` ++number_of_transpositions;` 486 ` }` 487 488 ` // generate the householder vector, store it below the diagonal` 489 ` RealScalar beta;` 490 ` qr.col(k).tail(rows-k).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);` 491 `// qr.col(k).tail(rows-k).makeStableHouseholderInPlace(hCoeffs.coeffRef(k), beta);` 492 493 ` // apply the householder transformation to the diagonal coefficient` 494 ` qr.coeffRef(k,k) = beta;` 495 496 ` // remember the maximum absolute value of diagonal coefficients` 497 ` if(abs(beta) > maxpivot) maxpivot = abs(beta);` 498 499 ` // apply the householder transformation` 500 ` qr.bottomRightCorner(rows-k, cols-k-1)` 501 ` .applyHouseholderOnTheLeft(qr.col(k).tail(rows-k-1), hCoeffs.coeffRef(k), &temp.coeffRef(k+1));` 502 503 ` // update our table of squared norms of the columns` 504 ` colSqNorms.tail(cols-k-1) -= qr.row(k).tail(cols-k-1).cwiseAbs2();` 505 ` }` 506 507 ` colsPermutation.setIdentity(PermIndexType(cols));` 508 ` for(PermIndexType k = 0; k < nonzero_pivots; ++k)` 509 ` colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(colsTranspositions.coeff(k)));` 510 ` ` 511 ` det_pq = (number_of_transpositions%2) ? -1 : 1;` 512 `}` 415 `/** Performs the QR factorization of the given matrix \a matrix. The result of` 513 `/** Performs the QR factorization of the given matrix \a matrix. The result of` 416 ` * the factorization is stored into \c *this, and a reference to \c *this` 514 ` * the factorization is stored into \c *this, and a reference to \c *this` 417 ` * is returned.` 515 ` * is returned.` 418 ` *` 516 ` *` 419 ` * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)` 517 ` * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)` 420 ` */` 518 ` */` 421 `template` 519 `template` 422 `ColPivHouseholderQR& ColPivHouseholderQR::compute(const MatrixType& matrix)` 520 `ColPivHouseholderQR& ColPivHouseholderQR::compute(const MatrixType& matrix)` 423 `{` 521 `{` 424 ` using std::abs;` 522 ` m_qr = matrix;` 425 ` Index rows = matrix.rows();` 523 ` colpivhouseholder_qr_inplace(m_qr, m_hCoeffs, m_colSqNorms, m_colsTranspositions, m_colsPermutation, ` 426 ` Index cols = matrix.cols();` 524 ` m_temp, m_maxpivot, m_nonzero_pivots, m_det_pq);` 427 ` Index size = matrix.diagonalSize();` 428 ` ` 525 ` ` 429 ` // the column permutation is stored as int indices, so just to be sure:` 430 ` eigen_assert(cols<=NumTraits::highest());` 431 432 ` m_qr = matrix;` 433 ` m_hCoeffs.resize(size);` 434 435 ` m_temp.resize(cols);` 436 437 ` m_colsTranspositions.resize(matrix.cols());` 438 ` Index number_of_transpositions = 0;` 439 440 ` m_colSqNorms.resize(cols);` 441 ` for(Index k = 0; k < cols; ++k)` 442 ` m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();` 443 444 ` RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits::epsilon()) / RealScalar(rows);` 445 446 ` m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)` 447 ` m_maxpivot = RealScalar(0);` 448 449 ` for(Index k = 0; k < size; ++k)` 450 ` {` 451 ` // first, we look up in our table m_colSqNorms which column has the biggest squared norm` 452 ` Index biggest_col_index;` 453 ` RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);` 454 ` biggest_col_index += k;` 455 456 ` // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute` 457 ` // the actual squared norm of the selected column.` 458 ` // Note that not doing so does result in solve() sometimes returning inf/nan values` 459 ` // when running the unit test with 1000 repetitions.` 460 ` biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();` 461 462 ` // we store that back into our table: it can't hurt to correct our table.` 463 ` m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;` 464 465 ` // if the current biggest column is smaller than epsilon times the initial biggest column,` 466 ` // terminate to avoid generating nan/inf values.` 467 ` // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)` 468 ` // repetitions of the unit test, with the result of solve() filled with large values of the order` 469 ` // of 1/(size*epsilon).` 470 ` if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))` 471 ` {` 472 ` m_nonzero_pivots = k;` 473 ` m_hCoeffs.tail(size-k).setZero();` 474 ` m_qr.bottomRightCorner(rows-k,cols-k)` 475 ` .template triangularView()` 476 ` .setZero();` 477 ` break;` 478 ` }` 479 480 ` // apply the transposition to the columns` 481 ` m_colsTranspositions.coeffRef(k) = biggest_col_index;` 482 ` if(k != biggest_col_index) {` 483 ` m_qr.col(k).swap(m_qr.col(biggest_col_index));` 484 ` std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));` 485 ` ++number_of_transpositions;` 486 ` }` 487 488 ` // generate the householder vector, store it below the diagonal` 489 ` RealScalar beta;` 490 ` m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);` 491 492 ` // apply the householder transformation to the diagonal coefficient` 493 ` m_qr.coeffRef(k,k) = beta;` 494 495 ` // remember the maximum absolute value of diagonal coefficients` 496 ` if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);` 497 498 ` // apply the householder transformation` 499 ` m_qr.bottomRightCorner(rows-k, cols-k-1)` 500 ` .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));` 501 502 ` // update our table of squared norms of the columns` 503 ` m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();` 504 ` }` 505 506 ` m_colsPermutation.setIdentity(PermIndexType(cols));` 507 ` for(PermIndexType k = 0; k < m_nonzero_pivots; ++k)` 508 ` m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));` 509 510 ` m_det_pq = (number_of_transpositions%2) ? -1 : 1;` 511 ` m_isInitialized = true;` 526 ` m_isInitialized = true;` 512 527 513 ` return *this;` 528 ` return *this;` 514 `}` 529 `}` 515 530 516 `namespace internal {` 531 `namespace internal {` 517 532 518 `template` 533 `template`
 Lines 10-17 (-)a/lapack/complex_double.cpp (+1 lines) 10 `#define SCALAR std::complex` 10 `#define SCALAR std::complex` 11 `#define SCALAR_SUFFIX z` 11 `#define SCALAR_SUFFIX z` 12 `#define SCALAR_SUFFIX_UP "Z"` 12 `#define SCALAR_SUFFIX_UP "Z"` 13 `#define REAL_SCALAR_SUFFIX d` 13 `#define REAL_SCALAR_SUFFIX d` 14 `#define ISCOMPLEX 1` 14 `#define ISCOMPLEX 1` 15 15 16 `#include "cholesky.cpp"` 16 `#include "cholesky.cpp"` 17 `#include "lu.cpp"` 17 `#include "lu.cpp"` 18 `#include "qr.cpp"`
 Lines 10-17 (-)a/lapack/complex_single.cpp (+1 lines) 10 `#define SCALAR std::complex` 10 `#define SCALAR std::complex` 11 `#define SCALAR_SUFFIX c` 11 `#define SCALAR_SUFFIX c` 12 `#define SCALAR_SUFFIX_UP "C"` 12 `#define SCALAR_SUFFIX_UP "C"` 13 `#define REAL_SCALAR_SUFFIX s` 13 `#define REAL_SCALAR_SUFFIX s` 14 `#define ISCOMPLEX 1` 14 `#define ISCOMPLEX 1` 15 15 16 `#include "cholesky.cpp"` 16 `#include "cholesky.cpp"` 17 `#include "lu.cpp"` 17 `#include "lu.cpp"` 18 `#include "qr.cpp"`
 Lines 10-17 (-)a/lapack/double.cpp (+1 lines) 10 `#define SCALAR double` 10 `#define SCALAR double` 11 `#define SCALAR_SUFFIX d` 11 `#define SCALAR_SUFFIX d` 12 `#define SCALAR_SUFFIX_UP "D"` 12 `#define SCALAR_SUFFIX_UP "D"` 13 `#define ISCOMPLEX 0` 13 `#define ISCOMPLEX 0` 14 14 15 `#include "cholesky.cpp"` 15 `#include "cholesky.cpp"` 16 `#include "lu.cpp"` 16 `#include "lu.cpp"` 17 `#include "eigenvalues.cpp"` 17 `#include "eigenvalues.cpp"` 18 `#include "qr.cpp"`
 Lines 1-21 Lines 72-79 EIGEN_LAPACK_FUNC(syev,(char *jobz, char (-)a/lapack/eigenvalues.cpp (-1 / +58 lines) 1 `// This file is part of Eigen, a lightweight C++ template library` 1 `// This file is part of Eigen, a lightweight C++ template library` 2 `// for linear algebra.` 2 `// for linear algebra.` 3 `//` 3 `//` 4 `// Copyright (C) 2011 Gael Guennebaud ` 4 `// Copyright (C) 2011 Gael Guennebaud ` 5 `// Copyright (C) 2013 Désiré Nuentsa-Wakam ` 5 `//` 6 `//` 6 `// This Source Code Form is subject to the terms of the Mozilla` 7 `// This Source Code Form is subject to the terms of the Mozilla` 7 `// Public License v. 2.0. If a copy of the MPL was not distributed` 8 `// Public License v. 2.0. If a copy of the MPL was not distributed` 8 `// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.` 9 `// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.` 9 10 10 `#include "common.h"` 11 `#include "common.h"` 11 `#include ` 12 `#include ` 12 13 13 `// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges` 14 `// computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A` 14 `EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info))` 15 `EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info))` 15 `{` 16 `{` 16 ` // TODO exploit the work buffer` 17 ` // TODO exploit the work buffer` 17 ` bool query_size = *lwork==-1;` 18 ` bool query_size = *lwork==-1;` 18 ` ` 19 ` ` 19 ` *info = 0;` 20 ` *info = 0;` 20 ` if(*jobz!='N' && *jobz!='V') *info = -1;` 21 ` if(*jobz!='N' && *jobz!='V') *info = -1;` 21 ` else if(UPLO(*uplo)==INVALID) *info = -2;` 22 ` else if(UPLO(*uplo)==INVALID) *info = -2;` 72 ` }` 73 ` }` 73 ` ` 74 ` ` 74 ` vector(w,*n) = eig.eigenvalues();` 75 ` vector(w,*n) = eig.eigenvalues();` 75 ` if(computeVectors)` 76 ` if(computeVectors)` 76 ` matrix(a,*n,*n,*lda) = eig.eigenvectors();` 77 ` matrix(a,*n,*n,*lda) = eig.eigenvectors();` 77 ` ` 78 ` ` 78 ` return 0;` 79 ` return 0;` 79 `}` 80 `}` 81 82 `// computes all the eigenvalues, and optionally, the eigenvectors of a real generalized symmetric-definite eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x` 83 84 `EIGEN_LAPACK_FUNC(sygv, (int *itype, char *jobz, char *uplo, int *n, Scalar *a, int *lda, Scalar *b, int *ldb, Scalar *w, Scalar *work, int *lwork, int *info))` 85 `{` 86 ` bool query_size = (*lwork == -1);` 87 ` *info = 0;` 88 89 ` if (*itype < 1 || *itype > 3) *info = -1;` 90 ` else if (!(*jobz=='V' || *jobz == 'N')) *info = -2;` 91 ` else if (UPLO(*uplo)==INVALID) *info = -3;` 92 ` else if (*n < 0) *info = -4;` 93 ` else if (*lda == (std::max)(1, *n)) *info = -6;` 94 ` else if (*ldb == (std::max)(1, *n)) *info = -8;` 95 96 ` // TODO Estimate workspace and return if asked by query_size` 97 98 ` if (*info != 0)` 99 ` {` 100 ` int e = -*info;` 101 ` return xerbla_(SCALAR_SUFFIX_UP"SYGV", &e, 5);` 102 ` }` 103 104 ` if(*n == 0) return 0;` 105 106 ` // Call the GeneralizedSelfAdjointEigenSolver class` 107 ` PlainMatrixType matA(*n, *n), matB(*n, *n);` 108 ` int option;` 109 ` if (*jobz=='V' || *jobz=='v') option = ComputeEigenvectors;` 110 ` if (*itype == 1) option |= Ax_lBx;` 111 ` else if (*itype == 2) option |= ABx_lx;` 112 ` else option |= BAx_lx;` 113 114 ` if(UPLO(*uplo)==UP) matA = matrix(a,*n,*n,*lda).adjoint();` 115 ` else matA = matrix(a,*n,*n,*lda);` 116 117 ` if(UPLO(*uplo)==UP) matB = matrix(b,*n,*n,*ldb).adjoint();` 118 ` else matB = matrix(b,*n,*n,*ldb);` 119 120 ` GeneralizedSelfAdjointEigenSolver es(matA, matB, option);` 121 ` if(es.info() == NoConvergence)` 122 ` {` 123 ` vector(w,*n).setZero();` 124 ` if((option&ComputeEigenvectors) == ComputeEigenvectors)` 125 ` matrix(a, *n,*n, *lda).setIdentity();` 126 ` //FIXME Set the correct return code in *info` 127 ` return 0;` 128 ` }` 129 130 ` // Get the eigenvalues` 131 ` vector(w, *n) = es.eigenvalues();` 132 ` //Get the eigenvectors` 133 ` matrix(a, *n, *n, *lda) = es.eigenvectors();` 134 135 ` return 0;` 136 `}`
 Lines 32-47 EIGEN_LAPACK_FUNC(getrf,(int *m, int *n, (-)a/lapack/lu.cpp (+1 lines) 32 ` ::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions));` 32 ` ::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions));` 33 33 34 ` for(int i=0; i=0)` 37 ` if(ret>=0)` 38 ` *info = ret+1;` 38 ` *info = ret+1;` 39 39 40 ` ` 40 ` return 0;` 41 ` return 0;` 41 `}` 42 `}` 42 43 43 `//GETRS solves a system of linear equations` 44 `//GETRS solves a system of linear equations` 44 `// A * X = B or A' * X = B` 45 `// A * X = B or A' * X = B` 45 `// with a general N-by-N matrix A using the LU factorization computed by GETRF` 46 `// with a general N-by-N matrix A using the LU factorization computed by GETRF` 46 `EIGEN_LAPACK_FUNC(getrs,(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info))` 47 `EIGEN_LAPACK_FUNC(getrs,(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info))` 47 `{` 48 `{`
 Line 0 (-)a/lapack/qr.cpp (+227 lines) 1 `// This file is part of Eigen, a lightweight C++ template library` 2 `// for linear algebra.` 3 `//` 4 `// Copyright (C) 2013 Desire Nuentsa ` 5 `//` 6 `// This Source Code Form is subject to the terms of the Mozilla` 7 `// Public License v. 2.0. If a copy of the MPL was not distributed` 8 `// with this file, You can obtain one at http://mozilla.org/MPL/2.0/` 9 10 `#include "lapack_common.h"` 11 `#include ` 12 13 `// computes a QR factorization of a general M-by-N matrix A` 14 `EIGEN_LAPACK_FUNC(geqrf, (int *m, int *n, Scalar *a, int *lda, Scalar *tau, Scalar *work,` 15 ` int *lwork, int *info))` 16 `{` 17 ` bool query_size = (*lwork == -1);` 18 ` *info = 0;` 19 20 ` if(*m < 0) *info = -1;` 21 ` else if(*n < 0) *info = -2;` 22 ` else if(*lda < (std::max)(1` 23 ` , *m)) *info = -4;` 24 ` else if((*lwork < (std::max)(1, *n)) && !query_size) *info = -7;` 25 ` if (*info != 0)` 26 ` {` 27 ` int e = -*info;` 28 ` return xerbla_(SCALAR_SUFFIX_UP"GEQRF", &e, 6);` 29 ` }` 30 `// if (query_size)` 31 `// {` 32 `// *lwork = 0;` 33 `// return 0;` 34 `// }` 35 ` if (*m == 0 || *n == 0)` 36 ` {` 37 ` work[0] = 1;` 38 ` return 0;` 39 ` }` 40 ` MatrixType mat(a, *m, *n, OuterStride<>(*lda));` 41 ` CompactVectorType hcoeffs(tau, (std::min)(*m,*n));` 42 ` householder_qr_inplace_blocked(mat, hcoeffs, 32, work);` 43 ` hcoeffs = hcoeffs.conjugate();` 44 ` return 0;` 45 `}` 46 47 ` /* Compute the QR factorization with column pivoting */` 48 ` /* FIXME xgeqpf is deprecated in LAPACK, We should define xgeqp3 instead` 49 ` * However, the implementation of column pivoting QR in Eigen is more similar to xgeqpf than xgeqp3, in terms of performance (BLAS3 operations)` 50 ` */` 51 ` EIGEN_LAPACK_FUNC(geqpf, (int *m, int *n, Scalar *a, int *lda, int *jpvt, Scalar *tau,` 52 ` Scalar *work, int *info))` 53 ` {` 54 ` *info = 0;` 55 ` if(*m < 0) *info = -1;` 56 ` else if(*n < 0) *info = -2;` 57 ` else if(*lda < (std::max)(1, *m)) *info = -4;` 58 59 ` if (*info != 0)` 60 ` {` 61 ` int e = -*info;` 62 ` return xerbla_(SCALAR_SUFFIX_UP"GEQPF", &e, 6);` 63 ` }` 64 ` if((std::min)(*m, *n) == 0) return 0;` 65 66 ` MatrixType qr(a, *m, *n, OuterStride<>(*lda));` 67 ` CompactVectorType hcoeffs(tau, (std::min)(*m, *n));` 68 69 ` Matrix colSqNorms;` 70 ` Matrix colsTranspositions;` 71 ` CompactVectorType temp(work, *n);` 72 ` RealScalar maxpivot;` 73 ` typename MatrixType::Index nonzero_pivots, det_pq;` 74 75 ` //FIXME Find a way to deal with 'free' columns specified in jpvt.` 76 ` PermutationMatrix colsPerm(*n);` 77 78 ` // Call ColPivHouseholder kernel` 79 ` colpivhouseholder_qr_inplace(qr, hcoeffs, colSqNorms, colsTranspositions,` 80 ` colsPerm, temp, maxpivot, nonzero_pivots, det_pq);` 81 ` hcoeffs = hcoeffs.conjugate();` 82 ` for(int i = 0; i < *n; i++) jpvt[i] = colsPerm.indices()(i)+1;` 83 84 ` return 0;` 85 ` }` 86 ` template` 87 ` int apply_reflectors(const char *method, char *side, char *trans, int *m, int *n, int *k,` 88 ` T *a, int *lda, T *tau, T *c, int *ldc, T *work, int *lwork, int *info)` 89 ` {` 90 91 ` *info = 0;` 92 ` bool query_size = (*lwork==-1);` 93 ` int nq; // order of Q` 94 ` int nw; // minimum dimension of Work` 95 ` if(*side == 'L')` 96 ` {` 97 ` nq = *m;` 98 ` nw = *n;` 99 ` } else {` 100 ` nq = *n;` 101 ` nw = *m;` 102 ` }` 103 ` if( (*side != 'L') && (*side != 'R') ) *info = -1;` 104 ` else if( (*trans != 'N') && (*trans != 'T') && (*trans != 'C') ) *info = -2;` 105 ` else if (*m < 0) *info = -3;` 106 ` else if (*n < 0) *info = -4;` 107 ` else if (*k < 0 || *k > nq) *info = -5;` 108 ` else if (*lda < (std::max)(1, nq)) *info = -7;` 109 ` else if (*ldc < (std::max)(1, *m)) *info = -10;` 110 ` else if (*lwork < (std::max)(1, nw) && !query_size) *info = -12;` 111 112 ` if(*info != 0)` 113 ` {` 114 ` int e = -*info;` 115 ` return xerbla_(method, &e, 6);` 116 ` }` 117 ` if(query_size) return 0;` 118 ` if ( *m == 0 || *n == 0 || *k == 0)` 119 ` {` 120 ` work[0] = 1;` 121 ` return 0;` 122 ` }` 123 124 ` //Map input arrays to Eigen matrices` 125 ` //NOTE Here, k should be the the number of nonzero pivots as returned by xgeqpf` 126 ` int row = (*side == 'L') ? *m : *n;` 127 ` Map, 0, OuterStride<> > qr(a, row, *k, OuterStride<>(*lda));` 128 ` Map, 0, OuterStride<> > dst(c, *m, *n, OuterStride<>(*ldc));` 129 ` Map > hcoeffs(tau, *k);` 130 131 ` Matrix matrixQ(*m, *n);` 132 ` matrixQ.setIdentity();` 133 ` matrixQ = householderSequence(qr, hcoeffs);` 134 135 ` // Apply the Householder reflectors` 136 ` if (*side == 'L')` 137 ` {` 138 ` if(*trans == 'N') dst = matrixQ * dst;` 139 ` //dst.applyOnTheLeft(householderSequence(qr, hcoeffs));` 140 ` else if(*trans =='T') dst = matrixQ.transpose() * dst;` 141 ` //dst.applyOnTheLeft(householderSequence(qr, hcoeffs).transpose());` 142 ` else dst = matrixQ.adjoint() * dst;` 143 ` //dst.applyOnTheLeft(householderSequence(qr, hcoeffs.conjugate())/*.adjoint()*/);` 144 ` }` 145 ` else` 146 ` {` 147 ` if(*trans == 'N') dst = dst * matrixQ;` 148 ` //dst.applyOnTheRight(householderSequence(qr, hcoeffs));` 149 ` else if(*trans =='T') dst = dst * matrixQ.transpose();` 150 ` //dst.applyOnTheRight(householderSequence(qr, hcoeffs).transpose());` 151 ` else dst = dst * matrixQ.adjoint();` 152 ` //dst.applyOnTheRight(householderSequence(qr, hcoeffs.conjugate())/*.adjoint()*/);` 153 ` }` 154 ` return 0;` 155 ` }` 156 157 ` /* Apply the product of k elementary reflectors to a real matrix` 158 ` */` 159 ` EIGEN_LAPACK_FUNC(ormqr, (char *side, char *trans, int *m, int *n, int *k, RealScalar *a,` 160 ` int *lda, RealScalar *tau, RealScalar *c, int *ldc, RealScalar *work,` 161 ` int *lwork, int *info))` 162 ` {` 163 ` int err = apply_reflectors(SCALAR_SUFFIX_UP"ORMQR", side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info);` 164 ` return err;` 165 ` }` 166 167 ` /* Apply the product of k elementary reflectors to a complex matrix */` 168 ` EIGEN_LAPACK_FUNC(unmqr, (char *side, char *trans, int *m, int *n, int *k, Complex *a,` 169 ` int *lda, Complex *tau, Complex *c, int *ldc, Complex *work,` 170 ` int *lwork, int *info))` 171 ` {` 172 ` int err = apply_reflectors(SCALAR_SUFFIX_UP"UNMQR", side, trans, m, n, k, a, lda, tau, c, ldc,` 173 ` work, lwork, info);` 174 ` return err;` 175 ` }` 176 177 `template` 178 `int generate_Q(const char *method, int *m, int *n, int *k, T *a, int *lda, T *tau,` 179 ` T *work, int *lwork, int *info)` 180 `{` 181 ` bool query_size = (*lwork == -1);` 182 ` *info = 0;` 183 ` if (*m < 0) *info = -1;` 184 ` else if (*n < 0 || *n > *m) *info = -2;` 185 ` else if (*k < 0 || *k > *n) *info = -3;` 186 ` else if (*lda < (std::max)(1, *m)) *info = -5;` 187 ` else if (*lwork < (std::max)(1, *n) && !query_size) *info = -8;` 188 ` if (*info != 0)` 189 ` {` 190 ` int e = -*info;` 191 ` return xerbla_(method, &e, 6);` 192 ` }` 193 ` if (query_size)` 194 ` {` 195 ` work[0] = 0;` 196 ` return 0;` 197 ` }` 198 ` if (*n <= 0)` 199 ` {` 200 ` work[0] = 1;` 201 ` return 0;` 202 ` }` 203 ` ` 204 ` Map, 0, OuterStride<> > qr(a, *m, *n, OuterStride<>(*lda));` 205 ` Matrix tA(*m, *k);` 206 ` tA = matrix(a, *m, *k, *lda);` 207 ` Map > hcoeffs(tau, *k);` 208 ` qr.setIdentity();` 209 ` qr = householderSequence(tA, hcoeffs) * qr;` 210 ` ` 211 ` return 0;` 212 `}` 213 `/* Generates a m-by-n real matrix Q from the k elementary reflectors */` 214 `EIGEN_LAPACK_FUNC(orgqr, (int *m, int *n, int *k, RealScalar *a, int *lda, RealScalar *tau,` 215 ` RealScalar *work, int *lwork, int *info))` 216 `{` 217 ` int err = generate_Q(SCALAR_SUFFIX_UP"ORGQR", m, n, k, a, lda, tau, work, lwork, info);` 218 ` return err;` 219 `}` 220 221 `/* Generates a m-by-n complex matrix Q from the k elementary reflectors */` 222 `EIGEN_LAPACK_FUNC(ungqr, (int *m, int *n, int *k, Complex *a, int *lda, Complex *tau,` 223 ` Complex *work, int *lwork, int *info))` 224 `{` 225 ` int err = generate_Q(SCALAR_SUFFIX_UP"UNGQR", m, n, k, a, lda, tau, work, lwork, info);` 226 ` return err;` 227 `}`
 Lines 10-17 (-)a/lapack/single.cpp (+1 lines) 10 `#define SCALAR float` 10 `#define SCALAR float` 11 `#define SCALAR_SUFFIX s` 11 `#define SCALAR_SUFFIX s` 12 `#define SCALAR_SUFFIX_UP "S"` 12 `#define SCALAR_SUFFIX_UP "S"` 13 `#define ISCOMPLEX 0` 13 `#define ISCOMPLEX 0` 14 14 15 `#include "cholesky.cpp"` 15 `#include "cholesky.cpp"` 16 `#include "lu.cpp"` 16 `#include "lu.cpp"` 17 `#include "eigenvalues.cpp"` 17 `#include "eigenvalues.cpp"` 18 `#include "qr.cpp"`