template<typename number> class Vector;
-class iVector;
-
/**
*/
void diagadd (const number s);
- /**
- * w+=part(A)*v. Conditional partial
- * Matrix-vector-multiplication <p>
- * (used elements of v determined by x)
- */
- template<typename number2>
- void gsmult (Vector<number2>& w, const Vector<number2>& v, const iVector& x) const;
-
-
/**
* Output of the matrix in user-defined format.
*/
// $Id$
#include <lac/vector.h>
-#include <lac/ivector.h>
#include <lac/fullmatrix.h>
+#include <vector>
#include <cmath>
#include <cstdlib>
#include <cstdio>
template <typename number>
template <typename number2>
-void FullMatrix<number>::gsmult (Vector<number2>& dst, const Vector<number2>& src, const iVector& gl) const
-{
- Assert(n() == m(), ExcNotQuadratic());
- Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
- Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
- Assert(gl.n() == n(), ExcDimensionMismatch(gl.n(), n()));
-
- number2 s;
- if ((n()==3) && (m()==3))
- {
- number2 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))
- {
- number2 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))
- {
- number2 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
- {
- number* e = val;
- const unsigned int size_m = m(),
- size_n = n();
- for (unsigned int i=0; i<size_m; ++i)
- {
- s = 0.;
- for (unsigned int j=0; j<size_n; ++j)
- if(gl(i)<gl(j)) s += src(j) * *(e++);
- dst(i) += s;
- }
- }
-}
-
-
-
-template <typename number>
-template <typename number2>
-void FullMatrix<number>::Tvmult (Vector<number2>& dst,
- const Vector<number2>& src,
- const bool adding) const
+void FullMatrix<number>::Tvmult (Vector<number2> &dst,
+ const Vector<number2> &src,
+ const bool adding) const
{
Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m()));
FullMatrix<number>::gauss_jordan()
{
Assert (dim_range == dim_image, ExcNotQuadratic());
- iVector p(n());
+ vector<unsigned int> p(n());
unsigned int i,j,k,r;
number max, hr;
- for (i=0; i<n(); ++i) p(i) = i;
+ for (i=0; i<n(); ++i)
+ p[i] = i;
for (j=0; j<n(); ++j)
{
{
hr = el(j,k) ; el(j,k) = el(r,k) ; el(r,k) = hr;
}
- i = p(j) ; p(j) = p(r) ; p(r) = i;
+ i = p[j];
+ p[j] = p[r];
+ p[r] = i;
}
// transformation
Vector<number> hv(n());
for (i=0; i<n(); ++i)
{
- for (k=0; k<n(); ++k) hv(p(k)) = el(i,k);
+ for (k=0; k<n(); ++k) hv(p[k]) = el(i,k);
for (k=0; k<n(); ++k) el(i,k) = hv(k);
}
}
+++ /dev/null
-/*---------------------------- ivector.h ---------------------------*/
-/* $Id$ */
-#ifndef __ivector_H
-#define __ivector_H
-/*---------------------------- ivector.h ---------------------------*/
-
-
-// This file is part of the DEAL Library
-// DEAL is Copyright(1995) by
-// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
-
-#include <base/exceptions.h>
-
-
-/**
- * Integer Vector.
- * Memory for Components is supplied explicitly <p>
- * ( ! Amount of memory needs not to comply with actual dimension due to reinitializations ! ) <p>
- * - all defined methods for iVectors are supplied <p>
- * - operators available are `=` and `( )` <p>
- * CONVENTIONS for used `equations` : <p>
- * - THIS vector is always named `U` <p>
- * - vectors are always uppercase , scalars are lowercase
- */
-class iVector
-{
- friend class dFMatrix;
-
- protected:
-
- /// Dimension. Actual number of components
- unsigned int dim;
-
- /// Dimension. Determines amount of reserved memory , evtl. >DIM !
- unsigned int maxdim;
-
- /// Component-array.
- int *val;
-
- public:
-
- /**@name 1: Basic Object-handling */
- //@{
- /**
- * Dummy-Constructor. Dimension=1
- */
- iVector ();
-
- /**
- * Copy-Constructor. Dimension set to that of V , <p>
- * all components are copied from V
- */
- iVector (const iVector& V);
-
- /**
- * Constructor. Dimension = N (>0)
- */
- iVector (const unsigned int N);
-
- /**
- * Destructor. Clears memory
- */
- ~iVector ();
-
- /**
- * U(0-N) = s . Fill all components
- */
- iVector& operator = (const int s);
-
- /**
- * U = V . Copy all components
- */
- iVector& operator = (const iVector& V);
-
- /**
- * Change Dimension. <p>
- * Set dimension to N <p>
- * ! reserved memory for This remains unchanged ! <p>
- * on fast==false vector is filled by 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==false vector is filled by 0.
- */
- void reinit (const iVector& V, const bool fast=false);
-
- /**
- * Set all elements to zero.
- */
- void clear ();
-
- /**
- * Inquire Dimension. returns Dimension ,
- * INLINE
- */
- unsigned int n () const;
- //@}
-
-
- /**@name 2: Data-Access
- */
- //@{
- /**
- * Access Components. returns U(i) ,
- * INLINE
- */
- int operator() (const unsigned int i) const;
-
- /**
- * Access Components. returns U(i) ,
- * INLINE
- */
- int& operator() (const unsigned int i);
- //@}
-
-
- /**@name 3: Modification of vectors
- */
- //@{
- /**
- * U+=V . Simple addition
- */
- void add (const iVector& V);
-
- /**
- * U+=a*V . Simple addition
- */
- void add (const unsigned int a, const iVector& V);
-
- /**
- * U=a*V . Replacing
- */
- void equ (const unsigned int a, const iVector& V);
- //@}
-
- /**
- * Exception
- */
- DeclException2 (ExcInvalidIndex,
- int, int,
- << "The given index " << arg1
- << " should be less than " << arg2 << ".");
- /**
- * Exception
- */
- DeclException1 (ExcInvalidNumber,
- int,
- << "The provided number is invalid here: " << arg1);
- /**
- * Exception
- */
- DeclException0 (ExcOutOfMemory);
-};
-
-
-/*-------------------------Inline functions -------------------------------*/
-
-
-inline unsigned int iVector::n() const
-{
- return dim;
-}
-
-
-
-inline int iVector::operator() (unsigned int i) const
-{
- Assert (i<dim, ExcInvalidIndex(i,dim));
- return val[i];
-}
-
-
-
-inline int& iVector::operator() (unsigned int i)
-{
- Assert (i<dim, ExcInvalidIndex(i,dim));
- return val[i];
-}
-
-
-
-
-/*---------------------------- ivector.h ---------------------------*/
-/* end of #ifndef __ivector_H */
-#endif
-/*---------------------------- ivector.h ---------------------------*/
//forward declarations
template <typename number> class Vector;
template <typename number> class SparseMatrix;
-
-class iVector;
class ostream;
//////////
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& rows, const iVector& cols);
/**
* Print the sparsity of the matrix
template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
-template void FullMatrix<TYPEMAT>::gsmult(Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const iVector&) const;
#undef TYPEVEC
#define TYPEVEC float
template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
-template void FullMatrix<TYPEMAT>::gsmult(Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const iVector&) const;
#undef TYPERES
#define TYPERES float
template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
-template void FullMatrix<TYPEMAT>::gsmult(Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const iVector&) const;
#undef TYPEVEC
#define TYPEVEC float
template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
-template void FullMatrix<TYPEMAT>::gsmult(Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const iVector&) const;
#undef TYPERES
+++ /dev/null
-// $Id$
-
-// This file is part of the DEAL Library
-// DEAL is Copyright(1995) by
-// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
-
-#include <lac/ivector.h>
-
-
-
-iVector::iVector () :
- dim(1),
- maxdim(1),
- val(new int[1])
-{
- Assert (val != 0, ExcOutOfMemory());
- clear ();
-}
-
-
-iVector::iVector (const unsigned int n) :
- dim(n),
- maxdim(n),
- val(new int[n])
-{
- Assert (n>0, ExcInvalidNumber(n));
-
- Assert (val != 0, ExcOutOfMemory());
- clear ();
-}
-
-
-iVector::iVector (const iVector& v) :
- dim(v.dim),
- maxdim(v.maxdim),
- val(0)
-{
- reinit(v.dim,1);
- for (unsigned int i=0; i<dim; i++)
- val[i] = v.val[i];
-}
-
-
-
-void iVector::reinit (const unsigned int n, const bool fast)
-{
- Assert (n>0, ExcInvalidNumber(n));
- if (n>maxdim)
- {
- if (val != 0) delete[] val;
- val = new int[n];
- Assert (val != 0, ExcOutOfMemory());
- maxdim = n;
- }
- dim = n;
- if (!fast)
- clear ();
-}
-
-
-
-void iVector::reinit (const iVector& v, const bool fast)
-{
- const unsigned int n = v.n();
- if (n>maxdim)
- {
- delete[] val;
- val = new int[n];
- Assert (val != 0, ExcOutOfMemory());
- maxdim = n;
- }
- dim = n;
- if (!fast)
- clear ();
-}
-
-
-
-iVector::~iVector ()
-{
- if (val)
- delete[] val;
-}
-
-
-
-void iVector::clear () {
- for (unsigned int i=0; i<dim; ++i)
- val[i] = 0;
-}
-
-
-void iVector::add (const iVector& v)
-{
- for (unsigned int i=0; i<dim; ++i) val[i] += v(i);
-}
-
-
-
-void iVector::add (const unsigned int a, const iVector& v)
-{
- for (unsigned int i=0; i<dim; ++i) val[i] += a * v(i);
-}
-
-
-
-void iVector::equ (const unsigned int a, const iVector& u)
-{
- for (unsigned int i=0; i<dim; ++i) val[i] = a*u.val[i];
-}
-
-
-
-iVector& iVector::operator = (const int s)
-{
-// if (s==0.) memset(val,sizeof(*val) * dim, 0);
- for (unsigned int i=0; i<dim; ++i) val[i] = s;
- return *this;
-}
-
-
-iVector& iVector::operator = (const iVector& v)
-{
- if (v.dim != dim) reinit(v,1);
-
- for (unsigned int i=0; i<dim; ++i) val[i] = v.val[i];
- return *this;
-}
#include <lac/sparsematrix.h>
-#include <lac/ivector.h>
#include <iostream>
#include <iomanip>
-void
-SparseMatrixStruct::add_matrix (const iVector& rowcols)
-{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
- unsigned int i,j;
- for (i=0;i<rowcols.n();i++)
- for (j=0;j<rowcols.n();j++)
- add(rowcols(i), rowcols(j));
-}
-
-
-
-void
-SparseMatrixStruct::add_matrix (const iVector& rows, const iVector& cols)
-{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
- unsigned int i,j;
- for (i=0;i<rows.n();i++)
- for (j=0;j<cols.n();j++)
- add(rows(i), cols(j));
-}
-
-
-
void
SparseMatrixStruct::print_gnuplot (ostream &out) const
{