00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceLinearSolveDriver.hpp"
00032 #include "PlayaLinearSolverDecl.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "PlayaTabs.hpp"
00035
00036
00037 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00038 #include "PlayaLinearSolverImpl.hpp"
00039 #endif
00040
00041
00042
00043
00044 using namespace Sundance;
00045 using namespace Sundance;
00046 using namespace Sundance;
00047 using namespace Sundance;
00048 using namespace Sundance;
00049 using namespace Sundance;
00050 using namespace Sundance;
00051 using namespace Teuchos;
00052 using namespace Playa;
00053 using namespace std;
00054
00055
00056 SolverState<double>
00057 LinearSolveDriver::solve(const LinearSolver<double>& solver,
00058 const LinearOperator<double>& A,
00059 const Array<Vector<double> >& rhs,
00060 const Array<RCP<DiscreteSpace> >& solutionSpace,
00061 const Array<Array<string> >& names,
00062 int verb,
00063 Expr& soln) const
00064 {
00065 Tabs tab(0);
00066 Array<Vector<double> > solnVec(rhs.size());
00067 SolverState<double> state;
00068
00069 for (int i=0; i<rhs.size(); i++)
00070 {
00071 Tabs tab1;
00072
00073 solnVec[i] = rhs[i].copy();
00074
00075 SUNDANCE_MSG2(verb, tab1 << "solving with RHS #" << i
00076 << " of " << rhs.size());
00077
00078 state = solver.solve(A, rhs[i], solnVec[i]);
00079
00080 SUNDANCE_MSG2(verb, tab1 << "solve completed with status="
00081 << state.stateDescription());
00082
00083
00084 if (state.finalState() != SolveConverged)
00085 {
00086 TeuchosOStringStream ss;
00087 ss << "Solve failed! state = "
00088 << state.stateDescription()
00089 << "\nmessage=" << state.finalMsg()
00090 << "\niters taken = " << state.finalIters()
00091 << "\nfinal residual = " << state.finalResid();
00092
00093
00094 if (dumpBadMatrix())
00095 {
00096 if (A.ptr().get() != 0)
00097 {
00098 ofstream osA(badMatrixFilename().c_str());
00099 A.print(osA);
00100 ss << "\nmatrix written to " << badMatrixFilename();
00101 }
00102 else
00103 {
00104 ss << "\nthe matrix is null! Evil is afoot in your code...";
00105 }
00106 if (rhs[i].ptr().get() != 0)
00107 {
00108 ofstream osb(badVectorFilename().c_str());
00109 rhs[i].print(osb);
00110 ss << "\nRHS vector written to " << badVectorFilename();
00111 }
00112 else
00113 {
00114 ss << "\nthe RHS vector is null! Evil is afoot in your code...";
00115 }
00116 }
00117
00118
00119 TEUCHOS_TEST_FOR_EXCEPTION(solveFailureIsFatal(),
00120 std::runtime_error, TEUCHOS_OSTRINGSTREAM_GET_C_STR(ss));
00121
00122
00123 return state;
00124 }
00125 }
00126
00127
00128 if (soln.ptr().get()==0)
00129 {
00130 soln = formSolutionExpr(solnVec, solutionSpace, names, verb);
00131 }
00132 else
00133 {
00134 writeIntoSolutionExpr(solnVec, soln, verb);
00135 }
00136
00137 return state;
00138
00139 }
00140
00141
00142
00143 Expr LinearSolveDriver::formSolutionExpr(
00144 const Array<Vector<double> >& solnVector,
00145 const Array<RCP<DiscreteSpace> >& solutionSpace,
00146 const Array<Array<string> >& names,
00147 int verb) const
00148 {
00149 Array<Expr> cols(solnVector.size());
00150
00151 for (int m=0; m<solnVector.size(); m++)
00152 {
00153 Array<Expr> col(solutionSpace.size());
00154
00155 for (int i=0; i<col.size(); i++)
00156 {
00157 std::string name = names[m][i];
00158 if (col.size() > 1) name += "[" + Teuchos::toString(i) + "]";
00159 col[i] = new DiscreteFunction(*(solutionSpace[i]),
00160 solnVector[m].getBlock(i), name);
00161 }
00162 if (col.size() > 1)
00163 {
00164 cols[m] = new ListExpr(col);
00165 }
00166 else
00167 {
00168 cols[m] = col[0];
00169 }
00170 }
00171
00172 if (cols.size() > 1)
00173 {
00174 return new ListExpr(cols);;
00175 }
00176 else
00177 {
00178 return cols[0];
00179 }
00180 }
00181
00182
00183 void LinearSolveDriver::writeIntoSolutionExpr(
00184 const Array<Vector<double> >& solnVector,
00185 Expr soln,
00186 int verb) const
00187 {
00188 #ifdef BLAH
00189 TEUCHOS_TEST_FOR_EXCEPTION(soln.size() != solnVector.size(),
00190 std::runtime_error,
00191 "soln=" << soln << " soln.size()=" << soln.size() << " while solnVector.size()=" << solnVector.size());
00192 #endif
00193 for (int i=0; i<solnVector.size(); i++)
00194 {
00195 Expr u = soln[i];
00196 setDiscreteFunctionVector(u, solnVector[i]);
00197 }
00198 }