From b3b93f92038fe37d2f8f2091997eb8dcda859954 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 4 Mar 1998 16:59:32 +0000 Subject: [PATCH] Initial revision git-svn-id: https://svn.dealii.org/trunk@11 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/types.h | 37 ++ deal.II/lac/include/lac/dfmatrix.h | 158 +++++ deal.II/lac/include/lac/dsmatrix.h | 140 +++++ deal.II/lac/include/lac/dvector.h | 100 ++++ deal.II/lac/include/lac/ivector.h | 83 +++ deal.II/lac/include/lac/vectorbase.h | 37 ++ deal.II/lac/source/dfmatrix.cc | 862 +++++++++++++++++++++++++++ deal.II/lac/source/dsmatrix.cc | 361 +++++++++++ deal.II/lac/source/dvector.cc | 259 ++++++++ deal.II/lac/source/ivector.cc | 103 ++++ 10 files changed, 2140 insertions(+) create mode 100644 deal.II/base/include/base/types.h create mode 100644 deal.II/lac/include/lac/dfmatrix.h create mode 100644 deal.II/lac/include/lac/dsmatrix.h create mode 100644 deal.II/lac/include/lac/dvector.h create mode 100644 deal.II/lac/include/lac/ivector.h create mode 100644 deal.II/lac/include/lac/vectorbase.h create mode 100644 deal.II/lac/source/dfmatrix.cc create mode 100644 deal.II/lac/source/dsmatrix.cc create mode 100644 deal.II/lac/source/dvector.cc create mode 100644 deal.II/lac/source/ivector.cc diff --git a/deal.II/base/include/base/types.h b/deal.II/base/include/base/types.h new file mode 100644 index 0000000000..fe77a49a17 --- /dev/null +++ b/deal.II/base/include/base/types.h @@ -0,0 +1,37 @@ +// $Id$ + +#ifndef __base_types_h +#define __base_types_h + +#ifndef _stdlib_h +#include +#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 index 0000000000..1682fa36c9 --- /dev/null +++ b/deal.II/lac/include/lac/dfmatrix.h @@ -0,0 +1,158 @@ +// $Id$ + +#ifndef __lac_dfmatrix_h +#define __lac_dfmatrix_h + +#ifndef __stdio_h +#include +#endif +#ifndef __lac_dvector_h +#include +#endif +#ifndef __lac_ivector_h +#include +#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 diff --git a/deal.II/lac/include/lac/dsmatrix.h b/deal.II/lac/include/lac/dsmatrix.h new file mode 100644 index 0000000000..351d552f9b --- /dev/null +++ b/deal.II/lac/include/lac/dsmatrix.h @@ -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 +#ifndef __lac_dvector_h +#include +#endif +#ifndef __lac_ivector_h +#include +#endif + + +class ostream; +template 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 index 0000000000..43cb079ab4 --- /dev/null +++ b/deal.II/lac/include/lac/dvector.h @@ -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 +#ifndef __base_types_h +#include +#endif +//#ifndef __base_errors_h +//#include +//#endif +#ifndef __lac_vectorbase_h +#include +#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 index 0000000000..c5586793b5 --- /dev/null +++ b/deal.II/lac/include/lac/ivector.h @@ -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 +#endif +//#ifndef __base_errors_h +//#include +//#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 index 0000000000..a619fd133d --- /dev/null +++ b/deal.II/lac/include/lac/vectorbase.h @@ -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 index 0000000000..c7dd222200 --- /dev/null +++ b/deal.II/lac/source/dfmatrix.cc @@ -0,0 +1,862 @@ +// $Id$ + +static const char* OBJFILE = "DEAL $RCSfile$ $Revision$"; + +#include +#include +#include +#include +#include + +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=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=0;i--) + { + s = src(i); + for (j=i+1;j=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 +#include +#include + +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= 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=rows), IntError(IntError::Range,i)); + //THROW1((j<0) || (j>=cols), IntError(IntError::Range,j)); + int k; + for (k=rowstart[i]; k=0 ; i--) val[i] = 0; +} + +void +dSMatrixStruct::print_gnuplot (ostream &out) const +{ + for (int i=0; i=0) + out << i << " " << -colnums[j] << endl; +} + + +int +dSMatrixStruct::bandwidth () const +{ + int b=0; + for (int i=0; i=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=0;i--) + { + for (j=cols.rowstart[i];ji) 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=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=0) + { + if (i +#include + +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;imaxdim) + { + 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;idim), 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 + +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;imaxdim) + { + 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