/**
* Dimension. Actual number of Columns
*/
- int dim_range;
+ unsigned int dim_range;
/**
* Dimension. Actual number of Rows
*/
- int dim_image;
+ unsigned int dim_image;
/**
* Dimension. Determines amount of reserved memory
*/
- int val_size;
+ unsigned int val_size;
/**
* Initialization . initialize memory for Matrix <p>
* ( m rows , n columns )
*/
- void init(int m, int n);
+ void init (const unsigned int m, const unsigned int n);
/**
* Access Elements. returns A(i,j)
*/
- double& el(int i, int j) { return val[i*dim_range+j]; }
+ double& el (const unsigned int i, const unsigned int j) {
+ return val[i*dim_range+j];
+ };
/**
* Access Elements. returns A(i,j)
*/
- double el(int i, int j) const { return val[i*dim_range+j]; }
+ double el (const unsigned int i, const unsigned int j) const {
+ return val[i*dim_range+j];
+ };
public:
* Constructor. Dimension = (n,n) <p>
* -> quadratic matrix (n rows , n columns)
*/
- dFMatrix(int n = 1) { init(n,n); }
+ dFMatrix (const unsigned int n = 1) { init(n,n); }
/**
* Constructor. Dimension = (m,n) <p>
* -> rectangular matrix (m rows , n columns)
*/
- dFMatrix(int m,int n) { init(m,n); }
+ dFMatrix (const unsigned int m,unsigned int n) { init(m,n); }
/**
* Copy constructor. Be very careful with
* huge amount of computing time for large
* matrices!!
*/
- dFMatrix(const dFMatrix&);
+ dFMatrix (const dFMatrix&);
/**
* Destructor. Clears memory
/**
* U(0-m,0-n) = s . Fill all elements
*/
- void fill(const dFMatrix& src, int i = 0, int j = 0);
+ void fill (const dFMatrix& src,
+ const unsigned int i=0, const unsigned int j=0);
/**
* Change Dimension.
* Set dimension to (m,n) <p>
* ( reinit rectangular matrix )
*/
- void reinit(int m, int n);
+ void reinit (const unsigned int m, const unsigned int n);
/**
* Change Dimension.
* Set dimension to (n,n) <p>
* ( reinit quadratic matrix )
*/
- void reinit(int n) { reinit(n,n); }
+ void reinit (const unsigned int n) { reinit(n,n); }
/**
* Adjust Dimension.
* Set dimension to ( m(B),n(B) ) <p>
* ( adjust to dimensions of another matrix B )
*/
- void reinit(const dFMatrix& B) { reinit(B.m(), B.n()); }
+ void reinit (const dFMatrix& B) { reinit(B.m(), B.n()); }
/**
* Inquire Dimension (Row) . returns Number of Rows
*/
- int m() const { return dim_image; }
+ unsigned int m() const { return dim_image; }
/**
* Inquire Dimension (Col) . returns Number of Columns
*/
- int n() const { return dim_range; }
+ unsigned int n () const { return dim_range; }
//@}
* Access Elements. returns element at relative 'address' i <p>
* ( -> access to A(i/n , i mod n) )
*/
- double el(int i) const { return val[i]; }
+ double el (const unsigned int i) const { return val[i]; }
/**
* Access Elements. returns A(i,j)
*/
- double operator() (int i, int j) const
+ double operator() (const unsigned int i, const unsigned int j) const
{
- Assert ((i>=0) && (i<dim_image), ExcInvalidIndex (i, dim_image));
- Assert ((j>=0) && (j<dim_range), ExcInvalidIndex (i, dim_range));
+ Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
+ Assert (j<dim_range, ExcInvalidIndex (i, dim_range));
return el(i,j);
}
/**
* Access Elements. returns A(i,j)
*/
- double& operator() (int i, int j)
+ double& operator() (const unsigned int i, const unsigned int j)
{
- Assert ((i>=0) && (i<dim_image), ExcInvalidIndex (i, dim_image));
- Assert ((j>=0) && (j<dim_range), ExcInvalidIndex (i, dim_range));
+ Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
+ Assert (j<dim_range, ExcInvalidIndex (i, dim_range));
return el(i,j);
}
/**
* A+=B . Simple addition
*/
- void add(double s, const dFMatrix& B);
+ void add (const double s, const dFMatrix& B);
/**
* A+=Transp(B).
* Simple addition of the transpose of B to this
*/
- void Tadd(double s, const dFMatrix& B);
+ void Tadd (const double s, const dFMatrix& B);
/**
* C=A*B.
* Matrix-matrix-multiplication
*/
- void mmult(dFMatrix& C, const dFMatrix& B) const;
+ void mmult (dFMatrix& C, const dFMatrix& B) const;
/**
* C=Transp(A)*B.
* Matrix-matrix-multiplication using
* transpose of this
*/
- void Tmmult(dFMatrix& C, const dFMatrix& B) const;
+ void Tmmult (dFMatrix& C, const dFMatrix& B) const;
/**
* w (+)= A*v.
* ( application of this to a vector v )
* flag adding=true : w=+A*v
*/
- void vmult(dVector& w, const dVector& v,const int adding = 0) const;
-
+ void vmult (dVector& w, const dVector& v, const bool adding=false) const;
+
/**
* w (+)= Transp(A)*v.
* Matrix-vector-multiplication ; <p>
* (application of transpose of this to a vector v)
* flag adding=true : w=+A*v
*/
- void Tvmult(dVector& w,const dVector& v,const int adding=0) const;
+ void Tvmult (dVector& w, const dVector& v, const bool adding=false) const;
/**
* A=Inverse(A). Inversion of this by
* Gauss-Jordan-algorithm
*/
- void gauss_jordan();
+ void gauss_jordan ();
/**
* Computes the determinant of a matrix.
* A(i,1-n)+=s*A(j,1-n).
* Simple addition of rows of this
*/
- void add_row(int i, double s, int j);
+ void add_row (const unsigned int i, const double s, const unsigned int j);
/**
* A(i,1-n)+=s*A(j,1-n)+t*A(k,1-n).
* Multiple addition of rows of this
*/
- void add_row(int i, double s, int j, double t, int k);
+ void add_row (const unsigned int i,
+ const double s, const unsigned int j,
+ const double t, const unsigned int k);
/**
* A(1-n,i)+=s*A(1-n,j).
* Simple addition of columns of this
*/
- void add_col(int i, double s, int j);
+ void add_col (const unsigned int i, const double s, const unsigned int j);
/**
* A(1-n,i)+=s*A(1-n,j)+t*A(1-n,k).
* Multiple addition of columns of this
*/
- void add_col(int i, double s, int j, double t, int k);
+ void add_col (const unsigned int i,
+ const double s, const unsigned int j,
+ const double t, const unsigned int k);
/**
* Swap A(i,1-n) <-> A(j,1-n).
* Swap rows i and j of this
*/
- void swap_row(int i, int j);
+ void swap_row (const unsigned int i, const unsigned int j);
/**
* Swap A(1-n,i) <-> A(1-n,j).
* Swap columns i and j of this
*/
- void swap_col(int i, int j);
+ void swap_col (const unsigned int i, const unsigned int j);
//@}
* w=b-A*v.
* Residual calculation , returns |w|
*/
- double residual(dVector& w, const dVector& v, const dVector& b) const;
+ double residual (dVector& w, const dVector& v, const dVector& b) const;
/**
* Inversion of lower triangle .
*/
- void forward(dVector& dst, const dVector& src) const;
+ void forward (dVector& dst, const dVector& src) const;
/**
* Inversion of upper triangle .
*/
- void backward(dVector& dst, const dVector& src) const;
+ void backward (dVector& dst, const dVector& src) const;
/**
* QR - factorization of a matrix.
* triangle contains the resulting matrix R, <p>
* the lower the incomplete factorization matrices.
*/
- void householder(dVector& y);
+ void householder (dVector& y);
/**
* Least - Squares - Approximation by QR-factorization.
*/
- double least_squares(dVector& dst, dVector& src);
+ double least_squares (dVector& dst, dVector& src);
/**
* A(i,i)+=B(i,1-n). Addition of complete
* rows of B to diagonal-elements of this ; <p>
* ( i = 1 ... m )
*/
- void add_diag(double s, const dFMatrix& B);
+ void add_diag (const double s, const dFMatrix& B);
/**
* A(i,i)+=s i=1-m.
* Add constant to diagonal elements of this
*/
- void diagadd(const double& src);
+ void diagadd (const double& src);
/**
* w+=part(A)*v. Conditional partial
* Matrix-vector-multiplication <p>
* (used elements of v determined by x)
*/
- void gsmult(dVector& w, const dVector& v,const iVector& x) const;
+ void gsmult (dVector& w, const dVector& v, const iVector& x) const;
/**
* Output of the matrix in user-defined format.
*/
- void print(FILE* fp, const char* format = 0) const;
+ void print (FILE* fp, const char* format = 0) const;
//@}
/**
class dSMatrixStruct
{
friend class dSMatrix;
- int max_dim;
- int rows, cols;
- int vec_len, max_vec_len;
- int max_row_len;
- int* rowstart;
+ unsigned int max_dim;
+ unsigned int rows, cols;
+ unsigned int vec_len, max_vec_len;
+ unsigned int max_row_len;
+ unsigned int* rowstart;
int* colnums;
- int compressed;
+ bool compressed;
public:
//////////
- void reinit(int m, int n, int max_per_row);
+ void reinit (const unsigned int m, const unsigned int n,
+ const unsigned int max_per_row);
//////////
- dSMatrixStruct(int m, int n, int max_per_row);
+ dSMatrixStruct (const unsigned int m, const unsigned int n,
+ const unsigned int max_per_row);
//////////
- dSMatrixStruct(int n, int max_per_row);
+ dSMatrixStruct (const unsigned int n, const unsigned int max_per_row);
//////////
- ~dSMatrixStruct();
+ ~dSMatrixStruct ();
//////////
- void compress();
+ void compress ();
//////////
- int operator () (int i, int j);
+ int operator() (const unsigned int i, const unsigned int j);
//////////
- void add(int i, int j);
+ void add (const unsigned int i, const unsigned int j);
//////////
- void add_matrix(int n, int* rowcols);
+ void add_matrix (const unsigned int n, const int* rowcols);
//////////
- void add_matrix(int m, int n, int* rows, int* cols);
+ void add_matrix (const unsigned int m, const unsigned int n,
+ const int* rows, const int* cols);
//////////
- void add_matrix(const iVector& rowcols);
+ void add_matrix (const iVector& rowcols);
//////////
- void add_matrix(const iVector& rows, const iVector& cols);
+ void add_matrix (const iVector& rows, const iVector& cols);
void print_gnuplot (ostream &) const;
- int n_rows () const;
- int n_cols () const;
- int bandwidth () const;
+ unsigned int n_rows () const;
+ unsigned int n_cols () const;
+ unsigned int bandwidth () const;
/**
* Return whether the structure is
* avoid programs relying on outdated
* information!
*/
- const int * get_rowstart_indices () const;
+ const unsigned int * get_rowstart_indices () const;
/**
* This is kind of an expert mode: get
{
dSMatrixStruct * cols;
double* val;
- int max_len;
+ unsigned int max_len;
public:
/**
dSMatrix ();
//
- dSMatrix(dSMatrixStruct& c);
+ dSMatrix (dSMatrixStruct& c);
//
- ~dSMatrix();
+ ~dSMatrix ();
//
- void reinit();
+ void reinit ();
//
void reinit (dSMatrixStruct &);
void clear ();
//
- int m() const;
+ unsigned int m () const;
//
- int n() const;
+ unsigned int n () const;
//
- void set (int i, int j, double value);
+ void set (const unsigned int i, const unsigned int j,
+ const double value);
//
- void add (int i, int j, double value);
+ void add (const unsigned int i, const unsigned int j,
+ const double value);
/**
* Return the value of the entry (i,j).
* throws an exception if the wanted
* element does not exist in the matrix.
*/
- double operator () (const int i, const int j) const;
+ double operator () (const unsigned int i, const unsigned int j) const;
/**
* Return the main diagonal element in
* involve searching for the right column
* number.
*/
- double diag_element (const int i) const;
+ double diag_element (const unsigned int i) const;
/**
* This is kind of an expert mode: get
* avoid programs relying on outdated
* information!
*/
- double global_entry (const int i) const;
+ double global_entry (const unsigned int i) const;
/**
* Same as above, but with write access.
* You certainly know what you do?
*/
- double & global_entry (const int i);
+ double & global_entry (const unsigned int i);
//
void vmult (dVector& dst,const dVector& src) const;
//
- void Tvmult(dVector& dst,const dVector& src) const;
+ void Tvmult (dVector& dst,const dVector& src) const;
//
- double residual (dVector& dst,const dVector& x,const dVector& b);
+ double residual (dVector& dst, const dVector& x, const dVector& b);
//
- void Jacobi_precond(dVector& dst,const dVector& src,double om = 1.);
+ void Jacobi_precond (dVector& dst, const dVector& src,
+ const double om = 1.);
//
- void SSOR_precond(dVector& dst,const dVector& src,double om = 1.);
+ void SSOR_precond (dVector& dst, const dVector& src,
+ const double om = 1.);
//
- void SOR_precond(dVector& dst,const dVector& src,double om = 1.);
+ void SOR_precond (dVector& dst, const dVector& src,
+ const double om = 1.);
//
- void SSOR(dVector& dst,double om = 1.);
+ void SSOR (dVector& dst, const double om = 1.);
//
- void SOR(dVector& dst,double om = 1.);
+ void SOR (dVector& dst, const double om = 1.);
//
- void precondition(dVector& dst,const dVector& src) { dst=src; }
+ void precondition (dVector& dst, const dVector& src) { dst=src; }
const dSMatrixStruct & get_sparsity_pattern () const;
/*---------------------- Inline functions -----------------------------------*/
inline
-int dSMatrixStruct::n_rows () const {
+unsigned int dSMatrixStruct::n_rows () const {
return rows;
};
inline
-int dSMatrixStruct::n_cols () const {
+unsigned int dSMatrixStruct::n_cols () const {
return cols;
};
inline
-const int * dSMatrixStruct::get_rowstart_indices () const {
+const unsigned int * dSMatrixStruct::get_rowstart_indices () const {
return rowstart;
};
inline
-int dSMatrix::m () const
+unsigned int dSMatrix::m () const
{
return cols->rows;
};
inline
-int dSMatrix::n () const
+unsigned int dSMatrix::n () const
{
return cols->cols;
};
inline
-void dSMatrix::set (int i, int j, double value) {
+void dSMatrix::set (const unsigned int i, const unsigned int j,
+ const double value) {
Assert (cols->operator()(i,j) != -1,
ExcInvalidIndex(i,j));
val[cols->operator()(i,j)] = value;
inline
-void dSMatrix::add (int i, int j, double value) {
+void dSMatrix::add (const unsigned int i, const unsigned int j,
+ const double value) {
Assert (cols->operator()(i,j) != -1,
ExcInvalidIndex(i,j));
val[cols->operator()(i,j)] += value;
inline
-double dSMatrix::operator () (const int i, const int j) const {
+double dSMatrix::operator () (const unsigned int i, const unsigned int j) const {
Assert (cols->operator()(i,j) != -1,
ExcInvalidIndex(i,j));
return val[cols->operator()(i,j)];
inline
-double dSMatrix::diag_element (const int i) const {
+double dSMatrix::diag_element (const unsigned int i) const {
Assert (m() == n(), ExcMatrixNotSquare());
- Assert ((0<=i) && (i<max_len), ExcInvalidIndex1(i));
+ Assert (i<max_len, ExcInvalidIndex1(i));
// Use that the first element in each
// row of a square matrix is the main
inline
-double dSMatrix::global_entry (const int j) const {
+double dSMatrix::global_entry (const unsigned int j) const {
return val[j];
};
inline
-double & dSMatrix::global_entry (const int j) {
+double & dSMatrix::global_entry (const unsigned int j) {
return val[j];
};
protected:
/// Dimension. Actual number of components
- int dim;
+ unsigned int dim;
/// Dimension. Determines amount of reserved memory , evtl. >DIM !
- int maxdim;
+ unsigned int maxdim;
/// Component-array.
double *val;
/**
* Dummy-Constructor. Dimension=0
*/
- dVector();
+ dVector ();
/**
* Copy-Constructor. Dimension set to that of V , <p>
* all components are copied from V
*/
- dVector(const dVector& V);
+ dVector (const dVector& V);
/**
* Constructor. Dimension = N (>0)
*/
- dVector(int N);
+ dVector (const unsigned int N);
/**
* Destructor. Clears memory
*/
- ~dVector();
+ ~dVector ();
/**
* Set all entries to zero. Equivalent to
/**
* U(0-N) = s . Fill all components
*/
- dVector& operator=(double s);
+ dVector& operator= (const double s);
/**
* U = V . Copy all components
*/
- dVector& operator=(const dVector& V);
+ dVector& operator= (const dVector& V);
/**
* U = U * V . Scalar Produkt
*/
- double operator*(const dVector& V) const;
+ double operator* (const dVector& V) const;
/**
* Return square of the norm.
* Change Dimension. <p>
* Set dimension to N <p>
* ! reserved memory for This remains unchanged ! <p>
- * on fast=0 vector is filled by 0.
+ * on fast==false vector is filled by 0.
*/
- void reinit(int N, int fast = 0);
+ void reinit (const unsigned int N, const bool fast=false);
/**
* Adjust Dimension. <p>
* Set dimension to n(V) <p>
* ! reserved memory for This remains unchanged ! <p>
* ! components of V are not copied in any case ! <p>
- * on fast=0 vector is filled by 0.
+ * on fast==false vector is filled by 0.
*/
- void reinit(const dVector& V, int fast = 0);
+ void reinit (const dVector& V, const bool fast=false);
/**
* Inquire Dimension. returns Dimension ,
* INLINE
*/
- int n() const;
+ unsigned int n () const;
//@}
* Access Components. returns U(i) ,
* INLINE
*/
- double operator()(int i) const;
+ double operator() (const unsigned int i) const;
/**
* Access Components. returns U(i) ,
* INLINE
*/
- double& operator()(int i);
+ double& operator() (const unsigned int i);
//@}
/**
* U(0-DIM)+=s . Addition of S to all components
*/
- void add(const double s);
+ void add (const double s);
/**
* U+=V . Simple addition
*/
- void add(const dVector& V);
+ void add (const dVector& V);
/**
* U+=a*V . Simple addition
*/
- void add(double a, const dVector& V);
+ void add (const double a, const dVector& V);
/**
* U+=a*V+b*W . Multiple addition
*/
- void add(double a, const dVector& V, double b, const dVector& W);
+ void add (const double a, const dVector& V,
+ const double b, const dVector& W);
/**
* U=s*U+V . Scaling + simple addition
*/
- void sadd(double s, const dVector& V);
+ void sadd (const double s, const dVector& V);
/**
* U=s*U+a*V . Scaling + simple addition
*/
- void sadd(double s, double a, const dVector& V);
+ void sadd (const double s, const double a, const dVector& V);
/**
* U=s*U+a*V+b*W . Scaling + multiple addition
*/
- void sadd(double s, double a, const dVector& V, double b, const dVector& W);
+ void sadd (const double s, const double a,
+ const dVector& V, const double b, const dVector& W);
/**
* U=s*U+a*V+b*W+c*X. Scaling + multiple addition
*/
- void sadd(double s, double a, const dVector& V, double b, const dVector& W,
- double c, const dVector& X);
+ void sadd (const double s, const double a,
+ const dVector& V, const double b, const dVector& W,
+ const double c, const dVector& X);
/**
* U=s*U . Scaling
*/
- void equ(double s);
+ void equ (const double s);
/**
* U=a*V . Replacing
*/
- void equ(double a, const dVector& V);
+ void equ (const double a, const dVector& V);
/**
* U=a*V+b*W . Replacing by sum
*/
- void equ(double a, const dVector& V, double b, const dVector& W);
+ void equ (const double a, const dVector& V,
+ const double b, const dVector& W);
//@}
/**
* U(i)=0 . ONE Component only
*/
- void czero(int i);
+ void czero (const unsigned int i);
/**
* U(i)=a*V(j) . Replacing
*/
- void cequ(int i, const VectorBase& V, double a, int j);
+ void cequ(unsigned int i, const VectorBase& V,
+ double a, unsigned int j);
/**
* U(i)=a*V(j)+b*V(k) . Replacing by sum
*/
- void cequ(int i, const VectorBase& V, double a, int j, double b, int k);
+ void cequ (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k);
/**
* U(i)=a*V(j)+b*V(k)+c*V(l)+d*V(m) . Replacing by sum
*/
- void cequ(int i, const VectorBase& V, double a, int j, double b,
- int k, double c, int l, double d, int m);
+ void cequ (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k,
+ const double c, const unsigned int l,
+ const double d, const unsigned int m);
/**
* U(i)+=a*V(j) . Simple addition
*/
- void cadd(int i, const VectorBase& V, double a, int j);
+ void cadd (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j);
/**
* U(i)+=a*V(j)+b*V(k). Multiple addition
*/
- void cadd(int i, const VectorBase& V, double a, int j, double b, int k);
+ void cadd (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k);
/**
* U(i)+=a*V(j)+b*V(k)+c*V(l)+d*V(m) . Multiple addition
*/
- void cadd(int i, const VectorBase& V, double a, int j, double b,
- int k, double c, int l, double d, int m);
+ void cadd (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k,
+ const double c, const unsigned int l,
+ const double d, const unsigned int m);
//@}
*/
//@{
///
- virtual const char* name() const;
+ virtual const char* name () const;
/**
* Output of vector in user-defined format.
*/
- void print(FILE* fp, const char* format = 0) const;
+ void print (FILE* fp, const char* format = 0) const;
/**
* Output of vector in user-defined format.
*/
- void print(const char* format = 0) const;
+ void print (const char* format = 0) const;
/**
* Print to given stream, one element per line.
-inline int dVector::n() const
+inline unsigned int dVector::n () const
{
return dim;
}
-inline double dVector::operator() (int i) const
+inline double dVector::operator() (const unsigned int i) const
{
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
return val[i];
}
-inline double& dVector::operator() (int i)
+inline double& dVector::operator() (const unsigned int i)
{
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
return val[i];
}
protected:
/// Dimension. Actual number of components
- int dim;
+ unsigned int dim;
/// Dimension. Determines amount of reserved memory , evtl. >DIM !
- int maxdim;
+ unsigned int maxdim;
/// Component-array.
int *val;
/**
* Dummy-Constructor. Dimension=1
*/
- iVector();
+ iVector ();
/**
* Copy-Constructor. Dimension set to that of V , <p>
* all components are copied from V
*/
- iVector(const iVector& V);
+ iVector (const iVector& V);
/**
* Constructor. Dimension = N (>0)
*/
- iVector(int N);
+ iVector (const unsigned int N);
/**
* Destructor. Clears memory
*/
- ~iVector();
+ ~iVector ();
/**
* U(0-N) = s . Fill all components
*/
- iVector& operator=(int s);
+ iVector& operator = (const int s);
/**
* U = V . Copy all components
*/
- iVector& operator=(const iVector& V);
+ iVector& operator = (const iVector& V);
/**
* Change Dimension. <p>
* Set dimension to N <p>
* ! reserved memory for This remains unchanged ! <p>
- * on fast=0 vector is filled by 0.
+ * on fast==false vector is filled by 0.
*/
- void reinit(int N, int fast = 0);
+ void reinit (const unsigned int N, const bool fast=false);
/**
* Adjust Dimension. <p>
* Set dimension to n(V) <p>
* ! reserved memory for This remains unchanged ! <p>
* ! components of V are not copied in any case ! <p>
- * on fast=0 vector is filled by 0.
+ * on fast==false vector is filled by 0.
*/
- void reinit(const iVector& V, int fast = 0);
+ void reinit (const iVector& V, const bool fast=false);
/**
* Set all elements to zero.
* Inquire Dimension. returns Dimension ,
* INLINE
*/
- int n() const;
+ unsigned int n () const;
//@}
* Access Components. returns U(i) ,
* INLINE
*/
- int operator()(int i) const;
+ int operator() (const unsigned int i) const;
/**
* Access Components. returns U(i) ,
* INLINE
*/
- int& operator()(int i);
+ int& operator() (const unsigned int i);
//@}
/**
* U+=V . Simple addition
*/
- void add(const iVector& V);
+ void add (const iVector& V);
/**
* U+=a*V . Simple addition
*/
- void add(int a, const iVector& V);
+ void add (const unsigned int a, const iVector& V);
/**
* U=a*V . Replacing
*/
- void equ(int a, const iVector& V);
+ void equ (const unsigned int a, const iVector& V);
//@}
/**
-inline int iVector::n() const
+inline unsigned int iVector::n() const
{
return dim;
}
-inline int iVector::operator() (int i) const
+inline int iVector::operator() (unsigned int i) const
{
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
return val[i];
}
-inline int& iVector::operator() (int i)
+inline int& iVector::operator() (unsigned int i)
{
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
return val[i];
}
/**
* U(i)=0 . ONE Component only
*/
- virtual void czero(int i) = 0;
+ virtual void czero (const unsigned int i) = 0;
/**
* U(i)=a*V(j) . Replacing
*/
- virtual void cequ(int i, const VectorBase& V, double a, int j) = 0;
+ virtual void cequ (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j) = 0;
/**
* U(i)=a*V(j)+b*V(k) . Replacing by sum
*/
- virtual void cequ(int i, const VectorBase& V, double a, int j, double b, int k) = 0;
+ virtual void cequ (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k) = 0;
/**
* U(i)=a*V(j)+b*V(k)+c*V(l)+d*V(m) . Replacing by sum
*/
- virtual void cequ(int i, const VectorBase& V, double a, int j, double b,
- int k, double c, int l, double d, int m) = 0;
+ virtual void cequ (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k,
+ const double c, const unsigned int l,
+ const double d, const unsigned int m) = 0;
/**
* U(i)+=a*V(j) . Simple addition
*/
- virtual void cadd(int i, const VectorBase& V, double a, int j) = 0;
+ virtual void cadd (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j) = 0;
/**
* U(i)+=a*V(j)+b*V(k). Multiple addition
*/
- virtual void cadd(int i, const VectorBase& V, double a, int j, double b, int k) = 0;
+ virtual void cadd (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k) = 0;
/**
* U(i)+=a*V(j)+b*V(k)+c*V(l)+d*V(m) . Multiple addition
*/
- virtual void cadd(int i, const VectorBase& V, double a, int j, double b,
- int k, double c, int l, double d, int m) = 0;
+ virtual void cadd (const unsigned int i, const VectorBase& V,
+ const double a, const unsigned int j,
+ const double b, const unsigned int k,
+ const double c, const unsigned int l,
+ const double d, const unsigned int m) = 0;
//@}
CXX = c++
INCLUDE = -I../include -I../../base/include
-CXXFLAGS.g= -DDEBUG -g -Wall -pedantic -Wconversion \
- -Winline -Woverloaded-virtual $(INCLUDE)
+CXXFLAGS.g= -DDEBUG -g -Wall -W -pedantic -Wconversion \
+ -Winline -Woverloaded-virtual -fno-rtti -fno-exceptions \
+ $(INCLUDE)
CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \
- -funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE)
+ -funroll-loops -felide-constructors -fnonnull-objects \
+ -fno-rtti -fno-exceptions $(INCLUDE)
+
+ifeq ($(shell uname),Linux)
+CXX = /home/wolf/bin/gcc-2.8.1/bin/c++
+endif
ifeq ($(shell uname),SunOS)
-#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS)
-#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
-INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \
- -I/usr/local/source/libstdc++-2.8.0/libio
-CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8)
-CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8)
+INCLUDE2.8= -I/usr/local/source/libstdc++-2.8.0/libstdc++ \
+ -I/usr/local/source/libstdc++-2.8.0/libio
+CXXFLAGS := $(CXXFLAGS) $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) $(INCLUDE2.8)
endif
+
cc-files = $(wildcard *.cc)
o-files = $(cc-files:.cc=.o)
go-files = $(cc-files:.cc=.go)
lib.g.a: ../lib/liblac.g.a($(go-files))
clean:
- rm -f *.o *.go *~ Makefile.dep
+ rm -f *.o *.go *~ Makefile.dep ../lib/lib*
dFMatrix::dFMatrix (const dFMatrix &m)
{
init (m.dim_image, m.dim_range);
- for (int i=0; i!=dim_image*dim_range; ++i)
+ for (unsigned int i=0; i!=dim_image*dim_range; ++i)
val[i] = m.val[i];
};
-void dFMatrix::init(int mm, int nn)
+void dFMatrix::init (const unsigned int mm, const unsigned int nn)
{
val_size = nn*mm;
val = new double[val_size];
clear ();
}
-dFMatrix::~dFMatrix()
+dFMatrix::~dFMatrix ()
{
delete[] val;
}
-void dFMatrix::reinit(int mm, int nn)
+void dFMatrix::reinit (const unsigned int mm, const unsigned int nn)
{
if (val_size<nn*mm)
{
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;
+ for (unsigned int i=0; i<nn*mm; ++i) val[i] = 0;
}
}
-void dFMatrix::vmult(dVector& dst, const dVector& src,const int adding) const
+void dFMatrix::vmult (dVector& dst, const dVector& src,
+ const bool adding) const
{
Assert(dst.n() == m(), ExcDimensionMismatch(dst.n(), m()));
Assert(src.n() == n(), ExcDimensionMismatch(src.n(), n()));
else
{
double* e = val;
- for (int i=0;i<m();i++)
+ const unsigned int size_m = m(),
+ size_n = n();
+ for (unsigned int i=0; i<size_m; ++i)
{
s = 0.;
- for (int j=0;j<n();j++) s += src(j) * *(e++);
+ for (unsigned int j=0; j<size_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
+void dFMatrix::gsmult (dVector& dst, const dVector& src, const iVector& gl) const
{
Assert(n() == m(), ExcNotQuadratic());
Assert(dst.n() == n(), ExcDimensionMismatch(dst.n(), n()));
else
{
double* e = val;
- for (int i=0;i<m();i++)
+ const unsigned int size_m = m(),
+ size_n = n();
+ for (unsigned int i=0; i<size_m; ++i)
{
s = 0.;
- for (int j=0;j<n();j++) if(gl(i)<gl(j)) s += src(j) * *(e++);
+ for (unsigned int j=0; j<size_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
+void dFMatrix::Tvmult (dVector& dst, const dVector& src, const bool adding) const
{
Assert(dst.n() == n(), ExcDimensionMismatch(dst.n(), n()));
Assert(src.n() == m(), ExcDimensionMismatch(src.n(), m()));
- int i,j;
+ unsigned int i,j;
double s;
- for (i=0;i<m();i++)
+ const unsigned int size_m = m(),
+ size_n = n();
+ for (i=0; i<size_m; ++i)
{
s = 0.;
- for (j=0;j<n();j++) s += src(j) * el(j,i);
+ for (j=0; j<size_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,
+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;
+ unsigned int i,j;
double s, res = 0.;
- for (i=0;i<n();i++)
+ const unsigned int size_m = m(),
+ size_n = n();
+ for (i=0; i<size_n; ++i)
{
s = right(i);
- for (j=0;j<m();j++) s -= src(j) * el(i,j);
+ for (j=0; j<size_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
+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());
+ unsigned int i,j;
+ unsigned int nu = MIN(m(),n());
double s;
- for (i=0;i<nu;i++)
+ for (i=0; i<nu; ++i)
{
s = src(i);
- for (j=0;j<i;j++) s -= dst(j) * el(i,j);
+ 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
+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());
+ unsigned int j;
+ unsigned int nu = MIN(m(),n());
double s;
- for (i=nu-1;i>=0;i--)
+ for (int i=nu-1; i>=0; --i)
{
s = src(i);
- for (j=i+1;j<nu;j++) s -= dst(j) * el(i,j);
+ 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];
+ unsigned int nn = n()*m();
+ for (unsigned int i=0; i<nn; ++i)
+ val[i] = M.val[i];
return *this;
}
-void dFMatrix::fill(const dFMatrix& src, int i, int j)
+void dFMatrix::fill (const dFMatrix& src,
+ const unsigned int i, const unsigned int j)
{
Assert (n() >= src.n() + j, ExcInvalidDestination(n(), src.n(), j));
Assert (m() >= src.m() + i, ExcInvalidDestination(m(), src.m(), i));
- for (int ii=0; ii<src.m() ; ii++)
- for (int jj=0; jj<src.n() ; jj++)
+ for (unsigned int ii=0; ii<src.m() ; ++ii)
+ for (unsigned 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)
+void dFMatrix::add_row (const unsigned int i,
+ const double s, const unsigned int j)
{
- int k;
- for (k=0;k<m();k++) el(i,k) += s*el(j,k);
+ for (unsigned int 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)
+void dFMatrix::add_row (const unsigned int i, const double s,
+ const unsigned int j, const double t,
+ const unsigned int k)
{
- int l;
- for (l=0;l<m();l++) el(i,l) += s*el(j,l) + t*el(k,l);
+ const unsigned int size_m = m();
+ for (unsigned l=0; l<size_m; ++l)
+ el(i,l) += s*el(j,l) + t*el(k,l);
}
-void dFMatrix::add_col(int i, double s, int j)
+void dFMatrix::add_col (const unsigned int i, const double s,
+ const unsigned int j)
{
- int k;
- for (k=0;k<n();k++) el(k,i) += s*el(k,j);
+ for (unsigned int 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)
+void dFMatrix::add_col (const unsigned int i, const double s,
+ const unsigned int j, const double t,
+ const unsigned int k)
{
- int l;
- for (l=0;l<n();l++) el(l,i) += s*el(l,j) + t*el(l,k);
+ for (unsigned int l=0; l<n(); ++l)
+ el(l,i) += s*el(l,j) + t*el(l,k);
}
-void dFMatrix::swap_row(int i, int j)
+void dFMatrix::swap_row (const unsigned int i, const unsigned int j)
{
- int k;
double s;
- for (k=0;k<m();k++)
+ for (unsigned int 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)
+void dFMatrix::swap_col (const unsigned int i, const unsigned int j)
{
- int k;
double s;
- for (k=0;k<n();k++)
+ for (unsigned int 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)
+void dFMatrix::diagadd (const double& src)
{
Assert (m() == n(), ExcDimensionMismatch(m(),n()));
- for (int i=0;i<n();i++)
+ for (unsigned int i=0; i<n(); ++i)
el(i,i) += src;
}
-void dFMatrix::mmult(dFMatrix& dst, const dFMatrix& src) const
+void dFMatrix::mmult (dFMatrix& dst, const dFMatrix& src) const
{
Assert (n() == src.m(), ExcDimensionMismatch(n(), src.m()));
- int i,j,k;
+ unsigned 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++)
+ for (j=0; j<src.n(); ++j)
{
s = 0.;
for (k=0;k<n();k++) s+= el(i,k) * src.el(k,j);
}
}
-/*void dFMatrix::mmult(dFMatrix& dst, const dFMatrix& src) const
+/*void dFMatrix::mmult (dFMatrix& dst, const dFMatrix& src) const
{
Assert (m() == src.n(), ExcDimensionMismatch(m(), src.n()));
- int i,j,k;
+ unsigned int i,j,k;
double s = 1.;
dst.reinit(n(), src.m());
}
}*/
-void dFMatrix::Tmmult(dFMatrix& dst, const dFMatrix& src) const
+void dFMatrix::Tmmult (dFMatrix& dst, const dFMatrix& src) const
{
Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
- int i,j,k;
+ unsigned int i,j,k;
double s = 1.;
dst.reinit(m(), src.m());
{
Assert (m() == src.n(), ExcDimensionMismatch(m(), src.n()));
- int i,j,k;
+ unsigned int i,j,k;
double s = 1.;
dst.reinit(n(), src.m());
}
}*/
-void dFMatrix::print(FILE* f, const char* format) const
+void dFMatrix::print (FILE* f, const char* format) const
{
if (!format) format = " %5.2f";
- int i,j;
+ unsigned int i,j;
for (i=0;i<m();i++)
{
for (j=0;j<n();j++) fprintf(f, format, el(i,j));
}
}
-void dFMatrix::add(double s,const dFMatrix& src)
+void dFMatrix::add (const double s,const dFMatrix& src)
{
Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m()));
Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
}
else
{
- for (int i = n()*m()-1; i>=0 ; i--) val[i] += s * src.el(i);
+ const unsigned int size = n()*m();
+ for (unsigned int i=0; i<size; i++)
+ val[i] += s * src.el(i);
}
}
-void dFMatrix::add_diag(double s,const dFMatrix& src)
+void dFMatrix::add_diag (const double s, const dFMatrix& src)
{
Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m()));
Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
}
else
{
- for (int i = n()*m()-1; i>=0 ; i--) val[i] += s * src.el(i);
+ const unsigned int size = n()*m();
+ for (unsigned int i=0; i<size; i++)
+ val[i] += s * src.el(i);
}
}
-void dFMatrix::Tadd(double s,const dFMatrix& src)
+void dFMatrix::Tadd (const double s, const dFMatrix& src)
{
Assert (m() == n(), ExcNotQuadratic());
Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m()));
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)
+ for (unsigned int i=0; i<dim_image; ++i)
+ for (unsigned int j=0; j<dim_range; ++j)
if (el(i,j) != m.el(i,j)) return false;
return true;
};
void dFMatrix::clear () {
- for (int i=0; i<val_size; ++i)
+ for (unsigned int i=0; i<val_size; ++i)
val[i] = 0.;
};
//////////
template<class T>
-inline void swap(T* a, T* b)
+inline void swap (T* a, T* b)
{
T x = *a;
*a = *b;
//////////
template<class T>
-inline void simple_sort(long n, T* field)
+inline void simple_sort (const long n, T* field)
{
long i,j;
for (i=1;i<n;i++)
}
template<class T>
-inline void heapsort_sift(T* a, long l, long r)
+inline void heapsort_sift (T* a, const long l, const long r)
{
long i = l;
long j = 2*i;
//////////
template<class T>
-inline void heapsort(int n, T* field)
+inline void heapsort (const int n, T* field)
{
field--;
//////////
template<class T>
-inline void _quicksort(long r, T* a, long l)
+inline void _quicksort (const long r, T* a, const long l)
{
long i = l;
long j = r;
}
template<class T>
-inline void quicksort(long r, T* a)
+inline void quicksort (const long r, T* a)
{
_quicksort(r,a,1);
}
-dSMatrixStruct::dSMatrixStruct(int m, int n, int max_per_row)
+dSMatrixStruct::dSMatrixStruct (const unsigned int m, const unsigned int n,
+ const unsigned int max_per_row)
: max_dim(0),
max_vec_len(0),
rowstart(0),
-dSMatrixStruct::dSMatrixStruct(int n, int max_per_row)
+dSMatrixStruct::dSMatrixStruct (const unsigned int n,
+ const unsigned int max_per_row)
: max_dim(0),
max_vec_len(0),
rowstart(0),
void
-dSMatrixStruct::reinit(int m, int n, int max_per_row)
+dSMatrixStruct::reinit (const unsigned int m, const unsigned int n,
+ const unsigned int max_per_row)
{
Assert (m>0, ExcInvalidNumber(m));
Assert (n>0, ExcInvalidNumber(n));
{
if (rowstart) delete[] rowstart;
max_dim = rows;
- rowstart = new int[max_dim+1];
+ rowstart = new unsigned int[max_dim+1];
}
if (vec_len > max_vec_len)
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;
+
+ for (unsigned int i=0; i<=rows; i++)
+ rowstart[i] = i * max_per_row;
+ for (int i = vec_len-1; i>=0; i--)
+ colnums[i] = -1;
if (rows == cols)
- for (i=0;i<rows;i++)
+ for (unsigned int i=0;i<rows;i++)
colnums[rowstart[i]] = i;
compressed = 0;
void
-dSMatrixStruct::compress()
+dSMatrixStruct::compress ()
{
if (compressed) return;
- int i,j,k,l,s;
+ unsigned 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++)
+ 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)
+ for (j=rowstart[i],l=0; j<rowstart[i+1]; ++j,++l)
entries[l] = colnums[j];
heapsort(max_row_len, entries);
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];
- }
+ if (entries[l]>=0)
+ if (((signed int)i!=entries[l]) || (rows!=cols))
+ colnums[k++] = entries[l];
+
rowstart[i] = s;
s = k;
}
vec_len = rowstart[i] = s;
- compressed = 1;
+ compressed = true;
delete[] entries;
}
int
-dSMatrixStruct::operator () (int i, int j)
+dSMatrixStruct::operator () (const unsigned int i, const unsigned int j)
{
- for (int k=rowstart[i] ; k<rowstart[i+1] ; k++)
- if (colnums[k] == j) return k;
+ for (unsigned int k=rowstart[i] ; k<rowstart[i+1] ; k++)
+ if (colnums[k] == (signed int)j) return k;
return -1;
}
void
-dSMatrixStruct::add (int i, int j)
+dSMatrixStruct::add (const unsigned int i, const unsigned int j)
{
- Assert ((i>=0) && (i<rows), ExcInvalidIndex(i,rows));
- Assert ((j>=0) && (j<cols), ExcInvalidIndex(j,cols));
+ Assert (i<rows, ExcInvalidIndex(i,rows));
+ Assert (j<cols, ExcInvalidIndex(j,cols));
- for (int k=rowstart[i]; k<rowstart[i+1]; k++)
+ for (unsigned int k=rowstart[i]; k<rowstart[i+1]; k++)
{
// entry already exists
- if (colnums[k] == j) return;
+ if (colnums[k] == (signed int)j) return;
// empty entry found, put new
// entry here
if (colnums[k] == -1)
void
-dSMatrixStruct::add_matrix(int n, int* rowcols)
+dSMatrixStruct::add_matrix (const unsigned int n, const int* rowcols)
{
- int i,j;
+ unsigned 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)
+dSMatrixStruct::add_matrix (const unsigned int m, const unsigned int n,
+ const int* rows, const int* cols)
{
- int i,j;
+ unsigned 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)
+dSMatrixStruct::add_matrix (const iVector& rowcols)
{
- int i,j;
+ unsigned 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)
+dSMatrixStruct::add_matrix (const iVector& rows, const iVector& cols)
{
- int i,j;
+ unsigned int i,j;
for (i=0;i<rows.n();i++)
for (j=0;j<cols.n();j++)
add(rows(i), cols(j));
void
dSMatrixStruct::print_gnuplot (ostream &out) const
{
- for (int i=0; i<rows; ++i)
- for (int j=rowstart[i]; j<rowstart[i+1]; ++j)
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
if (colnums[j]>=0)
out << i << " " << -colnums[j] << endl;
}
-int
+unsigned int
dSMatrixStruct::bandwidth () const
{
- int b=0;
- for (int i=0; i<rows; ++i)
- for (int j=rowstart[i]; j<rowstart[i+1]; ++j)
+ unsigned int b=0;
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
if (colnums[j]>=0)
{
- if (abs(i-colnums[j]) > b)
+ if ((unsigned int)abs(i-colnums[j]) > b)
b = abs(i-colnums[j]);
}
else
void
-dSMatrix::reinit()
+dSMatrix::reinit ()
{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert (cols->compressed, ExcNotCompressed());
void
-dSMatrix::vmult(dVector& dst,const dVector& src) const
+dSMatrix::vmult (dVector& dst, const dVector& src) const
{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert(m() == dst.n(), ExcDimensionsDontMatch(m(),dst.n()));
Assert(n() == src.n(), ExcDimensionsDontMatch(n(),src.n()));
- for (int i=0;i<m();i++)
+ for (unsigned int i=0;i<m();i++)
{
double s = 0.;
- for (int j=cols->rowstart[i]; j < cols->rowstart[i+1] ;j++)
+ for (unsigned int j=cols->rowstart[i]; j < cols->rowstart[i+1] ;j++)
{
int p = cols->colnums[j];
s += val[j] * src(p);
void
-dSMatrix::Tvmult(dVector& dst,const dVector& src) const
+dSMatrix::Tvmult (dVector& dst, const dVector& src) const
{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert(n() == dst.n(), ExcDimensionsDontMatch(n(),dst.n()));
Assert(m() == src.n(), ExcDimensionsDontMatch(m(),src.n()));
- int i;
+ unsigned 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++)
+ for (unsigned 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)
+dSMatrix::residual (dVector& dst, const dVector& u, const dVector& b)
{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert(m() == dst.n(), ExcDimensionsDontMatch(m(),dst.n()));
double s,norm=0.;
- for (int i=0;i<m();i++)
+ for (unsigned int i=0;i<m();i++)
{
s = b(i);
- for (int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
+ for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
int p = cols->colnums[j];
s -= val[j] * u(p);
}
void
-dSMatrix::Jacobi_precond(dVector& dst,const dVector& src,double om)
+dSMatrix::Jacobi_precond (dVector& dst, const dVector& src, const double om)
{
Assert (cols != 0, ExcMatrixNotInitialized());
- int n = src.n();
+ const unsigned int n = src.n();
- for (int i=0;i<n;++i)
+ for (unsigned 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)
+dSMatrix::SSOR_precond (dVector& dst, const dVector& src, const double om)
{
Assert (cols != 0, ExcMatrixNotInitialized());
- int p,n = src.n();
- int i,j;
+ int p;
+ unsigned int n = src.n();
+ unsigned int j;
- for (i=0;i<n;i++)
+ for (unsigned 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);
+ if (p<(signed int)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 (unsigned int i=0; i<n; i++)
+ dst(i) *= (2.-om)*val[cols->rowstart[i]];
- for (i=n-1;i>=0;i--)
+ for (int i=n-1; i>=0; i--)
{
for (j=cols->rowstart[i];j<cols->rowstart[i+1];j++)
{
}
void
-dSMatrix::SOR_precond(dVector& dst,const dVector& src,double om)
+dSMatrix::SOR_precond (dVector& dst, const dVector& src, const double om)
{
Assert (cols != 0, ExcMatrixNotInitialized());
dst = src;
}
void
-dSMatrix::SOR(dVector& dst, double om)
+dSMatrix::SOR (dVector& dst, const double om)
{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert(n() == m(), ExcDimensionsDontMatch(n(),m()));
Assert(m() == dst.n(), ExcDimensionsDontMatch(m(),dst.n()));
- for (int i=0;i<m();i++)
+ for (unsigned int i=0;i<m();i++)
{
double s = dst(i);
- for (int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
+ for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
int p = cols->colnums[j];
- if (p<i)
+ if (p<(signed int)i)
s -= val[j] * dst(p);
}
dst(i) = s * om / val[cols->rowstart[i]];
}
void
-dSMatrix::SSOR(dVector& dst, double om)
+dSMatrix::SSOR (dVector& dst, const double om)
{
Assert (cols != 0, ExcMatrixNotInitialized());
- int p,n = dst.n();
- int i,j;
+ int p;
+ unsigned int n = dst.n();
+ unsigned int j;
double s;
- for (i=0;i<n;i++)
+ for (unsigned int i=0; i<n; i++)
{
s = 0.;
for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
dst(i) /= val[cols->rowstart[i]];
}
- for (i=n-1;i>=0;i--)
+ for (int i=n-1; i>=0; i--) // this time, i is signed, but alsways positive!
{
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);
+ if ((unsigned int)i<j) s += val[j] * dst(p);
}
}
dst(i) -= s * om / val[cols->rowstart[i]];
void dSMatrix::print (ostream &out) const {
Assert (cols != 0, ExcMatrixNotInitialized());
- for (int i=0; i<cols->rows; ++i)
- for (int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
+ for (unsigned int i=0; i<cols->rows; ++i)
+ for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << endl;
};
#include <lac/dvector.h>
#include <math.h>
-dVector::dVector() :
+dVector::dVector () :
dim(0),
maxdim(0),
val(0)
{}
-dVector::dVector(int n) :
+dVector::dVector (const unsigned int n) :
dim(n),
maxdim(n),
val(0)
}
-dVector::dVector(const dVector& v) :
+dVector::dVector (const dVector& v) :
VectorBase(v),
dim(v.n()),
maxdim(v.n()),
{
val = new double[maxdim];
Assert (val != 0, ExcOutOfMemory());
- for (int i=0; i<dim; i++)
+ for (unsigned int i=0; i<dim; i++)
val[i] = v.val[i];
}
}
-void dVector::reinit(int n, int fast)
+void dVector::reinit (const unsigned int n, const bool fast)
{
Assert (n>0, ExcInvalidNumber(n));
if (n>maxdim)
-void dVector::reinit(const dVector& v, int fast)
+void dVector::reinit (const dVector& v, const bool fast)
{
- int n = v.n();
+ const unsigned int n = v.n();
if (n>maxdim)
{
if (val) delete[] val;
-dVector::~dVector()
+dVector::~dVector ()
{
if (val) delete[] val;
}
void dVector::clear () {
- for (int i=0; i<dim; ++i)
+ for (unsigned int i=0; i<dim; ++i)
val[i] = 0.;
}
// use modern processors better by
// allowing pipelined commands to be
// executed in parallel
- for (int i=0; i<(dim/4); ++i)
+ for (unsigned int i=0; i<(dim/4); ++i)
{
sum0 += val[4*i] * v.val[4*i];
sum1 += val[4*i+1] * v.val[4*i+1];
sum3 += val[4*i+3] * v.val[4*i+3];
};
// add up remaining elements
- for (int i=(dim/4)*4; i<dim; ++i)
+ for (unsigned int i=(dim/4)*4; i<dim; ++i)
sum0 += val[i] * v.val[i];
return sum0+sum1+sum2+sum3;
// use modern processors better by
// allowing pipelined commands to be
// executed in parallel
- for (int i=0; i<(dim/4); ++i)
+ for (unsigned int i=0; i<(dim/4); ++i)
{
sum0 += val[4*i] * val[4*i];
sum1 += val[4*i+1] * val[4*i+1];
sum3 += val[4*i+3] * val[4*i+3];
};
// add up remaining elements
- for (int i=(dim/4)*4; i<dim; ++i)
+ for (unsigned int i=(dim/4)*4; i<dim; ++i)
sum0 += val[i] * val[i];
return sum0+sum1+sum2+sum3;
-void dVector::add(const double v)
+void dVector::add (const double v)
{
- for (int i = 0; i < dim; i++) val[i] += v;
+ for (unsigned int i = 0; i < dim; i++) val[i] += v;
}
-void dVector::add(const dVector& v)
+void dVector::add (const dVector& v)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
- for (int i = 0; i < dim; i++) val[i] += v(i);
+ for (unsigned int i = 0; i < dim; i++) val[i] += v(i);
}
-void dVector::add(double a, const dVector& v)
+void dVector::add (const double a, const dVector& v)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
- for (int i = 0; i < dim; i++) val[i] += a * v(i);
+ for (unsigned int i = 0; i < dim; i++) val[i] += a * v(i);
}
-void dVector::add(double a, const dVector& v, double b, const dVector& w)
+void dVector::add (const double a, const dVector& v,
+ const double b, const dVector& w)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
- for (int i = 0; i < dim; i++) val[i] += a * v.val[i] + b * w.val[i];
+ for (unsigned int i = 0; i < dim; i++)
+ val[i] += a * v.val[i] + b * w.val[i];
}
-void dVector::sadd(double x, const dVector& v)
+void dVector::sadd (const double x, const dVector& v)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
- for (int i = 0; i < dim; i++) val[i] = x * val[i] + v.val[i];
+ for (unsigned int i = 0; i < dim; i++)
+ val[i] = x * val[i] + v.val[i];
}
-void dVector::sadd(double x, double a, const dVector& v)
+void dVector::sadd (const double x, const double a, const dVector& v)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
- for (int i = 0; i < dim; i++) val[i] = x * val[i] + a * v.val[i];
+ for (unsigned int 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)
+void dVector::sadd (const double x, const double a,
+ const dVector& v, const double b, const dVector& w)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
- for (int i = 0; i < dim; i++) val[i] = x * val[i] + a * v.val[i] + b * w.val[i];
+ for (unsigned int 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)
+void dVector::sadd (const double x, const double a,
+ const dVector& v, const double b,
+ const dVector& w, const double c, const dVector& y)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
Assert (dim == y.dim, ExcDimensionsDontMatch(dim, y.dim));
- for (int i = 0; i < dim; i++)
+ for (unsigned int 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)
+void dVector::equ (const double a, const dVector& u,
+ const double b, const dVector& v)
{
Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
- for (int i=0;i<dim;i++) val[i] = a*u.val[i] + b*v.val[i];
+ for (unsigned int i=0; i<dim; i++)
+ val[i] = a*u.val[i] + b*v.val[i];
}
-void dVector::equ(double a)
+void dVector::equ (const double a)
{
- for (int i=0;i<dim;i++) val[i] *= a;
+ for (unsigned int i=0; i<dim; i++)
+ val[i] *= a;
}
-void dVector::equ(double a, const dVector& u)
+void dVector::equ (const double a, const dVector& u)
{
Assert (dim == u.dim, ExcDimensionsDontMatch(dim, u.dim));
- for (int i=0;i<dim;i++) val[i] = a*u.val[i];
+ for (unsigned int i=0; i<dim; i++)
+ val[i] = a*u.val[i];
}
-dVector& dVector::operator = (double s)
+dVector& dVector::operator = (const double s)
{
- for (int i=0;i<dim;i++) val[i] = s;
+ for (unsigned int i=0; i<dim; i++)
+ val[i] = s;
return *this;
}
{
if (v.dim != dim) reinit(v,1);
- for (int i=0;i<dim;i++) val[i] = v.val[i];
+ for (unsigned int i=0; i<dim; i++)
+ val[i] = v.val[i];
return *this;
}
-void dVector::cadd(int i, const VectorBase& V, double s, int j)
+void dVector::cadd (const unsigned int i, const VectorBase& V,
+ const double s, const unsigned int j)
{
const dVector& v = (const dVector&) V;
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
- Assert ((j>=0) || (j<v.dim), ExcInvalidIndex(j,v.dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ Assert (j<v.dim, ExcInvalidIndex(j,v.dim));
val[i] += s*v.val[j];
}
-void dVector::cadd(int i, const VectorBase& V, double s, int j, double t, int k)
+void dVector::cadd(unsigned int i, const VectorBase& V,
+ double s, unsigned int j,
+ double t, unsigned int k)
{
const dVector& v = (const dVector&) V;
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
- Assert ((j>=0) || (j<v.dim), ExcInvalidIndex(j,v.dim));
- Assert ((k>=0) || (k<v.dim), ExcInvalidIndex(k,v.dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ Assert (j<v.dim, ExcInvalidIndex(j,v.dim));
+ Assert (k<v.dim, ExcInvalidIndex(k,v.dim));
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)
+void dVector::cadd (const unsigned int i, const VectorBase& V,
+ const double s, const unsigned int j,
+ const double t, const unsigned int k,
+ const double q, const unsigned int l,
+ const double r, const unsigned int m)
{
const dVector& v = (const dVector&) V;
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
- Assert ((j>=0) || (j<v.dim), ExcInvalidIndex(j,v.dim));
- Assert ((k>=0) || (k<v.dim), ExcInvalidIndex(k,v.dim));
- Assert ((l>=0) || (l<v.dim), ExcInvalidIndex(l,v.dim));
- Assert ((m>=0) || (m<v.dim), ExcInvalidIndex(m,v.dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ Assert (j<v.dim, ExcInvalidIndex(j,v.dim));
+ Assert (k<v.dim, ExcInvalidIndex(k,v.dim));
+ Assert (l<v.dim, ExcInvalidIndex(l,v.dim));
+ Assert (m<v.dim, ExcInvalidIndex(m,v.dim));
val[i] += s*v.val[j] + t*v.val[k] + q*v.val[l] + r*v.val[m];
}
-void dVector::czero(int i)
+void dVector::czero (const unsigned int i)
{
- Assert ((i>=0) && (i<dim), ExcInvalidIndex(i,dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
val[i] = 0.;
}
-void dVector::cequ(int i, const VectorBase& V, double s, int j)
+void dVector::cequ (const unsigned int i, const VectorBase& V,
+ const double s, const unsigned int j)
{
const dVector& v = (const dVector&) V;
- Assert ((i>=0) && (i<dim), ExcInvalidIndex(i,dim));
- Assert ((j>=0) && (j<v.dim), ExcInvalidIndex(j,v.dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ Assert (j<v.dim, ExcInvalidIndex(j,v.dim));
val[i] = s*v.val[j];
}
-void dVector::cequ(int i, const VectorBase& V, double s, int j, double t, int k)
+void dVector::cequ (const unsigned int i, const VectorBase& V,
+ const double s, const unsigned int j,
+ const double t, const unsigned int k)
{
const dVector& v = (const dVector&) V;
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
- Assert ((j>=0) || (j<v.dim), ExcInvalidIndex(j,v.dim));
- Assert ((k>=0) || (k<v.dim), ExcInvalidIndex(k,v.dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ Assert (j<v.dim, ExcInvalidIndex(j,v.dim));
+ Assert (k<v.dim, ExcInvalidIndex(k,v.dim));
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)
+void dVector::cequ (const unsigned int i, const VectorBase& V,
+ const double s, const unsigned int j,
+ const double t, const unsigned int k,
+ const double q, const unsigned int l,
+ const double r, const unsigned int m)
{
const dVector& v = (const dVector&) V;
- Assert ((i>=0) || (i<dim), ExcInvalidIndex(i,dim));
- Assert ((j>=0) || (j<v.dim), ExcInvalidIndex(j,v.dim));
- Assert ((k>=0) || (k<v.dim), ExcInvalidIndex(k,v.dim));
- Assert ((l>=0) || (l<v.dim), ExcInvalidIndex(l,v.dim));
- Assert ((m>=0) || (m<v.dim), ExcInvalidIndex(m,v.dim));
+ Assert (i<dim, ExcInvalidIndex(i,dim));
+ Assert (j<v.dim, ExcInvalidIndex(j,v.dim));
+ Assert (k<v.dim, ExcInvalidIndex(k,v.dim));
+ Assert (l<v.dim, ExcInvalidIndex(l,v.dim));
+ Assert (m<v.dim, ExcInvalidIndex(m,v.dim));
val[i] = s*v.val[j] + t*v.val[k] + q*v.val[l] + r*v.val[m];
}
-const char* dVector::name() const
+const char* dVector::name () const
{
return "dVector";
}
-void dVector::print(FILE* f, const char* format) const
+void dVector::print (FILE* f, const char* format) const
{
if (!format) format = " %5.2f";
- for (int j=0;j<n();j++)
+ for (unsigned int j=0;j<n();j++)
fprintf(f, format, val[j]);
fputc('\n',f);
}
-void dVector::print(const char* format) const
+void dVector::print (const char* format) const
{
if (!format) format = " %5.2f";
- for (int j=0;j<n();j++)
+ for (unsigned int j=0;j<n();j++)
printf (format, val[j]);
printf ("\n");
}
void dVector::print (ostream &out) const {
- for (int i=0; i<n(); ++i)
+ for (unsigned int i=0; i<n(); ++i)
out << val[i] << endl;
};
-iVector::iVector() :
+iVector::iVector () :
dim(1),
maxdim(1),
val(new int[1])
}
-iVector::iVector(int n) :
+iVector::iVector (const unsigned int n) :
dim(n),
maxdim(n),
val(new int[n])
}
-iVector::iVector(const iVector& v) :
+iVector::iVector (const iVector& v) :
dim(v.dim),
maxdim(v.maxdim),
val(0)
{
reinit(v.dim,1);
- int i;
- for (i=0; i<dim; i++)
+ for (unsigned int i=0; i<dim; i++)
val[i] = v.val[i];
}
-void iVector::reinit(int n, int fast)
+void iVector::reinit (const unsigned int n, const bool fast)
{
Assert (n>0, ExcInvalidNumber(n));
if (n>maxdim)
-void iVector::reinit(const iVector& v, int fast)
+void iVector::reinit (const iVector& v, const bool fast)
{
- int n = v.n();
+ const unsigned int n = v.n();
if (n>maxdim)
{
delete[] val;
-iVector::~iVector()
+iVector::~iVector ()
{
if (val)
delete[] val;
void iVector::clear () {
- for (int i=0; i<dim; ++i)
+ for (unsigned int i=0; i<dim; ++i)
val[i] = 0;
}
-void iVector::add(const iVector& v)
+void iVector::add (const iVector& v)
{
- for (int i=0; i<dim; ++i) val[i] += v(i);
+ for (unsigned int i=0; i<dim; ++i) val[i] += v(i);
}
-void iVector::add(int a, const iVector& v)
+void iVector::add (const unsigned int a, const iVector& v)
{
- for (int i=0; i<dim; ++i) val[i] += a * v(i);
+ for (unsigned int i=0; i<dim; ++i) val[i] += a * v(i);
}
-void iVector::equ(int a, const iVector& u)
+void iVector::equ (const unsigned int a, const iVector& u)
{
- for (int i=0; i<dim; ++i) val[i] = a*u.val[i];
+ for (unsigned int i=0; i<dim; ++i) val[i] = a*u.val[i];
}
-iVector& iVector::operator = (int s)
+iVector& iVector::operator = (const int s)
{
// if (s==0.) memset(val,sizeof(*val) * dim, 0);
- for (int i=0; i<dim; ++i) val[i] = s;
+ for (unsigned int i=0; i<dim; ++i) val[i] = s;
return *this;
}
{
if (v.dim != dim) reinit(v,1);
- for (int i=0; i<dim; ++i) val[i] = v.val[i];
+ for (unsigned int i=0; i<dim; ++i) val[i] = v.val[i];
return *this;
}