implementation.h

高效的c++科学算法库

/*                                                           -*- c++ -*-
 ***********************************************************************
 * Scientific Library (GNU Public Licence)
 *
 * module sparse
 *
 * Author: Pierre Aubert [e-mail]
 *
 * $Id: implementation.h,v 1.2 1998/11/16 17:22:13 paubert Exp $
 *
 * Suggestions: [e-mail]
 * Bugs:   [e-mail] 
 *
 * For more information, please see the sl++ Home Page:
 * http://wwwinfo.cern.ch/~ldeniau/sl.html 
 *
 ***********************************************************************
 */

#ifndef SL_MATRIX_STRUCTURES_SSS_IMPLEMENTATION_H
#define SL_MATRIX_STRUCTURES_SSS_IMPLEMENTATION_H

#ifdef HAVE_NAMESPACE
namespace sl 
{
#endif

      // ----------------------------------------------------------------------
      //                                                            CONSTRUCTOR
      // ----------------------------------------------------------------------
    SymmetricSkyline::SymmetricSkyline  ( )
        : my_n(0),
          my_data_length(0),
          my_profil_is_set(false),
          my_matrix_is_factored(false)
    {}
    
  
      // ----------------------------------------------------------------------
      //                                                            CONSTRUCTOR
      // ----------------------------------------------------------------------
    SymmetricSkyline::SymmetricSkyline  ( Index const& n )
        : my_n(n),                              
          my_data_length(n),                    
          my_profil_is_set(false),              
          my_matrix_is_factored(false)
    {
          // creation of vectors
        my_haut = VectorUInt(n);
        my_diag = VectorUInt(n);
          // init data
        my_haut = 1;
        for( Index i=1 ; i<=n ; ++i )
            my_diag(i) = i;
    }
    
      // ----------------------------------------------------------------------
      //                                                            CONSTRUCTOR
      // ----------------------------------------------------------------------
    SymmetricSkyline::SymmetricSkyline  ( Index const& n, Index const& p )
        : my_n(n),
          my_data_length(n),
          my_profil_is_set(false),
          my_matrix_is_factored(false) 
    {
        assert( n == p );

          // creation of vectors
        my_haut = VectorUInt(n);
        my_diag = VectorUInt(n);
      
          // init data
        my_haut = 1;
        for( Index i=1 ; i<=n ; ++i )
            my_diag(i) = i;
    }
  
      // ----------------------------------------------------------------------
      //                                                            CONSTRUCTOR
      // ----------------------------------------------------------------------
    SymmetricSkyline::SymmetricSkyline  ( Index const& n,
                                          Index const& p, 
                                          Index const &nb_store )
        : my_n(n),
          my_data_length(nb_store),
          my_profil_is_set(false),
          my_matrix_is_factored(false) 
    {
        assert( n == p );
        assert( nb_store >= n);
            
          // creation of vectors
        my_haut = VectorUInt(n);
        my_diag = VectorUInt(n);
            
          // init data
        my_haut = 1;
        for( Index i=1 ; i<=my_diag.rows() ; ++i )
            my_diag(i) = i;
    }
    
      // ----------------------------------------------------------------------
      //                                                             DESTRUCTOR
      // ----------------------------------------------------------------------
    SymmetricSkyline::~SymmetricSkyline()
    {
          // ~my_haut;
          // ~my_diag;
    }

      // ----------------------------------------------------------------------
      //                                                               resize()
      // ----------------------------------------------------------------------
    void SymmetricSkyline::resize( Index const new_size )
    {
        if( new_size != my_n || new_size != my_n )
            throw eNotYetImplemented;
    }

      // ----------------------------------------------------------------------
      //                                                               resize()
      // ----------------------------------------------------------------------
    void SymmetricSkyline::resize( Index const n, Index const p )
    {
        assert( n == p );
        resize(n);
    }

      // ----------------------------------------------------------------------
      //                                                              display()
      // ----------------------------------------------------------------------
    ostream& SymmetricSkyline::display(ostream& os, bool const verbose) const
    {
        const Index block_size = 24;
      
        Index i,j;
      
        Index mm = rows();
        Index mn = cols();
        Index mp = size();
        
        os << "Sparse matrix: "
           << mm << "x" << mn
           << " with "  << mp << " stored elements,"
           << " density = "   << mp*100/(mm*mn) << "%,";
        if( isFactored() )
            os << " matrix is factored," << endl;
        
        if( verbose )
            {
                os << "  my_diag = " << diag() << endl;
                os << "  my_haut = " << haut() << endl;
            }

        if( mm >= block_size )
            os << " reductions factors are " << mm/block_size << " and " << mn/block_size;

        os << '.' << endl;

        if( mm < block_size )
            {
                for (i = 1; i <= mm ; ++i)
                    {
                        os << "[ ";
                        for (j = 1; j <= mn ; ++j){
                              // are we in the stored part ?
                            if( j >= i ){
                                  // no element
                                if( j-i >= haut(j) )
                                    os << " ";
                                else
                                    os << "*";
                            } else {
                                  // use symmetrie
                                  // no element
                                if( i-j >= haut(i) )
                                    os << " ";
                                else
                                    os << "*";
                            }
                        } // end for j
                        os << " ]" << endl;
                    } // end for i
            }
        else
            {
                Index k1,l1,k2,l2;
                Index r,s;

                r = mm/block_size;
                s = mn/block_size;

                  // remark: the last rows and cols are not ploted
                for (k1 = 1; k1 < mm/r ; ++k1){
                      // start a block of r-lines
                    os << "[ ";
                    for (l1 = 1; l1 < mn/s ; ++l1){
                        bool nonnul = false;
                        for (k2 = 1; k2 <= r ; ++k2){
                              // compute indice (line)
                            i = k1*r+k2;
                                
                            for (l2 = 1; l2 <= s ; ++l2){
                                  // compute indice (column)
                                j = l1*s+l2;
                                            
                                  // are we in the stored part ?
                                if( j >= i ){
                                      // no element
                                    if( j-i < haut(j) )
                                        nonnul = true;
                                } else {
                                      // no element
                                    if( i-j < haut(i) )
                                        nonnul = true;
                                }
                            } // end l2 
                        } // end k2
                          // if theyre is something in the small r x s rectangle
                        if( nonnul )
                            os << "*";
                        else
                            os << " ";
                          // end a block of r-line x s-column
                    } // end l1 
                    os << " ]" << endl;
                } // end k1
            } // end if mm large
        
        
        return os;
    }

      // ----------------------------------------------------------------------
      //                                                                 diag()
      // ----------------------------------------------------------------------
    Index SymmetricSkyline::diag(Index const& i) const
    {
        return my_diag(i);
    }

      // ----------------------------------------------------------------------
      //                                                           isFactored()
      // ----------------------------------------------------------------------
    Bool SymmetricSkyline::isFactored() const
    {
        return my_matrix_is_factored;
    }

      // ----------------------------------------------------------------------
      //                                                     setNotIsFactored()
      // ----------------------------------------------------------------------
    void SymmetricSkyline::setIsNotFactored()
    {
        my_matrix_is_factored = false;
    }

      // ----------------------------------------------------------------------
      //                                                        setIsFactored()
      // ----------------------------------------------------------------------
    void SymmetricSkyline::setIsFactored() 
    {
        my_matrix_is_factored = true;
    }

      // ----------------------------------------------------------------------
      //                                                                 haut()
      // ----------------------------------------------------------------------
    Index SymmetricSkyline::haut(Index const& i) const
    {
        return my_haut(i);
    }
  
      // ----------------------------------------------------------------------
      //                                                                 rows()
      // ----------------------------------------------------------------------
    Index SymmetricSkyline::rows() const
    {
        assert( my_n == my_diag.rows() );
        assert( my_n == my_haut.rows() );
        return my_n;
    }

      // ----------------------------------------------------------------------
      //                                                                 cols()
      // ----------------------------------------------------------------------
    Index SymmetricSkyline::cols() const
    {
        assert( my_n == my_diag.rows() );
        assert( my_n == my_haut.rows() );
        return my_n;
    }

      // ----------------------------------------------------------------------
      //                                                                 size()
      // ----------------------------------------------------------------------
    size_t SymmetricSkyline::size() const
    {
        size_t sz = my_diag(my_n)+my_haut(my_n)-1;

        assert( sz <= my_data_length );
      
          // return the max 
        return sz;
    }

      // ----------------------------------------------------------------------
      //                                                             map_fast()
      // ----------------------------------------------------------------------