functions_matvect.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_FUNCTIONS_MATVECT_CXX

/*
  Functions defined in this file:

  A*U -> U
  Mlt(A, U)

  alpha.A*U + beta.V -> V
  MltAdd(alpha, A, U, beta, V)
*/

namespace Seldon
{


  /////////
  // Mlt //


  template <class T0, class Prop0, class Storage0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const Matrix<T0, Prop0, Storage0, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& B)
  {
    int ma = A.GetM();
    int na = A.GetN();

#ifdef SELDON_CHECK_BOUNDARIES
    if (na != B.GetM() || ma != B.GetM())
      throw WrongDim("Mlt(const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> C(B);

    T1 zero(0);
    T1 temp;

    for (int i = 0; i < ma; i++)
      {
	temp = zero;
	for (int j = 0; j < na; j++)
	  temp += A(i, j) * C(j);
	B(i) = temp;
      }
  }


  /*** Sparse matrices ***/


  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const Matrix<T0, Prop0, RowSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDARIES
    int na = A.GetN();

    if (na != C.GetM())
      throw WrongDim("Mlt(const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);

    T1 zero(0);
    T1 temp;

    int* ptr = A.GetPtr();
    int* ind = A.GetInd();
    typename Matrix<T0, Prop0, RowSparse, Allocator0>::pointer
      data = A.GetData();

    for (int i = 0; i < ma; i++)
      {
	temp = zero;
	for (int j = ptr[i]; j < ptr[i+1]; j++)
	  temp += data[j] * B(ind[j]);
	C(i) = temp;
      }
  }


  /*** Complex sparse matrices ***/


  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const Matrix<T0, Prop0, RowComplexSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int i, j;

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDARIES
    int na = A.GetN();

    if (na != C.GetM())
      throw WrongDim("Mlt(const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);

    T1 zero(0);
    T1 temp;

    int* real_ptr = A.GetRealPtr();
    int* imag_ptr = A.GetImagPtr();
    int* real_ind = A.GetRealInd();
    int* imag_ind = A.GetImagInd();
    typename Matrix<T0, Prop0, RowComplexSparse, Allocator0>::pointer
      real_data = A.GetRealData();
    typename Matrix<T0, Prop0, RowComplexSparse, Allocator0>::pointer
      imag_data = A.GetImagData();

    for (i = 0; i < ma; i++)
      {
	temp = zero;
	for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
	  temp += real_data[j] * B(real_ind[j]);
	C(i) = temp;
      }

    for (i = 0; i < ma; i++)
      {
	temp = zero;
	for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
	  temp += complex<T0>(T0(0), imag_data[j]) * B(imag_ind[j]);
	C(i) += temp;
      }
  }


  /*** Symmetric sparse matrices ***/


  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const Matrix<T0, Prop0, RowSymSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDARIES
    int na = A.GetN();

    if (na != C.GetM())
      throw WrongDim("Mlt(const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);

    int i, j;
    T1 zero(0);
    T1 temp;

    int* ptr = A.GetPtr();
    int* ind = A.GetInd();
    typename Matrix<T0, Prop0, RowSymSparse, Allocator0>::pointer
      data = A.GetData();

    for (i = 0; i < ma; i++)
      {
	temp = zero;
	for (j = ptr[i]; j < ptr[i + 1]; j++)
	  temp += data[j] * B(ind[j]);
	C(i) = temp;
      }
    for (i = 0; i < ma-1; i++)
      for (j = ptr[i]; j < ptr[i + 1]; j++)
	if (ind[j] != i)
	  C(ind[j]) += data[j] * B(i);
  }


  /*** Symmetric complex sparse matrices ***/


  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const Matrix<T0, Prop0, RowSymComplexSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDARIES
    if (A.GetN() != C.GetM())
      throw WrongDim("Mlt(const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);

    int i, j;
    T1 zero(0);
    T1 temp;

    int* real_ptr = A.GetRealPtr();
    int* imag_ptr = A.GetImagPtr();
    int* real_ind = A.GetRealInd();
    int* imag_ind = A.GetImagInd();
    typename Matrix<T0, Prop0, RowSymComplexSparse, Allocator0>::pointer
      real_data = A.GetRealData();
    typename Matrix<T0, Prop0, RowSymComplexSparse, Allocator0>::pointer
      imag_data = A.GetImagData();

    for (i = 0; i<ma; i++)
      {
	temp = zero;
	for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
	  temp += real_data[j] * B(real_ind[j]);
	C(i) = temp;
      }
    for (i = 0; i<ma-1; i++)
      for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
	if (real_ind[j] != i)
	  C(real_ind[j]) += real_data[j] * B(i);

    for (i = 0; i < ma; i++)
      {
	temp = zero;
	for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
	  temp += complex<T0>(T0(0), imag_data[j]) * B(imag_ind[j]);
	C(i) += temp;
      }
    for (i = 0; i<ma-1; i++)
      for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
	if (imag_ind[j] != i)
	  C(imag_ind[j]) += complex<T0>(T0(0), imag_data[j]) * B(i);
  }


  /*** Sparse matrices, *Trans ***/


  // NoTrans.
  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const class_SeldonNoTrans& Trans,
	   const Matrix<T0, Prop0, RowSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    Mlt(A, C);
  }


  // Trans.
  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const class_SeldonTrans& Trans,
	   const Matrix<T0, Prop0, RowSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int i, j;

    int ma = A.GetM();
    int na = A.GetN();

#ifdef SELDON_CHECK_BOUNDARIES
    if (ma != C.GetM())
      throw WrongDim("Mlt(class_SeldonTrans, const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);
    C.Zero();

    int* ptr = A.GetPtr();
    int* ind = A.GetInd();
    typename Matrix<T0, Prop0, RowSparse, Allocator0>::pointer
      data = A.GetData();

    for (i = 0; i < ma; i++)
      for (j = ptr[i]; j < ptr[i + 1]; j++)
	C(ind[j]) += data[j] * B(i);
  }


  // ConjTrans.
  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const class_SeldonConjTrans& Trans,
	   const Matrix<T0, Prop0, RowSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int i, j;

    int ma = A.GetM();
    int na = A.GetN();

#ifdef SELDON_CHECK_BOUNDARIES
    if (ma != C.GetM())
      throw WrongDim("Mlt(class_SeldonConjTrans, const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);
    C.Zero();

    int* ptr = A.GetPtr();
    int* ind = A.GetInd();
    typename Matrix<complex<T0>, Prop0, RowSparse, Allocator0>::pointer
      data = A.GetData();

    for (i = 0; i < ma; i++)
      for (j = ptr[i]; j < ptr[i + 1]; j++)
	C(ind[j]) += conj(data[j]) * B(i);
  }


  /*** Complex sparse matrices, *Trans ***/


  // NoTrans.
  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const class_SeldonNoTrans& Trans,
	   const Matrix<T0, Prop0, RowComplexSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    Mlt(A, C);
  }


  // Trans.
  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const class_SeldonTrans& Trans,
	   const Matrix<T0, Prop0, RowComplexSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int i, j;

    int ma = A.GetM();
    int na = A.GetN();

#ifdef SELDON_CHECK_BOUNDARIES
    if (na != C.GetM())
      throw WrongDim("Mlt(class_SeldonTrans, const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);
    C.Zero();

    int* real_ptr = A.GetRealPtr();
    int* imag_ptr = A.GetImagPtr();
    int* real_ind = A.GetRealInd();
    int* imag_ind = A.GetImagInd();
    typename Matrix<T0, Prop0, RowComplexSparse, Allocator0>::pointer
      real_data = A.GetRealData();
    typename Matrix<T0, Prop0, RowComplexSparse, Allocator0>::pointer
      imag_data = A.GetImagData();

    for (i = 0; i < ma; i++)
      for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
	C(real_ind[j]) += real_data[j] * B(i);

    for (i = 0; i < ma; i++)
      for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
	C(imag_ind[j]) += complex<T0>(T0(0), imag_data[j]) * B(i);
  }


  // ConjTrans.
  template <class T0, class Prop0, class Allocator0,
	    class T1, class Storage1, class Allocator1>
  void Mlt(const class_SeldonConjTrans& Trans,
	   const Matrix<T0, Prop0, RowComplexSparse, Allocator0>& A,
	   Vector<T1, Storage1, Allocator1>& C)
  {
    int i, j;

    int ma = A.GetM();
    int na = A.GetN();

#ifdef SELDON_CHECK_BOUNDARIES
    if (na != C.GetM())
      throw WrongDim("Mlt(class_SeldonConjTrans, const Matrix, Vector)",
		     "Unable to multiply matrix and vector.");
#endif

    Vector<T1, Storage1, Allocator1> B(C);
    C.Zero();

    int* real_ptr = A.GetRealPtr();
    int* imag_ptr = A.GetImagPtr();
    int* real_ind = A.GetRealInd();