--- /dev/null
+// $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
--- /dev/null
+// $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ß-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
--- /dev/null
+// $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
--- /dev/null
+// $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
--- /dev/null
+// $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
--- /dev/null
+// $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
--- /dev/null
+// $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;
+};
--- /dev/null
+// $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;
+};
--- /dev/null
+// $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");
+}
--- /dev/null
+// $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;
+}