1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@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/.
#ifndef EIGEN_SUITESPARSEQRSUPPORT_H
#define EIGEN_SUITESPARSEQRSUPPORT_H
namespace Eigen {
template<typename MatrixType> class SPQR;
template<typename SPQRType> struct SPQRMatrixQReturnType;
template<typename SPQRType> struct SPQRMatrixQTransposeReturnType;
template <typename SPQRType, typename Derived> struct SPQR_QProduct;
namespace internal {
template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> >
{
typedef typename SPQRType::MatrixType ReturnType;
};
template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> >
{
typedef typename SPQRType::MatrixType ReturnType;
};
template <typename SPQRType, typename Derived> struct traits<SPQR_QProduct<SPQRType, Derived> >
{
typedef typename Derived::PlainObject ReturnType;
};
} // End namespace internal
/**
* \ingroup SPQRSupport_Module
* \class SPQR
* \brief Sparse QR factorization based on SuiteSparseQR library
*
* This class is used to perform a multithreaded and multifrontal QR decomposition
* of sparse matrices. The result is then used to solve linear leasts_square systems.
* Clearly, a QR factorization is returned such that A*P = Q*R where :
*
* P is the column permutation. Use colsPermutation() to get it.
*
* Q is the orthogonal matrix represented as Householder reflectors.
* Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
* You can then apply it to a vector.
*
* R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix.
* NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index
*
* \tparam _MatrixType The type of the sparse matrix A, must be a SparseMatrix<>, either row-major or column-major.
* NOTE
*
*/
template<typename _MatrixType>
class SPQR
{
public:
typedef typename _MatrixType::Scalar Scalar;
typedef typename _MatrixType::RealScalar RealScalar;
typedef UF_long Index ;
typedef SparseMatrix<Scalar, _MatrixType::Flags, Index> MatrixType;
public:
SPQR()
: m_ordering(SPQR_ORDERING_DEFAULT),
m_allow_tol(SPQR_DEFAULT_TOL),
m_tolerance (NumTraits<Scalar>::epsilon())
{
cholmod_l_start(&m_cc);
}
SPQR(const _MatrixType& matrix) : SPQR()
{
compute(matrix);
}
~SPQR()
{
// Calls SuiteSparseQR_free()
cholmod_free_sparse(&m_H, &m_cc);
cholmod_free_dense(&m_HTau, &m_cc);
delete[] m_E;
delete[] m_HPinv;
}
void compute(const MatrixType& matrix)
{
MatrixType mat(matrix);
cholmod_sparse A;
A = viewAsCholmod(mat);
Index col = matrix.cols();
m_rank = SuiteSparseQR<Scalar>(m_ordering, m_tolerance, col, &A,
&m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
if (!m_cR)
{
m_info = NumericalIssue;
m_isInitialized = false;
return;
}
m_info = Success;
m_isInitialized = true;
}
template<typename Rhs, typename Dest>
void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
{
eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
eigen_assert(b.cols()==1 && "This method is for vectors only");
//Compute Q^T * b
// NOTE : We may have called directly the corresponding routines in SPQR codes.
// This version is used to test directly the corresponding part of the code
dest = matrixQ().transpose() * b;
// Solves with the triangular matrix R
Dest y;
y = this->matrixQR().template triangularView<Upper>().solve(dest.derived());
// Apply the column permutation //TODO Check the terminology behind the permutation
for (int j = 0; j < y.size(); j++) dest(m_E[j]) = y(j);
m_info = Success;
}
/// Get the sparse triangular matrix R. It is a sparse matrix
MatrixType matrixQR() const
{
MatrixType R;
R = viewAsEigen<Scalar, MatrixType::Flags, Index>(*m_cR);
return R;
}
/// Get an expression of the matrix Q
SPQRMatrixQReturnType<SPQR> matrixQ() const
{
return SPQRMatrixQReturnType<SPQR>(*this);
}
/// Get the permutation that was applied to columns of A
Index *colsPermutation() { return m_E; }
/// Set the fill-reducing ordering method to be used
void setOrdering(int ord) { m_ordering = ord;}
/// Set the tolerance tol to treat columns with 2-norm < =tol as zero
void setTolerance(RealScalar tol) { m_tolerance = tol; }
/// Return a pointer to SPQR workspace
cholmod_common *cc() const { return &m_cc; }
cholmod_sparse * H() const { return m_H; }
Index *HPinv() const { return m_HPinv; }
cholmod_dense* HTau() const { return m_HTau; }
/** \brief Reports whether previous computation was successful.
*
* \returns \c Success if computation was succesful,
* \c NumericalIssue if the sparse QR can not be computed
*/
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "Decomposition is not initialized.");
return m_info;
}
protected:
bool m_isInitialized;
bool m_analysisIsOk;
bool m_factorizationIsOk;
mutable ComputationInfo m_info;
int m_ordering; // Ordering method to use, see SPQR's manual
int m_allow_tol; // Allow to use some tolerance during numerical factorization.
RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
mutable Index *m_E; // The permutation applied to columns
mutable cholmod_sparse *m_H; //The householder vectors
mutable Index *m_HPinv; // The row permutation of H
mutable cholmod_dense *m_HTau; // The Householder coefficients
mutable Index m_rank; // The rank of the matrix
mutable cholmod_common m_cc; // Workspace and parameters
};
template <typename SPQRType, typename Derived>
struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
{
typedef typename SPQRType::Scalar Scalar;
//Define the constructor to get reference to argument types
SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
// Assign to a vector
template<typename ResType>
void evalTo(ResType& res) const
{
cholmod_dense y_cd;
cholmod_dense *x_cd;
int method = m_transpose ? SPQR_QTX : SPQR_QX;
cholmod_common *cc = m_spqr.cc();
y_cd = viewAsCholmod(m_other.const_cast_derived());
x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.H(), m_spqr.HTau(), m_spqr.HPinv(), &y_cd, cc);
res = Matrix<Scalar,ResType::RowsAtCompileTime,ResType::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol);
cholmod_free_dense(&x_cd, cc);
}
const SPQRType& m_spqr;
const Derived& m_other;
bool m_transpose;
};
template<typename SPQRType>
struct SPQRMatrixQReturnType{
SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
template<typename Derived>
SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other)
{
return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(),false);
}
// To use for operations with the transpose of Q
SPQRMatrixQTransposeReturnType<SPQRType> transpose() const
{
return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr);
}
const SPQRType& m_spqr;
};
template<typename SPQRType>
struct SPQRMatrixQTransposeReturnType{
SPQRMatrixQTransposeReturnType(const SPQRType& spqr) : m_spqr(spqr) {}
template<typename Derived>
SPQR_QProduct<SPQRType,Derived> operator*(const MatrixBase<Derived>& other)
{
return SPQR_QProduct<SPQRType,Derived>(m_spqr,other.derived(), true);
}
const SPQRType& m_spqr;
};
}// End namespace Eigen
#endif
|
