template <typename T>
void PetscLinearSolver<T>::init ( PetscMatrix<T>* matrix )
{
  // Initialize the data structures if not done so already.
  if (!this->initialized())
    {
      this->_is_initialized = true;
      
      int ierr=0;

// 2.1.x & earlier style      
#if PETSC_VERSION_LESS_THAN(2,2,0)
      
      // Create the linear solver context
      ierr = SLESCreate (libMesh::COMM_WORLD, &_sles);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Create the Krylov subspace & preconditioner contexts
      ierr = SLESGetKSP       (_sles, &_ksp);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = SLESGetPC        (_sles, &_pc);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Set user-specified  solver and preconditioner types
      this->set_petsc_solver_type();
      
      // Set the options from user-input
      // Set runtime options, e.g.,
      //      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
      //  These options will override those specified above as long as
      //  SLESSetFromOptions() is called _after_ any other customization
      //  routines.
      
      ierr = SLESSetFromOptions (_sles);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);

// 2.2.0 & newer style
#else
// ATTENTION      
      // Create the linear solver context
      ierr = KSPCreate (libMesh::COMM_WORLD, &_ksp);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
    
      //ierr = PCCreate (libMesh::COMM_WORLD, &_pc);
      //     CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Create the preconditioner context
      ierr = KSPGetPC        (_ksp, &_pc);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Set operators. The input matrix works as the preconditioning matrix
      ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),DIFFERENT_NONZERO_PATTERN);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);       

      // Set user-specified  solver and preconditioner types
      this->set_petsc_solver_type();

      // Set the options from user-input
      // Set runtime options, e.g.,
      //      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
      //  These options will override those specified above as long as
      //  KSPSetFromOptions() is called _after_ any other customization
      //  routines.
      
      ierr = KSPSetFromOptions (_ksp);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
      // NOT NECESSARY!!!!
      //ierr = PCSetFromOptions (_pc);
      //CHKERRABORT(libMesh::COMM_WORLD,ierr);

	       
#endif
	     
      // Have the Krylov subspace method use our good initial guess
      // rather than 0, unless the user requested a KSPType of
      // preonly, which complains if asked to use initial guesses.
#if PETSC_VERSION_LESS_THAN(3,0,0)
      KSPType ksp_type;
#else
      const KSPType ksp_type;
#endif

      ierr = KSPGetType (_ksp, &ksp_type);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);

      if (strcmp(ksp_type, "preonly"))
        {
          ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
                 CHKERRABORT(libMesh::COMM_WORLD,ierr);
        }
      
      // Notify PETSc of location to store residual history.
      // This needs to be called before any solves, since
      // it sets the residual history length to zero.  The default
      // behavior is for PETSc to allocate (internally) an array
      // of size 1000 to hold the residual norm history.
      ierr = KSPSetResidualHistory(_ksp,
				   PETSC_NULL,   // pointer to the array which holds the history
				   PETSC_DECIDE, // size of the array holding the history
				   PETSC_TRUE);  // Whether or not to reset the history for each solve. 
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);      
      if(this->_preconditioner)
      {
        this->_preconditioner->set_matrix(*matrix);
        PCShellSetContext(_pc,(void*)this->_preconditioner);
        PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
        PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
      }
    }
}


template <typename T>
std::pair<unsigned int, Real> 
PetscLinearSolver<T>::solve (SparseMatrix<T>&  matrix_in,
			     SparseMatrix<T>&  precond_in,
			     NumericVector<T>& solution_in,
			     NumericVector<T>& rhs_in,
			     const double tol,
			     const unsigned int m_its)
{
  START_LOG("solve()", "PetscLinearSolver");
  
  // Make sure the data passed in are really of Petsc types
  PetscMatrix<T>* matrix   = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
  PetscMatrix<T>* precond  = libmesh_cast_ptr<PetscMatrix<T>*>(&precond_in);
  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
  PetscVector<T>* rhs      = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);

  this->init (matrix);

  int ierr=0;
  int its=0, max_its = static_cast<int>(m_its);
  PetscReal final_resid=0.;

  // Close the matrices and vectors in case this wasn't already done.
  matrix->close ();
  precond->close ();
  solution->close ();
  rhs->close ();

//   // If matrix != precond, then this means we have specified a
//   // special preconditioner, so reset preconditioner type to PCMAT.
//   if (matrix != precond)
//     {
//       this->_preconditioner_type = USER_PRECOND;
//       this->set_petsc_preconditioner_type ();
//     }

// 2.1.x & earlier style      
#if PETSC_VERSION_LESS_THAN(2,2,0)

  if(_restrict_solve_to_is!=NULL)
    {
      libmesh_not_implemented();
    }
      
  // Set operators. The input matrix works as the preconditioning matrix
  ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
			  DIFFERENT_NONZERO_PATTERN);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  // Set the tolerances for the iterative solver.  Use the user-supplied
  // tolerance for the relative residual & leave the others at default values.
  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
 			   PETSC_DEFAULT, max_its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);


  // Solve the linear system
  ierr = SLESSolve (_sles, rhs->vec(), solution->vec(), &its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);


  // Get the norm of the final residual to return to the user.
  ierr = KSPGetResidualNorm (_ksp, &final_resid);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

// 2.2.0
#elif PETSC_VERSION_LESS_THAN(2,2,1)
      
  if(_restrict_solve_to_is!=NULL)
    {
      libmesh_not_implemented();
    }
      
  // Set operators. The input matrix works as the preconditioning matrix
  //ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
	 //			 SAME_NONZERO_PATTERN);
	 //       CHKERRABORT(libMesh::COMM_WORLD,ierr);


  // Set the tolerances for the iterative solver.  Use the user-supplied
  // tolerance for the relative residual & leave the others at default values.
  // Convergence is detected at iteration k if
  // ||r_k||_2 < max(rtol*||b||_2 , abstol)
  // where r_k is the residual vector and b is the right-hand side.  Note that
  // it is the *maximum* of the two values, the larger of which will almost
  // always be rtol*||b||_2. 
  ierr = KSPSetTolerances (_ksp,
			   tol,           // rtol   = relative decrease in residual  (1.e-5)
			   PETSC_DEFAULT, // abstol = absolute convergence tolerance (1.e-50)
 			   PETSC_DEFAULT, // dtol   = divergence tolerance           (1.e+5)
			   max_its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);


  // Set the solution vector to use
  ierr = KSPSetSolution (_ksp, solution->vec());
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Set the RHS vector to use
  ierr = KSPSetRhs (_ksp, rhs->vec());
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
   
  // Solve the linear system
  ierr = KSPSolve (_ksp);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Get the number of iterations required for convergence
  ierr = KSPGetIterationNumber (_ksp, &its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Get the norm of the final residual to return to the user.
  ierr = KSPGetResidualNorm (_ksp, &final_resid);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
// 2.2.1 & newer style
#else
      
  Mat submat = NULL;
  Mat subprecond = NULL;
  Vec subrhs = NULL;
  Vec subsolution = NULL;
  VecScatter scatter = NULL;
  PetscMatrix<Number>* subprecond_matrix = NULL;

  // Set operators.  Also restrict rhs and solution vector to
  // subdomain if neccessary.
  if(_restrict_solve_to_is!=NULL) //ATTENTION: _restrict_solve_to_is is false
    {
      size_t is_local_size = this->_restrict_solve_to_is_local_size();

      ierr = VecCreate(libMesh::COMM_WORLD,&subrhs);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecSetFromOptions(subrhs);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      ierr = VecCreate(libMesh::COMM_WORLD,&subsolution);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecSetFromOptions(subsolution);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

#if PETSC_VERSION_LESS_THAN(3,1,0)
      ierr = MatGetSubMatrix(matrix->mat(),
			     _restrict_solve_to_is,_restrict_solve_to_is,
			     PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = MatGetSubMatrix(precond->mat(),
			     _restrict_solve_to_is,_restrict_solve_to_is,
			     PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
      ierr = MatGetSubMatrix(matrix->mat(),
			     _restrict_solve_to_is,_restrict_solve_to_is,
			     MAT_INITIAL_MATRIX,&submat);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = MatGetSubMatrix(precond->mat(),
			     _restrict_solve_to_is,_restrict_solve_to_is,
			     MAT_INITIAL_MATRIX,&subprecond);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif

      /* Since removing columns of the matrix changes the equation
	 system, we will now change the right hand side to compensate
	 for this.  Note that this is not necessary if \p SUBSET_ZERO
	 has been selected.  */
      if(_subset_solve_mode!=SUBSET_ZERO)
	{
	  _create_complement_is(rhs_in);
	  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;

	  Vec subvec1 = NULL;
	  Mat submat1 = NULL;
	  VecScatter scatter1 = NULL;

	  ierr = VecCreate(libMesh::COMM_WORLD,&subvec1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  ierr = VecSetFromOptions(subvec1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);

	  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);

	  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);

	  ierr = VecScale(subvec1,-1.0);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);


#if PETSC_VERSION_LESS_THAN(3,1,0)
	  ierr = MatGetSubMatrix(matrix->mat(),
				 _restrict_solve_to_is,_restrict_solve_to_is_complement,
				 PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
	  ierr = MatGetSubMatrix(matrix->mat(),
				 _restrict_solve_to_is,_restrict_solve_to_is_complement,
				 MAT_INITIAL_MATRIX,&submat1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif

	  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);

	  ierr = VecScatterDestroy(scatter1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  ierr = VecDestroy(subvec1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  ierr = MatDestroy(submat1);
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	}

      ierr = KSPSetOperators(_ksp, submat, subprecond,
			     this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      if(this->_preconditioner)
	{
	  subprecond_matrix = new PetscMatrix<Number>(subprecond);
	  this->_preconditioner->set_matrix(*subprecond_matrix);
	}
    }
  else
    {/*
      ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
			     this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);


      if(this->_preconditioner)
	this->_preconditioner->set_matrix(matrix_in);*/
	
	//ATTENTION
	if (this->same_preconditioner)
	{
		ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(), SAME_PRECONDITIONER);
		std::cout << "***** DEBUG: SAME_PRECONDITIONER is applied. "<< std::endl;
	}
	else
	{
		// ATTENTION: my system matrix depends on the frequency, for each frequency I have then multiple RHS as well.
		//            for each frequency I call this solve() function with same_preconditioner=false at first, for the rest of the
		//            RHS the function is call with same_preconditioner=true. weired thing is if the solve() has been called, then 
		//						I cannot use the DIFFERENT_NONZERO_PATTERN any longer, however, the option SAME_NONZERO_PATTERN works fine.
		ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(), SAME_NONZERO_PATTERN); //DIFFERENT_NONZERO_PATTERN
// 		ierr = PCSetFromOptions (_pc);
// 		CHKERRABORT(libMesh::COMM_WORLD,ierr);
		std::cout << "***** DEBUG: DIFFERENT_NONZERO_PATTERN is applied. "<< std::endl;
		
	}
	CHKERRABORT(libMesh::COMM_WORLD,ierr);


	if(this->_preconditioner) //ATTENTION: the variable _preconditioner is NULL in this case
		this->_preconditioner->set_matrix(matrix_in);
    }

  // Set the tolerances for the iterative solver.  Use the user-supplied
  // tolerance for the relative residual & leave the others at default values.
  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
 			   PETSC_DEFAULT, max_its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
		 
//   ierr = KSPSetFromOptions (_ksp);
//       CHKERRABORT(libMesh::COMM_WORLD,ierr);
  // Solve the linear system
  if(_restrict_solve_to_is!=NULL)
    {
      ierr = KSPSolve (_ksp, subrhs, subsolution);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }
  else
    {
      ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }
	 
  // Get the number of iterations required for convergence
  ierr = KSPGetIterationNumber (_ksp, &its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Get the norm of the final residual to return to the user.
  ierr = KSPGetResidualNorm (_ksp, &final_resid);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  if(_restrict_solve_to_is!=NULL) 
    {
      switch(_subset_solve_mode)
	{
	case SUBSET_ZERO:
	  ierr = VecZeroEntries(solution->vec());
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  break;

	case SUBSET_COPY_RHS:
	  ierr = VecCopy(rhs->vec(),solution->vec());
	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  break;

	case SUBSET_DONT_TOUCH:
	  /* Nothing to do here.  */
	  break;

	}
      ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      ierr = VecScatterDestroy(scatter);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      if(this->_preconditioner)
	{
	  /* Before we delete subprecond_matrix, we should give the
	     _preconditioner a different matrix.  */
	  this->_preconditioner->set_matrix(matrix_in);
	  delete subprecond_matrix;
	  subprecond_matrix = NULL;
	}

      ierr = VecDestroy(subsolution);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecDestroy(subrhs);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = MatDestroy(submat);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = MatDestroy(subprecond);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }


#endif

  STOP_LOG("solve()", "PetscLinearSolver");
  // return the # of its. and the final residual norm.
  return std::make_pair(its, final_resid);
}