std::vector<Vector<number> > & eigenvectors,
const int itype = 1);
+ /**
+ * Compute the singular value
+ * decomposition of the
+ * matrix using LAPACK function
+ * Xgesdd.
+ *
+ * Requires that the #state is
+ * LAPACKSupport::matrix, fills
+ * the data members #wr, #svd_u,
+ * and #svd_vt, and leaves the
+ * object in the #state
+ * LAPACKSupport::svd.
+ */
+ void compute_svd();
+ /**
+ * Compute the inverse of the
+ * matrix by singular value
+ * decomposition.
+ *
+ * Requires that #state is either
+ * LAPACKSupport::matrix or
+ * LAPACKSupport::svd. In the
+ * first case, this function
+ * calls compute_svd(). After
+ * this function, the object will
+ * have the #state
+ * LAPACKSupport::inverse_svd.
+ *
+ * For a singular value
+ * decomposition, the inverse is
+ * simply computed by replacing
+ * all singular values by their
+ * reciprocal values. If the
+ * matrix does not have maximal
+ * rank, singular values 0 are
+ * not touched, thus computing
+ * the minimal norm right inverse
+ * of the matrix.
+ *
+ * The parameter
+ * <tt>threshold</tt> determines,
+ * when a singular value should
+ * be considered zero. It is the
+ * ratio of the smallest to the
+ * largest nonzero singular
+ * value
+ * <i>s</i><sub>max</sub>. Thus,
+ * the inverses of all singular
+ * values less than
+ * <i>s</i><sub>max</sub>/<tt>threshold</tt>
+ * will be set to zero.
+ */
+ void compute_inverse_svd (const double threshold = 0.);
+
/**
* Retrieve eigenvalue after
* compute_eigenvalues() was
std::complex<number>
eigenvalue (const unsigned int i) const;
+ /**
+ * Retrieve singular values after
+ * compute_svd() or
+ * compute_inverse_svd() was
+ * called.
+ */
+ number
+ singular_value (const unsigned int i) const;
+
/**
* Print the matrix and allow
* formatting of entries.
* permutations applied for
* pivoting in the
* LU-factorization.
+ *
+ * Also used as the scratch array
+ * IWORK for LAPACK functions
+ * needing it.
*/
std::vector<int> ipiv;
*/
std::vector<number> vr;
- /**
- * The matrix <i>U</i> in the singular value decomposition
- * <i>USV<sup>T</sup></i>.
- */
+ /**
+ * The matrix <i>U</i> in the
+ * singular value decomposition
+ * <i>USV<sup>T</sup></i>.
+ */
boost::shared_ptr<LAPACKFullMatrix<number> > svd_u;
- /**
- * The matrix <i>V<sup>T</sup></i> in the singular value decomposition
- * <i>USV<sup>T</sup></i>.
- */
+ /**
+ * The matrix
+ * <i>V<sup>T</sup></i> in the
+ * singular value decomposition
+ * <i>USV<sup>T</sup></i>.
+ */
boost::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
};
}
+template <typename number>
+number
+LAPACKFullMatrix<number>::singular_value (const unsigned int i) const
+{
+ Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state));
+ AssertIndexRange(i,wr.size());
+
+ return wr[i];
+}
+
+
DEAL_II_NAMESPACE_CLOSE
return *this;
}
-#ifdef HAVE_LIBBLAS
template <typename number>
void
Vector<number> &w,
const Vector<number> &v,
const bool adding) const
-{
- Assert ((state == matrix) || (state == inverse_matrix),
- ExcState(state));
-
+{
const int mm = this->n_rows();
const int nn = this->n_cols();
const number alpha = 1.;
const number beta = (adding ? 1. : 0.);
-
- gemv("N", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one);
+ const number null = 0.;
+
+ switch (state)
+ {
+ case matrix:
+ case inverse_matrix:
+ {
+ AssertDimension(v.size(), this->n_cols());
+ AssertDimension(w.size(), this->n_rows());
+
+ gemv("N", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one);
+ break;
+ }
+ case svd:
+ {
+ AssertDimension(v.size(), this->n_cols());
+ AssertDimension(w.size(), this->n_rows());
+ // Compute V^T v
+ work.resize(std::max(mm,nn));
+ gemv("N", &nn, &nn, &alpha, svd_vt->data(), &nn, v.val, &one, &null, &work[0], &one);
+ // Multiply by singular values
+ for (unsigned int i=0;i<wr.size();++i)
+ work[i] *= wr[i];
+ // Multiply with U
+ gemv("N", &mm, &mm, &alpha, svd_u->data(), &mm, &work[0], &one, &beta, w.val, &one);
+ break;
+ }
+ case inverse_svd:
+ {
+ AssertDimension(w.size(), this->n_cols());
+ AssertDimension(v.size(), this->n_rows());
+ // Compute U^T v
+ work.resize(std::max(mm,nn));
+ gemv("T", &mm, &mm, &alpha, svd_u->data(), &mm, v.val, &one, &null, &work[0], &one);
+ // Multiply by singular values
+ for (unsigned int i=0;i<wr.size();++i)
+ work[i] *= wr[i];
+ // Multiply with V
+ gemv("T", &nn, &nn, &alpha, svd_vt->data(), &nn, &work[0], &one, &beta, w.val, &one);
+ break;
+ }
+ default:
+ Assert (false, ExcState(state));
+ }
}
const Vector<number> &v,
const bool adding) const
{
- Assert (state == matrix, ExcState(state));
-
const int mm = this->n_rows();
const int nn = this->n_cols();
const number alpha = 1.;
const number beta = (adding ? 1. : 0.);
-
- gemv("T", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one);
+ const number null = 0.;
+
+ switch (state)
+ {
+ case matrix:
+ case inverse_matrix:
+ {
+ AssertDimension(w.size(), this->n_cols());
+ AssertDimension(v.size(), this->n_rows());
+
+ gemv("T", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one);
+ break;
+ }
+ case svd:
+ {
+ AssertDimension(w.size(), this->n_cols());
+ AssertDimension(v.size(), this->n_rows());
+
+ // Compute U^T v
+ work.resize(std::max(mm,nn));
+ gemv("T", &mm, &mm, &alpha, svd_u->data(), &mm, v.val, &one, &null, &work[0], &one);
+ // Multiply by singular values
+ for (unsigned int i=0;i<wr.size();++i)
+ work[i] *= wr[i];
+ // Multiply with V
+ gemv("T", &nn, &nn, &alpha, svd_vt->data(), &nn, &work[0], &one, &beta, w.val, &one);
+ break;
+ case inverse_svd:
+ {
+ AssertDimension(v.size(), this->n_cols());
+ AssertDimension(w.size(), this->n_rows());
+
+ // Compute V^T v
+ work.resize(std::max(mm,nn));
+ gemv("N", &nn, &nn, &alpha, svd_vt->data(), &nn, v.val, &one, &null, &work[0], &one);
+ // Multiply by singular values
+ for (unsigned int i=0;i<wr.size();++i)
+ work[i] *= wr[i];
+ // Multiply with U
+ gemv("N", &mm, &mm, &alpha, svd_u->data(), &mm, &work[0], &one, &beta, w.val, &one);
+ break;
+ }
+ }
+ default:
+ Assert (false, ExcState(state));
+ }
}
-#else
-
template <typename number>
void
-LAPACKFullMatrix<number>::vmult (
- Vector<number> &,
- const Vector<number> &,
- const bool ) const
+LAPACKFullMatrix<number>::compute_lu_factorization()
{
- Assert(false, ExcNeedsBLAS());
-}
+ Assert(state == matrix, ExcState(state));
+ const int mm = this->n_rows();
+ const int nn = this->n_cols();
+ number* values = const_cast<number*> (this->data());
+ ipiv.resize(mm);
+ int info;
+ getrf(&mm, &nn, values, &mm, &ipiv[0], &info);
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
-template <typename number>
-void
-LAPACKFullMatrix<number>::Tvmult (
- Vector<number> &,
- const Vector<number> &,
- const bool ) const
-{
- Assert(false, ExcNeedsBLAS());
+ state = lu;
}
-#endif
-
-#ifdef HAVE_LIBLAPACK
-
template <typename number>
void
-LAPACKFullMatrix<number>::compute_lu_factorization()
+LAPACKFullMatrix<number>::compute_svd()
{
Assert(state == matrix, ExcState(state));
+ state = LAPACKSupport::unusable;
+
const int mm = this->n_rows();
const int nn = this->n_cols();
number* values = const_cast<number*> (this->data());
- ipiv.resize(mm);
+ wr.resize(std::max(mm,nn));
+ std::fill(wr.begin(), wr.end(), 0.);
+ ipiv.resize(8*mm);
+
+ svd_u = boost::shared_ptr<LAPACKFullMatrix<number> >(new LAPACKFullMatrix<number>(mm,mm));
+ svd_vt = boost::shared_ptr<LAPACKFullMatrix<number> >(new LAPACKFullMatrix<number>(nn,nn));
+ number* mu = const_cast<number*> (svd_u->data());
+ number* mvt = const_cast<number*> (svd_vt->data());
int info;
- getrf(&mm, &nn, values, &mm, &ipiv[0], &info);
- Assert(info >= 0, ExcInternalError());
- Assert(info == 0, LACExceptions::ExcSingular());
+ // see comment on this #if
+ // below. Another reason to love Petsc
+#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE
- state = lu;
+ // First determine optimal
+ // workspace size
+ work.resize(1);
+ int lwork = -1;
+ gesdd(&LAPACKSupport::A, &mm, &nn, values, &mm,
+ &wr[0], mu, &mm, mvt, &nn,
+ &work[0], &lwork, &ipiv[0], &info);
+ Assert (info==0, LAPACKSupport::ExcErrorCode("gesdd", info));
+ // Now resize work array and
+ lwork = work[0] + .5;
+#else
+ int lwork = 3*std::min(mm,nn) +
+ std::max(std::max(mm,nn),4*std::min(mm,nn)*std::min(mm,nn)+4*std::min(mm,nn));
+#endif
+ work.resize(lwork);
+ // Do the actual SVD.
+ gesdd(&LAPACKSupport::A, &mm, &nn, values, &mm,
+ &wr[0], mu, &mm, mvt, &nn,
+ &work[0], &lwork, &ipiv[0], &info);
+ Assert (info==0, LAPACKSupport::ExcErrorCode("gesdd", info));
+
+ work.resize(0);
+ ipiv.resize(0);
+
+ state = LAPACKSupport::svd;
+}
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_inverse_svd(const double threshold)
+{
+ if (state == LAPACKSupport::matrix)
+ compute_svd();
+
+ Assert (state==LAPACKSupport::svd, ExcState(state));
+
+ const double lim = wr[0]*threshold;
+ for (unsigned int i=0;i<wr.size();++i)
+ {
+ if (wr[i] > lim)
+ wr[i] = 1./wr[i];
+ else
+ wr[i] = 0.;
+ }
+ state = LAPACKSupport::inverse_svd;
}
state = LAPACKSupport::State(eigenvalues | unusable);
}
-#else
-
-template <typename number>
-void
-LAPACKFullMatrix<number>::compute_lu_factorization()
-{
-Assert(false, ExcNeedsLAPACK());
-
-}
-
-
-template <typename number>
-void
-LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>&, bool) const
-{
- Assert(false, ExcNeedsLAPACK());
-}
-
-
-template <typename number>
-void
-LAPACKFullMatrix<number>::compute_eigenvalues(const bool /*right*/,
- const bool /*left*/)
-{
- Assert(false, ExcNeedsLAPACK());
-}
-
-
-#endif
-
-
-// template <typename number>
-// LAPACKFullMatrix<number>::()
-// {}
-
// template <typename number>
// LAPACKFullMatrix<number>::()
projects / dealii.git / commitdiff