RE: Re: RE: lots of problems VC7.1

//
//
//     $COPYRIGHT$
//                                                                           
//

#ifndef ITL_SSOR_H
#define ITL_SSOR_H

#include <vector>
#include <functional>
#include "itl/preconditioner/detail/preconditioner.h"
#include "mtl/meta_if.h"
#include "mtl/meta_equal.h"
#include "mtl/matrix.h"

namespace itl {

  using mtl::ROW_MAJOR; using mtl::COL_MAJOR;
  using mtl::generators_error;

  //:   SSOR preconditioner.
  //<codeblock>
  // Usage:
  //     Matrix A;
  //     SSOR<Matrix> precond(A);
  //     qmr(A, x, b, precond.left(), precond.right(), iter);
  //     cg(A, x, b, precond(), iter);
  //</codeblock>
  // Matrix has to be unsymmetric. 
  // For symmetric one, use incomlete cholesky.
  // Notes: The idea under a concrete Preconditioner such 
  // as Incomplete Cholesky is to create a Preconditioner
  // object to use in iterative methods. 
  //
  //!definition: ssor.h
  //!category: itl,functors
  //!component: type
  //!tparam: Matrix - Matrix 
  //!example: ssor.cc
  //
  template <class Matrix>
  class SSOR {
    typedef typename Matrix::value_type T;
    typedef typename Matrix::orientation Orien;
    typedef typename mtl::matrix< T, mtl::rectangle<>, 
				  mtl::compressed<int, mtl::external>,
				  Orien >::type LUMatrix;
        enum { Orien_id = Orien::id };
  public:

    typedef typename mtl::IF< EQUAL < Orien_id, ROW_MAJOR >::RET, 
                         preconditioner<LUMatrix, LUMatrix, mtl::lower, 
                                        mtl::unit_upper>,
            typename mtl::IF< EQUAL < Orien_id, COL_MAJOR >::RET, 
                         preconditioner<LUMatrix, LUMatrix, mtl::unit_lower,
                                        mtl::upper>,
                         generators_error
    >::RET
    >::RET  Precond;

    typedef typename mtl::IF< EQUAL < Orien_id, ROW_MAJOR >::RET, 
                         preconditioner1<LUMatrix, LUMatrix, mtl::lower, 
                                         mtl::unit_upper>,
            typename mtl::IF< EQUAL < Orien_id, COL_MAJOR >::RET, 
                         preconditioner1<LUMatrix, LUMatrix, mtl::unit_lower, 
                                         mtl::upper>,
                         generators_error
    >::RET
    >::RET  Left;

    typedef typename mtl::IF< EQUAL < Orien_id, ROW_MAJOR >::RET, 
                         preconditioner2<LUMatrix, LUMatrix, mtl::lower, 
                                         mtl::unit_upper>,
            typename mtl::IF< EQUAL < Orien_id, COL_MAJOR >::RET, 
                         preconditioner2<LUMatrix, LUMatrix, mtl::unit_lower, 
                                         mtl::upper>,
                         generators_error
    >::RET
    >::RET  Right;

  
    SSOR(const Matrix& A) 
      :L_val(A.nnz()),     U_val(A.nnz()),
       L_ind(A.nnz()),     U_ind(A.nnz()),
       L_ptr(A.nrows()+1), U_ptr(A.nrows()+1)
    {
      do_ssor(A, Orien());
    }

  private:
    void
    do_ssor(const Matrix& A, mtl::row_tag)
    {

      using std::upper_bound;
      using std::bind2nd;
      using std::multiplies;
      using std::transform;
      int L_loc=0, U_loc=0;
      L_ptr[0] = 0; 
      U_ptr[0] = 0;
#if 1
      typename Matrix::const_iterator A_i = A.begin();
      if (mtl::not_at(A_i,A.end())) do {
        typename Matrix::OneD A_row = *A_i;
	typedef typename Matrix::OneD::IndexArray::value_type IA_T;
        IA_T i=A_i.index(); 
        typename Matrix::OneD::IndexArray::const_iterator 
          diag_location = 
          upper_bound(A_row.nz_struct().begin(), A_row.nz_struct().end(), i);

        int d = diag_location - A_row.nz_struct().begin();

        std::copy(A_row.begin(), A_row.begin()+d, L_val.begin()+L_loc);
        std::copy(A_row.nz_struct().begin(), A_row.nz_struct().begin()+d, 
             L_ind.begin()+L_loc);

        L_loc += d;
        L_ptr[i+1] = L_loc; 

        T A_ii = *(A_row.begin()+d-1);
        A_ii = T(1)/A_ii;

        transform(A_row.begin()+d, A_row.end(), U_val.begin()+U_loc, 
                  bind2nd(multiplies<T>(), A_ii));

        std::copy(A_row.nz_struct().begin()+d, A_row.nz_struct().end(), 
             U_ind.begin()+U_loc);

        U_loc += A_row.nnz() - d;
        U_ptr[i+1] = U_loc;

        ++A_i;
      } while (mtl::not_at(A_i,A.end()));
#else
      //This is better way in term of general, however, It need to 
      // measure and compare the performance
      typename Matrix::const_iterator A_i = A.begin();
      if (mtl::not_at(A_i,A.end())) do {
	typename Matrix::size_type i = A_i.index();

	T A_ii;
	
        typename Matrix::OneD::const_iterator A_ij = (*A_i).begin();
	if (mtl::not_at(A_ij,(*A_i).end())) do {
	  typename Matrix::size_type j = A_ij.index();
	  if ( j > i ) {
	    U_val[U_loc] = *A_ij;
	    U_ind[U_loc] = j;
	    ++U_loc;
	  } else {
	    L_val[L_loc] = *A_ij;
	    L_ind[L_loc] = j;
	    ++L_loc;
	    if ( i == j ) A_ii = T(1) / *A_ij;
	  }
	  ++A_ij;
	} while (mtl::not_at(A_ij, (*A_i).end()));
	
	U_ptr[i+1] = U_loc;
	L_ptr[i+1] = L_loc;

	//scale upper parts
        transform(U_val.begin()+U_ptr[i], U_val.begin()+U_loc, 
		  U_val.begin()+U_ptr[i], bind2nd(multiplies<T>(), A_ii));

        ++A_i;
      } while (mtl::not_at(A_i,A.end()));
#endif
      L_val.resize(L_loc);
      U_val.resize(U_loc);
      L_ind.resize(L_loc);
      U_ind.resize(U_loc);

      L = LUMatrix(A.nrows(), A.ncols(), 
                   L_loc, &(*L_val.begin()), &(*L_ptr.begin()), 
                   &(*L_ind.begin()));

      U = LUMatrix(A.nrows(), A.ncols(), 
                   U_loc, &(*U_val.begin()), &(*U_ptr.begin()), 
                   &(*U_ind.begin()));
    }

    void
    do_ssor(const Matrix& A, mtl::column_tag)
    {
      using std::upper_bound;
      using std::bind2nd;
      using std::multiplies;
      using std::transform;
      int L_loc=0, U_loc=0;
      L_ptr[0] = 0; 
      U_ptr[0] = 0;

      typename Matrix::const_iterator A_i = A.begin();
      if (mtl::not_at(A_i,A.end())) do {
        const typename Matrix::OneD& A_col = *A_i;

        int i=A_i.index();
        
        typedef typename Matrix::OneD::IndexArray::const_iterator index_iter;
        
        index_iter diag_location = 
          upper_bound(A_col.nz_struct().begin(), A_col.nz_struct().end(), i);

        int d = diag_location - A_col.nz_struct().begin();

        std::copy(A_col.begin(), A_col.begin()+d, U_val.begin()+U_loc);
        std::copy(A_col.nz_struct().begin(), A_col.nz_struct().begin()+d, 
             U_ind.begin()+U_loc);

        U_loc += d;
        U_ptr[i+1] = U_loc; 
        T A_ii = *(A_col.begin()+d-1);
        A_ii = 1./A_ii;

        transform(A_col.begin()+d, A_col.end(), L_val.begin()+L_loc, 
                  bind2nd(multiplies<T>(), A_ii));

        std::copy(A_col.nz_struct().begin()+d, A_col.nz_struct().end(), 
             L_ind.begin()+L_loc);

        L_loc += A_col.nnz() - d;
        L_ptr[i+1] = L_loc;

        ++A_i;
      } while (mtl::not_at(A_i,A.end()));

      L_val.resize(L_loc);
      U_val.resize(U_loc);
      L_ind.resize(L_loc);
      U_ind.resize(U_loc);


      L = LUMatrix(A.nrows(), A.ncols(), 
                   L_loc, &(*L_val.begin()), &(*L_ptr.begin()), 
                   &(*L_ind.begin()));

      U = LUMatrix(A.nrows(), A.ncols(), 
                   U_loc, &(*U_val.begin()), &(*U_ptr.begin()), 
                   &(*U_ind.begin()));
    }
  
  public:
    //: return  a right or Left Preconditioner object.
    Precond operator()() { return Precond(L, U); }
    //: return the Left part of a Split Preconditioner
    Left left() { return Left(L, U); }
    //: return the Right part of a Split Preconditioner
    Right right() { return Right(L, U); }

  private:

    LUMatrix L, U;

    std::vector<T> L_val;
    std::vector<T> U_val;
    std::vector<int> L_ind;
    std::vector<int> U_ind;
    std::vector<int> L_ptr;
    std::vector<int> U_ptr;

  };


}

#endif

//
//     $COPYRIGHT$
//
//

#ifndef ITL_CHOLESKY_H
#define ITL_CHOLESKY_H

#include <vector>
#include <assert.h>
#include <algorithm>

#include "itl/itl.h"
#include "itl/preconditioner/detail/preconditioner.h"
#include "mtl/meta_if.h"



namespace itl {

  /* VC++ work around for ambiguous function */
   inline void check_symm(mtl::symmetric_tag) { }

    template < class Shape >
    inline void check_symm(Shape)
    {
      std::cout << "Matrix is not symmetric. Abort." << std::endl;
      assert(0);
    }


  //: Incomplete Cholesky Preconditioner.
  //  For use with symmetric matrices.
  //
  //<codeblock>
  // Usage:
  //    SymMatrix A;
  //    cholesky< SymMartix > precond(A);
  //    qmr(A, x, b, precond.left(), precond.right(), iter);
  //    cg(A, x, b, precond(), iter);
  //</codeblock>
  //
  //Notes: The idea under a concrete Preconditioner such 
  //as Incomplete Cholesky is to create a Preconditioner
  //object to use in iterative methods. 
  //
  //
  //!definition: cholesky.h
  //!example: cholesky.cc
  //!category: itl,functors
  //!component: type
  //!tparam: Matrix - A symmetric Matrix 
  //
  template < class Matrix >
  class cholesky {
    typedef typename Matrix::value_type T;
    typedef typename Matrix::orientation Orien;
    typedef Matrix SymMatrix;
    enum { Orien_id = Orien::id };

    typedef typename mtl::matrix< T, mtl::rectangle<>, 
                         mtl::compressed<int, mtl::external>, mtl::column_major >::type Matrix1; 

    typedef typename mtl::matrix< T, mtl::rectangle<>, 
                         mtl::compressed<int, mtl::external>, mtl::row_major >::type Matrix2; 

    typedef typename mtl::IF< EQUAL < Orien_id, mtl::ROW_MAJOR >::RET, 
                         Matrix2,
            typename mtl::IF< EQUAL < Orien_id, mtl::COL_MAJOR >::RET, 
                         Matrix1,
			 mtl::generators_error
    >::RET
    >::RET TriMatrix;


  public:

    typedef preconditioner < Matrix1, Matrix2, mtl::lower, mtl::upper> Precond;
    typedef preconditioner1< Matrix1, Matrix2, mtl::lower, mtl::upper> Left;
    typedef preconditioner2< Matrix1, Matrix2, mtl::lower, mtl::upper> Right;
  
    cholesky(const SymMatrix& A) 
      : Tri_val(A.nnz()), Tri_ind(A.nnz()), Tri_ptr(A.nrows()+1)
    {
          typedef typename mtl::matrix_traits<SymMatrix>::shape Shape;
      check_symm(Shape()); // JGS change to compile-time test
      do_cholesky(A, Orien()); 
    }
  private:
 
    void 
    do_cholesky(const SymMatrix& A, mtl::row_tag)
    {
      using std::copy;
      if ( A.is_upper() ) { 
        int Tri_loc= 0;
        Tri_ptr[0] = 0;
        typename SymMatrix::const_iterator A_i = A.begin();
        if (mtl::not_at(A_i,A.end())) do {
          typename SymMatrix::OneD A_row = *A_i;
          
          int i=A_i.index(); 
          
          std::copy(A_row.begin(), A_row.end(), Tri_val.begin()+Tri_loc);
          std::copy(A_row.nz_struct().begin(), A_row.nz_struct().end(), 
		    Tri_ind.begin()+Tri_loc);
          
          Tri_loc += A_row.nnz();
          Tri_ptr[i+1] = Tri_loc;
          
          ++A_i;
        } while (mtl::not_at(A_i,A.end()));
        
        int d, g, h, i, j, k, n = A.nrows();
        T z;
        
        for (k = 0; k < n - 1; k++) {
          d = Tri_ptr[k];
          z = Tri_val[d] = sqrt(Tri_val[d]);
          
          for (i = d + 1; i < Tri_ptr[k+1]; i++)
            Tri_val[i] /= z;
          
          for (i = d + 1; i < Tri_ptr[k+1]; i++) {
            z = Tri_val[i];
            h = Tri_ind[i];
            g = i;
            
            for (j = Tri_ptr[h] ; j < Tri_ptr[h+1]; j++)
              for ( ; g < Tri_ptr[k+1] && Tri_ind[g+1] <= Tri_ind[j]; g++)
                if (Tri_ind[g] == Tri_ind[j])
                  Tri_val[j] -= z * Tri_val[g];
          }
        }
        d = Tri_ptr[n-1];
        Tri_val[d] = sqrt(Tri_val[d]);
        

        Tri = TriMatrix(A.nrows(), A.ncols(), 
                        Tri_loc, &Tri_val[0], &Tri_ptr[0], 
                        &Tri_ind[0]);
      } else { 
        std::cout << "Warning: It is not so efficient as symmetric row-wise upper Matrix" << std::endl;
        assert(0);
      }
    }
    void 
    do_cholesky(const SymMatrix& A, mtl::column_tag)
    {
      using std::copy;
      if ( A.is_upper() ) {
        std::cout << "Warning: It is not so efficient as symmetric column-wise lower Matrix" << std::endl;
        assert(0);
      } else {
        int Tri_loc= 0;
        Tri_ptr[0] = 0;
        typename SymMatrix::const_iterator A_i = A.begin();
        if (mtl::not_at(A_i,A.end())) do {
          typename SymMatrix::OneD A_row = *A_i;

          int i=A_i.index(); 

          std::copy(A_row.begin(), A_row.end(), Tri_val.begin()+Tri_loc);
          std::copy(A_row.nz_struct().begin(), A_row.nz_struct().end(), 
		    Tri_ind.begin()+Tri_loc);

          Tri_loc += A_row.nnz();
          Tri_ptr[i+1] = Tri_loc;

          ++A_i;
        } while (mtl::not_at(A_i,A.end()));

        int d, g, h, i, j, k, n = A.nrows();
        T z;

        for (k = 0; k < n - 1; k++) {
          d = Tri_ptr[k];
          z = Tri_val[d] = sqrt(Tri_val[d]);

          for (i = d + 1; i < Tri_ptr[k+1]; i++)
            Tri_val[i] /= z;

          for (i = d + 1; i < Tri_ptr[k+1]; i++) {
            z = Tri_val[i];
            h = Tri_ind[i];
            g = i;

            for (j = Tri_ptr[h] ; j < Tri_ptr[h+1]; j++)
              for ( ; g < Tri_ptr[k+1] && Tri_ind[g+1] <= Tri_ind[j]; g++)
                if (Tri_ind[g] == Tri_ind[j])
                  Tri_val[j] -= z * Tri_val[g];
          }
        }
        d = Tri_ptr[n-1];
        Tri_val[d] = sqrt(Tri_val[d]);

        Tri = TriMatrix(A.nrows(), A.ncols(), 
                        Tri_loc, &Tri_val[0], &Tri_ptr[0], 
                        &Tri_ind[0]);
      }
    }

    inline Precond pre0(mtl::row_tag) {
      return  Precond(mtl::trans(Tri), Tri);
    }
    inline Precond pre0(mtl::column_tag) {
      return  Precond(Tri, mtl::trans(Tri));
    }

    inline Left pre1(mtl::row_tag) {
      return  Left(mtl::trans(Tri), Tri);
    }
    inline Left pre1(mtl::column_tag) {
      return  Left(Tri, mtl::trans(Tri));
    }
    inline Right pre2(mtl::row_tag) {
      return  Right(mtl::trans(Tri), Tri);
    }
    inline Right pre2(mtl::column_tag) {
      return  Right(Tri, mtl::trans(Tri));
    }


  public:
    //:return a left or right Preconditioner object
    inline Precond operator()() { 
      return pre0(Orien());
    }
    //: return a left part of Split Preconditioner object
    inline Left left() { return pre1(Orien()); }
    //: return a right part of Split Preconditioner object
    inline Right right() { return pre2(Orien()); }

    void print() {
    }

  private:
    TriMatrix Tri;
    std::vector<T> Tri_val;
    std::vector<int> Tri_ind;
    std::vector<int> Tri_ptr;
  };


}

#endif

// -*- c++ -*-
//
//     $COPYRIGHT$
//

#ifndef ITL_ILU_H
#define ITL_ILU_H


#include <algorithm>
#include <vector>
#include "mtl/entry.h"
#include "itl/preconditioner/detail/preconditioner.h"
#include "mtl/norm.h"


namespace itl {

  //:  Incomplete LU without fill-in Preconditioner.
  //<codeblock>
  //  Usage:
  //    Matrix A;
  //    ILU<Matrix> precond(A);
  //    qmr(A, x, b, precond.left(), precond.right(), iter);
  //</codeblock>
  // Matrix has to be unsymmetric. 
  //For symmetric one, use incomlete cholesky.
  //Notes: The idea under a concrete Preconditioner such 
  //as Incomplete Cholesky is to create a Preconditioner
  //object to use in iterative methods.
  //
  //!definition: ilu.h
  //!example: ilu.cc
  //!category: itl,functors
  //!component: type
  //!tparam: Matrix -  An MTL Matrix 
  //
  template <class Matrix>
  class ILU  {
    typedef typename Matrix::value_type T;
    typedef typename Matrix::size_type sizeT;
    typedef typename Matrix::orientation Orien;
    typedef typename mtl::matrix< T, mtl::rectangle<>, 
			 mtl::compressed<sizeT, mtl::external>, Orien >::type LUMatrix; 
  public:
    //: The preconditioner type
    typedef  preconditioner<LUMatrix, LUMatrix,
			    mtl::unit_lower, mtl::upper> Precond;
    //: The left preconditioner type
    typedef preconditioner1<LUMatrix, LUMatrix, 
			    mtl::unit_lower, mtl::upper> Left;
    //: The right preconditioner type
    typedef preconditioner2<LUMatrix, LUMatrix,
			    mtl::unit_lower, mtl::upper> Right;

    //: Default Constructor  
    ILU() {}
  
    //: Construct from Matrix A
    ILU(const Matrix& A_) 
      :L_val(A_.nnz()),     U_val(A_.nnz()),
       L_ind(A_.nnz()),     U_ind(A_.nnz()),
       L_ptr(A_.nrows()+1), U_ptr(A_.nrows()+1)
    {
      do_ilu(A_, Orien());
    }
  private:
    void 
    do_ilu(const Matrix& A_, mtl::column_tag)
    {
      using std::upper_bound;
      sizeT L_loc=0, U_loc=0;
      L_ptr[0] = 0; 
      U_ptr[0] = 0;

      typename Matrix::const_iterator A_i = A_.begin();
      if (mtl::not_at(A_i,A_.end())) do {
	typename Matrix::OneD A_row = *A_i;
	sizeT i = A_i.index();
	typename Matrix::OneD::IndexArray::const_iterator diag_location = 
	  std::upper_bound(A_row.nz_struct().begin(), 
			   A_row.nz_struct().end(), i);

	int d = diag_location - A_row.nz_struct().begin();

	std::copy(A_row.begin(), A_row.begin()+d, U_val.begin()+U_loc);
	std::copy(A_row.nz_struct().begin(), A_row.nz_struct().begin()+d, 
	     U_ind.begin()+U_loc);

	U_loc += d;
	U_ptr[i+1] = U_loc; 

	std::copy(A_row.begin()+d, A_row.end(), L_val.begin()+L_loc);
	std::copy(A_row.nz_struct().begin()+d, A_row.nz_struct().end(), 
	     L_ind.begin()+L_loc);

	L_loc += A_row.nnz() - d;
	L_ptr[i+1] = L_loc;

	++A_i;
      } while (mtl::not_at(A_i,A_.end()));

      L_val.resize(L_loc);
      U_val.resize(U_loc);
      L_ind.resize(L_loc);
      U_ind.resize(U_loc);

      sizeT i, j, qn, pn, rn; 
      for (i = 0; i < A_.nrows() - 1; i++) {
	T multiplier = U_val[U_ptr[i+1]-1];
    
	for (j = L_ptr[i]; j < L_ptr[i+1]; j++)
	  L_val[j] /= multiplier;
    
	for (j = U_ptr[i+1]; j < U_ptr[i+2]-1; j++) {
	  multiplier = U_val[j];
	  qn = j + 1;
	  rn = L_ptr[i+1];
	  for (pn = L_ptr[U_ind[j]]; 
	       L_ind[pn] <= i + 1 && pn < L_ptr[U_ind[j]+1]; 
	       pn++) {
	    while (U_ind[qn] < L_ind[pn] && qn < U_ptr[i+2])
	      qn++;
	    if (L_ind[pn] == U_ind[qn] && qn < U_ptr[i+2])
	      U_val[qn] -= multiplier * L_val[pn];
	  }
	  for (; pn < L_ptr[U_ind[j]+1]; pn++) {
	    while (L_ind[rn] < L_ind[pn] && rn < L_ptr[i+2])
	      rn++;
	    if (L_ind[pn] == L_ind[rn] && rn < L_ptr[i+2])
	      L_val[rn] -= multiplier * L_val[pn];
	  }
	}
      }

      L = LUMatrix(A_.nrows(), A_.ncols(), 
		   L_loc, &L_val[0], &L_ptr[0], 
		   &L_ind[0]);
      
      U = LUMatrix(A_.nrows(), A_.ncols(), 
		   U_loc, &U_val[0], &U_ptr[0], 
		   &U_ind[0]);
    }

    void
    do_ilu(const Matrix A_, mtl::row_tag)
    {
      sizeT L_loc=0, U_loc=0;
      L_ptr[0] = 0; 
      U_ptr[0] = 0;

      typename Matrix::const_iterator A_i = A_.begin();
      if (mtl::not_at(A_i,A_.end())) do {
	typename Matrix::OneD A_row = *A_i;
	typedef typename Matrix::OneD::IndexArray::value_type IA_T;
	IA_T i = A_i.index();
	typename Matrix::OneD::IndexArray::const_iterator diag_location = 
	  std::lower_bound(A_row.nz_struct().begin(), 
			   A_row.nz_struct().end(), i);

	int d = diag_location - A_row.nz_struct().begin();

	std::copy(A_row.begin(), A_row.begin()+d, L_val.begin()+L_loc);
	std::copy(A_row.nz_struct().begin(), A_row.nz_struct().begin()+d, 
	     L_ind.begin()+L_loc);

	L_loc += d;
	L_ptr[i+1] = L_loc; 

	std::copy(A_row.begin()+d, A_row.end(), U_val.begin()+U_loc);
	std::copy(A_row.nz_struct().begin()+d, A_row.nz_struct().end(), 
	     U_ind.begin()+U_loc);

	U_loc += A_row.nnz() - d;
	U_ptr[i+1] = U_loc;

	++A_i;
      } while (mtl::not_at(A_i,A_.end()));

      L_val.resize(L_loc);
      U_val.resize(U_loc);
      L_ind.resize(L_loc);
      U_ind.resize(U_loc);

      sizeT i, j, qn, pn, rn; 
      for (i = 1; i < A_.nrows(); i++) {
	for (j = L_ptr[i]; j < L_ptr[i+1]; j++) {
	  pn = U_ptr[L_ind[j]];

	  T multiplier = (L_val[j] /= U_val[pn]);

	  qn = j + 1;
	  rn = U_ptr[i];

	  for (pn++; U_ind[pn] < i && pn < U_ptr[L_ind[j]+1]; pn++) {
	    while (L_ind[qn] < U_ind[pn] && qn < L_ptr[i+1])
	      qn++;
	    if (U_ind[pn] == L_ind[qn] && qn < L_ptr[i+1])
	      L_val[qn] -= multiplier * U_val[pn];
	  }
	  for (; pn < U_ptr[L_ind[j]+1]; pn++) {
	    while (U_ind[rn] < U_ind[pn] && rn < U_ptr[i+1])
	      rn++;
	    if (U_ind[pn] == U_ind[rn] && rn < U_ptr[i+1])
	      U_val[rn] -= multiplier * U_val[pn];
	  }
	}
      }


      L = LUMatrix(A_.nrows(), A_.ncols(), 
		   L_loc, &(*L_val.begin()), &(*L_ptr.begin()), 
		   &(*L_ind.begin()));

      U = LUMatrix(A_.nrows(), A_.ncols(), 
		   U_loc, &(*U_val.begin()), &(*U_ptr.begin()), 
		   &(*U_ind.begin()));

    }

  public:
  

    //:return a left or right Preconditioner object.
    inline Precond operator()()  { return Precond(L, U); }
    //: return the left part of a Split Preconditioner object
    inline Left left() { return Left(L, U); }
    //: return the right part of a Split Preconditioner object
    inline Right right() { return Right(L, U); }
    void print() {
      print_all_matrix(L);
      print_all_matrix(U);
    }

  private:
    LUMatrix L; 
    LUMatrix U;
    std::vector<T> L_val;
    std::vector<T> U_val;
    std::vector<sizeT> L_ind;
    std::vector<sizeT> U_ind;
    std::vector<sizeT> L_ptr;
    std::vector<sizeT> U_ptr;
  };


}

#endif