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_BorderedSolver_LowerTriangularBlockElimination.H"
00043 #include "LOCA_BorderedSolver_AbstractOperator.H"
00044 #include "LOCA_GlobalData.H"
00045 #include "LOCA_ErrorCheck.H"
00046 #include "LOCA_MultiContinuation_ConstraintInterface.H"
00047 #include "LOCA_MultiContinuation_MultiVecConstraint.H"
00048 #include "Teuchos_LAPACK.hpp"
00049
00050 LOCA::BorderedSolver::LowerTriangularBlockElimination::
00051 LowerTriangularBlockElimination(
00052 const Teuchos::RCP<LOCA::GlobalData>& global_data) :
00053 globalData(global_data)
00054 {
00055 }
00056
00057 LOCA::BorderedSolver::LowerTriangularBlockElimination::
00058 ~LowerTriangularBlockElimination()
00059 {
00060 }
00061
00062
00063 NOX::Abstract::Group::ReturnType
00064 LOCA::BorderedSolver::LowerTriangularBlockElimination::
00065 solve(Teuchos::ParameterList& params,
00066 const LOCA::BorderedSolver::AbstractOperator& op,
00067 const LOCA::MultiContinuation::ConstraintInterface& B,
00068 const NOX::Abstract::MultiVector::DenseMatrix& C,
00069 const NOX::Abstract::MultiVector* F,
00070 const NOX::Abstract::MultiVector::DenseMatrix* G,
00071 NOX::Abstract::MultiVector& X,
00072 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00073 {
00074 string callingFunction =
00075 "LOCA::BorderedSolver::LowerTriangularBlockElimination::solve()";
00076 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00077 NOX::Abstract::Group::ReturnType status;
00078
00079
00080 bool isZeroF = (F == NULL);
00081 bool isZeroG = (G == NULL);
00082 bool isZeroB = B.isDXZero();
00083 bool isZeroX = isZeroF;
00084 bool isZeroY = isZeroG && (isZeroB || isZeroX);
00085
00086
00087 if (isZeroX)
00088 X.init(0.0);
00089 else {
00090
00091 status = op.applyInverse(params, *F, X);
00092 finalStatus =
00093 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00094 finalStatus,
00095 callingFunction);
00096 }
00097
00098
00099 if (isZeroY)
00100 Y.putScalar(0.0);
00101 else {
00102
00103 if (isZeroG)
00104 B.multiplyDX(-1.0, X, Y);
00105 else {
00106 Y.assign(*G);
00107 if (!isZeroB && !isZeroX) {
00108 NOX::Abstract::MultiVector::DenseMatrix T(Y.numRows(),Y.numCols());
00109 B.multiplyDX(1.0, X, T);
00110 Y -= T;
00111 }
00112 }
00113
00114
00115 NOX::Abstract::MultiVector::DenseMatrix M(C);
00116 int *ipiv = new int[M.numRows()];
00117 Teuchos::LAPACK<int,double> L;
00118 int info;
00119 L.GETRF(M.numRows(), M.numCols(), M.values(), M.stride(), ipiv, &info);
00120 if (info != 0) {
00121 status = NOX::Abstract::Group::Failed;
00122 finalStatus =
00123 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00124 status,
00125 finalStatus,
00126 callingFunction);
00127 }
00128 L.GETRS('N', M.numRows(), Y.numCols(), M.values(), M.stride(), ipiv,
00129 Y.values(), Y.stride(), &info);
00130 delete [] ipiv;
00131 if (info != 0) {
00132 status = NOX::Abstract::Group::Failed;
00133 finalStatus =
00134 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00135 status,
00136 finalStatus,
00137 callingFunction);
00138 }
00139 }
00140
00141 return finalStatus;
00142 }
00143
00144 NOX::Abstract::Group::ReturnType
00145 LOCA::BorderedSolver::LowerTriangularBlockElimination::
00146 solve(Teuchos::ParameterList& params,
00147 const LOCA::BorderedSolver::AbstractOperator& op,
00148 const NOX::Abstract::MultiVector& B,
00149 const NOX::Abstract::MultiVector::DenseMatrix& C,
00150 const NOX::Abstract::MultiVector* F,
00151 const NOX::Abstract::MultiVector::DenseMatrix* G,
00152 NOX::Abstract::MultiVector& X,
00153 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00154 {
00155
00156 LOCA::MultiContinuation::MultiVecConstraint cB(Teuchos::rcp(&B,false));
00157
00158
00159 return solve(params, op, cB, C, F, G, X, Y);
00160 }
00161
00162 NOX::Abstract::Group::ReturnType
00163 LOCA::BorderedSolver::LowerTriangularBlockElimination::
00164 solveTranspose(Teuchos::ParameterList& params,
00165 const LOCA::BorderedSolver::AbstractOperator& op,
00166 const LOCA::MultiContinuation::ConstraintInterface& B,
00167 const NOX::Abstract::MultiVector::DenseMatrix& C,
00168 const NOX::Abstract::MultiVector* F,
00169 const NOX::Abstract::MultiVector::DenseMatrix* G,
00170 NOX::Abstract::MultiVector& X,
00171 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00172 {
00173 string callingFunction =
00174 "LOCA::BorderedSolver::LowerTriangularBlockElimination::solveTranspose()";
00175 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00176 NOX::Abstract::Group::ReturnType status;
00177
00178
00179 bool isZeroF = (F == NULL);
00180 bool isZeroG = (G == NULL);
00181 bool isZeroB = B.isDXZero();
00182 bool isZeroX = isZeroF;
00183 bool isZeroY = isZeroG && (isZeroB || isZeroX);
00184
00185
00186 if (isZeroX)
00187 X.init(0.0);
00188 else {
00189
00190 status = op.applyInverseTranspose(params, *F, X);
00191 finalStatus =
00192 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00193 finalStatus,
00194 callingFunction);
00195 }
00196
00197
00198 if (isZeroY)
00199 Y.putScalar(0.0);
00200 else {
00201
00202 if (isZeroG)
00203 B.multiplyDX(-1.0, X, Y);
00204 else {
00205 Y.assign(*G);
00206 if (!isZeroB && !isZeroX) {
00207 NOX::Abstract::MultiVector::DenseMatrix T(Y.numRows(),Y.numCols());
00208 B.multiplyDX(1.0, X, T);
00209 Y -= T;
00210 }
00211 }
00212
00213
00214 NOX::Abstract::MultiVector::DenseMatrix M(C);
00215 int *ipiv = new int[M.numRows()];
00216 Teuchos::LAPACK<int,double> L;
00217 int info;
00218 L.GETRF(M.numRows(), M.numCols(), M.values(), M.stride(), ipiv, &info);
00219 if (info != 0) {
00220 status = NOX::Abstract::Group::Failed;
00221 finalStatus =
00222 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00223 status,
00224 finalStatus,
00225 callingFunction);
00226 }
00227 L.GETRS('T', M.numRows(), Y.numCols(), M.values(), M.stride(), ipiv,
00228 Y.values(), Y.stride(), &info);
00229 delete [] ipiv;
00230 if (info != 0) {
00231 status = NOX::Abstract::Group::Failed;
00232 finalStatus =
00233 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00234 status,
00235 finalStatus,
00236 callingFunction);
00237 }
00238 }
00239
00240 return finalStatus;
00241 }
00242
00243 NOX::Abstract::Group::ReturnType
00244 LOCA::BorderedSolver::LowerTriangularBlockElimination::
00245 solveTranspose(Teuchos::ParameterList& params,
00246 const LOCA::BorderedSolver::AbstractOperator& op,
00247 const NOX::Abstract::MultiVector& B,
00248 const NOX::Abstract::MultiVector::DenseMatrix& C,
00249 const NOX::Abstract::MultiVector* F,
00250 const NOX::Abstract::MultiVector::DenseMatrix* G,
00251 NOX::Abstract::MultiVector& X,
00252 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00253 {
00254
00255 LOCA::MultiContinuation::MultiVecConstraint cB(Teuchos::rcp(&B,false));
00256
00257
00258 return solveTranspose(params, op, cB, C, F, G, X, Y);
00259 }