// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* where all allocated objects
* should have been released. Since
* this exception is thrown, some
- * where still allocated.
+ * were still allocated.
*/
DeclException1 (ExcMemoryLeak, int,
<< "Destroying memory handler while " << arg1
<< " objects are still allocated");
/**
- * An error occured reading or
+ * An error occurred reading or
* writing a file.
*/
DeclException0 (ExcIO);
/**
- * An error occured opening the named file.
+ * An error occurred opening the named file.
*
* The constructor takes a single
* argument of type <tt>char*</tt>
<h3>lac</h3>
<ol>
+ <li>
+ <p>
+ New: There are new functions FullMatrix::cholesky and
+ FullMatrix::outer_product. FullMatrix::cholesky finds the Cholesky
+ decomposition of a matrix in lower triangular form.
+ FullMatrix::outer_product calculates <tt>*this</tt> $= VW^T$ where $V$
+ and $W$ are vectors.
+ <br>
+ (Jean Marie Linhart 2009/07/27)
+ </p>
+ </li>
+
+ <li>
+ <p>
+ Fixed: The TrilinosWrappers::MPI::BlockVector class declares an assignment
+ operator from the non-Trilinos BlockVector class but it could not be
+ compiled due to an oversight. This is now fixed.
+ <br>
+ (WB 2009/06/29)
+ </p>
+ </li>
+
<li>
<p>
New: Based on work by Francisco Alvaro, the existing SLEPcWrappers now
</p>
</li>
- <li>
- <p>
- Fixed: The TrilinosWrappers::MPI::BlockVector class declares an assignment
- operator from the non-Trilinos BlockVector class but it could not be
- compiled due to an oversight. This is now fixed.
- <br>
- (WB 2009/06/29)
- </p>
- </li>
-
<li>
<p>
Fixed: The TrilinosWrappers::BlockVector class declares an assignment
const FullMatrix<number2> &B);
/**
- * Add transose of a rectangular block.
+ * Add transpose of a rectangular block.
*
* A rectangular block of the
* matrix <tt>src</tt> is
template <typename number2>
void invert (const FullMatrix<number2> &M);
+ /**
+ * Assign the Cholesky decomposition
+ * of the given matrix to <tt>*this</tt>.
+ * The given matrix must be symmetric
+ * positive definite.
+ *
+ * ExcMatrixNotPositiveDefinite
+ * will be thrown in the case that the
+ * matrix is not positive definite.
+ */
+ template <typename number2>
+ void cholesky (const FullMatrix<number2> &A);
+
+ /**
+ * <tt>*this</tt> = $V W^T$ where $V,W$
+ * are vectors of the same length.
+ */
+ template <typename number2>
+ void outer_product (const Vector<number2> &V,
+ const Vector<number2> &W);
+
/**
* Assign the left_inverse of the given matrix
* to <tt>*this</tt>. The calculation being
* Exception
*/
DeclException0 (ExcSourceEqualsDestination);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcMatrixNotPositiveDefinite);
//@}
friend class Accessor;
break;
}
+
default:
// if no inversion is
};
}
+
+template <typename number>
+template <typename number2>
+void
+FullMatrix<number>::cholesky (const FullMatrix<number2> &A)
+{
+ Assert (!A.empty(), ExcEmptyMatrix());
+ Assert (A.n() == A.m(),
+ ExcNotQuadratic());
+ // Matrix must be symmetric.
+ Assert(A.relative_symmetry_norm2() < 1.0e-10, ExcMessage("A must be symmetric."));
+
+ if (PointerComparison::equal(&A, this))
+ {
+ // avoid overwriting source
+ // by destination matrix:
+ FullMatrix<number2> A2 = A;
+ cholesky(A2);
+ }
+ else
+ {
+ /* reinit *this to 0 */
+ this->reinit(A.m(), A.n());
+
+ double SLik2 = 0.0, SLikLjk = 0.0;
+ for (unsigned int i=0; i< this->n_cols(); i++){
+ SLik2 = 0.0;
+ for (unsigned int j = 0; j < i; j++){
+ SLik2 += (*this)(i,j)*(*this)(i,j);
+ SLikLjk = 0.0;
+ for (unsigned int k =0; k<j; k++)
+ {
+ SLikLjk = (*this)(i,k)*(*this)(j,k);
+ };
+ (*this)(i,j) = (1./(*this)(j,j))*(A(i,j) - SLikLjk);
+ }
+ AssertThrow (A(i,i) - SLik2 >= 0,
+ ExcMatrixNotPositiveDefinite());
+
+ (*this)(i,i) = std::sqrt(A(i,i) - SLik2);
+ }
+ }
+ }
+
+
+template <typename number>
+template <typename number2>
+void
+FullMatrix<number>::outer_product (const Vector<number2> &V,
+ const Vector<number2> &W)
+{
+ Assert (V.size() == W.size(), ExcMessage("Vectors V, W must be the same size."));
+ this->reinit(V.size(), V.size());
+
+ for(unsigned int i = 0; i<this->n(); i++)
+ {
+ for(unsigned int j = 0; j< this->n(); j++)
+ {
+ (*this)(i,j) = V(i)*W(j);
+ }
+ }
+}
+
+
template <typename number>
template <typename number2>
void
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template void FullMatrix<S>::copy_from<1>(
Tensor<2,1>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned);
-
+
template void FullMatrix<S>::copy_from<2>(
Tensor<2,2>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned);
-
+
template void FullMatrix<S>::copy_from<3>(
Tensor<2,3>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned);
-
+
template void FullMatrix<S>::copy_to<1>(
Tensor<2,1>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned);
template
void FullMatrix<S1>::precondition_Jacobi<S2> (
Vector<S2> &, const Vector<S2> &, const S1) const;
+
+ template
+ void FullMatrix<S1>::cholesky<S2> (const FullMatrix<S2>&);
+
+ template
+ void FullMatrix<S1>::outer_product<S2> (const Vector<S2>&,
+ const Vector<S2>&);
}
for (S1, S2, S3 : REAL_SCALARS)