Method Test

comp/preconditioner.cpp:87–179

  
  void Preconditioner :: Test () const
  {
    cout << IM(1) << "Compute eigenvalues" << endl;
    const BaseMatrix & amat = GetAMatrix();
    const BaseMatrix & pre = GetMatrix();
    int eigenretval;

    if ( !this->uselapack )
      {
        EigenSystem eigen (amat, pre);
        eigen.SetPrecision(1e-30);
        eigen.SetMaxSteps(1000); 
        
        eigen.SetPrecision(1e-15);
        eigenretval = eigen.Calc();
        eigen.PrintEigenValues (*testout);
        cout << IM(1) << " Min Eigenvalue : "  << eigen.EigenValue(1) << endl; 
        cout << IM(1) << " Max Eigenvalue : " << eigen.MaxEigenValue() << endl; 
        cout << IM(1) << " Condition   " << eigen.MaxEigenValue()/eigen.EigenValue(1) << endl; 
        (*testout) << " Min Eigenvalue : "  << eigen.EigenValue(1) << endl; 
        (*testout) << " Max Eigenvalue : " << eigen.MaxEigenValue() << endl; 
        
        if(testresult_ok) *testresult_ok = eigenretval;
        if(testresult_min) *testresult_min = eigen.EigenValue(1);
        if(testresult_max) *testresult_max = eigen.MaxEigenValue();
        
        
        //    (*testout) << " Condition   " << eigen.MaxEigenValue()/eigen.EigenValue(1) << endl; 
        //    for (int i = 1; i < min2 (10, eigen.NumEigenValues()); i++)
        //      cout << "cond(i) = " << eigen.MaxEigenValue() / eigen.EigenValue(i) << endl;
        (*testout) << " Condition   " << eigen.MaxEigenValue()/eigen.EigenValue(1) << endl;
        
      }

    else
      {
#ifdef LAPACK   
        int n = amat.Height();
	int n_elim = 0;
	BitArray internaldofs (n);
	internaldofs.Clear();

	for ( int i = 0; i < n; i++ )
	  {
	    FlatArray rowindices = 
	      dynamic_cast<const BaseSparseMatrix&>(amat).GetRowIndices(i);
	    if ( rowindices.Size() <= 1 )
	      internaldofs.SetBit(i);
	    else
	      n_elim++;
	  }

        Matrix<Complex> mat(n_elim), mat2(n_elim), ev(n_elim);
        BaseVector & v1 = *amat.CreateColVector();
        BaseVector & v2 = *amat.CreateColVector();
        FlatVector<Complex> fv1 = v1.FVComplex();
        // FlatVector<Complex> fv2 = v2.FVComplex();