]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Initial revision
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 4 Mar 1998 16:59:32 +0000 (16:59 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 4 Mar 1998 16:59:32 +0000 (16:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@11 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/types.h [new file with mode: 0644]
deal.II/lac/include/lac/dfmatrix.h [new file with mode: 0644]
deal.II/lac/include/lac/dsmatrix.h [new file with mode: 0644]
deal.II/lac/include/lac/dvector.h [new file with mode: 0644]
deal.II/lac/include/lac/ivector.h [new file with mode: 0644]
deal.II/lac/include/lac/vectorbase.h [new file with mode: 0644]
deal.II/lac/source/dfmatrix.cc [new file with mode: 0644]
deal.II/lac/source/dsmatrix.cc [new file with mode: 0644]
deal.II/lac/source/dvector.cc [new file with mode: 0644]
deal.II/lac/source/ivector.cc [new file with mode: 0644]

diff --git a/deal.II/base/include/base/types.h b/deal.II/base/include/base/types.h
new file mode 100644 (file)
index 0000000..fe77a49
--- /dev/null
@@ -0,0 +1,37 @@
+// $Id$
+
+#ifndef __base_types_h
+#define __base_types_h
+
+#ifndef _stdlib_h
+#include <stdlib.h>
+#endif
+
+/// Maximum Number of boundary conditions (divided by 8 )
+
+#define BOUNDARY_CONDITION_BYTES 1
+
+/// Maximum number of multi purpose flags per cell ( / 8 )
+
+#define CELL_FLAG_BYTES 1
+
+/// Maximum number of multi purpose flags per vertex ( / 8 )
+
+#define VERTEX_FLAG_BYTES 1
+
+#define UNSIGNED unsigned long
+
+typedef double COORD;
+#define DOUBLECOORD
+
+typedef int SHORT;
+
+#define ABS(x) ( ((x)>=0) ? (x) : -(x) )
+#define MAX(a,b) ( ((a)>(b)) ? (a) : (b) )
+#define MIN(a,b) ( ((a)<(b)) ? (a) : (b) )
+#define SIGN(a) ( ((a)<0) ? (-1) : 1 )
+#ifndef NULL
+#define NULL 0
+#endif
+
+#endif
diff --git a/deal.II/lac/include/lac/dfmatrix.h b/deal.II/lac/include/lac/dfmatrix.h
new file mode 100644 (file)
index 0000000..1682fa3
--- /dev/null
@@ -0,0 +1,158 @@
+// $Id$
+
+#ifndef __lac_dfmatrix_h
+#define __lac_dfmatrix_h
+
+#ifndef __stdio_h
+#include <stdio.h>
+#endif
+#ifndef __lac_dvector_h
+#include <lac/dvector.h>
+#endif
+#ifndef __lac_ivector_h
+#include <lac/ivector.h>
+#endif
+
+/// Double precision full matrix
+class dFMatrix
+{
+  double* val;
+  int dim_range, dim_image, val_size;
+  void init(int n, int m);
+  dFMatrix(const dFMatrix&);
+
+  double& el(int i, int j)  { return val[i*dim_range+j]; }
+  double el(int i, int j) const { return val[i*dim_range+j]; }
+  double el(int i) const { return val[i]; }
+ public:
+  /// Number of rows
+  int m() const { return dim_image; }
+  /// Number of columns
+  int n() const { return dim_range; }
+
+
+  /// copy constructor  
+  dFMatrix(const dFMatrix&);
+  /// Constructor for quadratic n x n matrices
+  dFMatrix(int n = 1) { init(n,n); }
+  /// Constructor for rectangular matrices with m rows and n columns.
+  dFMatrix(int m,int n) { init(m,n); }
+  /// Destructor.
+  ~dFMatrix();
+
+  /// Reinitialization of rectangular matrix.
+  void reinit(int m, int n);
+  /// Reinitialization of quadratic matrix
+  void reinit(int n) { reinit(n,n); }
+  /// Reinitialization to the same dimensions of another matrix.
+  void reinit(const dFMatrix& A) { reinit(A.m(), A.n()); }
+
+  /// Read access to matrix elements.
+  double operator() (int i, int j) const
+  {
+    //THROW2((i<0) || (i>=dim_image), IntError(IntError::Range,i));
+    //THROW2((j<0) || (j>=dim_range), IntError(IntError::Range,j));
+    return el(i,j);
+  }
+
+  /// Read-Write access to matrix elements
+  double& operator() (int i, int j)
+  {
+    //THROW2((i<0) || (i>=dim_image), IntError(IntError::Range,i));
+    //THROW2((j<0) || (j>=dim_range), IntError(IntError::Range,j));
+    return el(i,j);
+  }
+
+  /// Assignment operator
+  dFMatrix& operator = (const dFMatrix&);
+  /** Copying matrix entries.
+   * Copy the elements of matrix src into this beginning at
+   * element (i,j)
+   */
+
+  void fill(const dFMatrix& src, int i = 0, int j = 0);
+  void add(double s, const dFMatrix& src);
+  void add_diag(double s, const dFMatrix& src);
+  void Tadd(double s, const dFMatrix& src);
+
+  ///
+    /* Add different rows of a matrix
+     * a(i,.) += s * a(j,.)
+     */
+  void add_row(int i, double s, int j);
+  ///
+    /**
+     * Add different rows of a matrix
+     * a(i,.) += s * a(j,.) + t * a(k,.)
+     */
+  void add_row(int i, double s, int j, double t, int k);
+  ///
+    /**
+     * Add different columns of a matrix
+     * a(.,i) += s * a(.,j)
+     */
+  void add_col(int i, double s, int j);
+  ///
+    /**
+     * Add different columns of a matrix
+     * a(.,i) += s * a(.,j) + t * a(.,k)
+     */
+  void add_col(int i, double s, int j, double t, int k);
+  /// Exchange contents of rows i and j
+  void swap_row(int i, int j);
+  /// Exchange contents of columns i and j
+  void swap_col(int i, int j);
+  /// Adding a scalar value on the diagonal
+  void diagadd(const double& src);
+
+  /// Matrix-matrix-multiplication dst = this * src
+  void mmult(dFMatrix& dst, const dFMatrix& src) const;
+  ///
+  void Tmmult(dFMatrix& dst, const dFMatrix& src) const;
+  ///
+    /*
+     * Application of a matrix to a vector.
+     *
+     * [flag] adding determines if the result is copied to dst or added.
+     *
+     * dst (+)= this * src
+     */
+  void vmult(dVector& dst, const dVector& src,const int adding = 0) const;
+  ///
+  void gsmult(dVector& dst, const dVector& src,const iVector& gl) const;
+  /// Application of the transpose matrix to a vector.
+  void Tvmult(dVector& dst,const dVector& src,const int adding=0) const;
+  /// 
+  double residual(dVector& dst, const dVector& src, const dVector& right) const;
+  /// Inversion of lower triangle
+  void forward(dVector& dst, const dVector& src) const;
+  /// Inversion of upper triangle
+  void backward(dVector& dst, const dVector& src) const;
+  /// Replace this by its inverse matrix calculated with Gau&szlig;-Jordan algorithm.
+  void gauss_jordan();
+
+  ///
+  /*
+   * QR - factorization of a matrix.
+   * The orthogonal transformation Q is applied to the vector y and the
+   * matrix. After execution of householder, the upper triangle contains
+   * the resulting matrix R, the lower the incomplete factorization matrices.
+   */
+  void householder(dVector& y);
+
+  /// Least - Squares - Approximation by QR-factorization.
+  double least_squares(dVector& dst, dVector& src);
+
+  /// Output of the matrix in user-defined format.
+  void print(FILE* fp, const char* format = 0) const;
+
+
+                                    /**
+                                     * Comparison operator. Be careful with
+                                     * this thing, it may eat up huge amounts
+                                     * of computing time!
+                                     */
+    bool operator == (const dFMatrix &) const;
+    
+};
+#endif
diff --git a/deal.II/lac/include/lac/dsmatrix.h b/deal.II/lac/include/lac/dsmatrix.h
new file mode 100644 (file)
index 0000000..351d552
--- /dev/null
@@ -0,0 +1,140 @@
+// $Id$
+
+// This file is part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+
+
+#ifndef __lac_dsmatrix_h
+#define __lac_dsmatrix_h
+
+#include <math.h>
+#ifndef __lac_dvector_h
+#include <lac/dvector.h>
+#endif
+#ifndef __lac_ivector_h
+#include <lac/ivector.h>
+#endif
+
+
+class ostream;
+template <int dim> class DoFHandler;
+
+/*
+CLASS
+   dSMatrixStruct
+   */
+class dSMatrixStruct
+{
+  friend class dSMatrix;
+  int max_dim;
+  int rows, cols;
+  int vec_len, max_vec_len;
+  int max_row_len;
+  int* rowstart;
+  int* colnums;
+  int compressed;
+
+public:
+  //////////
+  void reinit(int m, int n, int max_per_row);
+  //////////
+  dSMatrixStruct(int m, int n, int max_per_row) 
+    : max_dim(0), max_vec_len(0), rowstart(0), colnums(0)
+      {
+       reinit(m,n,max_per_row);
+      }
+  //////////
+  dSMatrixStruct(int n, int max_per_row)
+    : max_dim(0), max_vec_len(0), rowstart(0), colnums(0)
+      {
+       reinit(n,n,max_per_row);
+      }
+  //////////
+  ~dSMatrixStruct()
+    {
+      delete[] rowstart;
+      delete[] colnums;
+    }
+  //////////
+  void compress();
+  //////////
+  int operator () (int i, int j);
+  //////////
+  void add(int i, int j);
+  //////////
+  void add_matrix(int n, int* rowcols);
+  //////////
+  void add_matrix(int m, int n, int* rows, int* cols);
+  //////////
+  void add_matrix(const iVector& rowcols);
+  //////////
+  void add_matrix(const iVector& rows, const iVector& cols);
+
+  void print_gnuplot (ostream &) const;
+  int n_rows () const {  return rows;   };
+  int n_cols () const {  return cols;   };
+  int bandwidth () const;
+
+    friend class ConstraintMatrix;
+    friend class DoFHandler<1>;
+    friend class DoFHandler<2>;
+};
+
+/*
+CLASS
+   dSMatrix
+   */
+class dSMatrix
+{
+  dSMatrixStruct& cols;
+  double* val;
+  int max_len;
+ public:
+  //////////
+  void reinit();
+  //////////
+  dSMatrix(dSMatrixStruct& c)
+    : cols(c), val(0), max_len(0)
+    {
+      reinit();
+    }
+  //////////
+  ~dSMatrix()
+    {
+      delete[] val;
+    }
+
+  //////////
+  int m() const { return cols.rows; }
+  //////////
+  int n() const { return cols.cols; }
+
+  //////////
+  void set(int i,int j,double value) { val[cols(i,j)] = value; }
+  //////////
+  void add(int i,int j,double value) { val[cols(i,j)]+= value; }
+  
+  //////////
+  void vmult (dVector& dst,const dVector& src);
+  //////////
+  void Tvmult(dVector& dst,const dVector& src);
+  
+  //////////
+  double residual (dVector& dst,const dVector& x,const dVector& b);
+  //////////
+  void Jacobi_precond(dVector& dst,const dVector& src,double om = 1.);
+  //////////
+  void SSOR_precond(dVector& dst,const dVector& src,double om = 1.);
+  //////////
+  void SOR_precond(dVector& dst,const dVector& src,double om = 1.);
+  //////////
+  void SSOR(dVector& dst,double om = 1.);
+  //////////
+  void SOR(dVector& dst,double om = 1.);
+  //////////
+  void precondition(dVector& dst,const dVector& src) { dst=src; }
+
+    void print (ostream &) const;
+};
+#endif
diff --git a/deal.II/lac/include/lac/dvector.h b/deal.II/lac/include/lac/dvector.h
new file mode 100644 (file)
index 0000000..43cb079
--- /dev/null
@@ -0,0 +1,100 @@
+// $Id$
+
+// This file is part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+
+
+#ifndef __lac_dvector_h
+#define __lac_dvector_h
+#include <stdio.h>
+#ifndef __base_types_h
+#include <base/types.h>
+#endif
+//#ifndef __base_errors_h
+//#include <base/errors.h>
+//#endif
+#ifndef __lac_vectorbase_h
+#include <lac/vectorbase.h>
+#endif
+
+/*
+CLASS
+   dVector
+   */
+class dVector : public VectorBase
+{
+  friend class dFMatrix;
+protected:
+  int dim, maxdim;
+  double *val;
+public:
+  dVector();
+  dVector(const dVector& v);
+  dVector(int n);
+  ~dVector();
+  void reinit(int n, int fast = 0);
+  void reinit(const dVector&, int fast = 0);
+
+  int n() const; // Abfrage der Dimension
+
+  double operator()(int i) const; // read-only Zugriff
+  double& operator()(int i); //Zugriff auf die Komponenten
+
+  double operator*(const dVector& v) const; //Skalarprodukt
+
+  dVector& operator=(double s);
+  dVector& operator=(const dVector& v);
+
+  // GROUP: Addition
+
+  void add(const double);
+  void add(const dVector&);
+  void add(double, const dVector&);
+  void add(double, const dVector&, double, const dVector&);
+
+  void sadd(double, const dVector&);
+  void sadd(double, double, const dVector&);
+  void sadd(double, double, const dVector&, double, const dVector&);
+  void sadd(double, double, const dVector&, double, const dVector&, 
+           double, const dVector&);
+
+  void equ(double);
+  void equ(double, const dVector&);
+  void equ(double, const dVector&, double, const dVector&);
+
+  void czero(int);
+  void cequ(int, const VectorBase&, double, int);
+  void cequ(int, const VectorBase&, double, int, double, int);
+  void cequ(int, const VectorBase&, double, int, double, int, double, int, double, int);
+
+  void cadd(int, const VectorBase&, double, int);
+  void cadd(int, const VectorBase&, double, int, double, int);
+  void cadd(int, const VectorBase&, double, int, double, int, double, int, double, int);
+
+  virtual const char* name() const;
+  //
+  // Output of the vector in user-defined format.
+  //
+  void print(FILE* fp, const char* format = 0) const;
+  void print(const char* format = 0) const;
+};
+
+inline int dVector::n() const
+{
+  return dim;
+}
+
+inline double dVector::operator() (int i) const
+{
+  THROW2( (i<0) || (i>=dim), IntError(IntError::Range,i));
+  return val[i];
+}
+
+inline double& dVector::operator() (int i)
+{
+  THROW2( (i<0) || (i>=dim), IntError(IntError::Range,i));
+  return val[i];
+}
+
+#endif
diff --git a/deal.II/lac/include/lac/ivector.h b/deal.II/lac/include/lac/ivector.h
new file mode 100644 (file)
index 0000000..c558679
--- /dev/null
@@ -0,0 +1,83 @@
+// $Id$
+
+// This file is part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+
+
+#ifndef __lac_ivector_h
+#define __lac_ivector_h
+#ifndef __base_types_h
+#include <base/types.h>
+#endif
+//#ifndef __base_errors_h
+//#include <deal/errors.h>
+//#endif
+
+/*
+CLASS
+   iVector
+   */
+class iVector
+{
+  friend class dFMatrix;
+protected:
+  //////////
+  int dim, maxdim;
+  //////////
+  int *val;
+public:
+  //////////
+  iVector();
+  //////////
+  iVector(const iVector& v);
+  //////////
+  iVector(int n);
+  //////////
+  ~iVector();
+  //////////
+  void reinit(int n, int fast = 0);
+  //////////
+  void reinit(const iVector&, int fast = 0);
+
+  //////////
+  int n() const; // Abfrage der Dimension
+
+  //////////
+  int operator()(int i) const; //read-only Zugriff
+  //////////
+  int& operator()(int i); //Zugriff auf die Komponenten
+
+  //////////
+  iVector& operator=(int i);
+  //////////
+  iVector& operator=(const iVector& v);
+
+  //////////
+  void add(const iVector&);
+  //////////
+  void add(int, const iVector&);
+
+  // Zuweisung
+
+  //////////
+  void equ(int, const iVector&);
+};
+
+inline int iVector::n() const
+{
+  return dim;
+}
+
+inline int iVector::operator() (int i) const
+{
+  THROW2( (i<0) || (i>=dim), IntError(IntError::Range,i));
+  return val[i];
+}
+
+inline int& iVector::operator() (int i)
+{
+  THROW2( (i<0) || (i>=dim), IntError(IntError::Range,i));
+  return val[i];
+}
+#endif
diff --git a/deal.II/lac/include/lac/vectorbase.h b/deal.II/lac/include/lac/vectorbase.h
new file mode 100644 (file)
index 0000000..a619fd1
--- /dev/null
@@ -0,0 +1,37 @@
+// $Id$
+
+// This file is part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+
+
+#ifndef __lac_vectorbase_h
+#define __lac_vectorbase_h
+
+/*
+CLASS
+   VectorBase
+   */
+class VectorBase
+{
+  
+  public:
+  virtual ~VectorBase() {}
+
+  // Komponentenweises Eintragen
+
+  virtual void czero(int) = 0;
+  virtual void cequ(int, const VectorBase&, double, int) = 0;
+  virtual void cequ(int, const VectorBase&, double, int, double, int) = 0;
+  virtual void cequ(int, const VectorBase&, double, int, double, int, double, int, double, int) = 0;
+
+  // Komponentenweise Addition
+
+  virtual void cadd(int, const VectorBase&, double, int) = 0;
+  virtual void cadd(int, const VectorBase&, double, int, double, int) = 0;
+  virtual void cadd(int, const VectorBase&, double, int, double, int, double, int, double, int) = 0;
+
+  virtual const char* name() const = 0;
+};
+
+#endif
diff --git a/deal.II/lac/source/dfmatrix.cc b/deal.II/lac/source/dfmatrix.cc
new file mode 100644 (file)
index 0000000..c7dd222
--- /dev/null
@@ -0,0 +1,862 @@
+// $Id$
+
+static const char* OBJFILE = "DEAL $RCSfile$ $Revision$";
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <lac/dfmatrix.h>
+#include <base/error.h>
+
+dFMatrix::dFMatrix (const dFMatrix &m)
+{
+  init (m.dim_image, m.dim_range);
+  for (int i=0; i!=dim_image*dim_range; ++i)
+    val[i] = m.val[i];
+};
+
+void dFMatrix::init(int mm, int nn)
+{
+  val_size = nn*mm;
+  val = new double[val_size];
+  dim_range = nn;
+  dim_image = mm;
+  for (int i = nn*mm-1; i>=0 ; i--) val[i] = 0;
+}
+
+dFMatrix::~dFMatrix()
+{
+  delete[] val;
+}
+
+void dFMatrix::reinit(int mm, int nn)
+{
+  if (val_size<nn*mm)
+    {
+      delete[] val;
+      init(mm, nn);
+    }
+  else
+    {
+      dim_range = nn;
+      dim_image = mm;
+//      memset(val, 0, sizeof(double)*nn*mm);
+      for (int i = nn*mm-1; i>=0 ; i--) val[i] = 0;
+    }
+}
+
+void dFMatrix::vmult(dVector& dst, const dVector& src,const int adding) const
+{
+  Assert(dst.n() == m(), ExcDimensionMismatch(dst.n(), m()));
+  Assert(src.n() == n(), ExcDimensionMismatch(src.n(), n()));
+
+  double s;
+  if ((n()==3) && (m()==3))
+  {
+    double s0,s1,s2;
+    s   = src(0);
+    s0  = s*val[0]; s1  = s*val[3]; s2  = s*val[6]; 
+    s   = src(1);
+    s0 += s*val[1]; s1 += s*val[4]; s2 += s*val[7];
+    s   = src(2);
+    s0 += s*val[2]; s1 += s*val[5]; s2 += s*val[8];
+
+    if (!adding)
+    {
+      dst(0) = s0;
+      dst(1) = s1;
+      dst(2) = s2;
+    }
+    else
+    {
+      dst(0) += s0;
+      dst(1) += s1;
+      dst(2) += s2;
+    }
+  }
+  else if ((n()==4) && (m()==4))
+  {
+    double s0,s1,s2,s3;
+    s = src(0);
+    s0 = s*val[0]; s1 = s*val[4]; s2 = s*val[8]; s3 = s*val[12];
+    s = src(1);
+    s0 += s*val[1]; s1 += s*val[5]; s2 += s*val[9]; s3 += s*val[13];
+    s = src(2);
+    s0 += s*val[2]; s1 += s*val[6]; s2 += s*val[10]; s3 += s*val[14];
+    s = src(3);
+    s0 += s*val[3]; s1 += s*val[7]; s2 += s*val[11]; s3 += s*val[15];
+    
+    if (!adding)
+    {
+      dst(0) = s0;
+      dst(1) = s1;
+      dst(2) = s2;
+      dst(3) = s3;
+    }
+    else
+    {
+      dst(0) += s0;
+      dst(1) += s1;
+      dst(2) += s2;
+      dst(3) += s3;
+    }
+  }
+  else if ((n()==8) && (m()==8))
+  {
+    double s0,s1,s2,s3,s4,s5,s6,s7;
+    s = src(0);
+    s0 = s*val[0]; s1 = s*val[8]; s2 = s*val[16]; s3 = s*val[24];
+    s4 = s*val[32]; s5 = s*val[40]; s6 = s*val[48]; s7 = s*val[56];
+    s = src(1);
+    s0 += s*val[1]; s1 += s*val[9]; s2 += s*val[17]; s3 += s*val[25];
+    s4 += s*val[33]; s5 += s*val[41]; s6 += s*val[49]; s7 += s*val[57];
+    s = src(2);
+    s0 += s*val[2]; s1 += s*val[10]; s2 += s*val[18]; s3 += s*val[26];
+    s4 += s*val[34]; s5 += s*val[42]; s6 += s*val[50]; s7 += s*val[58];
+    s = src(3);
+    s0 += s*val[3]; s1 += s*val[11]; s2 += s*val[19]; s3 += s*val[27];
+    s4 += s*val[35]; s5 += s*val[43]; s6 += s*val[51]; s7 += s*val[59];
+    s = src(4);
+    s0 += s*val[4]; s1 += s*val[12]; s2 += s*val[20]; s3 += s*val[28];
+    s4 += s*val[36]; s5 += s*val[44]; s6 += s*val[52]; s7 += s*val[60];
+    s = src(5);
+    s0 += s*val[5]; s1 += s*val[13]; s2 += s*val[21]; s3 += s*val[29];
+    s4 += s*val[37]; s5 += s*val[45]; s6 += s*val[53]; s7 += s*val[61];
+    s = src(6);
+    s0 += s*val[6]; s1 += s*val[14]; s2 += s*val[22]; s3 += s*val[30];
+    s4 += s*val[38]; s5 += s*val[46]; s6 += s*val[54]; s7 += s*val[62];
+    s = src(7);
+    s0 += s*val[7]; s1 += s*val[15]; s2 += s*val[23]; s3 += s*val[31];
+    s4 += s*val[39]; s5 += s*val[47]; s6 += s*val[55]; s7 += s*val[63];
+    
+    if (!adding)
+    {
+      dst(0) = s0;
+      dst(1) = s1;
+      dst(2) = s2;
+      dst(3) = s3;
+      dst(4) = s4;
+      dst(5) = s5;
+      dst(6) = s6;
+      dst(7) = s7;
+    }
+    else
+    {
+      dst(0) += s0;
+      dst(1) += s1;
+      dst(2) += s2;
+      dst(3) += s3;
+      dst(4) += s4;
+      dst(5) += s5;
+      dst(6) += s6;
+      dst(7) += s7;
+    }
+  }
+  else
+  {    
+    double* e = val;
+    for (int i=0;i<m();i++)
+    {
+      s = 0.;
+      for (int j=0;j<n();j++) s += src(j) * *(e++);
+      if (!adding) dst(i) = s;
+      else dst(i) += s;
+    }
+  }
+}
+
+void dFMatrix::gsmult(dVector& dst, const dVector& src,const iVector& gl) const
+{
+  Assert(n() == m(), ExcNotQuadratic());
+  Assert(dst.n() == n(), ExcDimensionMismatch(dst.n(), n()));
+  Assert(src.n() == n(), ExcDimensionMismatch(src.n(), n()));
+  Assert(gl.n() == n(), ExcDimensionMismatch(gl.n(), n()));
+  double s;
+  if ((n()==3) && (m()==3))
+  {
+    double s0=0.,s1=0.,s2=0.;
+    s   = src(0);
+    if(gl(1)<gl(0)) s1  = s*val[3]; if(gl(2)<gl(0))  s2  = s*val[6]; 
+    s   = src(1);
+    if(gl(0)<gl(1)) s0 += s*val[1]; if(gl(2)<gl(1))  s2 += s*val[7];
+    s   = src(2);
+    if(gl(0)<gl(2)) s0 += s*val[2]; if(gl(1)<gl(2))  s1 += s*val[5]; 
+
+    dst(0) += s0;
+    dst(1) += s1;
+    dst(2) += s2;
+  }
+  else if ((n()==4) && (m()==4))
+  {
+    double s0=0.,s1=0.,s2=0.,s3=0.;
+    s = src(0);
+    if(gl(1)<gl(0)) s1 = s*val[4];  if(gl(2)<gl(0)) s2 = s*val[8]; if(gl(3)<gl(0)) s3 = s*val[12];
+    s = src(1);
+    if(gl(0)<gl(1)) s0 += s*val[1]; if(gl(2)<gl(1)) s2 += s*val[9]; if(gl(3)<gl(1)) s3 += s*val[13];
+    s = src(2);
+    if(gl(0)<gl(2)) s0 += s*val[2]; if(gl(1)<gl(2)) s1 += s*val[6]; if(gl(3)<gl(2)) s3 += s*val[14];
+    s = src(3);
+    if(gl(0)<gl(3)) s0 += s*val[3]; if(gl(1)<gl(3)) s1 += s*val[7]; if(gl(2)<gl(3)) s2 += s*val[11];
+
+    dst(0) += s0;
+    dst(1) += s1;
+    dst(2) += s2;
+    dst(3) += s3;
+  }
+  else if ((n()==8) && (m()==8))
+  {
+    double s0=0.,s1=0.,s2=0.,s3=0.,s4=0.,s5=0.,s6=0.,s7=0.;
+    s = src(0);
+    if(gl(1)<gl(0)) s1 = s*val[8]; 
+    if(gl(2)<gl(0)) s2 = s*val[16]; 
+    if(gl(3)<gl(0)) s3 = s*val[24];
+    if(gl(4)<gl(0)) s4 = s*val[32]; 
+    if(gl(5)<gl(0)) s5 = s*val[40]; 
+    if(gl(6)<gl(0)) s6 = s*val[48]; 
+    if(gl(7)<gl(0)) s7 = s*val[56];
+    s = src(1);
+    if(gl(0)<gl(1)) s0 += s*val[1]; 
+    if(gl(2)<gl(1)) s2 += s*val[17]; 
+    if(gl(3)<gl(1)) s3 += s*val[25];
+    if(gl(4)<gl(1)) s4 += s*val[33]; 
+    if(gl(5)<gl(1)) s5 += s*val[41]; 
+    if(gl(6)<gl(1)) s6 += s*val[49]; 
+    if(gl(7)<gl(1)) s7 += s*val[57];
+    s = src(2);
+    if(gl(0)<gl(2)) s0 += s*val[2]; 
+    if(gl(1)<gl(2)) s1 += s*val[10]; 
+    if(gl(3)<gl(2)) s3 += s*val[26];
+    if(gl(4)<gl(2)) s4 += s*val[34]; 
+    if(gl(5)<gl(2)) s5 += s*val[42]; 
+    if(gl(6)<gl(2)) s6 += s*val[50]; 
+    if(gl(7)<gl(2)) s7 += s*val[58];
+    s = src(3);
+    if(gl(0)<gl(3)) s0 += s*val[3]; 
+    if(gl(1)<gl(3)) s1 += s*val[11]; 
+    if(gl(2)<gl(3)) s2 += s*val[19]; 
+    if(gl(4)<gl(3)) s4 += s*val[35]; 
+    if(gl(5)<gl(3)) s5 += s*val[43]; 
+    if(gl(6)<gl(3)) s6 += s*val[51]; 
+    if(gl(7)<gl(3)) s7 += s*val[59];
+    s = src(4);
+    if(gl(0)<gl(4)) s0 += s*val[4]; 
+    if(gl(1)<gl(4)) s1 += s*val[12]; 
+    if(gl(2)<gl(4)) s2 += s*val[20]; 
+    if(gl(3)<gl(4)) s3 += s*val[28];
+    if(gl(5)<gl(4)) s5 += s*val[44]; 
+    if(gl(6)<gl(4)) s6 += s*val[52]; 
+    if(gl(7)<gl(4)) s7 += s*val[60];
+    s = src(5);
+    if(gl(0)<gl(5)) s0 += s*val[5]; 
+    if(gl(1)<gl(5)) s1 += s*val[13]; 
+    if(gl(2)<gl(5)) s2 += s*val[21]; 
+    if(gl(3)<gl(5)) s3 += s*val[29];
+    if(gl(4)<gl(5)) s4 += s*val[37]; 
+    if(gl(6)<gl(5)) s6 += s*val[53]; 
+    if(gl(7)<gl(5)) s7 += s*val[61];
+    s = src(6);
+    if(gl(0)<gl(6)) s0 += s*val[6]; 
+    if(gl(1)<gl(6)) s1 += s*val[14]; 
+    if(gl(2)<gl(6)) s2 += s*val[22]; 
+    if(gl(3)<gl(6)) s3 += s*val[30];
+    if(gl(4)<gl(6)) s4 += s*val[38]; 
+    if(gl(5)<gl(6)) s5 += s*val[46]; 
+    if(gl(7)<gl(6)) s7 += s*val[62];
+    s = src(7);
+    if(gl(0)<gl(7)) s0 += s*val[7]; 
+    if(gl(1)<gl(7)) s1 += s*val[15]; 
+    if(gl(2)<gl(7)) s2 += s*val[23]; 
+    if(gl(3)<gl(7)) s3 += s*val[31];
+    if(gl(4)<gl(7)) s4 += s*val[39]; 
+    if(gl(5)<gl(7)) s5 += s*val[47]; 
+    if(gl(6)<gl(7)) s6 += s*val[55]; 
+    
+    dst(0) += s0;
+    dst(1) += s1;
+    dst(2) += s2;
+    dst(3) += s3;
+    dst(4) += s4;
+    dst(5) += s5;
+    dst(6) += s6;
+    dst(7) += s7;
+  }
+  else
+  {    
+    double* e = val;
+    for (int i=0;i<m();i++)
+    {
+      s = 0.;
+      for (int j=0;j<n();j++) if(gl(i)<gl(j)) s += src(j) * *(e++);
+      dst(i) += s;
+    }
+  }
+}
+
+void dFMatrix::Tvmult(dVector& dst, const dVector& src, const int adding) const
+{
+  Assert(dst.n() == n(), ExcDimensionMismatch(dst.n(), n()));
+  Assert(src.n() == m(), ExcDimensionMismatch(src.n(), m()));
+
+  int i,j;
+  double s;
+  for (i=0;i<m();i++)
+  {
+    s = 0.;
+    for (j=0;j<n();j++) s += src(j) * el(j,i);
+    if(!adding) dst(i) = s;
+    else dst(i) += s;
+  }
+}
+
+double dFMatrix::residual(dVector& dst, const dVector& src,
+                         const dVector& right) const
+{
+  Assert(dst.n() == m(), ExcDimensionMismatch(dst.n(), m()));
+  Assert(src.n() == n(), ExcDimensionMismatch(src.n(), n()));
+  Assert(right.n() == m(), ExcDimensionMismatch(right.n(), m()));
+  int i,j;
+  double s, res = 0.;
+  for (i=0;i<n();i++)
+    {
+      s = right(i);
+      for (j=0;j<m();j++) s -= src(j) * el(i,j);
+      dst(i) = s;
+      res += s*s;
+    }
+  return sqrt(res);
+}
+
+void dFMatrix::forward(dVector& dst, const dVector& src) const
+{
+  Assert(n() == m(), ExcNotQuadratic());
+  Assert(dst.n() == n(), ExcDimensionMismatch(dst.n(), n()));
+  Assert(src.n() == n(), ExcDimensionMismatch(src.n(), n()));
+  int i,j;
+  int nu = MIN(m(),n());
+  double s;
+  for (i=0;i<nu;i++)
+    {
+      s = src(i);
+      for (j=0;j<i;j++) s -= dst(j) * el(i,j);
+      dst(i) = s/el(i,i);
+    }
+}
+
+void dFMatrix::backward(dVector& dst, const dVector& src) const
+{
+  Assert(n() == m(), ExcNotQuadratic());
+  Assert(dst.n() == n(), ExcDimensionMismatch(dst.n(), n()));
+  Assert(src.n() == n(), ExcDimensionMismatch(src.n(), n()));
+  int i,j;
+  int nu = MIN(m(),n());
+  double s;
+  for (i=nu-1;i>=0;i--)
+    {
+      s = src(i);
+      for (j=i+1;j<nu;j++) s -= dst(j) * el(i,j);
+      dst(i) = s/el(i,i);
+    }
+}
+
+dFMatrix& dFMatrix::operator = (const dFMatrix& M)
+{
+  reinit(M);
+  int nn = n()*m();
+  for (int i=0;i<nn;i++) val[i] = M.val[i];
+  return *this;
+}
+
+void dFMatrix::fill(const dFMatrix& src, int i, int j)
+{
+  //THROW2 (n() < src.n() + j, Error(Error::Dimension,"fill",OBJFILE));
+  //THROW2 (m() < src.m() + i, Error(Error::Dimension,"fill",OBJFILE));
+
+  for (int ii=0; ii<src.m() ; ii++)
+    for (int jj=0; jj<src.n() ; jj++)
+      el(ii+i,jj+j) = src.el(ii,jj);
+}
+
+void dFMatrix::add_row(int i, double s, int j)
+{
+  int k;
+  for (k=0;k<m();k++) el(i,k) += s*el(j,k);
+}
+
+void dFMatrix::add_row(int i, double s, int j, double t, int k)
+{
+  int l;
+  for (l=0;l<m();l++) el(i,l) += s*el(j,l) + t*el(k,l);
+}
+
+void dFMatrix::add_col(int i, double s, int j)
+{
+  int k;
+  for (k=0;k<n();k++) el(k,i) += s*el(k,j);
+}
+
+void dFMatrix::add_col(int i, double s, int j, double t, int k)
+{
+  int l;
+  for (l=0;l<n();l++) el(l,i) += s*el(l,j) + t*el(l,k);
+}
+
+void dFMatrix::swap_row(int i, int j)
+{
+  int k;
+  double s;
+  for (k=0;k<m();k++)
+  {
+    s = el(i,k); el(i,k) = el(j,k); el(j,k) = s;
+  }
+}
+
+void dFMatrix::swap_col(int i, int j)
+{
+  int k;
+  double s;
+  for (k=0;k<n();k++)
+  {
+    s = el(k,i); el(k,i) = el(k,j); el(k,j) = s;
+  }
+}
+
+void dFMatrix::diagadd(const double& src)
+{
+  //THROW1 (m() != n(), Error(Error::Dimension,"diagadd"));
+  for (int i=0;i<n();i++)
+    el(i,i) += src;
+}
+
+void dFMatrix::mmult(dFMatrix& dst, const dFMatrix& src) const
+{
+  //THROW1 (n() != src.m(), Error(Error::Dimension,"mmult"));
+  int i,j,k;
+  double s = 1.;
+  dst.reinit(m(), src.n());
+
+  for (i=0;i<m();i++)
+    for (j=0;j<src.n();j++)
+      {
+       s = 0.;
+       for (k=0;k<n();k++) s+= el(i,k) * src.el(k,j);
+       dst.el(i,j) = s;
+      }
+}
+
+/*void dFMatrix::mmult(dFMatrix& dst, const dFMatrix& src) const
+{
+  int i,j,k;
+  double s = 1.;
+  THROW1 (m() != src.n(), Error(Error::Dimension,"mmult"));
+  dst.reinit(n(), src.m());
+
+  for (i=0;i<n();i++)
+    for (j=0;j<src.m();j++)
+      {
+       s = 0.;
+       for (k=0;k<m();k++) s+= el(i,k) * src.el(k,j);
+       dst.el(i,j) = s;
+      }
+}*/
+
+void dFMatrix::Tmmult(dFMatrix& dst, const dFMatrix& src) const
+{
+  //THROW1 (n() != src.n(), Error(Error::Dimension,"Tmmult"));
+  int i,j,k;
+  double s = 1.;
+  dst.reinit(m(), src.m());
+
+  for (i=0;i<m();i++)
+    for (j=0;j<src.m();j++)
+      {
+       s = 0.;
+       for (k=0;k<n();k++) s+= el(k,i) * src.el(k,j);
+       dst.el(i,j) = s;
+      }
+}
+
+/*void dFMatrix::Tmmult(dFMatrix& dst, const dFMatrix& src) const
+{
+  int i,j,k;
+  double s = 1.;
+  THROW1 (m() != src.n(), Error(Error::Dimension,"Tmmult"));
+  dst.reinit(n(), src.m());
+
+  for (i=0;i<n();i++)
+    for (j=0;j<src.m();j++)
+      {
+       s = 0.;
+       for (k=0;k<m();k++) s+= el(k,i) * src.el(k,j);
+       dst.el(i,j) = s;
+      }
+}*/
+
+void dFMatrix::print(FILE* f, const char* format) const
+{
+  if (!format) format = " %5.2f";
+  int i,j;
+  for (i=0;i<m();i++)
+    {
+      for (j=0;j<n();j++) fprintf(f, format, el(i,j));
+      fprintf(f, "\n");
+    }
+}
+
+void dFMatrix::add(double s,const dFMatrix& src)
+{
+  //THROW1 (m() != src.m(), Error(Error::Dimension,"add"));
+  //THROW1 (n() != src.n(), Error(Error::Dimension,"add"));
+  if ((n()==3) && (m()==3))
+  {
+    val[0] += s * src.el(0);
+    val[1] += s * src.el(1);
+    val[2] += s * src.el(2);
+    val[3] += s * src.el(3);
+    val[4] += s * src.el(4);
+    val[5] += s * src.el(5);
+    val[6] += s * src.el(6);
+    val[7] += s * src.el(7);
+    val[8] += s * src.el(8);
+  }
+  else if ((n()==4) && (m()==4))
+  {
+    val[0] += s * src.el(0);
+    val[1] += s * src.el(1);
+    val[2] += s * src.el(2);
+    val[3] += s * src.el(3);
+    val[4] += s * src.el(4);
+    val[5] += s * src.el(5);
+    val[6] += s * src.el(6);
+    val[7] += s * src.el(7);
+    val[8] += s * src.el(8);
+    val[9] += s * src.el(9);
+    val[10] += s * src.el(10);
+    val[11] += s * src.el(11);
+    val[12] += s * src.el(12);
+    val[13] += s * src.el(13);
+    val[14] += s * src.el(14);
+    val[15] += s * src.el(15);
+  }
+  else if ((n()==8) && (m()==8))
+  {
+    val[0] += s * src.el(0);
+    val[1] += s * src.el(1);
+    val[2] += s * src.el(2);
+    val[3] += s * src.el(3);
+    val[4] += s * src.el(4);
+    val[5] += s * src.el(5);
+    val[6] += s * src.el(6);
+    val[7] += s * src.el(7);
+    val[8] += s * src.el(8);
+    val[9] += s * src.el(9);
+    val[10] += s * src.el(10);
+    val[11] += s * src.el(11);
+    val[12] += s * src.el(12);
+    val[13] += s * src.el(13);
+    val[14] += s * src.el(14);
+    val[15] += s * src.el(15);
+    val[16] += s * src.el(16);
+    val[17] += s * src.el(17);
+    val[18] += s * src.el(18);
+    val[19] += s * src.el(19);
+
+    val[20] += s * src.el(20);
+    val[21] += s * src.el(21);
+    val[22] += s * src.el(22);
+    val[23] += s * src.el(23);
+    val[24] += s * src.el(24);
+    val[25] += s * src.el(25);
+    val[26] += s * src.el(26);
+    val[27] += s * src.el(27);
+    val[28] += s * src.el(28);
+    val[29] += s * src.el(29);
+
+    val[30] += s * src.el(30);
+    val[31] += s * src.el(31);
+    val[32] += s * src.el(32);
+    val[33] += s * src.el(33);
+    val[34] += s * src.el(34);
+    val[35] += s * src.el(35);
+    val[36] += s * src.el(36);
+    val[37] += s * src.el(37);
+    val[38] += s * src.el(38);
+    val[39] += s * src.el(39);
+
+    val[40] += s * src.el(40);
+    val[41] += s * src.el(41);
+    val[42] += s * src.el(42);
+    val[43] += s * src.el(43);
+    val[44] += s * src.el(44);
+    val[45] += s * src.el(45);
+    val[46] += s * src.el(46);
+    val[47] += s * src.el(47);
+    val[48] += s * src.el(48);
+    val[49] += s * src.el(49);
+
+    val[50] += s * src.el(50);
+    val[51] += s * src.el(51);
+    val[52] += s * src.el(52);
+    val[53] += s * src.el(53);
+    val[54] += s * src.el(54);
+    val[55] += s * src.el(55);
+    val[56] += s * src.el(56);
+    val[57] += s * src.el(57);
+    val[58] += s * src.el(58);
+    val[59] += s * src.el(59);
+
+    val[60] += s * src.el(60);
+    val[61] += s * src.el(61);
+    val[62] += s * src.el(62);
+    val[63] += s * src.el(63);
+  }
+  else
+  {
+    for (int i = n()*m()-1; i>=0 ; i--) val[i] += s * src.el(i);
+  }
+}
+
+void dFMatrix::add_diag(double s,const dFMatrix& src)
+{
+  //THROW1 (m() != src.m(), Error(Error::Dimension,"add_diag"));
+  //THROW1 (n() != src.n(), Error(Error::Dimension,"add_diag"));
+  if ((n()==3) && (m()==3))
+  {
+    val[0] += s * src.el(0);
+    val[0] += s * src.el(1);
+    val[0] += s * src.el(2);
+    val[3] += s * src.el(3);
+    val[3] += s * src.el(4);
+    val[3] += s * src.el(5);
+    val[6] += s * src.el(6);
+    val[6] += s * src.el(7);
+    val[6] += s * src.el(8);
+  }
+  else if ((n()==4) && (m()==4))
+  {
+    val[0] += s * src.el(0);
+    val[0] += s * src.el(1);
+    val[0] += s * src.el(2);
+    val[0] += s * src.el(3);
+    val[4] += s * src.el(4);
+    val[4] += s * src.el(5);
+    val[4] += s * src.el(6);
+    val[4] += s * src.el(7);
+    val[8] += s * src.el(8);
+    val[8] += s * src.el(9);
+    val[8] += s * src.el(10);
+    val[8] += s * src.el(11);
+    val[12] += s * src.el(12);
+    val[12] += s * src.el(13);
+    val[12] += s * src.el(14);
+    val[12] += s * src.el(15);
+  }
+  else if ((n()==8) && (m()==8))
+  {
+    val[0] += s * src.el(0);
+    val[0] += s * src.el(1);
+    val[0] += s * src.el(2);
+    val[0] += s * src.el(3);
+    val[0] += s * src.el(4);
+    val[0] += s * src.el(5);
+    val[0] += s * src.el(6);
+    val[0] += s * src.el(7);
+    val[8] += s * src.el(8);
+    val[8] += s * src.el(9);
+    val[8] += s * src.el(10);
+    val[8] += s * src.el(11);
+    val[8] += s * src.el(12);
+    val[8] += s * src.el(13);
+    val[8] += s * src.el(14);
+    val[8] += s * src.el(15);
+    val[16] += s * src.el(16);
+    val[16] += s * src.el(17);
+    val[16] += s * src.el(18);
+    val[16] += s * src.el(19);
+
+    val[16] += s * src.el(20);
+    val[16] += s * src.el(21);
+    val[16] += s * src.el(22);
+    val[16] += s * src.el(23);
+    val[24] += s * src.el(24);
+    val[24] += s * src.el(25);
+    val[24] += s * src.el(26);
+    val[24] += s * src.el(27);
+    val[24] += s * src.el(28);
+    val[24] += s * src.el(29);
+
+    val[24] += s * src.el(30);
+    val[24] += s * src.el(31);
+    val[32] += s * src.el(32);
+    val[32] += s * src.el(33);
+    val[32] += s * src.el(34);
+    val[32] += s * src.el(35);
+    val[32] += s * src.el(36);
+    val[32] += s * src.el(37);
+    val[32] += s * src.el(38);
+    val[32] += s * src.el(39);
+
+    val[40] += s * src.el(40);
+    val[40] += s * src.el(41);
+    val[40] += s * src.el(42);
+    val[40] += s * src.el(43);
+    val[40] += s * src.el(44);
+    val[40] += s * src.el(45);
+    val[40] += s * src.el(46);
+    val[40] += s * src.el(47);
+    val[48] += s * src.el(48);
+    val[48] += s * src.el(49);
+
+    val[48] += s * src.el(50);
+    val[48] += s * src.el(51);
+    val[48] += s * src.el(52);
+    val[48] += s * src.el(53);
+    val[48] += s * src.el(54);
+    val[48] += s * src.el(55);
+    val[56] += s * src.el(56);
+    val[56] += s * src.el(57);
+    val[56] += s * src.el(58);
+    val[56] += s * src.el(59);
+
+    val[56] += s * src.el(60);
+    val[56] += s * src.el(61);
+    val[56] += s * src.el(62);
+    val[56] += s * src.el(63);
+  }
+  else
+  {
+    for (int i = n()*m()-1; i>=0 ; i--) val[i] += s * src.el(i);
+  }
+}
+
+void dFMatrix::Tadd(double s,const dFMatrix& src)
+{
+  //THROW1 (m() != n(),     Error(Error::Dimension,"Tadd"));
+  //THROW1 (m() != src.m(), Error(Error::Dimension,"Tadd"));
+  //THROW1 (n() != src.n(), Error(Error::Dimension,"Tadd"));
+  if ((n()==3) && (m()==3))
+  {
+    val[0] += s * src.el(0);
+    val[1] += s * src.el(3);
+    val[2] += s * src.el(6);
+
+    val[3] += s * src.el(1);
+    val[4] += s * src.el(4);
+    val[5] += s * src.el(7);
+
+    val[6] += s * src.el(2);
+    val[7] += s * src.el(5);
+    val[8] += s * src.el(8);
+  }
+  else if ((n()==4) && (m()==4))
+  {
+    val[0] += s * src.el(0);
+    val[1] += s * src.el(4);
+    val[2] += s * src.el(8);
+    val[3] += s * src.el(12);
+
+    val[4] += s * src.el(1);
+    val[5] += s * src.el(5);
+    val[6] += s * src.el(9);
+    val[7] += s * src.el(13);
+
+    val[8] += s * src.el(2);
+    val[9] += s * src.el(6);
+    val[10] += s * src.el(10);
+    val[11] += s * src.el(14);
+
+    val[12] += s * src.el(3);
+    val[13] += s * src.el(7);
+    val[14] += s * src.el(11);
+    val[15] += s * src.el(15);
+  }
+  else if ((n()==8) && (m()==8))
+  {
+    val[0] += s * src.el(0);
+    val[1] += s * src.el(8);
+    val[2] += s * src.el(16);
+    val[3] += s * src.el(24);
+    val[4] += s * src.el(32);
+    val[5] += s * src.el(40);
+    val[6] += s * src.el(48);
+    val[7] += s * src.el(56);
+
+    val[8] += s * src.el(1);
+    val[9] += s * src.el(9);
+    val[10] += s * src.el(17);
+    val[11] += s * src.el(25);
+    val[12] += s * src.el(33);
+    val[13] += s * src.el(41);
+    val[14] += s * src.el(49);
+    val[15] += s * src.el(57);
+
+    val[16] += s * src.el(2);
+    val[17] += s * src.el(10);
+    val[18] += s * src.el(18);
+    val[19] += s * src.el(26);
+    val[20] += s * src.el(34);
+    val[21] += s * src.el(42);
+    val[22] += s * src.el(50);
+    val[23] += s * src.el(58);
+
+    val[24] += s * src.el(3);
+    val[25] += s * src.el(11);
+    val[26] += s * src.el(19);
+    val[27] += s * src.el(27);
+    val[28] += s * src.el(35);
+    val[29] += s * src.el(43);
+    val[30] += s * src.el(51);
+    val[31] += s * src.el(59);
+
+    val[32] += s * src.el(4);
+    val[33] += s * src.el(12);
+    val[34] += s * src.el(20);
+    val[35] += s * src.el(28);
+    val[36] += s * src.el(36);
+    val[37] += s * src.el(44);
+    val[38] += s * src.el(52);
+    val[39] += s * src.el(60);
+
+    val[40] += s * src.el(5);
+    val[41] += s * src.el(13);
+    val[42] += s * src.el(21);
+    val[43] += s * src.el(29);
+    val[44] += s * src.el(37);
+    val[45] += s * src.el(45);
+    val[46] += s * src.el(53);
+    val[47] += s * src.el(61);
+
+    val[48] += s * src.el(6);
+    val[49] += s * src.el(14);
+    val[50] += s * src.el(22);
+    val[51] += s * src.el(30);
+    val[52] += s * src.el(38);
+    val[53] += s * src.el(46);
+    val[54] += s * src.el(54);
+    val[55] += s * src.el(62);
+
+    val[56] += s * src.el(7);
+    val[57] += s * src.el(15);
+    val[58] += s * src.el(23);
+    val[59] += s * src.el(31);
+    val[60] += s * src.el(39);
+    val[61] += s * src.el(47);
+    val[62] += s * src.el(55);
+    val[63] += s * src.el(63);
+  }
+  else
+  {
+    THROWUNCOND(1, Error(Error::Dimension,"Tadd"));
+  }
+}
+
+bool
+dFMatrix::operator == (const dFMatrix &m) const
+{
+  bool q = (dim_range==m.dim_range) && (dim_image==m.dim_image);
+  if (!q) return false;
+
+  for (int i=0; i<dim_image; ++i)
+    for (int j=0; j<dim_range; ++j)
+      if (el(i,j) != m.el(i,j)) return false;
+  return true;
+};
diff --git a/deal.II/lac/source/dsmatrix.cc b/deal.II/lac/source/dsmatrix.cc
new file mode 100644 (file)
index 0000000..f55f02f
--- /dev/null
@@ -0,0 +1,361 @@
+// $Id$
+
+// This file is part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+
+static const char* OBJFILE = "DEAL $RCSfile$ $Revision$";
+
+#include <lac/dsmatrix.h>
+#include <ostream.h>
+#include <base/sort.h>
+
+void
+dSMatrixStruct::reinit(int m, int n, int max_per_row)
+{
+  THROW1(m<=0, IntError(IntError::IllegalDimension, m));
+  THROW1(n<=0, IntError(IntError::IllegalDimension, n));
+  THROW1(max_per_row<=0, IntError(IntError::IllegalDimension, max_per_row));
+  rows = m;
+  cols = n;
+  vec_len = m * max_per_row;
+  max_row_len = max_per_row;
+
+  if (rows > max_dim)
+    {
+      if (rowstart) delete[] rowstart;
+      max_dim = rows;
+      rowstart = new int[max_dim+1];
+    }
+
+  if (vec_len > max_vec_len)
+    {
+      if (colnums) delete[] colnums;
+      max_vec_len = vec_len;
+      colnums = new int[max_vec_len];
+    }
+  int i;
+  for (i=0;i<=rows;i++) rowstart[i] = i * max_per_row;
+  for (i = vec_len-1 ; i>=0 ; i--) colnums[i] = -1;
+
+  if (rows == cols)
+    for (i=0;i<rows;i++)
+      colnums[rowstart[i]] = i;
+
+  compressed = 0;
+}
+
+void
+dSMatrixStruct::compress()
+{
+  if (compressed) return;
+  int i,j,k,l,s;
+  int* entries;
+
+  entries = new int[max_row_len];
+
+  // Traverse all rows
+
+  for (i=0,k=0,s=0 ; i < rows ; i++)
+    {
+
+      // Sort entries in ascending order
+
+      for (j=rowstart[i],l=0; j < rowstart[i+1] ;++j,++l)
+       entries[l] = colnums[j];
+      heapsort(max_row_len, entries);
+
+      // Re-insert column numbers into the field
+
+      // Ensure diagonal entry first in quadratic matrix
+
+      if (cols == rows) colnums[k++] = i;
+
+      for (l=0; l < max_row_len ; ++l)
+       {
+         if (entries[l] >= 0)
+           if ((i!=entries[l]) || (rows!=cols))
+             colnums[k++] = entries[l];
+       }
+      rowstart[i] = s;
+      s = k;
+    }
+  vec_len = rowstart[i] = s;
+  compressed = 1;
+
+  delete[] entries;
+}
+
+int
+dSMatrixStruct::operator () (int i, int j)
+{
+  for (int k=rowstart[i] ; k<rowstart[i+1] ; k++)
+    if (colnums[k] == j) return k;
+  return -1;
+}
+
+void
+dSMatrixStruct::add (int i, int j)
+{
+  //THROW1((i<0) || (i>=rows), IntError(IntError::Range,i));
+  //THROW1((j<0) || (j>=cols), IntError(IntError::Range,j));
+  int k;
+  for (k=rowstart[i]; k<rowstart[i+1]; k++)
+    {
+      if (colnums[k] == j) return;
+      if (colnums[k] == -1)
+       {
+         colnums[k] = j;
+         return;
+       }
+    }
+}
+
+void
+dSMatrixStruct::add_matrix(int n, int* rowcols)
+{
+  int i,j;
+  for (i=0;i<n;i++)
+    for (j=0;j<n;j++)
+      add(rowcols[i], rowcols[j]);
+}
+
+
+void
+dSMatrixStruct::add_matrix(int m, int n, int* rows, int* cols)
+{
+  int i,j;
+  for (i=0;i<m;i++)
+    for (j=0;j<n;j++)
+      add(rows[i], cols[j]);
+}
+
+void
+dSMatrixStruct::add_matrix(const iVector& rowcols)
+{
+  int i,j;
+  for (i=0;i<rowcols.n();i++)
+    for (j=0;j<rowcols.n();j++)
+      add(rowcols(i), rowcols(j));
+}
+
+
+void
+dSMatrixStruct::add_matrix(const iVector& rows, const iVector& cols)
+{
+  int i,j;
+  for (i=0;i<rows.n();i++)
+    for (j=0;j<cols.n();j++)
+      add(rows(i), cols(j));
+}
+
+void
+dSMatrix::reinit()
+{
+  if(max_len<cols.vec_len)
+  {
+    if (val) delete[] val;
+    val = new double[cols.vec_len];
+    max_len = cols.vec_len;
+  }
+//  memset(val, 0, sizeof(*val) * cols.vec_len);
+  for (int i = cols.vec_len-1 ; i>=0 ; i--) val[i] = 0;
+}
+
+void
+dSMatrixStruct::print_gnuplot (ostream &out) const
+{
+  for (int i=0; i<rows; ++i)
+    for (int j=rowstart[i]; j<rowstart[i+1]; ++j)
+      if (colnums[j]>=0)
+       out << i << " " << -colnums[j] << endl;
+}
+
+
+int
+dSMatrixStruct::bandwidth () const
+{
+  int b=0;
+  for (int i=0; i<rows; ++i)
+    for (int j=rowstart[i]; j<rowstart[i+1]; ++j)
+      if (colnums[j]>=0) 
+       {
+         if (abs(i-colnums[j]) > b)
+           b = abs(i-colnums[j]);
+       }
+      else
+                                        // leave if at the end of
+                                        // the entries of this line
+       break;
+  return b;
+}
+
+void
+dSMatrix::vmult(dVector& dst,const dVector& src)
+{
+  // Assert(m() = dst.n())
+  // Assert(n() = src.n())
+
+  for (int i=0;i<m();i++)
+    {
+      double s = 0.;
+      for (int j=cols.rowstart[i]; j < cols.rowstart[i+1] ;j++) 
+       {
+         int p = cols.colnums[j];
+         s += val[j] * src(p);
+       }
+      dst(i) = s;
+    }
+}
+
+void
+dSMatrix::Tvmult(dVector& dst,const dVector& src) 
+{
+  // Assert(n() = dst.n())
+  // Assert(m() = src.n())
+
+  int i;
+  
+  for (i=0;i<n();i++) dst(i) = 0.;
+  
+  for (i=0;i<m();i++)
+    {
+      for (int j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+       {
+         int p = cols.colnums[j];
+         dst(p) += val[j] * src(i);
+       }
+    }
+}
+
+double
+dSMatrix::residual(dVector& dst,const dVector& u,const dVector& b)
+{
+  // Assert(m() = dst.n())
+  // Assert(n() = src.n())
+
+  double s,norm=0.;   
+  
+  for (int i=0;i<m();i++)
+    {
+      s = b(i);
+      for (int j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+       {
+         int p = cols.colnums[j];
+         s -= val[j] * u(p);
+       }
+      dst(i) = s;
+      norm += dst(i)*dst(i);
+    }
+  return sqrt(norm);
+}
+
+void
+dSMatrix::Jacobi_precond(dVector& dst,const dVector& src,double om)
+{
+  int n = src.n();
+
+  for (int i=0;i<n;++i)
+    {
+      dst(i) = om * src(i) * val[cols.rowstart[i]];
+    }
+}
+
+void
+dSMatrix::SSOR_precond(dVector& dst,const dVector& src,double om)
+{
+  int  p,n = src.n();
+  int  i,j;
+  
+  for (i=0;i<n;i++)
+    {
+      dst(i) = src(i);
+      for (j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+       {
+         p = cols.colnums[j];
+         if (p<i) dst(i) -= om* val[j] * dst(p);
+       }
+      dst(i) /= val[cols.rowstart[i]];
+    }
+  for (i=0;i<n;i++) dst(i) *= (2.-om)*val[cols.rowstart[i]];
+  
+  for (i=n-1;i>=0;i--)
+    {
+      for (j=cols.rowstart[i];j<cols.rowstart[i+1];j++)
+       {
+         p = cols.colnums[j];
+         if (p>i) dst(i) -= om* val[j] * dst(p);
+       }
+      dst(i) /= val[cols.rowstart[i]];
+    }
+}
+
+void
+dSMatrix::SOR_precond(dVector& dst,const dVector& src,double om)
+{
+  dst = src;
+  SOR(dst,om);
+}
+
+void
+dSMatrix::SOR(dVector& dst, double om)
+{
+  // Assert(m()==n())
+  // Assert(m()==dst.n())
+
+  for (int i=0;i<m();i++)
+    {
+      double s = dst(i);
+      for (int j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+       {
+         int p = cols.colnums[j];
+         if (p<i)
+           s -= val[j] * dst(p);
+       }
+      dst(i) = s * om / val[cols.rowstart[i]];
+    }
+}
+
+void
+dSMatrix::SSOR(dVector& dst, double om)
+{
+  int  p,n = dst.n();
+  int  i,j;
+  double s;
+  
+  for (i=0;i<n;i++)
+    {
+      s = 0.;
+      for (j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+       {
+         p = cols.colnums[j];
+         if (p>=0)
+           {
+             if (i>j) s += val[j] * dst(p);
+           }
+       }
+      dst(i) -= s * om;
+      dst(i) /= val[cols.rowstart[i]];
+    }
+
+  for (i=n-1;i>=0;i--)
+    {
+      s = 0.;
+      for (j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+       {
+         p = cols.colnums[j];
+         if (p>=0)
+           {
+             if (i<j) s += val[j] * dst(p);
+           }
+       }
+      dst(i) -= s * om / val[cols.rowstart[i]];
+    }
+}
+
+
+void dSMatrix::print (ostream &out) const {
+  for (int i=0; i<cols.rows; ++i)
+    for (int j=cols.rowstart[i]; j<cols.rowstart[i+1]; ++j)
+      out << "(" << i << "," << cols.colnums[j] << ") " << val[j] << endl;
+};
diff --git a/deal.II/lac/source/dvector.cc b/deal.II/lac/source/dvector.cc
new file mode 100644 (file)
index 0000000..1a7302e
--- /dev/null
@@ -0,0 +1,259 @@
+// $Id$
+
+// This file is part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+
+static const char* OBJFILE = "DEAL $RCSfile$ $Revision$";
+
+#include <lac/dvector.h>
+#include <math.h>
+
+dVector::dVector() : val(0)
+{
+  dim=maxdim=0;
+}
+
+dVector::dVector(int n) : val(0)
+{
+  dim = n;
+  maxdim = n;
+  THROW1(n<0, IntError(IntError::IllegalDimension, n, "dVector"));
+
+  if (n)
+    {
+      val = new double[maxdim];
+      THROWUNCOND(!val, Error(Error::NoMem,"dVector::dVector"));
+      (*this) = 0.;
+    }
+}
+
+dVector::dVector(const dVector& v) : val(0)
+{
+  dim = v.n();
+  maxdim = v.n();
+
+  if (dim)
+    {
+      val = new double[maxdim];
+      THROWUNCOND(!val, Error(Error::NoMem,"dVector::dVector"));
+      for (int i=0;i<dim;i++) val[i] = v.val[i];
+    }
+}
+
+void dVector::reinit(int n, int fast)
+{
+  THROW1(n<=0, IntError(IntError::IllegalDimension, n));
+  if (n>maxdim)
+  {
+    if (val) delete[] val;
+    val = new double[n];
+    THROWUNCOND(!val, Error(Error::NoMem,"dVector::reinit"));
+    maxdim = n;
+  }
+  dim = n;
+  if (!fast) (*this) = 0.;
+}
+
+void dVector::reinit(const dVector& v, int fast)
+{
+  int n = v.n();
+  if (n>maxdim)
+  {
+    if (val) delete[] val;
+    val = new double[n];
+    THROWUNCOND(!val, Error(Error::NoMem,"dVector::reinit"));
+    maxdim = n;
+  }
+  dim = n;
+  if (!fast) (*this) = 0.;
+}
+
+dVector::~dVector()
+{
+  if (val) delete[] val;
+}
+
+double dVector::operator * (const dVector& v) const
+{
+  int i;
+  double s;
+  for (i = 0, s = 0.; i < dim; i++)
+  {
+        s += val[i] * v.val[i];
+  }
+  return s;
+}
+
+void dVector::add(const double v)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] += v;
+}
+
+void dVector::add(const dVector& v)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] += v(i);
+}
+
+void dVector::add(double a, const dVector& v)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] += a * v(i);
+}
+
+void dVector::add(double a, const dVector& v, double b, const dVector& w)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] += a * v.val[i] + b * w.val[i];
+}
+
+void dVector::sadd(double x, const dVector& v)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] = x * val[i] + v.val[i];
+}
+
+void dVector::sadd(double x, double a, const dVector& v)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] = x * val[i] + a * v.val[i];
+}
+
+void dVector::sadd(double x, double a, const dVector& v, double b, const dVector& w)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] = x * val[i] + a * v.val[i] + b * w.val[i];
+}
+
+void dVector::sadd(double x, double a, const dVector& v, 
+                  double b, const dVector& w, double c, const dVector& y)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] = x * val[i] + a * v.val[i] + b * w.val[i] 
+                             + c * y.val[i];
+}
+
+void dVector::equ(double a, const dVector& u, double b, const dVector& v)
+{
+  int i;
+  for (i=0;i<dim;i++) val[i] = a*u.val[i] + b*v.val[i];
+}
+
+void dVector::equ(double a)
+{
+  int i;
+  for (i=0;i<dim;i++) val[i] *= a;
+}
+
+void dVector::equ(double a, const dVector& u)
+{
+  int i;
+  for (i=0;i<dim;i++) val[i] = a*u.val[i];
+}
+
+dVector& dVector::operator = (double s)
+{
+  int i;
+  for (i=0;i<dim;i++) val[i] = s;
+  return *this;
+}
+
+dVector& dVector::operator = (const dVector& v)
+{
+  if (v.dim != dim) reinit(v,1);
+
+  int i;
+  for (i=0;i<dim;i++) val[i] = v.val[i];
+  return *this;
+}
+
+void dVector::cadd(int i, const VectorBase& V, double s, int j)
+{
+  const dVector& v = (const dVector&) V;
+  THROW1((i<0) || (i>dim), IntError(IntError::Range,i,"cadd"));
+  THROW1((j<0) || (j>v.dim), IntError(IntError::Range,j,"cadd"));
+
+  val[i] += s*v.val[j];
+}
+
+void dVector::cadd(int i, const VectorBase& V, double s, int j, double t, int k)
+{
+  const dVector& v = (const dVector&) V;
+  THROW1((i<0) || (i>dim), IntError(IntError::Range,i,"cadd"));
+  THROW1((j<0) || (j>v.dim), IntError(IntError::Range,j,"cadd"));
+  THROW1((k<0) || (k>v.dim), IntError(IntError::Range,k,"cadd"));
+
+  val[i] += s*v.val[j] + t*v.val[k];
+}
+
+void dVector::cadd(int i, const VectorBase& V, double s, int j,
+                  double t, int k, double q, int l, double r, int m)
+{
+  const dVector& v = (const dVector&) V;
+  THROW1((i<0) || (i>dim), IntError(IntError::Range,i,"cadd"));
+  THROW1((j<0) || (j>v.dim), IntError(IntError::Range,j,"cadd"));
+  THROW1((k<0) || (k>v.dim), IntError(IntError::Range,k,"cadd"));
+  THROW1((l<0) || (l>v.dim), IntError(IntError::Range,l,"cadd"));
+  THROW1((m<0) || (m>v.dim), IntError(IntError::Range,m,"cadd"));
+
+  val[i] += s*v.val[j] + t*v.val[k] + q*v.val[l] + r*v.val[m];
+}
+
+void dVector::czero(int i)
+{
+    THROW1((i<0) || (i>dim), IntError(IntError::Range,i));
+    val[i] = 0.;
+}
+
+void dVector::cequ(int i, const VectorBase& V, double s, int j)
+{
+  const dVector& v = (const dVector&) V;
+  THROW1((i<0) || (i>dim), IntError(IntError::Range,i,"cequ"));
+  THROW1((j<0) || (j>v.dim), IntError(IntError::Range,j,"cequ"));
+
+  val[i] = s*v.val[j];
+}
+
+void dVector::cequ(int i, const VectorBase& V, double s, int j, double t, int k)
+{
+  const dVector& v = (const dVector&) V;
+  THROW1((i<0) || (i>dim), IntError(IntError::Range,i,"cequ"));
+  THROW1((j<0) || (j>v.dim), IntError(IntError::Range,j,"cequ"));
+  THROW1((k<0) || (k>v.dim), IntError(IntError::Range,k,"cequ"));
+
+  val[i] = s*v.val[j] + t*v.val[k];
+}
+
+void dVector::cequ(int i, const VectorBase& V, double s, int j,
+                  double t, int k, double q, int l, double r, int m)
+{
+  const dVector& v = (const dVector&) V;
+  THROW1((i<0) || (i>dim), IntError(IntError::Range,i,"cequ"));
+  THROW1((j<0) || (j>v.dim), IntError(IntError::Range,j,"cequ"));
+  THROW1((k<0) || (k>v.dim), IntError(IntError::Range,k,"cequ"));
+  THROW1((l<0) || (l>v.dim), IntError(IntError::Range,l,"cequ"));
+  THROW1((m<0) || (m>v.dim), IntError(IntError::Range,m,"cequ"));
+
+  val[i] = s*v.val[j] + t*v.val[k] + q*v.val[l] + r*v.val[m];
+}
+
+const char* dVector::name() const
+{
+  return "dVector";
+}
+
+void dVector::print(FILE* f, const char* format) const
+{
+  if (!format) format = " %5.2f";
+  for (int j=0;j<n();j++) fprintf(f, format, val[j]);
+  fputc('\n',f);
+}
+
+void dVector::print(const char* format) const
+{
+  if (!format) format = " %5.2f";
+  for (int j=0;j<n();j++) printf (format, val[j]);
+  printf ("\n");
+}
diff --git a/deal.II/lac/source/ivector.cc b/deal.II/lac/source/ivector.cc
new file mode 100644 (file)
index 0000000..4e72421
--- /dev/null
@@ -0,0 +1,103 @@
+// $Id$
+
+// This file is part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+
+static const char* OBJFILE = "DEAL $RCSfile$ $Revision$";
+
+#include <lac/ivector.h>
+
+iVector::iVector()
+{
+  dim=maxdim=1;
+  val = new int[1];
+  THROWUNCOND(!val, Error(Error::NoMem,"iVector::iVector"));
+  val[0] = 0;
+}
+
+iVector::iVector(int n)
+{
+  dim = n;
+  maxdim = n;
+  THROW1(n<=0, IntError(IntError::IllegalDimension, n, "iVector::iVector"));
+
+  val = new int[maxdim];
+  THROWUNCOND(!val, Error(Error::NoMem,"iVector::iVector"));
+  (*this) = 0;
+}
+
+iVector::iVector(const iVector& v)
+{
+  reinit(v.dim,1);
+  int i;
+  for (i=0;i<dim;i++) val[i] = v.val[i];
+}
+
+void iVector::reinit(int n, int fast)
+{
+  THROW1(n<=0, IntError(IntError::IllegalDimension, n, "iVector::reinit"));
+  if (n>maxdim)
+  {
+    delete[] val;
+    val = new int[n];
+    THROWUNCOND(!val, Error(Error::NoMem,"iVector::reinit"));
+    maxdim = n;
+  }
+  dim = n;
+  if (!fast) (*this) = 0;
+}
+
+void iVector::reinit(const iVector& v, int fast)
+{
+  int n = v.n();
+  if (n>maxdim)
+  {
+    delete[] val;
+    val = new int[n];
+    THROWUNCOND(!val, Error(Error::NoMem,"iVector::reinit"));
+    maxdim = n;
+  }
+  dim = n;
+  if (!fast) (*this) = 0;
+}
+
+iVector::~iVector()
+{
+  delete[] val;
+}
+
+void iVector::add(const iVector& v)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] += v(i);
+}
+
+void iVector::add(int a, const iVector& v)
+{
+  int i;
+  for (i = 0; i < dim; i++) val[i] += a * v(i);
+}
+
+void iVector::equ(int a, const iVector& u)
+{
+  int i;
+  for (i=0;i<dim;i++) val[i] = a*u.val[i];
+}
+
+iVector& iVector::operator = (int s)
+{
+  int i;
+//  if (s==0.) memset(val,sizeof(*val) * dim, 0);
+  for (i=0;i<dim;i++) val[i] = s;
+  return *this;
+}
+
+iVector& iVector::operator = (const iVector& v)
+{
+  if (v.dim != dim) reinit(v,1);
+
+  int i;
+  for (i=0;i<dim;i++) val[i] = v.val[i];
+  return *this;
+}

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.