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_LowRankUpdateOp.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::LowRankUpdateOp::LowRankUpdateOp(
00051 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00052 const Teuchos::RCP<Epetra_Operator>& jacOperator,
00053 const Teuchos::RCP<const Epetra_MultiVector>& U_multiVec,
00054 const Teuchos::RCP<const Epetra_MultiVector>& V_multiVec,
00055 bool setup_for_solve) :
00056 globalData(global_data),
00057 label("LOCA::Epetra::LowRankUpdateOp"),
00058 localMap(V_multiVec->NumVectors(), 0, jacOperator->Comm()),
00059 J(jacOperator),
00060 U(U_multiVec),
00061 V(V_multiVec),
00062 useTranspose(false),
00063 tmpMat(),
00064 JinvU(),
00065 lu(),
00066 ipiv(),
00067 lapack()
00068 {
00069 if (setup_for_solve) {
00070 int m = U->NumVectors();
00071
00072
00073 JinvU = Teuchos::rcp(new Epetra_MultiVector(U->Map(), m));
00074 J->ApplyInverse(*U, *JinvU);
00075
00076
00077 lu = Teuchos::rcp(new Epetra_MultiVector(localMap, m));
00078 lu->Multiply('T', 'N', 1.0, *V, *JinvU, 0.0);
00079 for (int i=0; i<m; i++)
00080 (*lu)[i][i] += 1.0;
00081
00082
00083 ipiv.resize(m);
00084 int info;
00085 lapack.GETRF(m, m, lu->Values(), m, &ipiv[0], &info);
00086 }
00087 }
00088
00089 LOCA::Epetra::LowRankUpdateOp::~LowRankUpdateOp()
00090 {
00091 }
00092
00093 int
00094 LOCA::Epetra::LowRankUpdateOp::SetUseTranspose(bool UseTranspose)
00095 {
00096 useTranspose = UseTranspose;
00097 return J->SetUseTranspose(UseTranspose);
00098 }
00099
00100 int
00101 LOCA::Epetra::LowRankUpdateOp::Apply(const Epetra_MultiVector& Input,
00102 Epetra_MultiVector& Result) const
00103 {
00104
00105 int m = Input.NumVectors();
00106 double n = Input.GlobalLength();
00107
00108
00109 int res = J->Apply(Input, Result);
00110
00111
00112 if (tmpMat == Teuchos::null || tmpMat->NumVectors() != m) {
00113 tmpMat = Teuchos::rcp(new Epetra_MultiVector(localMap, m, false));
00114 }
00115
00116 if (!useTranspose) {
00117
00118
00119 tmpMat->Multiply('T', 'N', 1.0, *V, Input, 0.0);
00120
00121
00122 Result.Multiply('N', 'N', 1.0, *U, *tmpMat, 1.0);
00123
00124 }
00125 else {
00126
00127
00128 tmpMat->Multiply('T', 'N', 1.0/n, *U, Input, 0.0);
00129
00130
00131 Result.Multiply('N', 'N', 1.0, *V, *tmpMat, 1.0);
00132
00133 }
00134
00135 return res;
00136 }
00137
00138 int
00139 LOCA::Epetra::LowRankUpdateOp::ApplyInverse(const Epetra_MultiVector& Input,
00140 Epetra_MultiVector& Result) const
00141 {
00142
00143 int k = Input.NumVectors();
00144
00145
00146 int m = U->NumVectors();
00147
00148
00149 int res = J->ApplyInverse(Input, Result);
00150
00151
00152 if (tmpMat == Teuchos::null || tmpMat->NumVectors() != k) {
00153 tmpMat = Teuchos::rcp(new Epetra_MultiVector(localMap, k, false));
00154 }
00155
00156 if (!useTranspose) {
00157
00158
00159 tmpMat->Multiply('T', 'N', 1.0, *V, Result, 0.0);
00160
00161
00162 int info;
00163 lapack.GETRS('N', m, k, lu->Values(), m, &ipiv[0], tmpMat->Values(), m,
00164 &info);
00165
00166
00167 Result.Multiply('N', 'N', -1.0, *JinvU, *tmpMat, 1.0);
00168 }
00169 else {
00170 globalData->locaErrorCheck->throwError(
00171 "LOCA::Epetra::LowRankUpdateOp::ApplyInverse",
00172 "Operator does not support transpose");
00173 return -1;
00174 }
00175
00176 return res;
00177 }
00178
00179 double
00180 LOCA::Epetra::LowRankUpdateOp::NormInf() const
00181 {
00182 double Jn;
00183 vector<double> un(U->NumVectors());
00184 vector<double> vn(V->NumVectors());
00185
00186 Jn = J->NormInf();
00187 U->NormInf(&un[0]);
00188 V->NormInf(&vn[0]);
00189
00190 for (unsigned int i=0; i<un.size(); i++)
00191 Jn += un[i];
00192 for (unsigned int i=0; i<vn.size(); i++)
00193 Jn += vn[i];
00194
00195 return Jn;
00196 }
00197
00198
00199 const char*
00200 LOCA::Epetra::LowRankUpdateOp::Label () const
00201 {
00202 return const_cast<char*>(label.c_str());
00203 }
00204
00205 bool
00206 LOCA::Epetra::LowRankUpdateOp::UseTranspose() const
00207 {
00208 return useTranspose;
00209 }
00210
00211 bool
00212 LOCA::Epetra::LowRankUpdateOp::HasNormInf() const
00213 {
00214 return J->HasNormInf();
00215 }
00216
00217 const Epetra_Comm &
00218 LOCA::Epetra::LowRankUpdateOp::Comm() const
00219 {
00220 return J->Comm();
00221 }
00222 const Epetra_Map&
00223 LOCA::Epetra::LowRankUpdateOp::OperatorDomainMap() const
00224 {
00225 return J->OperatorDomainMap();
00226 }
00227
00228 const Epetra_Map&
00229 LOCA::Epetra::LowRankUpdateOp::OperatorRangeMap() const
00230 {
00231 return J->OperatorRangeMap();
00232 }