matrix_triangpacked.cxx

Multivac 的Level set包

// Copyright (C) 2001-2004 Vivien Mallet
//
// This file is part of Seldon library.
// Seldon library provides matrices and vectors structures for
// linear algebra.
// 
// Seldon is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// Seldon is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License (file "license") for more details.
//
// For more information, please see the Seldon home page:
//     http://spacetown.free.fr/lib/seldon/

#ifndef SELDON_FILE_MATRIX_TRIANGPACKED_CXX

#include "Matrix_TriangPacked.hxx"

namespace Seldon
{


  /****************
   * CONSTRUCTORS *
   ****************/
  

  //! Default constructor.
  /*!
    On exit, the matrix is an empty 0x0 matrix.
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_TriangPacked<T, Prop, Storage, Allocator>
  ::Matrix_TriangPacked(): Matrix_Base<T, Allocator>()
  {
  }


  //! Main constructor.
  /*! Builds a i x j triangular matrix in packed form.
    \param i number of rows.
    \param j number of columns.
    \note 'j' is assumed to be equal to 'i' and is therefore discarded.
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_TriangPacked<T, Prop, Storage, Allocator>
  ::Matrix_TriangPacked(int i, int j): Matrix_Base<T, Allocator>(i, i)
  {

#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif

	this->data_ = this->allocator_.allocate((i * (i + 1)) / 2, this);

#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
	this->m_ = 0;
	this->n_ = 0;
	this->data_ = NULL;
	return;
      }
    if (this->data_ == NULL)
      {
	this->m_ = 0;
	this->n_ = 0;
	return;
      }
#endif

  }

  
  /**************
   * DESTRUCTOR *
   **************/


  //! Destructor.
  template <class T, class Prop, class Storage, class Allocator>
  inline Matrix_TriangPacked<T, Prop, Storage, Allocator>
  ::~Matrix_TriangPacked()
  {
    
#ifdef SELDON_CHECK_MEMORY
    try
      {
#endif
	
	if (this->data_ != NULL)
	  {
	    this->allocator_.deallocate(this->data_,
					(this->m_ * (this->m_ + 1)) / 2);
	    this->data_ = NULL;
	  }
	
#ifdef SELDON_CHECK_MEMORY
      }
    catch (...)
      {
	this->m_ = 0;
	this->n_ = 0;
	this->data_ = NULL;
      }
#endif

  }


  //! Clears the matrix.
  /*!
    Destructs the matrix.
    \warning On exit, the matrix is an empty 0x0 matrix.
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_TriangPacked<T, Prop, Storage, Allocator>::Clear()
  {
    this->~Matrix_TriangPacked();
  }


  /*******************
   * BASIC FUNCTIONS *
   *******************/


  //! Returns the number of elements stored in memory.
  /*!
    \return The number of elements stored in memory.
  */
  template <class T, class Prop, class Storage, class Allocator>
  int Matrix_TriangPacked<T, Prop, Storage, Allocator>::GetDataSize() const
  {
    return (this->m_ * (this->m_ + 1)) / 2;
  }


  /*********************
   * MEMORY MANAGEMENT *
   *********************/


  //! Reallocates memory to resize the matrix.
  /*!
    On exit, the matrix is a i x j matrix.
    \param i new number of rows.
    \param j new number of columns.
    \warning Depending on your allocator, data may be lost.
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_TriangPacked<T, Prop, Storage, Allocator>
  ::Reallocate(int i, int j)
  {
    if (i != this->m_)
      {
	this->m_ = i;
	this->n_ = i;

#ifdef SELDON_CHECK_MEMORY
	try
	  {
#endif

	    this->data_ =
	      reinterpret_cast<pointer>(this->allocator_.reallocate(this->data_,
								    (i*(i+1))/2,
								    this) );

#ifdef SELDON_CHECK_MEMORY
	  }
	catch (...)
	  {
	    this->m_ = 0;
	    this->n_ = 0;
	    this->data_ = NULL;
	    throw NoMemory("Matrix_TriangPacked::Reallocate(int, int)",
			   "Unable to reallocate memory for data_.");
	  }
	if (this->data_ == NULL)
	  {
	    this->m_ = 0;
	    this->n_ = 0;
	    throw NoMemory("Matrix_TriangPacked::Reallocate(int, int)",
			   "Unable to reallocate memory for data_.");
	  }
#endif

      }
  }


  //! Changes the size of the matrix and sets its data array
  //! (low level method).
  /*!
    The matrix is first cleared (memory is freed). The matrix is then resized
    to a i x j matrix, and the data array of the matrix is set to 'data'.
    'data' elements are not duplicated: the new data array of the matrix is
    the 'data' array. It is useful to create a matrix from pre-existing data.
    \param i new number of rows.
    \param j new number of columns.
    \param data new array storing elements.
    \warning 'data' has to be used carefully outside the object.
    Unless you use 'Nullify', 'data' will be freed by the destructor,
    which means that 'data' must have been allocated carefully. The matrix
    allocator should be compatible.
    \note This method should only be used by advanced users.
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline void Matrix_TriangPacked<T, Prop, Storage, Allocator>::
  SetData(int i, int j,
	  typename Matrix_TriangPacked<T, Prop, Storage, Allocator>
	  ::pointer data)
  {
    this->Clear();

    this->m_ = i;
    this->n_ = i;

    this->data_ = data;
  }


  /**********************************
   * ELEMENT ACCESS AND AFFECTATION *
   **********************************/


  //! Access operator.
  /*!
    Returns the value of element (i, j).
    \param i row index.
    \param j column index.
    \return Element (i, j) of the matrix.
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline typename Matrix_TriangPacked<T, Prop, Storage, Allocator>::value_type
  Matrix_TriangPacked<T, Prop, Storage, Allocator>::operator() (int i, int j)
  {

#ifdef SELDON_CHECK_BOUNDARIES
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_TriangPacked::operator()",
		     string("Index should be in [0, ") + to_str(this->m_-1)
		     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("Matrix_TriangPacked::operator()",
		     string("Index should be in [0, ") + to_str(this->n_-1)
		     + "], but is equal to " + to_str(j) + ".");
#endif

    if (Storage::UpLo())
      if (i > j)
	return 0;
      else
	return this->data_[Storage::GetFirst(i * this->n_
					     - (i * (i + 1)) / 2 + j,
					     (j * (j + 1)) / 2 + i)];
    else
      if (i < j)
	return 0;
      else
	return this->data_[Storage::GetFirst((i * (i + 1)) / 2 + j,
					     j * this->m_
					     - (j * (j + 1)) / 2 + i)];
  }

 
  //! Access operator.
  /*!
    Returns the value of element (i, j).
    \param i row index.
    \param j column index.
    \return Element (i, j) of the matrix.
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline typename Matrix_TriangPacked<T, Prop, Storage, Allocator>::value_type
  Matrix_TriangPacked<T, Prop, Storage, Allocator>
  ::operator() (int i, int j) const
  {
    
#ifdef SELDON_CHECK_BOUNDARIES
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_TriangPacked::operator()",
		     string("Index should be in [0, ") + to_str(this->m_-1)
		     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("Matrix_TriangPacked::operator()",
		     string("Index should be in [0, ") + to_str(this->n_-1)
		     + "], but is equal to " + to_str(j) + ".");
#endif
    
    if (Storage::UpLo())
      if (i > j)
	return 0;
      else
	return this->data_[Storage::GetFirst(i * this->n_
					     - (i * (i + 1)) / 2 + j,
					     (j * (j + 1)) / 2 + i)];
    else
      if (i < j)
	return 0;
      else
	return this->data_[Storage::GetFirst((i * (i + 1)) / 2 + j,
					     j * this->m_
					     - (j * (j + 1)) / 2 + i)];
  }


  //! Direct access method.
  /*!
    This method allows access to elements stored in memory, i.e. elements
    from the upper part. i <= j must be satisfied.
    \param i row index.
    \param j column index.
    \return The value of the matrix at (i, j).
  */
  template <class T, class Prop, class Storage, class Allocator>
  inline typename Matrix_TriangPacked<T, Prop, Storage, Allocator>::reference
  Matrix_TriangPacked<T, Prop, Storage, Allocator>::Val(int i, int j)
  {

#ifdef SELDON_CHECK_BOUNDARIES
    if (i < 0 || i >= this->m_)
      throw WrongRow("Matrix_TriangPacked::Val(int, int)",
		     string("Index should be in [0, ") + to_str(this->m_-1)
		     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("Matrix_TriangPacked::Val(int, int)",
		     string("Index should be in [0, ") + to_str(this->n_-1)
		     + "], but is equal to " + to_str(j) + ".");
    if (Storage::UpLo())
      {
	if (i > j)
	  throw WrongRow("Matrix_TriangPacked::Val(int, int)",
			 string("Attempted to access to element (")
			 + to_str(i) + ", " + to_str(j) + string(") but row")
			 + string(" index should not be strictly more")
			 + " than column index (upper triangular matrix).");
	return this->data_[Storage::GetFirst(i * this->n_
					     - (i * (i + 1)) / 2 + j,
					     (j * (j + 1)) / 2 + i)];
      }
    else
      {
	if (j > i)
	  throw WrongCol("Matrix_TriangPacked::Val(int, int)",
			 string("Attempted to access to element (")
			 + to_str(i) + ", " + to_str(j) + string(") but")
			 + string(" column index should not be strictly more")
			 + " than row index (lower triangular matrix).");
	return this->data_[Storage::GetFirst((i * (i + 1)) / 2 + j,
					     j * this->m_
					     - (j * (j + 1)) / 2 + i)];
      }
#endif

    if (Storage::UpLo())
      return this->data_[Storage::GetFirst(i * this->n_
					   - (i * (i + 1)) / 2 + j,
					   (j * (j + 1)) / 2 + i)];
    else
      return this->data_[Storage::GetFirst((i * (i + 1)) / 2 + j,
					   j * this->m_
					   - (j * (j + 1)) / 2 + i)];
  }

 
  //! Direct access method.