From b3b93f92038fe37d2f8f2091997eb8dcda859954 Mon Sep 17 00:00:00 2001 From: Guido Kanschat <dr.guido.kanschat@gmail.com> 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 <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 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 <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 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 <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 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 <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 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 <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 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 <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 index 0000000000..f55f02f804 --- /dev/null +++ b/deal.II/lac/source/dsmatrix.cc @@ -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 index 0000000000..1a7302e526 --- /dev/null +++ b/deal.II/lac/source/dvector.cc @@ -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 index 0000000000..4e7242151b --- /dev/null +++ b/deal.II/lac/source/ivector.cc @@ -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; +} -- 2.39.5