/home/vivien/public_html/.src_seldon/matrix_sparse/Relaxation_MatVect.cxx

// Copyright (C) 2003-2009 Marc Duruflé
// Copyright (C) 2001-2009 Vivien Mallet
//
// This file is part of the linear-algebra library Seldon,
// http://seldon.sourceforge.net/.
//
// Seldon is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 2.1 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 Lesser General Public License for
// more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Seldon. If not, see http://www.gnu.org/licenses/.


#ifndef SELDON_FILE_RELAXATION_MATVECT_CXX

/*
  Functions defined in this file

  SOR(A, X, B, omega, iter, type_ssor)

*/

namespace Seldon
{


  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, RowSparse, Allocator0>& A,
           Vector<T2, Storage2, Allocator2>& X,
           const Vector<T1, Storage1, Allocator1>& B,
           const T3& omega, int iter, int type_ssor = 2)
  {
    T1 temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif

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

    // Forward sweep.
    if (type_ssor % 2 == 0)
      for (int i = 0; i < iter; i++)
        for (int j = 0; j < ma; j++)
          {
            temp = T1(0);
            ajj = A(j, j);
            for (int k = ptr[j] ; k < ptr[j+1]; k++)
              temp += data[k] * X(ind[k]);

            temp = B(j) - temp + ajj * X(j);
            X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
          }

    // Backward sweep.
    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        for (int j = ma-1 ; j >= 0; j--)
          {
            temp = T1(0);
            ajj = A(j, j);
            for (int k = ptr[j]; k < ptr[j+1]; k++)
              temp += data[k] * X(ind[k]);

            temp = B(j) - temp + ajj * X(j);
            X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
          }
  }



  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, ArrayRowSparse, Allocator0>& A,
           Vector<T2, Storage2, Allocator2>& X,
           const Vector<T1, Storage1, Allocator1>& B,
           const T3& omega,
           int iter, int type_ssor = 2)
  {
    T1 temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif

    T0 ajj(1);

    // Forward sweep.
    if (type_ssor % 2 == 0)
      for (int i = 0; i < iter; i++)
        for (int j = 0; j < ma; j++)
          {
            temp = T1(0);
            ajj = T0(0);
            for (int k = 0; k < A.GetRowSize(j); k++)
              {
                temp += A.Value(j,k) * X(A.Index(j,k));
                if (A.Index(j,k) == j)
                  ajj += A.Value(j,k);
              }

            temp = B(j) - temp + ajj * X(j);
            X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
          }

    // Backward sweep.
    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        for (int j = ma-1; j >= 0; j--)
          {
            temp = T1(0);
            ajj = T0(0);
            for (int k = 0; k < A.GetRowSize(j); k++)
              {
                temp += A.Value(j,k) * X(A.Index(j,k));
                if (A.Index(j,k) == j)
                  ajj += A.Value(j,k);
              }

            temp = B(j) - temp + ajj * X(j);
            X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
          }
  }



  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, RowSymSparse, Allocator0>& A,
           Vector<T2, Storage2, Allocator2>& X,
           const Vector<T1, Storage1, Allocator1>& B,
           const T3& omega,int iter, int type_ssor = 2)
  {
    T1 temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif

    int* ptr = A.GetPtr();
    int* ind = A.GetInd();
    T0* data = A.GetData();

    Vector<T2, Storage2, Allocator2> Y(ma);
    Y.Zero();
    T0 ajj(1);
    int p;
    T0 val(0);
    // Let us consider the following splitting : A = D - L - U
    // D diagonal of A
    // L lower part of A
    // U upper part of A, A is symmetric, so L = U^t

    // Forward sweep
    // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
    if (type_ssor % 2 == 0)
      for (int i = 0; i < iter; i++)
        {
          for (int j = 0; j < ma; j++)
            {
              // First we do X = (U + (1-omega)/omega D) X + B
              temp = T1(0);
              ajj = T0(0);
              for (int k = ptr[j]; k < ptr[j+1]; k++)
                {
                  p = ind[k];
                  val = data[k];
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }

              temp = B(j) - temp;
              X(j) = (T2(1) - omega) / omega * ajj * X(j) + temp;
            }

          for (int j = 0; j < ma; j++)
            {
              ajj = T0(0);
              // Then we solve (D/omega - L) X = X
              for (int k = ptr[j]; k < ptr[j+1]; k++)
                {
                  p = ind[k];
                  val = data[k];
                  if (p == j)
                    ajj += val;
                }
              X(j) *= omega / ajj;
              for (int k = ptr[j]; k < ptr[j+1]; k++)
                {
                  p = ind[k];
                  val = data[k];
                  if (p != j)
                    X(p) -= val*X(j);
                }
            }
        }


    // Backward sweep.
    // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        {
          Y.Zero();
          for (int j = 0; j < ma; j++)
            {
              ajj = T0(0);
              // Then we compute X = (L + (1-omega)/omega D) X + B.
              for (int k = ptr[j]; k < ptr[j+1]; k++)
                {
                  p = ind[k];
                  val = data[k];
                  if (p == j)
                    ajj += val;
                  else
                    Y(p) += val*X(j);
                }

              X(j) = (T2(1) - omega) / omega * ajj * X(j) + B(j) - Y(j);
            }

          for (int j = ma-1; j >= 0; j--)
            {
              temp = T1(0);
              ajj = T0(0);
              // Then we solve (D/omega - U) X = X.
              for (int k = ptr[j]; k < ptr[j+1]; k++)
                {
                  p = ind[k];
                  val = data[k];
                  if (p == j)
                    ajj += val;
                  else
                    temp += val*X(p);
                }
              X(j) = (X(j) - temp) * omega / ajj;
            }
        }
  }



  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, ArrayRowSymSparse, Allocator0>& A,
           Vector<T2, Storage2, Allocator2>& X,
           const Vector<T1, Storage1, Allocator1>& B,
           const T3& omega,int iter, int type_ssor = 2)
  {
    T1 temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif


    Vector<T2,Storage2,Allocator2> Y(ma);
    Y.Zero();
    T0 ajj(1);
    int p; T0 val(0);
    // Let us consider the following splitting : A = D - L - U
    // D diagonal of A
    // L lower part of A
    // U upper part of A, A is symmetric, so L = U^t

    // Forward sweep.
    // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
    if (type_ssor % 2 == 0)
      for (int i = 0; i < iter; i++)
        {
          for (int j = 0; j < ma; j++)
            {
              // First we do X = (U + (1-omega)/omega D) X + B.
              temp = T1(0);
              ajj = T0(0);
              for (int k = 0; k < A.GetRowSize(j); k++)
                {
                  p = A.Index(j,k);
                  val = A.Value(j,k);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }

              temp = B(j) - temp;
              X(j) = (T2(1) - omega) / omega * ajj * X(j) + temp;
            }

          for (int j = 0; j < ma; j++)
            {
              ajj = T0(0);
              // Then we solve (D/omega - L) X = X.
              for (int k = 0; k < A.GetRowSize(j); k++)
                {
                  p = A.Index(j,k);
                  val = A.Value(j,k);
                  if (p == j)
                    ajj += val;
                }
              X(j) *= omega / ajj;
              for (int k = 0 ; k < A.GetRowSize(j); k++)
                {
                  p = A.Index(j,k);
                  val = A.Value(j,k);
                  if (p != j)
                    X(p) -= val*X(j);
                }
            }
        }


    // Backward sweep.
    // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        {
          Y.Zero();
          for (int j = 0; j < ma; j++)
            {
              ajj = T0(0);
              // Then we compute X = (L + (1-omega)/omega D) X + B.
              for (int k = 0; k < A.GetRowSize(j); k++)
                {
                  p = A.Index(j,k);
                  val = A.Value(j,k);
                  if (p == j)
                    ajj += val;
                  else
                    Y(p) += val*X(j);
                }

              X(j) = (T2(1) - omega) / omega * ajj * X(j) + B(j) - Y(j);
            }

          for (int j = ma-1; j >= 0; j--)
            {
              temp = T1(0);
              ajj = T0(0);
              // Then we solve (D/omega - U) X = X.
              for (int k = 0; k < A.GetRowSize(j); k++)
                {
                  p = A.Index(j,k);
                  val = A.Value(j,k);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }
              X(j) = (X(j) - temp) * omega / ajj;
            }
        }

  }



  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, RowComplexSparse, Allocator0>& A,
           Vector<complex<T2>, Storage2, Allocator2>& X,
           const Vector<complex<T1>, Storage1, Allocator1>& B,
           const T3& omega, int iter, int type_ssor = 2)
  {
    complex<T1> temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif

    int* ptr_real = A.GetRealPtr();
    int* ptr_imag = A.GetImagPtr();
    int* ind_real = A.GetRealInd();
    int* ind_imag = A.GetImagInd();
    typename Matrix<T0, Prop0, RowComplexSparse, Allocator0>::pointer
      data_real = A.GetRealData();
    typename Matrix<T0, Prop0, RowComplexSparse, Allocator0>::pointer
      data_imag = A.GetImagData();
    complex<T0> ajj(1);

    if (type_ssor % 2 == 0)
      for (int i = 0; i < iter; i++)
        for (int j = 0; j < ma; j++)
          {
            temp = complex<T1>(0);
            ajj = A(j,j);
            for (int k = ptr_real[j]; k < ptr_real[j+1]; k++)
              temp += data_real[k] * X(ind_real[k]);

            for (int k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
              temp += complex<T1>(0, data_imag[k]) * X(ind_imag[k]);

            temp = B(j) - temp + ajj*X(j);
            X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
          }

    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        for (int j = ma-1; j >= 0; j--)
          {
            temp = complex<T1>(0);
            ajj = A(j,j);
            for (int k = ptr_real[j]; k < ptr_real[j+1]; k++)
              temp += data_real[k] * X(ind_real[k]);

            for (int k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
              temp += complex<T1>(0, data_imag[k]) * X(ind_imag[k]);

            temp = B(j) - temp + ajj * X(j);
            X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
          }
  }



  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, ArrayRowComplexSparse, Allocator0>& A,
           Vector<complex<T2>, Storage2, Allocator2>& X,
           const Vector<complex<T1>, Storage1, Allocator1>& B,
           const T3& omega, int iter, int type_ssor = 2)
  {
    complex<T1> temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif

    complex<T0> ajj(1);

    if (type_ssor%2 == 0)
      for (int i = 0; i < iter; i++)
        {
          for (int j = 0; j < ma; j++)
            {
              temp = complex<T1>(0);
              ajj = complex<T0>(0);
              for (int k = 0; k < A.GetRealRowSize(j); k++)
                {
                  temp += A.ValueReal(j,k) * X(A.IndexReal(j,k));
                  if (A.IndexReal(j,k) == j)
                    ajj += complex<T0>(A.ValueReal(j,k), 0);
                }
              for (int k = 0; k < A.GetImagRowSize(j); k++)
                {
                  temp += complex<T0>(0,A.ValueImag(j,k))
                    * X(A.IndexImag(j,k));
                  if (A.IndexImag(j,k) == j)
                    ajj += complex<T0>(0, A.ValueImag(j,k));
                }

              temp = B(j) - temp + ajj * X(j);
              X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
            }
        }

    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        for (int j = ma-1; j >= 0; j--)
          {
            temp = complex<T1>(0);
            ajj = complex<T0>(0);
            for (int k = 0; k < A.GetRealRowSize(j); k++)
              {
                temp += A.ValueReal(j,k) * X(A.IndexReal(j,k));
                if (A.IndexReal(j,k) == j)
                  ajj += complex<T0>(A.ValueReal(j,k), 0);
              }
            for (int k = 0; k < A.GetImagRowSize(j); k++)
              {
                temp += complex<T0>(0, A.ValueImag(j,k))
                  * X(A.IndexImag(j,k));
                if (A.IndexImag(j,k) == j)
                  ajj += complex<T0>(0, A.ValueImag(j,k));
              }

            temp = B(j) - temp + ajj * X(j);
            X(j) = (T2(1) - omega) * X(j) + omega * temp / ajj;
          }
  }



  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, RowSymComplexSparse, Allocator0>& A,
           Vector<complex<T2>, Storage2, Allocator2>& X,
           const Vector<complex<T1>, Storage1, Allocator1>& B,
           const T3& omega, int iter, int type_ssor = 2)
  {
    complex<T1> temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif

    int* ptr_real = A.GetRealPtr();
    int* ptr_imag = A.GetImagPtr();
    int* ind_real = A.GetRealInd();
    int* ind_imag = A.GetImagInd();
    T0* data_real = A.GetRealData();
    T0* data_imag = A.GetImagData();

    Vector<complex<T2>, Storage2, Allocator2> Y(ma);
    Y.Zero();
    complex<T0> ajj(1);
    int p;
    complex<T0> val(0);

    // Let us consider the following splitting : A = D - L - U
    // D diagonal of A
    // L lower part of A
    // U upper part of A, A is symmetric, so L = U^t
    // forward sweep
    // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
    if (type_ssor % 2 == 0)
      for (int i = 0; i < iter; i++)
        {
          for (int j = 0; j < ma; j++)
            {
              // first we do X = (U + (1-omega)/omega D) X + B
              temp = complex<T1>(0);
              ajj = complex<T0>(0);
              for (int k = ptr_real[j]; k < ptr_real[j+1]; k++)
                {
                  p = ind_real[k];
                  val = complex<T0>(data_real[k], 0);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }
              for (int k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
                {
                  p = ind_imag[k];
                  val = complex<T0>(0, data_imag[k]);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }

              temp = B(j) - temp;
              X(j) = (T2(1) - omega) / omega * ajj * X(j) + temp;
            }

          for (int j = 0; j < ma; j++)
            {
              ajj = complex<T0>(0);
              // Then we solve (D/omega - L) X = X.
              for (int k = ptr_real[j]; k < ptr_real[j+1]; k++)
                {
                  p = ind_real[k];
                  val = complex<T0>(data_real[k], 0);
                  if (p == j)
                    ajj += val;
                }

              for (int k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
                {
                  p = ind_imag[k];
                  val = complex<T0>(0, data_imag[k]);
                  if (p == j)
                    ajj += val;
                }
              X(j) *= omega / ajj;

              for (int k = ptr_real[j]; k < ptr_real[j+1]; k++)
                {
                  p = ind_real[k];
                  val = complex<T0>(data_real[k], 0);
                  if (p != j)
                    X(p) -= val*X(j);
                }

              for (int k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
                {
                  p = ind_imag[k];
                  val = complex<T0>(0, data_imag[k]);
                  if (p != j)
                    X(p) -= val*X(j);
                }
            }
        }

    // Backward sweep.
    // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        {
          Y.Zero();

          for (int j = 0; j < ma; j++)
            {
              ajj = complex<T0>(0);
              // Then we compute X = (L + (1-omega)/omega D) X + B.
              for (int k = ptr_real[j]; k < ptr_real[j+1]; k++)
                {
                  p = ind_real[k];
                  val = complex<T0>(data_real[k], 0);
                  if (p == j)
                    ajj += val;
                  else
                    Y(p) += val * X(j);
                }

              for (int k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
                {
                  p = ind_imag[k];
                  val = complex<T0>(0, data_imag[k]);
                  if (p == j)
                    ajj += val;
                  else
                    Y(p) += val * X(j);
                }
              X(j) = (T2(1) - omega) / omega * ajj * X(j) + B(j) - Y(j);
            }

          for (int j = (ma-1); j >= 0; j--)
            {
              temp = complex<T1>(0);
              ajj = complex<T0>(0);
              // Then we solve (D/omega - U) X = X.
              for (int k = ptr_real[j]; k < ptr_real[j+1]; k++)
                {
                  p = ind_real[k];
                  val = complex<T0>(data_real[k], 0);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }

              for (int k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
                {
                  p = ind_imag[k];
                  val = complex<T0>(0, data_imag[k]);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }
              X(j) = (X(j) - temp) * omega / ajj;
            }
        }
  }



  template <class T0, class Prop0, class Allocator0,
            class T1, class Storage1, class Allocator1,
            class T2, class Storage2, class Allocator2, class T3>
  void SOR(const Matrix<T0, Prop0, ArrayRowSymComplexSparse, Allocator0>& A,
           Vector<complex<T2>, Storage2, Allocator2>& X,
           const Vector<complex<T1>, Storage1, Allocator1>& B,
           const T3& omega, int iter, int type_ssor = 2)
  {
    complex<T1> temp(0);

    int ma = A.GetM();

#ifdef SELDON_CHECK_BOUNDS
    int na = A.GetN();
    if (na != ma)
      throw WrongDim("SOR", "Matrix must be squared.");

    if (ma != X.GetLength() || ma != B.GetLength())
      throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
#endif

    Vector<complex<T2>, Storage2, Allocator2> Y(ma);
    Y.Zero();
    complex<T0> ajj(1);
    int p;
    complex<T0> val(0);
    // Let us consider the following splitting : A = D - L - U
    // D diagonal of A
    // L lower part of A
    // U upper part of A, A is symmetric, so L = U^t
    // forward sweep
    // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
    if (type_ssor % 2 == 0)
      for (int i = 0; i < iter; i++)
        {
          for (int j = 0; j < ma; j++)
            {
              // First we do X = (U + (1-omega)/omega D) X + B
              temp = complex<T1>(0);
              ajj = complex<T0>(0);
              for (int k = 0; k < A.GetRealRowSize(j); k++)
                {
                  p = A.IndexReal(j,k);
                  val = complex<T0>(A.ValueReal(j,k), 0);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }
              for (int k = 0; k < A.GetImagRowSize(j); k++)
                {
                  p = A.IndexImag(j,k);
                  val = complex<T0>(0, A.ValueImag(j,k));
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }

              temp = B(j) - temp;
              X(j) = (T2(1) - omega) / omega * ajj * X(j) + temp;
            }

          for (int j = 0; j < ma; j++)
            {
              ajj = complex<T0>(0);
              // Then we solve (D/omega - L) X = X.
              for (int k = 0; k < A.GetRealRowSize(j); k++)
                {
                  p = A.IndexReal(j,k);
                  val = complex<T0>(A.ValueReal(j,k), 0);
                  if (p == j)
                    ajj += val;
                }
              for (int k = 0; k < A.GetImagRowSize(j); k++)
                {
                  p = A.IndexImag(j,k);
                  val = complex<T0>(0, A.ValueImag(j,k));
                  if (p == j)
                    ajj += val;
                }
              X(j) *= omega / ajj;
              for (int k = 0; k < A.GetRealRowSize(j); k++)
                {
                  p = A.IndexReal(j,k);
                  val = complex<T0>(A.ValueReal(j,k), 0);
                  if (p != j)
                    X(p) -= val * X(j);
                }
              for (int k = 0; k < A.GetImagRowSize(j); k++)
                {
                  p = A.IndexImag(j,k);
                  val = complex<T0>(0, A.ValueImag(j,k));
                  if (p != j)
                    X(p) -= val*X(j);
                }
            }
        }

    // Backward sweep.
    // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
    if (type_ssor % 3 == 0)
      for (int i = 0; i < iter; i++)
        {
          Y.Zero();

          for (int j = 0; j < ma; j++)
            {
              ajj = complex<T0>(0);
              // Then we compute X = (L + (1-omega)/omega D) X + B.
              for (int k = 0; k < A.GetRealRowSize(j); k++)
                {
                  p = A.IndexReal(j,k);
                  val = complex<T0>(A.ValueReal(j,k), 0);
                  if (p == j)
                    ajj += val;
                  else
                    Y(p) += val * X(j);
                }
              for (int k = 0; k < A.GetImagRowSize(j); k++)
                {
                  p = A.IndexImag(j,k);
                  val = complex<T0>(0, A.ValueImag(j,k));
                  if (p == j)
                    ajj += val;
                  else
                    Y(p) += val * X(j);
                }
              X(j) = (T2(1) - omega) / omega * ajj * X(j) + B(j) - Y(j);
            }

          for (int j = ma-1; j >= 0; j--)
            {
              temp = complex<T1>(0);
              ajj = complex<T0>(0);
              // Then we solve (D/omega - U) X = X.
              for (int k = 0; k < A.GetRealRowSize(j); k++)
                {
                  p = A.IndexReal(j,k);
                  val = complex<T0>(A.ValueReal(j,k), 0);
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }
              for (int k = 0; k < A.GetImagRowSize(j); k++)
                {
                  p = A.IndexImag(j,k);
                  val = complex<T0>(0, A.ValueImag(j,k));
                  if (p == j)
                    ajj += val;
                  else
                    temp += val * X(p);
                }
              X(j) = (X(j) - temp) * omega / ajj;
            }
        }
  }


} // end namespace

#define SELDON_FILE_RELAXATION_MATVECT_CXX
#endif