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_AnasaziOperator_Cayley.H"
00043 #include "Teuchos_ParameterList.hpp"
00044 #include "LOCA_GlobalData.H"
00045 #include "LOCA_ErrorCheck.H"
00046
00047 LOCA::AnasaziOperator::Cayley::Cayley(
00048 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00049 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00050 const Teuchos::RCP<Teuchos::ParameterList>& eigenParams_,
00051 const Teuchos::RCP<Teuchos::ParameterList>& solverParams_,
00052 const Teuchos::RCP<LOCA::TimeDependent::AbstractGroup>& grp_)
00053 : globalData(global_data),
00054 myLabel("Cayley Transformation"),
00055 eigenParams(eigenParams_),
00056 solverParams(solverParams_),
00057 grp(grp_),
00058 tmp_r(),
00059 tmp_i(),
00060 sigma(0.0),
00061 mu(0.0)
00062 {
00063 sigma = eigenParams->get("Cayley Pole",0.0);
00064 mu = eigenParams->get("Cayley Zero",0.0);
00065 }
00066
00067 LOCA::AnasaziOperator::Cayley::~Cayley()
00068 {
00069 }
00070
00071 const string&
00072 LOCA::AnasaziOperator::Cayley::label() const
00073 {
00074 return myLabel;
00075 }
00076
00077 void
00078 LOCA::AnasaziOperator::Cayley::apply(const NOX::Abstract::MultiVector& input,
00079 NOX::Abstract::MultiVector& output) const
00080 {
00081 string callingFunction =
00082 "LOCA::AnasaziOperator::Cayley::apply()";
00083
00084 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00085 NOX::Abstract::Group::ReturnType status;
00086
00087
00088 if (tmp_r == Teuchos::null || tmp_r->numVectors() != input.numVectors())
00089 tmp_r = input.clone(NOX::ShapeCopy);
00090
00091
00092 status = grp->computeShiftedMatrix(1.0, -mu);
00093 finalStatus =
00094 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00095 finalStatus,
00096 callingFunction);
00097
00098
00099 status = grp->applyShiftedMatrixMultiVector(input, *tmp_r);
00100 finalStatus =
00101 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00102 finalStatus,
00103 callingFunction);
00104
00105
00106 status = grp->computeShiftedMatrix(1.0, -sigma);
00107 finalStatus =
00108 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00109 finalStatus,
00110 callingFunction);
00111
00112
00113 status = grp->applyShiftedMatrixInverseMultiVector(*solverParams, *tmp_r,
00114 output);
00115
00116 finalStatus =
00117 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00118 finalStatus,
00119 callingFunction);
00120 }
00121
00122 void
00123 LOCA::AnasaziOperator::Cayley::preProcessSeedVector(NOX::Abstract::MultiVector& ivec)
00124 {
00125
00126 string callingFunction =
00127 "LOCA::AnasaziOperator::Cayley::preProcessSeedVector()";
00128
00129 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00130 NOX::Abstract::Group::ReturnType status;
00131
00132
00133 if (tmp_r == Teuchos::null || tmp_r->numVectors() != ivec.numVectors())
00134 tmp_r = ivec.clone(NOX::ShapeCopy);
00135
00136
00137 status = grp->computeShiftedMatrix(0.0, 1.0);
00138 finalStatus =
00139 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00140 finalStatus,
00141 callingFunction);
00142
00143
00144 status = grp->applyShiftedMatrixMultiVector(ivec, *tmp_r);
00145 finalStatus =
00146 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00147 finalStatus,
00148 callingFunction);
00149
00150
00151 status = grp->computeShiftedMatrix(1.0, -sigma);
00152 finalStatus =
00153 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00154 finalStatus,
00155 callingFunction);
00156
00157
00158 status = grp->applyShiftedMatrixInverseMultiVector(*solverParams, *tmp_r,
00159 ivec);
00160
00161 finalStatus =
00162 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00163 finalStatus,
00164 callingFunction);
00165 }
00166
00167 void
00168 LOCA::AnasaziOperator::Cayley::transformEigenvalue(double& ev_r,
00169 double& ev_i) const
00170 {
00171
00172 double mag = (1.0 - ev_r)*(1.0 - ev_r) + ev_i*ev_i;
00173 ev_r = (sigma*(ev_r*ev_r + ev_i*ev_i) - (sigma+mu)*ev_r + mu) / mag;
00174 ev_i = (mu-sigma)*ev_i/mag;
00175 }
00176
00177 NOX::Abstract::Group::ReturnType
00178 LOCA::AnasaziOperator::Cayley::rayleighQuotient(
00179 const NOX::Abstract::Vector& evec_r,
00180 const NOX::Abstract::Vector& evec_i,
00181 double& rq_r, double& rq_i) const
00182 {
00183 string callingFunction =
00184 "LOCA::AnasaziOperator::Cayley::rayleighQuotient()";
00185
00186
00187 if (tmp_r == Teuchos::null)
00188 tmp_r = evec_r.createMultiVector(1, NOX::ShapeCopy);
00189 if (tmp_i == Teuchos::null)
00190 tmp_i = evec_i.createMultiVector(1, NOX::ShapeCopy);
00191
00192 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00193 NOX::Abstract::Group::ReturnType status;
00194
00195
00196 status = grp->computeJacobian();
00197 finalStatus =
00198 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00199 callingFunction);
00200
00201
00202 status = grp->applyJacobian(evec_r, (*tmp_r)[0]);
00203 finalStatus =
00204 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00205 callingFunction);
00206
00207 status = grp->applyJacobian(evec_i, (*tmp_i)[0]);
00208 finalStatus =
00209 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00210 callingFunction);
00211
00212 rq_r = evec_r.innerProduct((*tmp_r)[0]) + evec_i.innerProduct((*tmp_i)[0]);
00213 rq_i = evec_r.innerProduct((*tmp_i)[0]) - evec_i.innerProduct((*tmp_r)[0]);
00214
00215
00216 status = grp->computeShiftedMatrix(0.0, 1.0);
00217 finalStatus =
00218 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00219 callingFunction);
00220
00221
00222 status = grp->applyShiftedMatrix(evec_r, (*tmp_r)[0]);
00223 finalStatus =
00224 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00225 callingFunction);
00226
00227 status = grp->applyShiftedMatrix(evec_i, (*tmp_i)[0]);
00228 finalStatus =
00229 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00230 callingFunction);
00231
00232 double m_r =
00233 evec_r.innerProduct((*tmp_r)[0]) + evec_i.innerProduct((*tmp_i)[0]);
00234 double m_i =
00235 evec_r.innerProduct((*tmp_i)[0]) - evec_i.innerProduct((*tmp_r)[0]);
00236 double m = m_r*m_r + m_i*m_i;
00237
00238
00239 double tmp = (rq_r*m_r + rq_i*m_i) / m;
00240 rq_i = (rq_i*m_r - rq_r*m_i) / m;
00241 rq_r = tmp;
00242
00243 return finalStatus;
00244 }