]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
LU decomposition added to LAPACKFullMatrix
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Oct 2006 23:01:21 +0000 (23:01 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Oct 2006 23:01:21 +0000 (23:01 +0000)
some items in LAPACKFullMatrix fixed

git-svn-id: https://svn.dealii.org/trunk@13973 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/config.h.in
deal.II/configure
deal.II/configure.in
deal.II/contrib/blastemplates/lapack_templates.h.in
deal.II/lac/include/lac/exceptions.h
deal.II/lac/include/lac/lapack_full_matrix.h
deal.II/lac/include/lac/lapack_templates.h
deal.II/lac/source/lapack_full_matrix.cc

index e1016ea246dd1e9a104097e5729d776fbb82feb3..0f7d6ba18d23191bc095976b87817c58f37a2233 100644 (file)
 /* Define to 1 if you have the `dgesvd_' function. */
 #undef HAVE_DGESVD_
 
+/* Define to 1 if you have the `dgetrf_' function. */
+#undef HAVE_DGETRF_
+
+/* Define to 1 if you have the `dgetrs_' function. */
+#undef HAVE_DGETRS_
+
 /* Define to 1 if you have the `dstev_' function. */
 #undef HAVE_DSTEV_
 
 /* Define to 1 if you have the `sgesvd_' function. */
 #undef HAVE_SGESVD_
 
+/* Define to 1 if you have the `sgetrf_' function. */
+#undef HAVE_SGETRF_
+
+/* Define to 1 if you have the `sgetrs_' function. */
+#undef HAVE_SGETRS_
+
 /* Define to 1 if you have the `sstev_' function. */
 #undef HAVE_SSTEV_
 
index b8492dedbea66f2ffb6cd7ccafa725d0efc2b9bf..82b3612a3a8b330bf78e4dbf335484af5a303e4b 100755 (executable)
@@ -12067,9 +12067,114 @@ fi;
 
 
 
+for ac_func in dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_
+do
+as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
+echo "$as_me:$LINENO: checking for $ac_func" >&5
+echo $ECHO_N "checking for $ac_func... $ECHO_C" >&6
+if eval "test \"\${$as_ac_var+set}\" = set"; then
+  echo $ECHO_N "(cached) $ECHO_C" >&6
+else
+  cat >conftest.$ac_ext <<_ACEOF
+/* confdefs.h.  */
+_ACEOF
+cat confdefs.h >>conftest.$ac_ext
+cat >>conftest.$ac_ext <<_ACEOF
+/* end confdefs.h.  */
+/* Define $ac_func to an innocuous variant, in case <limits.h> declares $ac_func.
+   For example, HP-UX 11i <limits.h> declares gettimeofday.  */
+#define $ac_func innocuous_$ac_func
+
+/* System header to define __stub macros and hopefully few prototypes,
+    which can conflict with char $ac_func (); below.
+    Prefer <limits.h> to <assert.h> if __STDC__ is defined, since
+    <limits.h> exists even on freestanding compilers.  */
+
+#ifdef __STDC__
+# include <limits.h>
+#else
+# include <assert.h>
+#endif
+
+#undef $ac_func
+
+/* Override any gcc2 internal prototype to avoid an error.  */
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+/* We use char because int might match the return type of a gcc2
+   builtin and then its argument prototype would still apply.  */
+char $ac_func ();
+/* The GNU C library defines this for functions which it implements
+    to always fail with ENOSYS.  Some functions are actually named
+    something starting with __ and the normal name is an alias.  */
+#if defined (__stub_$ac_func) || defined (__stub___$ac_func)
+choke me
+#else
+char (*f) () = $ac_func;
+#endif
+#ifdef __cplusplus
+}
+#endif
+
+int
+main ()
+{
+return f != $ac_func;
+  ;
+  return 0;
+}
+_ACEOF
+rm -f conftest.$ac_objext conftest$ac_exeext
+if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
+  (eval $ac_link) 2>conftest.er1
+  ac_status=$?
+  grep -v '^ *+' conftest.er1 >conftest.err
+  rm -f conftest.er1
+  cat conftest.err >&5
+  echo "$as_me:$LINENO: \$? = $ac_status" >&5
+  (exit $ac_status); } &&
+        { ac_try='test -z "$ac_cxx_werror_flag"
+                        || test ! -s conftest.err'
+  { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
+  (eval $ac_try) 2>&5
+  ac_status=$?
+  echo "$as_me:$LINENO: \$? = $ac_status" >&5
+  (exit $ac_status); }; } &&
+        { ac_try='test -s conftest$ac_exeext'
+  { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
+  (eval $ac_try) 2>&5
+  ac_status=$?
+  echo "$as_me:$LINENO: \$? = $ac_status" >&5
+  (exit $ac_status); }; }; then
+  eval "$as_ac_var=yes"
+else
+  echo "$as_me: failed program was:" >&5
+sed 's/^/| /' conftest.$ac_ext >&5
+
+eval "$as_ac_var=no"
+fi
+rm -f conftest.err conftest.$ac_objext \
+      conftest$ac_exeext conftest.$ac_ext
+fi
+echo "$as_me:$LINENO: result: `eval echo '${'$as_ac_var'}'`" >&5
+echo "${ECHO_T}`eval echo '${'$as_ac_var'}'`" >&6
+if test `eval echo '${'$as_ac_var'}'` = yes; then
+  cat >>confdefs.h <<_ACEOF
+#define `echo "HAVE_$ac_func" | $as_tr_cpp` 1
+_ACEOF
+
+fi
+done
+
+
+
+
+
 
 
-for ac_func in dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ dstev_ sstev_
+for ac_func in dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_
 do
 as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
 echo "$as_me:$LINENO: checking for $ac_func" >&5
index ce8c0388234dc4b4ab673ee70256cc09a7bbc84a..1ee4c2c5cecb84b8c2c08f987e975070fe4b5361 100644 (file)
@@ -510,7 +510,8 @@ AC_ARG_WITH(lapack,
                           argument is given, use -llapack. Default is not to
                           use -llapack.],
   DEAL_II_WITH_LAPACK($withval))
-AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ dstev_ sstev_])
+AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_])
+AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_])
 
 dnl -------------------------------------------------------------
 dnl       set include paths of several libraries
index 1b23d0d64f443da55ecdd9b4d88e87be6f822dc9..89c8578f2b008806743debbd0e006359761b8dc8 100644 (file)
@@ -4,6 +4,15 @@ void dgemv_ (const char* trans, const int* m, const int* n,
             const double* alpha, const double* A, const int* lda,
             const double* x, const int* incx,
             const double* b, double* y, const int* incy);
+
+// Compute LU factorization
+void dgetrf_ (const int* m, const int* n, double* A,
+             const int* lda, int* ipiv, int* info);
+// Apply forward/backward substitution to LU factorization
+void dgetrs_ (const char* trans, const int* n, const int* nrhs,
+             const double* A, const int* lda, const int* ipiv,
+             double* b, const int* ldb, int* info);
+
 // Compute eigenvalues and vectors
 void dgeev_ (const char* jobvl, const char* jobvr,
             const int* n, double* A, const int* lda,
index 0e63732d8c7958ff3025e8a24e0cb0c7f6298cd4..d286ad19fe67707dfff2cfd168117af0ae139a9c 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -28,11 +28,17 @@ namespace LACExceptions
                                    */
   DeclException0 (ExcNotQuadratic);
 
+                                  /**
+                                   * The operation cannot be finished
+                                   * since the matrix is singular.
+                                   */
+  DeclException0 (ExcSingular);
+  
                                   /**
                                    * Block indices of two block
                                    * objects are different.
                                    */
-  DeclException0(ExcDifferentBlockIndices);
+  DeclException0 (ExcDifferentBlockIndices);
   
                                   /**
                                    * An error of a PETSc function was
index 945fe88b714d78ea1c53c99c0239cd89278bb6c1..6f00e2b46e0937fbadb09aea7f6eb0519f844858 100644 (file)
 
 
 #include <base/config.h>
+#include <base/smartpointer.h>
 #include <base/table.h>
 #include <lac/lapack_support.h>
+#include <lac/vector_memory.h>
 
 #include <vector>
 #include <complex>
 
 // forward declarations
 template<typename number> class Vector;
+template<typename number> class BlockVector;
 template<typename number> class FullMatrix;
-
-/*! @addtogroup Matrix1
- *@{
- */
+template<class VECTOR> class VectorMemory;
 
 
 /**
- * A variant of FullMatrix using LAPACK functions whereever possible.
+ * A variant of FullMatrix using LAPACK functions whereever
+ * possible. In order to do this, the matrix is stored in transposed
+ * order. The element access functions hide this fact by reverting the
+ * transposition.
  *
+ * @note In order to perform LAPACK functions, the class contains a lot of
+ * auxiliary data in the private section. The names of these data
+ * vectors are usually the names chosen for the arguments in the
+ * LAPACK documentation.
+ *
+ * @ingroup Matrix1
  * @author Guido Kanschat, 2005
  */
 template <typename number>
@@ -136,13 +145,19 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      * maximum size possible,
                                      * determined either by the size
                                      * of <tt>this</tt> or <tt>src</tt>.
+                                     *
+                                     * The final two arguments allow
+                                     * to enter a multiple of the
+                                     * source or its transpose.
                                      */
     template<class MATRIX>
     void fill (const MATRIX &src,
               const unsigned int dst_offset_i = 0,
               const unsigned int dst_offset_j = 0,
               const unsigned int src_offset_i = 0,
-              const unsigned int src_offset_j = 0);
+              const unsigned int src_offset_j = 0,
+              const number factor = 1.,
+              const bool transpose = false);
     
                                     /**
                                      * Matrix-vector-multiplication.
@@ -208,6 +223,25 @@ class LAPACKFullMatrix : public TransposeTable<number>
     void Tvmult_add (Vector<number>       &w,
                     const Vector<number> &v) const;
 
+                                    /**
+                                     * Compute the LU factorization
+                                     * of the matrix using LAPACK
+                                     * function Xgetrf.
+                                     */
+    void compute_lu_factorization ();
+
+                                    /**
+                                     * Solve the linear system with
+                                     * right hand side given by
+                                     * applying forward/backward
+                                     * substitution to the previously
+                                     * computed LU
+                                     * factorization. Uses LAPACK
+                                     * function Xgetrs.
+                                     */
+    void apply_lu_factorization (Vector<number>& v,
+                                const bool transposed) const;
+    
                                     /**
                                      * Compute eigenvalues of the
                                      * matrix. After this routine has
@@ -326,6 +360,14 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      * some LAPACK functions.
                                      */
     mutable std::vector<number> work;
+
+                                    /**
+                                     * The vector storing the
+                                     * permutations applied for
+                                     * pivoting in the
+                                     * LU-factorization.
+                                     */
+    std::vector<int> ipiv;
     
                                     /**
                                      * Real parts of
@@ -354,7 +396,35 @@ class LAPACKFullMatrix : public TransposeTable<number>
     std::vector<number> vr;
 };
 
-/**@}*/
+
+
+/**
+ * A preconditioner based on the LU-factorization of LAPACKFullMatrix.
+ *
+ * @ingroup Preconditioners
+ * @author Guido Kanschat, 2006
+ */
+template <typename number>
+class LUPrecondition
+  :
+  public Subscriptor
+{
+  public:
+    void initialize(const LAPACKFullMatrix<number>&);
+    void initialize(const LAPACKFullMatrix<number>&,
+                   VectorMemory<Vector<number> >&);
+    void vmult(Vector<number>&, const Vector<number>&) const;
+    void Tvmult(Vector<number>&, const Vector<number>&) const;
+    void vmult(BlockVector<number>&,
+              const BlockVector<number>&) const;
+    void Tvmult(BlockVector<number>&,
+               const BlockVector<number>&) const;
+  private:
+//    SmartPointer<const LAPACKFullMatrix<number> > matrix;
+//    SmartPointer<VectorMemory<Vector<number> > mem;
+};
+
+
 
 template <typename number>
 template <class MATRIX>
@@ -380,19 +450,21 @@ LAPACKFullMatrix<number>::fill (
   const unsigned int dst_offset_i,
   const unsigned int dst_offset_j,
   const unsigned int src_offset_i,
-  const unsigned int src_offset_j)
+  const unsigned int src_offset_j,
+  const number factor,
+  const bool transpose)
 {
   const typename MATRIX::const_iterator end = M.end();
   for (typename MATRIX::const_iterator entry = M.begin(src_offset_i);
        entry != end; ++entry)
     {
-      const unsigned int i = entry->row();
-      const unsigned int j = entry->column();
-
+      const unsigned int i = transpose ? entry->column() : entry->row();
+      const unsigned int j = transpose ? entry->row() : entry->column();
+      
       const unsigned int dst_i=dst_offset_i+i-src_offset_i;
       const unsigned int dst_j=dst_offset_j+j-src_offset_j;
       if (dst_i<this->n_rows() && dst_j<this->n_cols())
-       (*this)(dst_i, dst_j) = entry->value();
+       (*this)(dst_i, dst_j) = factor * entry->value();
     }
   
   state = LAPACKSupport::matrix;
index d6803b9eb01198613ee5bbae55bb316f7f564a55..3d772cbbc17f554be82904af6cd562074333c26a 100644 (file)
@@ -31,6 +31,18 @@ void sgemv_ (const char* trans, const int* m, const int* n,
             const float* alpha, const float* A, const int* lda,
             const float* x, const int* incx,
             const float* b, float* y, const int* incy);
+// Compute LU factorization
+void dgetrf_ (const int* m, const int* n, double* A,
+             const int* lda, int* ipiv, int* info);
+void sgetrf_ (const int* m, const int* n, float* A,
+             const int* lda, int* ipiv, int* info);
+// Apply forward/backward substitution to LU factorization
+void dgetrs_ (const char* trans, const int* n, const int* nrhs,
+             const double* A, const int* lda, const int* ipiv,
+             double* b, const int* ldb, int* info);
+void sgetrs_ (const char* trans, const int* n, const int* nrhs,
+             const float* A, const int* lda, const int* ipiv,
+             float* b, const int* ldb, int* info);
 // Compute eigenvalues and vectors
 void dgeev_ (const char* jobvl, const char* jobvr,
             const int* n, double* A, const int* lda,
@@ -121,13 +133,73 @@ gemv (const char* trans, const int* m, const int* n, const float* alpha, const f
 }
 #else
 inline void
-gemv (const char*, const int*, const int*, const float*, const float*, const int*, const float*, const int*, const float*, float*, const int*)
+gemv (const char*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*)
 {
   LAPACKSupport::ExcMissing("sgemv");
 }
 #endif
 
 
+#ifdef HAVE_DGETRF_
+inline void
+getrf (const int* m, const int* n, double* A, const int* lda, int* ipiv, int* info)
+{
+  dgetrf_ (m,n,A,lda,ipiv,info);
+}
+#else
+inline void
+getrf (const int*, const int*, double*, const int*, int*, int*)
+{
+  LAPACKSupport::ExcMissing("dgetrf");
+}
+#endif
+
+
+#ifdef HAVE_SGETRF_
+inline void
+getrf (const int* m, const int* n, float* A, const int* lda, int* ipiv, int* info)
+{
+  sgetrf_ (m,n,A,lda,ipiv,info);
+}
+#else
+inline void
+getrf (const int*, const int*, double*, const int*, int*, int*)
+{
+  LAPACKSupport::ExcMissing("sgetrf");
+}
+#endif
+
+
+#ifdef HAVE_DGETRS_
+inline void
+getrs (const char* trans, const int* n, const int* nrhs, const double* A, const int* lda, const int* ipiv, double* b, const int* ldb, int* info)
+{
+  dgetrs_ (trans,n,nrhs,A,lda,ipiv,b,ldb,info);
+}
+#else
+inline void
+getrs (const char*, const int*, const int*, const double*, const int*, const int*, double*, const int*, int*)
+{
+  LAPACKSupport::ExcMissing("dgetrs");
+}
+#endif
+
+
+#ifdef HAVE_SGETRS_
+inline void
+getrs (const char* trans, const int* n, const int* nrhs, const float* A, const int* lda, const int* ipiv, float* b, const int* ldb, int* info)
+{
+  sgetrs_ (trans,n,nrhs,A,lda,ipiv,b,ldb,info);
+}
+#else
+inline void
+getrs (const char*, const int*, const int*, const double*, const int*, const int*, double*, const int*, int*)
+{
+  LAPACKSupport::ExcMissing("sgetrs");
+}
+#endif
+
+
 #ifdef HAVE_DGEEV_
 inline void
 geev (const char* jobvl, const char* jobvr, const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, double* vl, const int* ldvl, double* vr, const int* ldva, double* work, const int* lwork, int* info)
@@ -151,7 +223,7 @@ geev (const char* jobvl, const char* jobvr, const int* n, float* A, const int* l
 }
 #else
 inline void
-geev (const char*, const char*, const int*, float*, const int*, float*, float*, float*, const int*, float*, const int*, float*, const int*, int*)
+geev (const char*, const char*, const int*, double*, const int*, double*, double*, double*, const int*, double*, const int*, double*, const int*, int*)
 {
   LAPACKSupport::ExcMissing("sgeev");
 }
@@ -181,7 +253,7 @@ geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sen
 }
 #else
 inline void
-geevx (const char*, const char*, const char*, const char*, const int*, float*, const int*, float*, float*, float*, const int*, float*, const int*, int*, int*, float*, float*, float*, float*, float*, const int*, int*, int*)
+geevx (const char*, const char*, const char*, const char*, const int*, double*, const int*, double*, double*, double*, const int*, double*, const int*, int*, int*, double*, double*, double*, double*, double*, const int*, int*, int*)
 {
   LAPACKSupport::ExcMissing("sgeevx");
 }
@@ -211,7 +283,7 @@ gesvd (int* jobu, int* jobvt, const int* n, const int* m, float* A, const int* l
 }
 #else
 inline void
-gesvd (int*, int*, const int*, const int*, float*, const int*, float*, float*, const int*, float*, const int*, float*, const int*, int*)
+gesvd (int*, int*, const int*, const int*, double*, const int*, double*, double*, const int*, double*, const int*, double*, const int*, int*)
 {
   LAPACKSupport::ExcMissing("sgesvd");
 }
@@ -241,7 +313,7 @@ stev (const char* jobz, const int* n, float* d, float* e, float* z, const int* l
 }
 #else
 inline void
-stev (const char*, const int*, float*, float*, float*, const int*, float*, int*)
+stev (const char*, const int*, double*, double*, double*, const int*, double*, int*)
 {
   LAPACKSupport::ExcMissing("sstev");
 }
index 599bd43d55d6a6c1d14617fc095f5040db8aa2f8..e65f1440cb8479387f5b8dee704ce4e68d0807c0 100644 (file)
@@ -25,7 +25,8 @@ using namespace LAPACKSupport;
 template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix(const unsigned int n)
                :
-               TransposeTable<number> (n,n)
+               TransposeTable<number> (n,n),
+               state(matrix)
 {}
 
 
@@ -34,17 +35,17 @@ LAPACKFullMatrix<number>::LAPACKFullMatrix(
   const unsigned int m,
   const unsigned int n)
                :
-               TransposeTable<number> (m,n)
+               TransposeTable<number> (m,n),
+               state(matrix)
 {}
 
 
 template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
                :
-               TransposeTable<number> (M)
-{
-  state = LAPACKSupport::matrix;
-}
+               TransposeTable<number> (M),
+               state(matrix)
+{}
 
 
 template <typename number>
@@ -95,7 +96,7 @@ LAPACKFullMatrix<number>::vmult (
   const Vector<number> &v,
   const bool            adding) const
 {
-  Assert (state == matrix, ExcInvalidState());
+  Assert (state == matrix, ExcState(state));
   
   const int mm = this->n_rows();
   const int nn = this->n_cols();
@@ -113,7 +114,7 @@ LAPACKFullMatrix<number>::Tvmult (
   const Vector<number> &v,
   const bool            adding) const
 {
-  Assert (state == matrix, ExcInvalidState());
+  Assert (state == matrix, ExcState(state));
   
   const int mm = this->n_rows();
   const int nn = this->n_cols();
@@ -152,12 +153,53 @@ LAPACKFullMatrix<number>::Tvmult (
 
 #ifdef HAVE_LIBLAPACK
 
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_lu_factorization()
+{
+  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());
+  
+  state = lu;
+}
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>& v,
+                                                const bool transposed) const
+{
+  Assert(state == lu, ExcState(state));
+  Assert(this->n_rows() == this->n_cols(),
+        LACExceptions::ExcNotQuadratic());
+  
+  const char* trans = transposed ? &T : &N;
+  const int nn = this->n_cols();
+  const number* values = this->data();
+  int info;
+
+  getrs(trans, &nn, &one, values, &nn, &ipiv[0],
+       v.begin(), &nn, &info);
+
+  Assert(info == 0, ExcInternalError());
+}
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_eigenvalues(
   const bool right,
   const bool left)
 {
+  Assert(state == matrix, ExcState(state));
   const int nn = this->n_cols();
   wr.resize(nn);
   wi.resize(nn);
@@ -232,6 +274,21 @@ LAPACKFullMatrix<number>::compute_eigenvalues(
 
 #else
 
+template <typename number>
+LAPACKFullMatrix<number>::compute_lu_factorization()
+{
+Assert(false, ExcNeedsLAPACK());
+  
+}
+
+
+template <typename number>
+LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>&, bool)
+{
+  Assert(false, ExcNeedsLAPACK());
+}
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_eigenvalues(const bool /*right*/,
@@ -243,6 +300,7 @@ LAPACKFullMatrix<number>::compute_eigenvalues(const bool /*right*/,
 
 #endif
 
+
 // template <typename number>
 // LAPACKFullMatrix<number>::()
 // {}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.