This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 206 | Differences between
and this patch

Collapse All | Expand All

(-)a/Eigen/src/SVD/JacobiSVD.h (-52 / +142 lines)
Lines 61-69 Link Here
61
> struct qr_preconditioner_impl {};
61
> struct qr_preconditioner_impl {};
62
62
63
template<typename MatrixType, int QRPreconditioner, int Case>
63
template<typename MatrixType, int QRPreconditioner, int Case>
64
struct qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
64
class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
65
{
65
{
66
  static bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
66
public:
67
  typedef typename MatrixType::Index Index;
68
  void resize(Index, Index) {}
69
  bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
67
  {
70
  {
68
    return false;
71
    return false;
69
  }
72
  }
Lines 72-205 Link Here
72
/*** preconditioner using FullPivHouseholderQR ***/
75
/*** preconditioner using FullPivHouseholderQR ***/
73
76
74
template<typename MatrixType>
77
template<typename MatrixType>
75
struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
78
class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
76
{
79
{
77
  static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
80
public:
81
  typedef typename MatrixType::Index Index;
82
83
  void resize(Index rows, Index cols)
84
  {
85
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = FullPivHouseholderQR<MatrixType>(rows, cols);
86
  }
87
88
  bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
78
  {
89
  {
79
    if(matrix.rows() > matrix.cols())
90
    if(matrix.rows() > matrix.cols())
80
    {
91
    {
81
      FullPivHouseholderQR<MatrixType> qr(matrix);
92
      m_qr.compute(matrix);
82
      svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
93
      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
83
      if(svd.m_computeFullU) svd.m_matrixU = qr.matrixQ();
94
      if(svd.m_computeFullU) svd.m_matrixU = m_qr.matrixQ();
84
      if(svd.computeV()) svd.m_matrixV = qr.colsPermutation();
95
      if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
85
      return true;
96
      return true;
86
    }
97
    }
87
    return false;
98
    return false;
88
  }
99
  }
100
private:
101
  FullPivHouseholderQR<MatrixType> m_qr;
89
};
102
};
90
103
91
template<typename MatrixType>
104
template<typename MatrixType>
92
struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
105
class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
93
{
106
{
94
  static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
107
public:
108
  typedef typename MatrixType::Index Index;
109
  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
110
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
111
  TransposeTypeWithSameStorageOrder;
112
113
  void resize(Index rows, Index cols)
114
  {
115
    if (cols != m_qr.rows() || rows != m_qr.cols())
116
    {
117
      m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
118
    }
119
    m_adjoint.resize(cols, rows);
120
  }
121
122
  bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
95
  {
123
  {
96
    if(matrix.cols() > matrix.rows())
124
    if(matrix.cols() > matrix.rows())
97
    {
125
    {
98
      typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
126
      m_adjoint = matrix.adjoint();
99
                     MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
127
      m_qr.compute(m_adjoint);
100
              TransposeTypeWithSameStorageOrder;
128
101
      FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
129
      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
102
      svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
130
      if(svd.m_computeFullV) svd.m_matrixV = m_qr.matrixQ();
103
      if(svd.m_computeFullV) svd.m_matrixV = qr.matrixQ();
131
      if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
104
      if(svd.computeU()) svd.m_matrixU = qr.colsPermutation();
105
      return true;
132
      return true;
106
    }
133
    }
107
    else return false;
134
    else return false;
108
  }
135
  }
136
private:
137
  FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
138
  TransposeTypeWithSameStorageOrder m_adjoint;
109
};
139
};
110
140
111
/*** preconditioner using ColPivHouseholderQR ***/
141
/*** preconditioner using ColPivHouseholderQR ***/
112
142
113
template<typename MatrixType>
143
template<typename MatrixType>
114
struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
144
class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
115
{
145
{
116
  static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
146
public:
147
  typedef typename MatrixType::Index Index;
148
149
  void resize(Index rows, Index cols)
150
  {
151
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = ColPivHouseholderQR<MatrixType>(rows, cols);
152
  }
153
154
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
117
  {
155
  {
118
    if(matrix.rows() > matrix.cols())
156
    if(matrix.rows() > matrix.cols())
119
    {
157
    {
120
      ColPivHouseholderQR<MatrixType> qr(matrix);
158
      m_qr.compute(matrix);
121
      svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
159
      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
122
      if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ();
160
      if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ();
123
      else if(svd.m_computeThinU) {
161
      else if(svd.m_computeThinU) {
124
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
162
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
125
        qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
163
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
126
      }
164
      }
127
      if(svd.computeV()) svd.m_matrixV = qr.colsPermutation();
165
      if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
128
      return true;
166
      return true;
129
    }
167
    }
130
    return false;
168
    return false;
131
  }
169
  }
170
171
private:
172
  ColPivHouseholderQR<MatrixType> m_qr;
132
};
173
};
133
174
134
template<typename MatrixType>
175
template<typename MatrixType>
135
struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
176
class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
136
{
177
{
137
  static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
178
public:
179
  typedef typename MatrixType::Index Index;
180
  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
181
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
182
  TransposeTypeWithSameStorageOrder;
183
184
  void resize(Index rows, Index cols)
185
  {
186
    if (cols != m_qr.rows() || rows != m_qr.cols())
187
    {
188
      m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
189
    }
190
    m_adjoint.resize(cols, rows);
191
  }
192
193
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
138
  {
194
  {
139
    if(matrix.cols() > matrix.rows())
195
    if(matrix.cols() > matrix.rows())
140
    {
196
    {
141
      typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
197
      m_adjoint = matrix.adjoint();
142
                     MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
198
      m_qr.compute(m_adjoint);
143
              TransposeTypeWithSameStorageOrder;
199
144
      ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
200
      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
145
      svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
201
      if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ();
146
      if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ();
147
      else if(svd.m_computeThinV) {
202
      else if(svd.m_computeThinV) {
148
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
203
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
149
        qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
204
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
150
      }
205
      }
151
      if(svd.computeU()) svd.m_matrixU = qr.colsPermutation();
206
      if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
152
      return true;
207
      return true;
153
    }
208
    }
154
    else return false;
209
    else return false;
155
  }
210
  }
211
212
private:
213
  ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
214
  TransposeTypeWithSameStorageOrder m_adjoint;
156
};
215
};
157
216
158
/*** preconditioner using HouseholderQR ***/
217
/*** preconditioner using HouseholderQR ***/
159
218
160
template<typename MatrixType>
219
template<typename MatrixType>
161
struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
220
class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
162
{
221
{
163
  static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
222
public:
223
  typedef typename MatrixType::Index Index;
224
225
  void resize(Index rows, Index cols)
226
  {
227
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = HouseholderQR<MatrixType>(rows, cols);
228
  }
229
230
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
164
  {
231
  {
165
    if(matrix.rows() > matrix.cols())
232
    if(matrix.rows() > matrix.cols())
166
    {
233
    {
167
      HouseholderQR<MatrixType> qr(matrix);
234
      m_qr.compute(matrix);
168
      svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
235
      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
169
      if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ();
236
      if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ();
170
      else if(svd.m_computeThinU) {
237
      else if(svd.m_computeThinU) {
171
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
238
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
172
        qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
239
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
173
      }
240
      }
174
      if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
241
      if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
175
      return true;
242
      return true;
176
    }
243
    }
177
    return false;
244
    return false;
178
  }
245
  }
246
private:
247
  HouseholderQR<MatrixType> m_qr;
179
};
248
};
180
249
181
template<typename MatrixType>
250
template<typename MatrixType>
182
struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
251
class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
183
{
252
{
184
  static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
253
public:
254
  typedef typename MatrixType::Index Index;
255
  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
256
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
257
  TransposeTypeWithSameStorageOrder;
258
259
  void resize(Index rows, Index cols)
260
  {
261
    if (cols != m_qr.rows() || rows != m_qr.cols())
262
    {
263
      m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
264
    }
265
    m_adjoint.resize(cols, rows);
266
  }
267
268
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
185
  {
269
  {
186
    if(matrix.cols() > matrix.rows())
270
    if(matrix.cols() > matrix.rows())
187
    {
271
    {
188
      typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
272
      m_adjoint = matrix.adjoint();
189
                     MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
273
      m_qr.compute(m_adjoint);
190
              TransposeTypeWithSameStorageOrder;
274
191
      HouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
275
      svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
192
      svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
276
      if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ();
193
      if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ();
194
      else if(svd.m_computeThinV) {
277
      else if(svd.m_computeThinV) {
195
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
278
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
196
        qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
279
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
197
      }
280
      }
198
      if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
281
      if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
199
      return true;
282
      return true;
200
    }
283
    }
201
    else return false;
284
    else return false;
202
  }
285
  }
286
287
private:
288
  HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
289
  TransposeTypeWithSameStorageOrder m_adjoint;
203
};
290
};
204
291
205
/*** 2x2 SVD implementation
292
/*** 2x2 SVD implementation
Lines 510-515 Link Here
510
    friend struct internal::svd_precondition_2x2_block_to_be_real;
597
    friend struct internal::svd_precondition_2x2_block_to_be_real;
511
    template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
598
    template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
512
    friend struct internal::qr_preconditioner_impl;
599
    friend struct internal::qr_preconditioner_impl;
600
601
    internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
602
    internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
513
};
603
};
514
604
515
template<typename MatrixType, int QRPreconditioner>
605
template<typename MatrixType, int QRPreconditioner>
Lines 541-546 Link Here
541
                          : m_computeThinV ? m_diagSize
631
                          : m_computeThinV ? m_diagSize
542
                          : 0);
632
                          : 0);
543
  m_workMatrix.resize(m_diagSize, m_diagSize);
633
  m_workMatrix.resize(m_diagSize, m_diagSize);
634
  m_qr_precond_morecols.resize(rows, cols);
635
  m_qr_precond_morerows.resize(rows, cols);
544
}
636
}
545
637
546
template<typename MatrixType, int QRPreconditioner>
638
template<typename MatrixType, int QRPreconditioner>
Lines 554-562 Link Here
554
  const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
646
  const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
555
647
556
  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
648
  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
557
649
  if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
558
  if(!internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>::run(*this, matrix)
559
  && !internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>::run(*this, matrix))
560
  {
650
  {
561
    m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
651
    m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
562
    if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
652
    if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
(-)a/Eigen/src/QR/FullPivHouseholderQR.h (-8 / +17 lines)
Lines 65-70 Link Here
65
    typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
65
    typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
66
    typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
66
    typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
67
    typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
67
    typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
68
    typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> TempVectorType;
68
    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
69
    typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
69
    typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
70
    typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
70
    typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
71
    typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
Lines 82-87 Link Here
82
        m_cols_transpositions(),
83
        m_cols_transpositions(),
83
        m_cols_permutation(),
84
        m_cols_permutation(),
84
        m_temp(),
85
        m_temp(),
86
        m_hhTemp(),
87
        m_q(),
85
        m_isInitialized(false),
88
        m_isInitialized(false),
86
        m_usePrescribedThreshold(false) {}
89
        m_usePrescribedThreshold(false) {}
87
90
Lines 98-103 Link Here
98
        m_cols_transpositions(cols),
101
        m_cols_transpositions(cols),
99
        m_cols_permutation(cols),
102
        m_cols_permutation(cols),
100
        m_temp(std::min(rows,cols)),
103
        m_temp(std::min(rows,cols)),
104
        m_hhTemp(rows),
105
        m_q(rows, rows),
101
        m_isInitialized(false),
106
        m_isInitialized(false),
102
        m_usePrescribedThreshold(false) {}
107
        m_usePrescribedThreshold(false) {}
103
108
Lines 108-113 Link Here
108
        m_cols_transpositions(matrix.cols()),
113
        m_cols_transpositions(matrix.cols()),
109
        m_cols_permutation(matrix.cols()),
114
        m_cols_permutation(matrix.cols()),
110
        m_temp(std::min(matrix.rows(), matrix.cols())),
115
        m_temp(std::min(matrix.rows(), matrix.cols())),
116
        m_hhTemp(matrix.rows()),
117
        m_q(matrix.rows(), matrix.rows()),
111
        m_isInitialized(false),
118
        m_isInitialized(false),
112
        m_usePrescribedThreshold(false)
119
        m_usePrescribedThreshold(false)
113
    {
120
    {
Lines 139-145 Link Here
139
      return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
146
      return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
140
    }
147
    }
141
148
142
    MatrixQType matrixQ(void) const;
149
    const MatrixQType& matrixQ(void);
143
150
144
    /** \returns a reference to the matrix where the Householder QR decomposition is stored
151
    /** \returns a reference to the matrix where the Householder QR decomposition is stored
145
      */
152
      */
Lines 351-356 Link Here
351
    IntRowVectorType m_cols_transpositions;
358
    IntRowVectorType m_cols_transpositions;
352
    PermutationType m_cols_permutation;
359
    PermutationType m_cols_permutation;
353
    RowVectorType m_temp;
360
    RowVectorType m_temp;
361
    TempVectorType m_hhTemp;
362
    MatrixQType m_q;
354
    bool m_isInitialized, m_usePrescribedThreshold;
363
    bool m_isInitialized, m_usePrescribedThreshold;
355
    RealScalar m_prescribedThreshold, m_maxpivot;
364
    RealScalar m_prescribedThreshold, m_maxpivot;
356
    Index m_nonzero_pivots;
365
    Index m_nonzero_pivots;
Lines 512-518 Link Here
512
521
513
/** \returns the matrix Q */
522
/** \returns the matrix Q */
514
template<typename MatrixType>
523
template<typename MatrixType>
515
typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<MatrixType>::matrixQ() const
524
const typename FullPivHouseholderQR<MatrixType>::MatrixQType& FullPivHouseholderQR<MatrixType>::matrixQ()
516
{
525
{
517
  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
526
  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
518
  // compute the product H'_0 H'_1 ... H'_n-1,
527
  // compute the product H'_0 H'_1 ... H'_n-1,
Lines 521-535 Link Here
521
  Index rows = m_qr.rows();
530
  Index rows = m_qr.rows();
522
  Index cols = m_qr.cols();
531
  Index cols = m_qr.cols();
523
  Index size = std::min(rows,cols);
532
  Index size = std::min(rows,cols);
524
  MatrixQType res = MatrixQType::Identity(rows, rows);
533
  m_q.setIdentity(rows, rows);
525
  Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
534
  m_hhTemp.resize(rows);
526
  for (Index k = size-1; k >= 0; k--)
535
  for (Index k = size-1; k >= 0; k--)
527
  {
536
  {
528
    res.block(k, k, rows-k, rows-k)
537
    m_q.block(k, k, rows-k, rows-k)
529
       .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
538
       .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &m_hhTemp.coeffRef(k));
530
    res.row(k).swap(res.row(m_rows_transpositions.coeff(k)));
539
    m_q.row(k).swap(m_q.row(m_rows_transpositions.coeff(k)));
531
  }
540
  }
532
  return res;
541
  return m_q;
533
}
542
}
534
543
535
/** \return the full-pivoting Householder QR decomposition of \c *this.
544
/** \return the full-pivoting Householder QR decomposition of \c *this.
(-)a/Eigen/src/Householder/HouseholderSequence.h (-4 / +22 lines)
Lines 285-308 Link Here
285
    /** \internal */
285
    /** \internal */
286
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
286
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
287
    {
287
    {
288
      Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
288
       Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> temp(dst.rows());
289
      applyThisOnTheRight(dst, temp);
290
    }
291
292
    /** \internal */
293
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst,
294
                                                            Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,
295
                                                                   1,Dest::MaxRowsAtCompileTime>& temp) const
296
    {
297
      temp.resize(dst.rows());
289
      for(Index k = 0; k < m_length; ++k)
298
      for(Index k = 0; k < m_length; ++k)
290
      {
299
      {
291
        Index actual_k = m_trans ? m_length-k-1 : k;
300
        Index actual_k = m_trans ? m_length-k-1 : k;
292
        dst.rightCols(rows()-m_shift-actual_k)
301
        dst.rightCols(rows()-m_shift-actual_k)
293
           .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
302
        .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
294
      }
303
      }
295
    }
304
    }
296
305
297
    /** \internal */
306
    /** \internal */
298
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
307
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
299
    {
308
    {
300
      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
309
      Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> temp(dst.cols());
310
      applyThisOnTheLeft(dst, temp);
311
    }
312
313
    /** \internal */
314
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst,
315
                                                           Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,
316
                                                           1,Dest::MaxColsAtCompileTime>& temp) const
317
    {
318
      temp.resize(dst.cols());
301
      for(Index k = 0; k < m_length; ++k)
319
      for(Index k = 0; k < m_length; ++k)
302
      {
320
      {
303
        Index actual_k = m_trans ? k : m_length-k-1;
321
        Index actual_k = m_trans ? k : m_length-k-1;
304
        dst.bottomRows(rows()-m_shift-actual_k)
322
        dst.bottomRows(rows()-m_shift-actual_k)
305
           .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
323
        .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
306
      }
324
      }
307
    }
325
    }
308
326
(-)a/Eigen/src/SVD/JacobiSVD.h (-6 / +32 lines)
Lines 145-154 Link Here
145
{
145
{
146
public:
146
public:
147
  typedef typename MatrixType::Index Index;
147
  typedef typename MatrixType::Index Index;
148
  typedef Matrix<typename MatrixType::Scalar,1, MatrixType::RowsAtCompileTime, RowMajor,
149
  1, MatrixType::MaxRowsAtCompileTime> RowVectorType;
148
150
149
  void resize(Index rows, Index cols)
151
  void resize(Index rows, Index cols)
150
  {
152
  {
151
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = ColPivHouseholderQR<MatrixType>(rows, cols);
153
    if (rows != m_qr.rows() || cols != m_qr.cols())
154
    {
155
      m_qr = ColPivHouseholderQR<MatrixType>(rows, cols);
156
    }
157
    m_vec.resize(rows);
152
  }
158
  }
153
159
154
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
160
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 160-166 Link Here
160
      if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ();
166
      if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ();
161
      else if(svd.m_computeThinU) {
167
      else if(svd.m_computeThinU) {
162
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
168
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
163
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
169
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_vec);
164
      }
170
      }
165
      if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
171
      if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
166
      return true;
172
      return true;
Lines 170-175 Link Here
170
176
171
private:
177
private:
172
  ColPivHouseholderQR<MatrixType> m_qr;
178
  ColPivHouseholderQR<MatrixType> m_qr;
179
  RowVectorType m_vec;
173
};
180
};
174
181
175
template<typename MatrixType>
182
template<typename MatrixType>
Lines 177-182 Link Here
177
{
184
{
178
public:
185
public:
179
  typedef typename MatrixType::Index Index;
186
  typedef typename MatrixType::Index Index;
187
188
  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::ColsAtCompileTime, RowMajor,
189
  1, MatrixType::MaxColsAtCompileTime> RowVectorType;
190
180
  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
191
  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
181
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
192
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
182
  TransposeTypeWithSameStorageOrder;
193
  TransposeTypeWithSameStorageOrder;
Lines 188-193 Link Here
188
      m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
199
      m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
189
    }
200
    }
190
    m_adjoint.resize(cols, rows);
201
    m_adjoint.resize(cols, rows);
202
    m_vec.resize(cols);
191
  }
203
  }
192
204
193
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
205
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 201-207 Link Here
201
      if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ();
213
      if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ();
202
      else if(svd.m_computeThinV) {
214
      else if(svd.m_computeThinV) {
203
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
215
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
204
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
216
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_vec);
205
      }
217
      }
206
      if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
218
      if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
207
      return true;
219
      return true;
Lines 212-217 Link Here
212
private:
224
private:
213
  ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
225
  ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
214
  TransposeTypeWithSameStorageOrder m_adjoint;
226
  TransposeTypeWithSameStorageOrder m_adjoint;
227
  RowVectorType m_vec;
215
};
228
};
216
229
217
/*** preconditioner using HouseholderQR ***/
230
/*** preconditioner using HouseholderQR ***/
Lines 221-230 Link Here
221
{
234
{
222
public:
235
public:
223
  typedef typename MatrixType::Index Index;
236
  typedef typename MatrixType::Index Index;
237
  typedef Matrix<typename MatrixType::Scalar,1, MatrixType::RowsAtCompileTime, RowMajor,
238
  1, MatrixType::MaxRowsAtCompileTime> RowVectorType;
224
239
225
  void resize(Index rows, Index cols)
240
  void resize(Index rows, Index cols)
226
  {
241
  {
227
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = HouseholderQR<MatrixType>(rows, cols);
242
    if (rows != m_qr.rows() || cols != m_qr.cols())
243
    {
244
      m_qr = HouseholderQR<MatrixType>(rows, cols);
245
    }
246
    m_vec.resize(rows);
228
  }
247
  }
229
248
230
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
249
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 236-242 Link Here
236
      if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ();
255
      if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ();
237
      else if(svd.m_computeThinU) {
256
      else if(svd.m_computeThinU) {
238
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
257
        svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
239
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
258
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_vec);
240
      }
259
      }
241
      if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
260
      if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
242
      return true;
261
      return true;
Lines 245-250 Link Here
245
  }
264
  }
246
private:
265
private:
247
  HouseholderQR<MatrixType> m_qr;
266
  HouseholderQR<MatrixType> m_qr;
267
  RowVectorType m_vec;
248
};
268
};
249
269
250
template<typename MatrixType>
270
template<typename MatrixType>
Lines 252-257 Link Here
252
{
272
{
253
public:
273
public:
254
  typedef typename MatrixType::Index Index;
274
  typedef typename MatrixType::Index Index;
275
276
  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::ColsAtCompileTime, RowMajor,
277
  1, MatrixType::MaxColsAtCompileTime> RowVectorType;
278
255
  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
279
  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
256
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
280
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
257
  TransposeTypeWithSameStorageOrder;
281
  TransposeTypeWithSameStorageOrder;
Lines 263-268 Link Here
263
      m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
287
      m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
264
    }
288
    }
265
    m_adjoint.resize(cols, rows);
289
    m_adjoint.resize(cols, rows);
290
    m_vec.resize(cols);
266
  }
291
  }
267
292
268
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
293
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 276-282 Link Here
276
      if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ();
301
      if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ();
277
      else if(svd.m_computeThinV) {
302
      else if(svd.m_computeThinV) {
278
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
303
        svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
279
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
304
        m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_vec);
280
      }
305
      }
281
      if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
306
      if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
282
      return true;
307
      return true;
Lines 287-292 Link Here
287
private:
312
private:
288
  HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
313
  HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
289
  TransposeTypeWithSameStorageOrder m_adjoint;
314
  TransposeTypeWithSameStorageOrder m_adjoint;
315
  RowVectorType m_vec;
290
};
316
};
291
317
292
/*** 2x2 SVD implementation
318
/*** 2x2 SVD implementation
(-)a/Eigen/src/SVD/JacobiSVD.h (-4 / +4 lines)
Lines 154-160 Link Here
154
    {
154
    {
155
      m_qr = ColPivHouseholderQR<MatrixType>(rows, cols);
155
      m_qr = ColPivHouseholderQR<MatrixType>(rows, cols);
156
    }
156
    }
157
    m_vec.resize(rows);
157
    m_vec.resize(cols);
158
  }
158
  }
159
159
160
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
160
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 199-205 Link Here
199
      m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
199
      m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
200
    }
200
    }
201
    m_adjoint.resize(cols, rows);
201
    m_adjoint.resize(cols, rows);
202
    m_vec.resize(cols);
202
    m_vec.resize(rows);
203
  }
203
  }
204
204
205
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
205
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 243-249 Link Here
243
    {
243
    {
244
      m_qr = HouseholderQR<MatrixType>(rows, cols);
244
      m_qr = HouseholderQR<MatrixType>(rows, cols);
245
    }
245
    }
246
    m_vec.resize(rows);
246
    m_vec.resize(cols);
247
  }
247
  }
248
248
249
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
249
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 287-293 Link Here
287
      m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
287
      m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
288
    }
288
    }
289
    m_adjoint.resize(cols, rows);
289
    m_adjoint.resize(cols, rows);
290
    m_vec.resize(cols);
290
    m_vec.resize(rows);
291
  }
291
  }
292
292
293
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
293
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
(-)a/Eigen/src/SVD/JacobiSVD.h (-31 / +19 lines)
Lines 65-71 Link Here
65
{
65
{
66
public:
66
public:
67
  typedef typename MatrixType::Index Index;
67
  typedef typename MatrixType::Index Index;
68
  void resize(Index, Index) {}
68
  void resize(Index, Index, bool) {}
69
  bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
69
  bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
70
  {
70
  {
71
    return false;
71
    return false;
Lines 80-86 Link Here
80
public:
80
public:
81
  typedef typename MatrixType::Index Index;
81
  typedef typename MatrixType::Index Index;
82
82
83
  void resize(Index rows, Index cols)
83
  void resize(Index rows, Index cols, bool)
84
  {
84
  {
85
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = FullPivHouseholderQR<MatrixType>(rows, cols);
85
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = FullPivHouseholderQR<MatrixType>(rows, cols);
86
  }
86
  }
Lines 110-116 Link Here
110
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
110
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
111
  TransposeTypeWithSameStorageOrder;
111
  TransposeTypeWithSameStorageOrder;
112
112
113
  void resize(Index rows, Index cols)
113
  void resize(Index rows, Index cols, bool)
114
  {
114
  {
115
    if (cols != m_qr.rows() || rows != m_qr.cols())
115
    if (cols != m_qr.rows() || rows != m_qr.cols())
116
    {
116
    {
Lines 148-160 Link Here
148
  typedef Matrix<typename MatrixType::Scalar,1, MatrixType::RowsAtCompileTime, RowMajor,
148
  typedef Matrix<typename MatrixType::Scalar,1, MatrixType::RowsAtCompileTime, RowMajor,
149
  1, MatrixType::MaxRowsAtCompileTime> RowVectorType;
149
  1, MatrixType::MaxRowsAtCompileTime> RowVectorType;
150
150
151
  void resize(Index rows, Index cols)
151
  void resize(Index rows, Index cols, bool computeThinU)
152
  {
152
  {
153
    if (rows != m_qr.rows() || cols != m_qr.cols())
153
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = ColPivHouseholderQR<MatrixType>(rows, cols);
154
    {
154
    if (computeThinU) m_vec.resize(cols);
155
      m_qr = ColPivHouseholderQR<MatrixType>(rows, cols);
156
    }
157
    m_vec.resize(cols);
158
  }
155
  }
159
156
160
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
157
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 192-205 Link Here
192
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
189
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
193
  TransposeTypeWithSameStorageOrder;
190
  TransposeTypeWithSameStorageOrder;
194
191
195
  void resize(Index rows, Index cols)
192
  void resize(Index rows, Index cols, bool computeThinV)
196
  {
193
  {
197
    if (cols != m_qr.rows() || rows != m_qr.cols())
194
    if (cols != m_qr.rows() || rows != m_qr.cols()) m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
198
    {
195
    if (computeThinV) m_vec.resize(rows);
199
      m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
200
    }
201
    m_adjoint.resize(cols, rows);
196
    m_adjoint.resize(cols, rows);
202
    m_vec.resize(rows);
203
  }
197
  }
204
198
205
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
199
  bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 223-230 Link Here
223
217
224
private:
218
private:
225
  ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
219
  ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
220
  RowVectorType m_vec;
226
  TransposeTypeWithSameStorageOrder m_adjoint;
221
  TransposeTypeWithSameStorageOrder m_adjoint;
227
  RowVectorType m_vec;
228
};
222
};
229
223
230
/*** preconditioner using HouseholderQR ***/
224
/*** preconditioner using HouseholderQR ***/
Lines 237-249 Link Here
237
  typedef Matrix<typename MatrixType::Scalar,1, MatrixType::RowsAtCompileTime, RowMajor,
231
  typedef Matrix<typename MatrixType::Scalar,1, MatrixType::RowsAtCompileTime, RowMajor,
238
  1, MatrixType::MaxRowsAtCompileTime> RowVectorType;
232
  1, MatrixType::MaxRowsAtCompileTime> RowVectorType;
239
233
240
  void resize(Index rows, Index cols)
234
  void resize(Index rows, Index cols, bool computeThinU)
241
  {
235
  {
242
    if (rows != m_qr.rows() || cols != m_qr.cols())
236
    if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = HouseholderQR<MatrixType>(rows, cols);
243
    {
237
    if (computeThinU) m_vec.resize(cols);
244
      m_qr = HouseholderQR<MatrixType>(rows, cols);
245
    }
246
    m_vec.resize(cols);
247
  }
238
  }
248
239
249
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
240
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 280-293 Link Here
280
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
271
  MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
281
  TransposeTypeWithSameStorageOrder;
272
  TransposeTypeWithSameStorageOrder;
282
273
283
  void resize(Index rows, Index cols)
274
  void resize(Index rows, Index cols, bool computeThinV)
284
  {
275
  {
285
    if (cols != m_qr.rows() || rows != m_qr.cols())
276
    if (cols != m_qr.rows() || rows != m_qr.cols()) m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
286
    {
277
    if (computeThinV) m_vec.resize(rows);
287
      m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows);
288
    }
289
    m_adjoint.resize(cols, rows);
278
    m_adjoint.resize(cols, rows);
290
    m_vec.resize(rows);
291
  }
279
  }
292
280
293
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
281
  bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
Lines 311-318 Link Here
311
299
312
private:
300
private:
313
  HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
301
  HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
302
  RowVectorType m_vec;
314
  TransposeTypeWithSameStorageOrder m_adjoint;
303
  TransposeTypeWithSameStorageOrder m_adjoint;
315
  RowVectorType m_vec;
316
};
304
};
317
305
318
/*** 2x2 SVD implementation
306
/*** 2x2 SVD implementation
Lines 657-664 Link Here
657
                          : m_computeThinV ? m_diagSize
645
                          : m_computeThinV ? m_diagSize
658
                          : 0);
646
                          : 0);
659
  m_workMatrix.resize(m_diagSize, m_diagSize);
647
  m_workMatrix.resize(m_diagSize, m_diagSize);
660
  m_qr_precond_morecols.resize(rows, cols);
648
  m_qr_precond_morecols.resize(rows, cols, m_computeThinV);
661
  m_qr_precond_morerows.resize(rows, cols);
649
  m_qr_precond_morerows.resize(rows, cols, m_computeThinU);
662
}
650
}
663
651
664
template<typename MatrixType, int QRPreconditioner>
652
template<typename MatrixType, int QRPreconditioner>

Return to bug 206