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
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "LOCA_Epetra_TransposeLinearSystem_ExplicitTranspose.H"
00043 #include "LOCA_GlobalData.H"
00044 #include "LOCA_ErrorCheck.H"
00045 #include "Teuchos_ParameterList.hpp"
00046 #include "NOX_Epetra_LinearSystem.H"
00047 #include "NOX_Epetra_Scaling.H"
00048 #include "Epetra_RowMatrix.h"
00049
00050 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00051 ExplicitTranspose(
00052 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00053 const Teuchos::RCP<Teuchos::ParameterList>& solverParams,
00054 const Teuchos::RCP<NOX::Epetra::LinearSystem>& linsys_) :
00055 globalData(global_data),
00056 linsys(linsys_),
00057 jac_trans(),
00058 prec_trans(),
00059 scaling_trans(),
00060 transposer(true)
00061 {
00062
00063 if (solverParams->isParameter("Transpose Scaling"))
00064 scaling_trans = (*solverParams).INVALID_TEMPLATE_QUALIFIER
00065 get< Teuchos::RCP<NOX::Epetra::Scaling> >("Transpose Scaling");
00066 }
00067
00068 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00069 ~ExplicitTranspose()
00070 {
00071 }
00072
00073 bool
00074 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00075 applyJacobianTransposeInverse(Teuchos::ParameterList ¶ms,
00076 const NOX::Epetra::Vector &input,
00077 NOX::Epetra::Vector &result)
00078 {
00079
00080 linsys->setJacobianOperatorForSolve(jac_trans);
00081 if (linsys->hasPreconditioner())
00082 linsys->setPrecOperatorForSolve(prec_trans);
00083
00084
00085 Teuchos::RCP<NOX::Epetra::Scaling> scaling_orig;
00086 if (scaling_trans != Teuchos::null) {
00087 scaling_orig = linsys->getScaling();
00088 linsys->resetScaling(scaling_trans);
00089 }
00090
00091
00092 bool res = linsys->applyJacobianInverse(params, input, result);
00093
00094
00095 if (scaling_trans != Teuchos::null)
00096 linsys->resetScaling(scaling_orig);
00097
00098 return res;
00099 }
00100
00101 bool
00102 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00103 createJacobianTranspose()
00104 {
00105
00106 Teuchos::RCP<Epetra_RowMatrix> jac =
00107 Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(linsys->getJacobianOperator());
00108
00109 if (jac == Teuchos::null)
00110 globalData->locaErrorCheck->throwError(
00111 string("LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::") +
00112 string("createJacobianTranspose()"),
00113 string("Jacobian operator must be of type Epetra_RowMatrix for ") +
00114 string("Explicit Transpose method"));
00115
00116
00117 if (jac_trans == Teuchos::null)
00118 jac_trans = Teuchos::rcp(&(transposer(*jac)), false);
00119 else
00120 transposer.fwd();
00121
00122 return true;
00123 }
00124
00125 bool
00126 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00127 createTransposePreconditioner(const NOX::Epetra::Vector& x,
00128 Teuchos::ParameterList& p)
00129 {
00130
00131 if (!linsys->hasPreconditioner())
00132 return true;
00133
00134
00135 bool res1 = linsys->destroyPreconditioner();
00136
00137
00138 Teuchos::RCP<NOX::Epetra::Scaling> scaling_orig;
00139 if (scaling_trans != Teuchos::null) {
00140 scaling_orig = linsys->getScaling();
00141 linsys->resetScaling(scaling_trans);
00142 }
00143
00144
00145 linsys->setJacobianOperatorForSolve(jac_trans);
00146 bool res2 = linsys->createPreconditioner(x, p, true);
00147 prec_trans = linsys->getGeneratedPrecOperator();
00148
00149
00150 if (scaling_trans != Teuchos::null)
00151 linsys->resetScaling(scaling_orig);
00152
00153 return res1 && res2;
00154 }
00155
00156 Teuchos::RCP<Epetra_Operator>
00157 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00158 getJacobianTransposeOperator()
00159 {
00160 return jac_trans;
00161 }
00162
00163 Teuchos::RCP<Epetra_Operator>
00164 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00165 getTransposePreconditioner()
00166 {
00167 return prec_trans;
00168 }
00169
00170 void
00171 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00172 setJacobianTransposeOperator(
00173 const Teuchos::RCP<Epetra_Operator>& new_jac_trans)
00174 {
00175 jac_trans = new_jac_trans;
00176 }
00177
00178 void
00179 LOCA::Epetra::TransposeLinearSystem::ExplicitTranspose::
00180 setTransposePreconditioner(
00181 const Teuchos::RCP<Epetra_Operator>& new_prec_trans)
00182 {
00183 prec_trans = new_prec_trans;
00184 }