Ifpack_Condest.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_Condest.h"
00003 #include "Ifpack_CondestType.h"
00004 #include "Ifpack_Preconditioner.h"
00005 #include "Epetra_Vector.h"
00006 #include "Epetra_LinearProblem.h"
00007 #include "Epetra_Map.h"
00008 #include "Epetra_RowMatrix.h"
00009 #ifdef HAVE_IFPACK_AZTECOO
00010 #include "AztecOO.h"
00011 #endif
00012 
00013 double Ifpack_Condest(const Ifpack_Preconditioner& IFP,
00014               const Ifpack_CondestType CT,
00015               const int MaxIters,
00016               const double Tol,
00017               Epetra_RowMatrix* Matrix)
00018 {
00019   double ConditionNumberEstimate = -1.0;
00020 
00021   if (CT == Ifpack_Cheap) {
00022 
00023     // Create a vector with all values equal to one
00024     Epetra_Vector Ones(IFP.OperatorDomainMap());
00025     Ones.PutScalar(1.0);
00026     // Create the vector of results
00027     Epetra_Vector OnesResult(IFP.OperatorRangeMap());
00028     // Compute the effect of the solve on the vector of ones
00029     IFPACK_CHK_ERR(IFP.ApplyInverse(Ones, OnesResult)); 
00030     // Make all values non-negative
00031     IFPACK_CHK_ERR(OnesResult.Abs(OnesResult)); 
00032     // Get the maximum value across all processors
00033     IFPACK_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); 
00034 
00035   }
00036   else if (CT == Ifpack_CG) {
00037 
00038 #ifdef HAVE_IFPACK_AZTECOO
00039     if (Matrix == 0)
00040       Matrix = (Epetra_RowMatrix*)&(IFP.Matrix());
00041 
00042     Epetra_Vector LHS(IFP.OperatorDomainMap());
00043     LHS.PutScalar(0.0);
00044     Epetra_Vector RHS(IFP.OperatorRangeMap());
00045     RHS.Random();
00046     Epetra_LinearProblem Problem;
00047     Problem.SetOperator(Matrix);
00048     Problem.SetLHS(&LHS);
00049     Problem.SetRHS(&RHS);
00050 
00051     AztecOO Solver(Problem);
00052     Solver.SetAztecOption(AZ_output,AZ_none);
00053     Solver.SetAztecOption(AZ_solver,AZ_cg_condnum);
00054     Solver.Iterate(MaxIters,Tol);
00055 
00056     const double* status = Solver.GetAztecStatus();
00057     ConditionNumberEstimate = status[AZ_condnum];
00058 #endif
00059 
00060   } else if (CT == Ifpack_GMRES) {
00061 
00062 #ifdef HAVE_IFPACK_AZTECOO
00063     if (Matrix == 0)
00064       Matrix = (Epetra_RowMatrix*)&(IFP.Matrix());
00065 
00066     Epetra_Vector LHS(IFP.OperatorDomainMap());
00067     LHS.PutScalar(0.0);
00068     Epetra_Vector RHS(IFP.OperatorRangeMap());
00069     RHS.Random();
00070     Epetra_LinearProblem Problem;
00071     Problem.SetOperator(Matrix);
00072     Problem.SetLHS(&LHS);
00073     Problem.SetRHS(&RHS);
00074 
00075     AztecOO Solver(Problem);
00076     Solver.SetAztecOption(AZ_solver,AZ_gmres_condnum);
00077     Solver.SetAztecOption(AZ_output,AZ_none);
00078     // the following can be problematic for large problems,
00079     // but any restart would destroy useful information about
00080     // the condition number.
00081     Solver.SetAztecOption(AZ_kspace,MaxIters);
00082     Solver.Iterate(MaxIters,Tol);
00083 
00084     const double* status = Solver.GetAztecStatus();
00085     ConditionNumberEstimate = status[AZ_condnum];
00086 #endif
00087   }
00088 
00089   return(ConditionNumberEstimate);
00090 
00091 }

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1