11 #ifndef EIGEN_BIDIAGONALIZATION_H
12 #define EIGEN_BIDIAGONALIZATION_H
92 template<
typename MatrixType>
108 tempData = tempVector.data();
111 for (
Index k = 0; ; ++k)
117 mat.col(k).tail(remainingRows)
118 .makeHouseholderInPlace(
mat.coeffRef(k,k),
diagonal[k]);
120 mat.bottomRightCorner(remainingRows, remainingCols)
121 .applyHouseholderOnTheLeft(
mat.col(k).tail(remainingRows-1),
mat.coeff(k,k), tempData);
123 if(k ==
cols-1)
break;
126 mat.row(k).tail(remainingCols)
127 .makeHouseholderInPlace(
mat.coeffRef(k,k+1), upper_diagonal[k]);
129 mat.bottomRightCorner(remainingRows-1, remainingCols)
130 .applyHouseholderOnTheRight(
mat.row(k).tail(remainingCols-1).adjoint(),
mat.coeff(k,k+1), tempData);
151 template<
typename MatrixType>
174 Scalar tau_u, tau_u_prev(0), tau_v;
176 for(
Index k = 0; k < bs; ++k)
178 Index remainingRows = brows - k;
179 Index remainingCols = bcols - k - 1;
181 SubMatType X_k1(
X.block(k,0, remainingRows,k) );
182 SubMatType V_k1(
A.block(k,0, remainingRows,k) );
185 SubColumnType v_k =
A.col(k).tail(remainingRows);
186 v_k -= V_k1 *
Y.row(k).head(k).adjoint();
187 if(k) v_k -= X_k1 *
A.col(k).head(k);
190 v_k.makeHouseholderInPlace(tau_v,
diagonal[k]);
194 SubMatType Y_k (
Y.block(k+1,0, remainingCols, k+1) );
195 SubMatType U_k1 (
A.block(0,k+1, k,remainingCols) );
203 SubColumnType y_k(
Y.col(k).tail(remainingCols) );
206 SubColumnType tmp(
Y.col(k).head(k) );
207 y_k.noalias() =
A.block(k,k+1, remainingRows,remainingCols).adjoint() * v_k;
208 tmp.noalias() = V_k1.adjoint() * v_k;
209 y_k.noalias() -= Y_k.leftCols(k) * tmp;
210 tmp.noalias() = X_k1.adjoint() * v_k;
211 y_k.noalias() -= U_k1.adjoint() * tmp;
216 SubRowType u_k(
A.row(k).tail(remainingCols) );
217 u_k = u_k.conjugate();
219 u_k -= Y_k *
A.row(k).head(k+1).adjoint();
220 if(k) u_k -= U_k1.adjoint() *
X.row(k).head(k).adjoint();
224 u_k.makeHouseholderInPlace(tau_u, upper_diagonal[k]);
232 SubColumnType x_k (
X.col(k).tail(remainingRows-1) );
236 SubColumnType tmp0 (
X.col(k).head(k) ),
237 tmp1 (
X.col(k).head(k+1) );
239 x_k.noalias() =
A.block(k+1,k+1, remainingRows-1,remainingCols) * u_k.transpose();
240 tmp0.noalias() = U_k1 * u_k.transpose();
241 x_k.noalias() -= X_k1.bottomRows(remainingRows-1) * tmp0;
242 tmp1.noalias() = Y_k.adjoint() * u_k.transpose();
243 x_k.noalias() -=
A.block(k+1,0, remainingRows-1,k+1) * tmp1;
246 u_k = u_k.conjugate();
249 if(k>0)
A.coeffRef(k-1,k) = tau_u_prev;
253 A.coeffRef(k-1,k) = tau_u_prev;
255 A.coeffRef(k,k) = tau_v;
259 A.coeffRef(bs-1,bs) = tau_u_prev;
262 if(bcols>bs && brows>bs)
264 SubMatType
A11(
A.bottomRightCorner(brows-bs,bcols-bs) );
265 SubMatType
A10(
A.block(bs,0, brows-bs,bs) );
266 SubMatType A01(
A.block(0,bs, bs,bcols-bs) );
268 A01(bs-1,0) = Literal(1);
269 A11.noalias() -=
A10 *
Y.topLeftCorner(bcols,bs).bottomRows(bcols-bs).adjoint();
270 A11.noalias() -=
X.topLeftCorner(brows,bs).bottomRows(brows-bs) * A01;
283 template<
typename MatrixType,
typename B
idiagType>
285 Index maxBlockSize=32,
298 MatrixType::RowsAtCompileTime,
301 MatrixType::MaxRowsAtCompileTime>
X(
rows,maxBlockSize);
303 MatrixType::ColsAtCompileTime,
306 MatrixType::MaxColsAtCompileTime>
Y(
cols,maxBlockSize);
310 for(k = 0; k <
size; k += blockSize)
330 BlockType
B =
A.block(k,k,brows,bcols);
336 if(k+bs==
cols || bcols<48)
339 &(bidiagonal.template diagonal<0>().coeffRef(k)),
340 &(bidiagonal.template diagonal<1>().coeffRef(k)),
347 upperbidiagonalization_blocked_helper<BlockType>(
B,
348 &(bidiagonal.template diagonal<0>().coeffRef(k)),
349 &(bidiagonal.template diagonal<1>().coeffRef(k)),
351 X.topLeftCorner(brows,bs),
352 Y.topLeftCorner(bcols,bs)
358 template<
typename _MatrixType>
365 eigen_assert(
rows >=
cols &&
"UpperBidiagonalization is only for Arices satisfying rows>=cols.");
372 &(m_bidiagonal.template diagonal<0>().coeffRef(0)),
373 &(m_bidiagonal.template diagonal<1>().coeffRef(0)),
376 m_isInitialized =
true;
380 template<
typename _MatrixType>
388 eigen_assert(
rows >=
cols &&
"UpperBidiagonalization is only for Arices satisfying rows>=cols.");
393 m_isInitialized =
true;
402 template<
typename Derived>
414 #endif // EIGEN_BIDIAGONALIZATION_H