symm_bdmatrix.h

Go to the documentation of this file.

// Version 2.0
//
//--------------------------------------------------------------------
// NOTE: The whole code in this files is based on the assumption that
// the main diagonal of an object of type Symm_BdMatrix has been
// allocated. This should be taken into account upon doing changes.
// 13th of march 1999
//--------------------------------------------------------------------
// Jose L. Tinoco,
// Feb. 1999
// $Id: symm_bdmatrix.h,v 1.15.2.10 2008/07/09 16:28:07 garloff Exp $
//--------------------------------------------------------------------

#ifndef _SYMM_BAND_MATRIX
#define _SYMM_BAND_MATRIX

#include "vector.h"
#include "matrix_sig.h"
//#include "solver/ilu0precond.h"

template <typename T> class Symm_BdMatrix;
// Avoid -fguiding-decls
#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "symm_bdmatrix_gd.h"
#endif

NAMESPACE_TBCI

// Forward - declarations
template <typename T> class Symm_BdMatrix;
//template <typename T, typename MatrixType = Symm_BdMatrix<T> > class ILU0_Symm_BdMatrixPreconditioner;
template <typename T> class ILU0_Symm_BdMatrixPreconditioner;


// Class declarations
//--------------------------------------------------------------------
template<typename T>
class Symm_BdMatrix : public Matrix_Sig<T>
{
        friend class ILU0_Symm_BdMatrixPreconditioner<T>;

        protected:
                T zero;
                T* elementPtr; // Allocation of all diagonals as a single vector
                unsigned int dimension;
                BVector<T*> rowPtr; // Elements point to entries in main diag
                BVector<unsigned int> rowOccupation; // which diags have been allocated
                BVector<unsigned int> conf; // ditto 

                void construct (const unsigned int N, const T&, T*&,
                                const BVector<unsigned int>&, BVector<T*>&, 
                                BVector<unsigned int>&);
                inline void create (const T& val, const unsigned int N);
                inline void create (const T& val, const unsigned int N, 
                                    const BVector<unsigned int>&);
                inline void destroy ();
        public:
                // Constructor, Destructor
                Symm_BdMatrix ()                                  { create(T(0),1); }
                Symm_BdMatrix (const unsigned int N)              { create(T(0),N); }
                Symm_BdMatrix (const T& val, const unsigned int N) { create(val,N);  }
                Symm_BdMatrix (const T& val, const unsigned int N, 
                               const BVector<unsigned int>& diagConf)  { create(val,N,diagConf); }
                //Symm_BdMatrix(const Matrix_Sig<T>&);

                inline void resize (const T&, const unsigned int);
                inline void resize (const unsigned int newDim)  { resize (0, newDim); }
                inline void resize (const T&, const unsigned int N, const BVector<unsigned int>&);
                
                ~Symm_BdMatrix ()  { destroy(); }

                // allow instantiation (Matrix_Sig)
                /*virtual*/ const char* mat_info () const  { return ("Symm_BdMatrix"); }

                // element access
                const T& operator() (const unsigned int i, const unsigned int k) const;
                const T& get (const unsigned int i, const unsigned int k) const
                { return (*this) (i, k); }
                friend STD__ ostream& operator<< FGD (STD__ ostream&, const Symm_BdMatrix<T>&);
                T& setval (const unsigned int i, const unsigned int k);
                void setval (const T& wert, const unsigned int i, const unsigned int k)
                { setval (i,k) = wert; }
                void autoinsert (const T& wert, const unsigned int i, const unsigned int k)
                { setval (i,k) = wert; }
                
                // Matrix-Vector-Multiplication
                TVector<T> operator * (const Vector<T>&) const;
                
                // Misc
                unsigned int rows () const  { return dimension; }
                unsigned int columns () const { return dimension; }
                void clear ();
                // typedefs
                typedef T value_type;
                typedef T element_type;
                typedef T aligned_value_type TALIGN(MIN_ALIGN2); 

};


/********************************************************************
 *******************     Memberdefinitions    ***********************
 ********************************************************************/


template<typename T>
inline void Symm_BdMatrix<T>::resize (const T& val, const unsigned int N)
{
        destroy ();
        create (val, N);
}

template<typename T>
inline void Symm_BdMatrix<T>::resize (const T& val, const unsigned int N,
                                      const BVector<unsigned int>& Diagconf)
{
        destroy ();
        create (val, N, Diagconf);
}



//--------------------------------------------------------------------
// Konstruktor: Allokiert nur die Hauptdiagonale, initialisiert alle
// Eintraege auf T(0) und setzt alle privaten Variablen
template<typename T>
inline void Symm_BdMatrix<T>::create (const T& wert, const unsigned int dim)
{
        dimension = dim;
        zero = T(0);
        conf.resize (0);
        rowPtr.resize (dim);
        rowOccupation.resize (dim);

        // check
        if(dim==0) 
                STD__ cerr << "void Symm_BdMatrix<T>::Symm_BdMatrix(const T&, const unsigned int)"  << STD__ endl
                << "Invalid constructor parameters: Matrix can not have dimension == 0" << STD__ endl;
        
        // erzeuge Speicherbereich
        construct (dim, wert, elementPtr, conf, rowPtr, rowOccupation);

}// End OF Symm_BdMatrix<T>::Symm_BdMatrix (const T&, const unsigned int)



//--------------------------------------------------------------------
// Constructor: Allokiert die Hauptdiagonale und die Nebendiagonalen
// mit den Nummern configVector(0), configVector(1), configVector(2),
// u.s.w., initiallisiert alle Eintraege auf T(0) und setzt alle privaten
// Variablen
template<typename T>
inline void Symm_BdMatrix<T>::create(const T& wert, const unsigned int dim, 
                                     const BVector<unsigned int>& configVector)
{
        dimension=dim;
        zero=T(0);
        conf.resize(configVector.size());
        rowPtr.resize(dim);
        rowOccupation.resize(dim);

        // check
        if(dim==0) 
                STD__ cerr << "void Symm_BdMatrix<T>::Symm_BdMatrix(const T&, const unsigned int, const BVector<unsigned int>&)"  << STD__ endl
                << "Invalid constructor parameters: Matrix can not have dimension == 0" << STD__ endl;
        if(configVector.size()>=dim) 
                STD__ cerr << "void Symm_BdMatrix<T>::Symm_BdMatrix(const T&, const unsigned int, const BVector<unsigned int>&)"  << STD__ endl
                << "Invalid constructor parameters: Number of off-diagonals larger than matrix dimension" << STD__ endl;
        for(unsigned i=0; i<configVector.size(); i++)
                if(configVector(i)>=dim) 
                        STD__ cerr << "void Symm_BdMatrix<T>::Symm_BdMatrix(const T&, const unsigned int, const BVector<unsigned int>&)"  << STD__ endl
                        << "Invalid constructor parameters: Off-diagonal is out of range" << STD__ endl;
        
        // erzeuge Speicherbereich
        construct(dim,wert,elementPtr,configVector,rowPtr,rowOccupation);
        conf = configVector;
}// Ende Symm_BdMatrix<T>::Symm_BdMatrix(const T&, const unsigned int, const unsigned int)



//--------------------------------------------------------------------
// Construct: Hilfefunktion. Erzeugt das interne "Datensystem", d.h.
// allokiert den noetigen Speicherbereich und besetzt die Vectoren
// rowPtr und rowOccupation
template<typename T>
inline void Symm_BdMatrix<T>::construct (const unsigned int dim, const T& value,
                                         T*& elemPtr, const BVector<unsigned int>& ConfVector,
                                         BVector<T*>& rPtr, BVector<unsigned int>& rowOccup)
{
        // berechne Groesse des Speicherbereiches
        unsigned int i, j, length = dim; // Laenge der Hauptdiagonale
        for (i=0; i<ConfVector.size(); i++) 
                length += dim - ConfVector(i); // addiere Laenge der Nebendiagonalen

        // allokiere Speicherbereich initialisiere alle Elemente auf value
        elemPtr = new T [length];
        for (i=0; i<length; i++) 
                elemPtr[i] = value;
        
        // besetze zPtr mit Zeigern auf die Zeilen (Eintraege der Hauptdiagonale)
        if (ConfVector.size()) {
                unsigned int next = 0;
                for (i=0; i<dim; i++) {
                        rPtr(i) = &(elemPtr[next]);
                        next += ConfVector.size()+1;
                        // verringere next gegen Ende des Speicherbereiches -> unten
                        // rechts in der Matrix
                        for (j = ConfVector.size(); j >= 1; j--) 
                                if(ConfVector(j-1) +i >=dim) next--;
                                else break;
                }
        }
        else 
                for (i=0; i<dim; i++) rPtr(i) = &(elemPtr[i]);

        // besetze rowOccup mit Null wenn die entspechende Diagonale nicht allokiert ist
        // sonst mit der Nummer der Stelle, an der angesprochene Eintrag im Zeilenblock im elemPtr
        // gespeichert ist
        for (i=0; i<dim; i++) 
                rowOccup(i) = 0;
        /*if(ConfVector.size())*/ 
        for (i=0; i<ConfVector.size(); i++)
                rowOccup(ConfVector(i)) = i+1;

}// Ende void Symm_BdMatrix<T>::construct(unsigned int,const T&,T*&,const BVector<unsigned int>&,Vector<T*>&,Vector<unsigned int>&)


#if 0
template <typename T>
Symm_BdMatrix<T>::Symm_BdMatrix(const Matrix_Sig<T>& m)
{
        create(T(0),1);
        for (unsigned i = 0; i < m.rows(); i++)
                for (unsigned j = 0; j < m.columns(); j++)
                        if (MATH__ fabs(m(i,j)) > 1e-15)
                                this->setval(i, j) = m (i,j);
}
#endif

//--------------------------------------------------------------------
// Destruktor: Gibt den belegten Speicherplatz wieder frei
template<typename T>
inline void Symm_BdMatrix<T>::destroy()
{
        delete[] elementPtr;
}// Ende Symm_BdMatrix<T>::~Symm_BdMatrix()



//-------------------------------------------------------------
// Setval: Wie operator() mit dem Unterschied, dass setval die zu
// (i,k) zugehoerige Nebendiagonale allokiert und alle Eintraege
// in ihr zu T(0) initiallisiert wenn diese nicht vorhanden ist,
// und gibt erst dann den Wert (i,k) zurueck, auf dem geschrieben
// werden kann 
template<typename T>
inline T& Symm_BdMatrix<T>::setval (const unsigned int i, const unsigned int k)
{
        // check
        if((i>=dimension) || (k>=dimension)) 
                STD__ cerr << "T& Symm_BdMatrix<T>::setval (const unsigned int, const unsigned int)" << STD__ endl 
                << "Invalid Index: Index out of range" << STD__ endl;
        
        // tausche k und i um, wenn i>k, also wenn das Element unterhalb der Hauptdiagonalen ist
        unsigned int l,m;
        if (i > k) { 
                l=k; m=i; 
        } else {
                l=i; m=k;
        }
        // (l,m) oberhalb der Hauptdiagonale

        // falls (l,m) in der Hauptdiagonale
        if (m == l)
                return *rowPtr(l);
        
        // falls (l,m) sonst irgendwo
        if(rowOccupation(m-l) != 0) {
                return *(rowPtr(l) + rowOccupation(m-l));
        }
        else {
                
                // erzeuge neuen Speicherbereich und schreibe neue Diagonale
                // zwischen den Diagonalen Nr. (found-1) und found
                T* newmem;
                BVector<unsigned int> newConfVector(conf.size() + 1);
                
                // schreibe newConfVector und finde Stelle found, in der
                // die neue Nebendiagonale geschrieben wird
                unsigned int j,g,h;
    int found(-1);
                for(g=j=0; j<conf.size(); j++,g++) {
                        // falls konf(f)-te Nebendiagonale weiter ausserhalb als die (m-l)-te
                        if( conf(j) > (m-l)) {newConfVector(g++)= m-l; found=j;break;}
                        newConfVector(g) = conf(j);
                }
                // schreibe rest ab
                for(; j<conf.size(); j++,g++) {
                        newConfVector(g) = conf(j);
                }
                // falls Stelle nicht gefunden -> neue Diagonale ist letzte
                if(found == -1) {newConfVector(conf.size()) = m-l; found=conf.size();}
 
                Vector<T*> newrowPtr(dimension);
                Vector<unsigned int> neuerowOccupation(dimension);
                
                // erzeuge neuen Speicherbereich
                construct(dimension,T(0),newmem,newConfVector,newrowPtr,neuerowOccupation);
                
                T* newPtr; 
                T* oldPtr;
                oldPtr = elementPtr;
                newPtr = newmem;
                
                // schreibe Werte auf neuen Speicherbereich ab
                for(h=1; h<dimension; h++){ // h springt von Zeile zu Zeile
                        j=0; // j laeuft innerhalb einer Zeile
                        while(newPtr != newrowPtr(h)) {
                                *(newPtr++) = *(oldPtr++);
                                if( ((int)j == found)&& (int(found+1) < int(newrowPtr(h) - newrowPtr(h-1))) )
                                        newPtr++;
                                j++;
                        }
                }
                *(newPtr) = *(oldPtr); // letzter Eintrag
                
                // zerstoere alten Speicherbereich
                delete[] elementPtr;
                
                // setze private Variablen auf neue Werte
                rowPtr = newrowPtr;
                rowOccupation = neuerowOccupation;              
                conf.resize(newConfVector.size())  = newConfVector;
                elementPtr = newmem;
                
                // endlich!!
                return *(rowPtr(l) + rowOccupation(m-l));
        }

        
}// Ende T& Symm_BdMatrix<T>::setval (const unsigned int, const unsigned int)



//-------------------------------------------------------------
// Zugriffoperator: Erlaubt Lesezugriff auf das Matrixelement
// mit den Indizes (i,k), wobei auf das erste Element oben links
// in der Matrix mit (0,0) zugegriffen wird
template<typename T>
inline const T& Symm_BdMatrix<T>::operator() (const unsigned int i, const unsigned int k) const
{
        // check
        if((i>=dimension) || (k>=dimension)) 
                STD__ cerr << "const T& Symm_BdMatrix<T>::operator() (const unsigned int, const unsigned int)const" << STD__ endl 
                << "Invalid Index: Index out of range" << STD__ endl;
        
        // tausche k und i um, wenn i>k, also wenn das Element unterhalb der Hauptdiagonalen ist
        unsigned int l,m;
        if(i>k) {l=k; m=i;}
        else {l=i; m=k;}
        // (l,m) oberhalb der Hauptdiagonale

        // falls (l,m) in der Hauptdiagonale
        if(m==l) return *rowPtr(l);
        
        // falls (l,m) sonst irgendwo
        if(rowOccupation(m-l) != 0) return  *(rowPtr(l) + rowOccupation(m-l));
        else return zero;

}// Ende const T& Symm_BdMatrix<T>::operator() (const unsigned int, const unsigned int)const



//-------------------------------------------------------------
// Ausgabenoperator: Schreibt alle Eintraege der Matrix matrix
// auf den ostream out in der natuerlichen Reihenfolge
template<typename T>
inline STD__ ostream& operator<< (STD__ ostream& out, const Symm_BdMatrix<T>& matrix)
{
        unsigned int i,k;
        for(i=0; i< matrix.dimension; i++){
                // schreibe Zeile
                for(k=0; k< matrix.dimension; k++) out << matrix(i,k) << " ";
                out << STD__ endl;
        }
        return out;
}// Ende std::ostream& operator<< (std::ostream&, const Symm_BdMatrix<T>&)



//-------------------------------------------------------------
// Multiplikationsoperator: Multipliziert *this mit Vector.
template<typename T>
inline TVector<T> Symm_BdMatrix<T>::operator* (const Vector<T>& Vector) const
{
        // check
        if (Vector.size() != dimension)
                STD__ cerr << "TVector<T> Symm_BdMatrix<T>::operator* (const Vector<T>&)const" << STD__ endl
                << "Matrix and Vector have different dimensions" << STD__ endl;

        TVector<T> result (Vector.size());

        // falls keine Nebendiagonalen allokiert
        unsigned int i;
        if (conf.size() == 0) {
                for (i=0; i<dimension; i++) 
                        result.setval(i) = ((*(rowPtr(i))) * Vector(i));
                return result;
        }

/*
        // multipliziere in der natuerlichen Reihenfolge, optimiert, 2 Schleifen
        unsigned int j, diagonal, confsize;
        confsize=conf.size();
        // falls mindestens eine Nebendiagonal allokiert
        T temp;
        unsigned max_off=conf(confsize-1);
        T* memPtr = elementPtr;
        unsigned int* konfPtr = conf.vecptr();

        // h haelt die Zeilennummer fest
        // Hauptdiagonale und obere Diagonanlen
        for(unsigned int h=0; h<dimension-1; h++)
        {
                // Hauptdiagonale nur einmal zaehlen
                temp = (*(memPtr)) * Vector(h);
                memPtr++;
                while(memPtr != rowPtr(h+1)) 
                {
                        temp += *(memPtr) * (Vector(h + *(konfPtr)));
                        konfPtr++;
                        memPtr++;
                }
                result.setval(h)=temp;
                // setze konfPtr zurueck
                konfPtr= conf.vecptr();
        }
        // letzter Summand
        result.setval(dimension-1) = (*(memPtr)) * Vector(h);
        
        //Unteren Nebendiagonalen
        for(i=0; i<max_off; i++)
        {
                temp=0;
                for(j=0; ((diagonal=conf(j))<=i) && (j<confsize); j++)
                        temp += *(rowPtr(i-diagonal) + j +1) * Vector(i-diagonal);
                result.setval(i)+=temp;
        }
        for(i=max_off; i<dimension; i++)
        {
                temp=0;
                for(j=0; j<confsize; j++)
                        temp += *(rowPtr(i-conf(j)) + j +1) * Vector(i-conf(j));
                result.setval(i)+=temp;
        }
*/


// /*
        // multipliziere in der natuerlichen Reihenfolge, optimiert (Nr. 1)
        unsigned int j, diagonal, confsize;
        confsize=conf.size();
        // falls mindestens eine Nebendiagonale allokiert
        T temp;
        unsigned max_off=conf(confsize-1);
        for(i=0; i<max_off; i++)
        {
                // Hauptdiagonale
                temp = ((*(rowPtr(i))) * Vector(i));
                // unterhalb der Hauptdiagonalen
                for(j=0; ((diagonal=conf(j))<=i) && (j<confsize); j++)
                        temp += *(rowPtr(i-diagonal) + j +1) * Vector(i-diagonal);
                // oberhalb der Hauptdiagonalen
                for(j=0; j<confsize; j++)
                        temp += *(rowPtr(i) + j +1) * Vector(i + conf(j));
                result.setval(i)=temp;
        }
        for(i=max_off; i<dimension-max_off; i++)
        {
                // Hauptdiagonale
                temp = ((*(rowPtr(i))) * Vector(i));
                for(j=0; j<confsize; j++)
                {
                        // unterhalb der Hauptdiagonalen
                        temp += *(rowPtr(i-conf(j)) + j +1) * Vector(i-conf(j));
                // oberhalb der Hauptdiagonalen
                        temp += *(rowPtr(i) + j +1) * Vector(i + conf(j));
                }               
                result.setval(i)=temp;
        }
        for(i=dimension-max_off; i<dimension; i++)
        {
                // Hauptdiagonale
                temp = ((*(rowPtr(i))) * Vector(i));
                // unterhalb der Hauptdiagonalen
                for(j=0; j<confsize; j++)
                        temp += *(rowPtr(i-conf(j)) + j +1) * Vector(i-conf(j));
                // oberhalb der Hauptdiagonalen
                for(j=0; (conf(j)<dimension-i) && (j<confsize); j++)
                        temp += *(rowPtr(i) + j +1) * Vector(i + conf(j));
                result.setval(i)=temp;
        }
// */

/*
        // multipliziere in der natuerliche Reihenfolge
        unsigned int j, diagonal, confsize;
        confsize=conf.size();
        // falls mindestens eine Nebendiagonale allokiert
        T temp;

        for(i=0; i<dimension; i++)
        {
                // Hauptdiagonale
                temp = ((*(rowPtr(i))) * Vector(i));
                
                // unterhalb der Hauptdiagonalen
                for(j=0; ((diagonal=conf(j))<=i) && (j<confsize); j++)
                        temp += *(rowPtr(i-diagonal) + j +1) * Vector(i-diagonal);
                
                // oberhalb der Hauptdiagonalen
                for(j=0; (conf(j)<dimension-i) && (j<confsize); j++)
                        temp += *(rowPtr(i) + j +1) * Vector(i + conf(j));
                
                result.setval(i)=temp;
        }
*/

/*
        // multipliziere in der Reihenfolge der Allokierung der Matrixelemente
        T* memPtr = elementPtr;
        unsigned int* konfPtr = conf.vecptr();
        result.clear();
        T temp;

        // h haelt die Zeilennummer fest
        for(unsigned int h=0; h<dimension-1; h++)
        {
                // Hauptdiagonale nur einmal zaehlen
                temp = (*(memPtr)) * Vector(h);
                memPtr++;
                
                // fuehre alle Operationen durch, in denen die Matrixelemente der
                // Zeile h (ohne Hauptdiagonale) enthalten sind
                while(memPtr != rowPtr(h+1)) 
                {
                        // oberhalb der Hauptdiagonale
                        temp += *(memPtr) * (Vector(h + *(konfPtr)));
                        // unterhalb der Hauptdiagonale
                        result.setval(h + *(konfPtr)) += *(memPtr) * (Vector(h));
                        
                        konfPtr++;
                        memPtr++;
                }
                result.setval(h)+=temp;
                
                // setze konfPtr zurueck
                konfPtr= conf.vecptr();
        }
        // letzter Summand
        result.setval(h) += (*(memPtr)) * Vector(h);
*/

/*
        // multipliziere in der Reihenfolge der Allokierung der Matrixelemente, (Nr. 2)
        // 2 Schleifen!!
        T* memPtr = elementPtr;
        unsigned int* konfPtr = conf.vecptr();
        T temp;

        // h haelt die Zeilennummer fest
        // Hauptdiagonale und obere Diagonanlen
        for(unsigned int h=0; h<dimension-1; h++)
        {
                // Hauptdiagonale nur einmal zaehlen
                temp = (*(memPtr)) * Vector(h);
                memPtr++;
                while(memPtr != rowPtr(h+1)) 
                {
                        temp += *(memPtr) * (Vector(h + *(konfPtr)));
                        konfPtr++;
                        memPtr++;
                }
                result.setval(h)=temp;
                // setze konfPtr zurueck
                konfPtr= conf.vecptr();
        }
        // letzter Summand
        result.setval(h) = (*(memPtr)) * Vector(h);
        // unterhalb der Hauptdiagonale
        memPtr = elementPtr;
        for(h=0; h<dimension-1; h++)
        {
                // Hauptdiagonale ueberspringen
                memPtr++;
                // setze konfPtr zurueck
                konfPtr= conf.vecptr();
                while(memPtr != rowPtr(h+1)) 
                {
                        result.setval(h + *(konfPtr)) += *(memPtr) * (Vector(h));
                        konfPtr++;
                        memPtr++;
                }
        }
*/
        
        return result;
}// Ende TVector<T> Symm_BdMatrix<T>::operator* (const Vector<T>&)const



//-------------------------------------------------------------
// Clear: Assigns 0 to all elements
template<typename T>
inline void Symm_BdMatrix<T>::clear()
{
        T* Ptr = elementPtr;
        while(Ptr != rowPtr(dimension-1)) (*Ptr++) = T(0);
        *Ptr = T(0); // last entry

}// end void Symm_BdMatrix<T>::clear()


NAMESPACE_END

#endif // !defined(_Symm_BdMatrix)

---