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 "Epetra_config.h"
00043 #include "LOCA_Epetra_CompactWYOp.H"
00044
00045 #include "Epetra_Map.h"
00046 #include "Epetra_Comm.h"
00047 #include "LOCA_GlobalData.H"
00048 #include "LOCA_ErrorCheck.H"
00049
00050 LOCA::Epetra::CompactWYOp::CompactWYOp(
00051 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00052 const Teuchos::RCP<const Epetra_Operator>& jacOperator,
00053 const Teuchos::RCP<const Epetra_MultiVector>& A_multiVec,
00054 const Teuchos::RCP<const Epetra_MultiVector>& Y_x_multiVec,
00055 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& Y_p_matrix,
00056 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& T_matrix) :
00057 globalData(global_data),
00058 label("LOCA::Epetra::CompactWYOp"),
00059 localMap(Y_x_multiVec->NumVectors(), 0, jacOperator->Comm()),
00060 J(jacOperator),
00061 A(A_multiVec),
00062 Y_x(Y_x_multiVec),
00063 Y_p(View, localMap, Y_p_matrix->values(), Y_p_matrix->stride(),
00064 Y_p_matrix->numCols()),
00065 T(View, localMap, T_matrix->values(), T_matrix->stride(),
00066 T_matrix->numCols()),
00067 tmpMat1(NULL),
00068 tmpMV(NULL),
00069 dblas()
00070 {
00071 }
00072
00073 LOCA::Epetra::CompactWYOp::~CompactWYOp()
00074 {
00075 if (tmpMat1 != NULL)
00076 delete tmpMat1;
00077 if (tmpMV != NULL)
00078 delete tmpMV;
00079 }
00080
00081 int
00082 LOCA::Epetra::CompactWYOp::SetUseTranspose(bool UseTranspose)
00083 {
00084 if (UseTranspose == false)
00085 return 0;
00086 else {
00087 globalData->locaErrorCheck->throwError(
00088 "LOCA::Epetra::CompactWYOp::SetUseTranspose",
00089 "Operator does not support transpose");
00090 return -1;
00091 }
00092 }
00093
00094 int
00095 LOCA::Epetra::CompactWYOp::Apply(const Epetra_MultiVector& Input,
00096 Epetra_MultiVector& Result) const
00097 {
00098 if (tmpMat1 == NULL || tmpMV == NULL) {
00099 globalData->locaErrorCheck->throwError(
00100 "LOCA::Epetra::CompactWYOp::Apply()",
00101 "Must call init() before Apply()!");
00102 return -1;
00103 }
00104
00105
00106 applyCompactWY(Input, *tmpMV, *tmpMat1);
00107
00108
00109 J->Apply(*tmpMV, Result);
00110
00111
00112 if (A.get() != NULL)
00113 Result.Multiply('N', 'N', 1.0, *A, *tmpMat1, 1.0);
00114
00115 return 0;
00116 }
00117
00118 int
00119 LOCA::Epetra::CompactWYOp::ApplyInverse(const Epetra_MultiVector& cInput,
00120 Epetra_MultiVector& Result) const
00121 {
00122 globalData->locaErrorCheck->throwError(
00123 "LOCA::Epetra::CompactWYOp::ApplyInverse",
00124 "Operator does not support ApplyInverse");
00125 return -1;
00126 }
00127
00128 double
00129 LOCA::Epetra::CompactWYOp::NormInf() const
00130 {
00131 double Jn;
00132 vector<double> an(A->NumVectors());
00133
00134 Jn = J->NormInf();
00135 A->NormInf(&an[0]);
00136
00137 for (unsigned int i=0; i<an.size(); i++)
00138 Jn += an[i];
00139
00140 return Jn;
00141 }
00142
00143
00144 const char*
00145 LOCA::Epetra::CompactWYOp::Label () const
00146 {
00147 return const_cast<char*>(label.c_str());
00148 }
00149
00150 bool
00151 LOCA::Epetra::CompactWYOp::UseTranspose() const
00152 {
00153 return false;
00154 }
00155
00156 bool
00157 LOCA::Epetra::CompactWYOp::HasNormInf() const
00158 {
00159 return J->HasNormInf();
00160 }
00161
00162 const Epetra_Comm &
00163 LOCA::Epetra::CompactWYOp::Comm() const
00164 {
00165 return J->Comm();
00166 }
00167 const Epetra_Map&
00168 LOCA::Epetra::CompactWYOp::OperatorDomainMap() const
00169 {
00170 return J->OperatorDomainMap();
00171 }
00172
00173 const Epetra_Map&
00174 LOCA::Epetra::CompactWYOp::OperatorRangeMap() const
00175 {
00176 return J->OperatorRangeMap();
00177 }
00178
00179 void
00180 LOCA::Epetra::CompactWYOp::init(const Epetra_MultiVector& x)
00181 {
00182 if (tmpMat1 != NULL)
00183 delete tmpMat1;
00184 if (tmpMV != NULL)
00185 delete tmpMV;
00186
00187 tmpMat1 = NULL;
00188 tmpMV = NULL;
00189
00190 tmpMat1 = new Epetra_MultiVector(localMap, 1, false);
00191 tmpMV = new Epetra_MultiVector(x.Map(), 1, false);
00192 }
00193
00194 void
00195 LOCA::Epetra::CompactWYOp::finish()
00196 {
00197 if (tmpMat1 != NULL)
00198 delete tmpMat1;
00199 if (tmpMV != NULL)
00200 delete tmpMV;
00201
00202 tmpMat1 = NULL;
00203 tmpMV = NULL;
00204 }
00205
00206 void
00207 LOCA::Epetra::CompactWYOp::applyCompactWY(const Epetra_MultiVector& x,
00208 Epetra_MultiVector& result_x,
00209 Epetra_MultiVector& result_p) const
00210 {
00211
00212 result_p.Multiply('T', 'N', 1.0, *Y_x, x, 0.0);
00213
00214
00215 dblas.TRMM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00216 Teuchos::NON_UNIT_DIAG, result_p.MyLength(),
00217 result_p.NumVectors(), 1.0, T.Values(), T.MyLength(),
00218 result_p.Values(), result_p.MyLength());
00219
00220
00221 result_x = x;
00222 result_x.Multiply('N', 'N', 1.0, *Y_x, result_p, 1.0);
00223
00224
00225 dblas.TRMM(Teuchos::LEFT_SIDE, Teuchos::LOWER_TRI, Teuchos::NO_TRANS,
00226 Teuchos::UNIT_DIAG, result_p.MyLength(),
00227 result_p.NumVectors(), 1.0, Y_p.Values(), Y_p.MyLength(),
00228 result_p.Values(), result_p.MyLength());
00229
00230 }