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_LAPACK_Group.H"
00043 #include "LOCA_GlobalData.H"
00044 #include "LOCA_ErrorCheck.H"
00045 #include "Teuchos_LAPACK.hpp"
00046
00047 LOCA::LAPACK::Group::Group(
00048 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00049 LOCA::LAPACK::Interface& interface) :
00050 NOX::LAPACK::Group(interface),
00051 LOCA::Abstract::Group(global_data),
00052 globalData(global_data),
00053 locaProblemInterface(interface),
00054 params(),
00055 shiftedSolver(jacSolver.getMatrix().numRows()),
00056 freq(0.0),
00057 isValidComplex(false),
00058 complexSolver(jacSolver.getMatrix().numRows())
00059 {
00060 }
00061
00062 LOCA::LAPACK::Group::Group(const LOCA::LAPACK::Group& source,
00063 NOX::CopyType type) :
00064 NOX::LAPACK::Group(source,type),
00065 LOCA::Abstract::Group(source,type),
00066 globalData(source.globalData),
00067 locaProblemInterface(source.locaProblemInterface),
00068 params(source.params),
00069 shiftedSolver(source.shiftedSolver),
00070 freq(source.freq),
00071 isValidComplex(source.isValidComplex),
00072 complexSolver(source.complexSolver)
00073 {
00074 }
00075
00076 LOCA::LAPACK::Group::~Group()
00077 {}
00078
00079 LOCA::LAPACK::Group&
00080 LOCA::LAPACK::Group::operator=(const LOCA::LAPACK::Group& source) {
00081
00082 NOX::LAPACK::Group::operator=(source);
00083 LOCA::Abstract::Group::copy(source);
00084
00085 globalData = source.globalData;
00086 params = source.params;
00087 shiftedSolver = source.shiftedSolver;
00088 freq = source.freq;
00089 isValidComplex = source.isValidComplex;
00090 complexSolver = source.complexSolver;
00091
00092 return *this;
00093 }
00094
00095 NOX::Abstract::Group&
00096 LOCA::LAPACK::Group::operator=(const NOX::Abstract::Group& source) {
00097 operator=(dynamic_cast<const LOCA::LAPACK::Group&>(source));
00098 return *this;
00099 }
00100
00101 NOX::LAPACK::Group&
00102 LOCA::LAPACK::Group::operator=(const NOX::LAPACK::Group& source) {
00103 operator=(dynamic_cast<const LOCA::LAPACK::Group&>(source));
00104 return *this;
00105 }
00106
00107 Teuchos::RCP<NOX::Abstract::Group>
00108 LOCA::LAPACK::Group::clone(NOX::CopyType type) const {
00109 return Teuchos::rcp(new LOCA::LAPACK::Group(*this, type));
00110 }
00111
00112 NOX::Abstract::Group::ReturnType
00113 LOCA::LAPACK::Group::computeF() {
00114 locaProblemInterface.setParams(params);
00115 return NOX::LAPACK::Group::computeF();
00116 }
00117
00118 NOX::Abstract::Group::ReturnType
00119 LOCA::LAPACK::Group::computeJacobian() {
00120 locaProblemInterface.setParams(params);
00121 return NOX::LAPACK::Group::computeJacobian();
00122 }
00123
00124 NOX::Abstract::Group::ReturnType
00125 LOCA::LAPACK::Group::applyJacobianTransposeInverse(
00126 Teuchos::ParameterList& p,
00127 const NOX::Abstract::Vector& input,
00128 NOX::Abstract::Vector& result) const
00129 {
00130
00131 if (!isJacobian()) {
00132 cerr << "ERROR: "
00133 << "LOCA::LAPACK::Group::applyJacobianTransposeInverse()"
00134 << " - invalid Jacobian" << endl;
00135 throw "NOX Error";
00136 }
00137
00138 const NOX::LAPACK::Vector& lapack_input =
00139 dynamic_cast<const NOX::LAPACK::Vector&> (input);
00140 NOX::LAPACK::Vector& lapack_result =
00141 dynamic_cast<NOX::LAPACK::Vector&> (result);
00142
00143
00144 lapack_result = lapack_input;
00145 bool res = jacSolver.solve(true, 1, &lapack_result(0));
00146
00147 return res ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00148 }
00149
00150 NOX::Abstract::Group::ReturnType
00151 LOCA::LAPACK::Group::applyJacobianTransposeInverseMultiVector(
00152 Teuchos::ParameterList& p,
00153 const NOX::Abstract::MultiVector& input,
00154 NOX::Abstract::MultiVector& result) const
00155 {
00156
00157 if (!isJacobian()) {
00158 cerr << "ERROR: "
00159 << "LOCA::LAPACK::Group::applyJacobianTransposeInverseMultiVector()"
00160 << " - invalid Jacobian" << endl;
00161 throw "NOX Error";
00162 }
00163
00164
00165 int nVecs = input.numVectors();
00166
00167 int m = jacSolver.getMatrix().numRows();
00168
00169
00170 NOX::LAPACK::Matrix<double> B(m,nVecs);
00171 const NOX::LAPACK::Vector* constVecPtr;
00172 for (int j=0; j<nVecs; j++) {
00173 constVecPtr = dynamic_cast<const NOX::LAPACK::Vector*>(&(input[j]));
00174 for (int i=0; i<m; i++)
00175 B(i,j) = (*constVecPtr)(i);
00176 }
00177
00178
00179 bool res = jacSolver.solve(true, nVecs, &B(0,0));
00180
00181 if (!res)
00182 return NOX::Abstract::Group::Failed;
00183
00184
00185 NOX::LAPACK::Vector* vecPtr;
00186 for (int j=0; j<nVecs; j++) {
00187 vecPtr = dynamic_cast<NOX::LAPACK::Vector*>(&(result[j]));
00188 for (int i=0; i<m; i++)
00189 (*vecPtr)(i) = B(i,j);
00190 }
00191
00192 return NOX::Abstract::Group::Ok;
00193 }
00194
00195 void
00196 LOCA::LAPACK::Group::copy(const NOX::Abstract::Group& source) {
00197 *this = source;
00198 }
00199
00200 void
00201 LOCA::LAPACK::Group::setParams(const LOCA::ParameterVector& p)
00202 {
00203 resetIsValid();
00204 params = p;
00205 }
00206
00207 void
00208 LOCA::LAPACK::Group::setParam(int paramID, double val)
00209 {
00210 resetIsValid();
00211 params.setValue(paramID, val);
00212 }
00213
00214 void
00215 LOCA::LAPACK::Group::setParam(string paramID, double val)
00216 {
00217 resetIsValid();
00218 params.setValue(paramID, val);
00219 }
00220
00221 const LOCA::ParameterVector&
00222 LOCA::LAPACK::Group::getParams() const
00223 {
00224 return params;
00225 }
00226
00227 double
00228 LOCA::LAPACK::Group::getParam(int paramID) const
00229 {
00230 return params.getValue(paramID);
00231 }
00232
00233 double
00234 LOCA::LAPACK::Group::getParam(string paramID) const
00235 {
00236 return params.getValue(paramID);
00237 }
00238
00239 void
00240 LOCA::LAPACK::Group::projectToDraw(const NOX::Abstract::Vector& x,
00241 double *px) const
00242 {
00243 const NOX::LAPACK::Vector& lx =
00244 dynamic_cast<const NOX::LAPACK::Vector&>(x);
00245 locaProblemInterface.projectToDraw(lx, px);
00246 }
00247
00248 int
00249 LOCA::LAPACK::Group::projectToDrawDimension() const
00250 {
00251 return locaProblemInterface.projectToDrawDimension();
00252 }
00253
00254 double
00255 LOCA::LAPACK::Group::computeScaledDotProduct(
00256 const NOX::Abstract::Vector& a,
00257 const NOX::Abstract::Vector& b) const
00258 {
00259 return a.innerProduct(b) / a.length();
00260 }
00261
00262 void
00263 LOCA::LAPACK::Group::printSolution(const double conParam) const
00264 {
00265 printSolution(xVector, conParam);
00266 }
00267
00268 void
00269 LOCA::LAPACK::Group::printSolution(const NOX::Abstract::Vector& x_,
00270 const double conParam) const
00271 {
00272 locaProblemInterface.printSolution(dynamic_cast<const NOX::LAPACK::Vector&>(x_), conParam);
00273 }
00274
00275 void
00276 LOCA::LAPACK::Group::scaleVector(NOX::Abstract::Vector& x) const
00277 {
00278 x.scale(1.0 / sqrt(static_cast<double>(x.length())));
00279 }
00280
00281 NOX::Abstract::Group::ReturnType
00282 LOCA::LAPACK::Group::computeShiftedMatrix(double alpha, double beta)
00283 {
00284
00285 bool res =
00286 locaProblemInterface.computeShiftedMatrix(alpha, beta, xVector,
00287 shiftedSolver.getMatrix());
00288
00289 if (res)
00290 return NOX::Abstract::Group::Ok;
00291 else
00292 return NOX::Abstract::Group::Failed;
00293 }
00294
00295 NOX::Abstract::Group::ReturnType
00296 LOCA::LAPACK::Group::applyShiftedMatrix(const NOX::Abstract::Vector& input,
00297 NOX::Abstract::Vector& result) const
00298 {
00299
00300 const NOX::LAPACK::Vector& lapack_input =
00301 dynamic_cast<const NOX::LAPACK::Vector&>(input);
00302 NOX::LAPACK::Vector& lapack_result =
00303 dynamic_cast<NOX::LAPACK::Vector&>(result);
00304
00305
00306 shiftedSolver.apply(false, 1, &lapack_input(0), &lapack_result(0));
00307
00308 return NOX::Abstract::Group::Ok;
00309 }
00310
00311 NOX::Abstract::Group::ReturnType
00312 LOCA::LAPACK::Group::applyShiftedMatrixMultiVector(
00313 const NOX::Abstract::MultiVector& input,
00314 NOX::Abstract::MultiVector& result) const
00315 {
00316
00317 int nVecs = input.numVectors();
00318
00319 int m = shiftedSolver.getMatrix().numRows();
00320
00321
00322 NOX::LAPACK::Matrix<double> B(m,nVecs);
00323 NOX::LAPACK::Matrix<double> C(m,nVecs);
00324 const NOX::LAPACK::Vector* constVecPtr;
00325 for (int j=0; j<nVecs; j++) {
00326 constVecPtr = dynamic_cast<const NOX::LAPACK::Vector*>(&(input[j]));
00327 for (int i=0; i<m; i++)
00328 B(i,j) = (*constVecPtr)(i);
00329 }
00330
00331
00332 shiftedSolver.apply(false, nVecs, &B(0,0), &C(0,0));
00333
00334
00335 NOX::LAPACK::Vector* vecPtr;
00336 for (int j=0; j<nVecs; j++) {
00337 vecPtr = dynamic_cast<NOX::LAPACK::Vector*>(&(result[j]));
00338 for (int i=0; i<m; i++)
00339 (*vecPtr)(i) = C(i,j);
00340 }
00341
00342 return NOX::Abstract::Group::Ok;
00343 }
00344
00345 NOX::Abstract::Group::ReturnType
00346 LOCA::LAPACK::Group::applyShiftedMatrixInverseMultiVector(
00347 Teuchos::ParameterList& params,
00348 const NOX::Abstract::MultiVector& input,
00349 NOX::Abstract::MultiVector& result) const
00350 {
00351
00352 int nVecs = input.numVectors();
00353
00354 int m = shiftedSolver.getMatrix().numRows();
00355
00356
00357 NOX::LAPACK::Matrix<double> B(m,nVecs);
00358 const NOX::LAPACK::Vector* constVecPtr;
00359 for (int j=0; j<nVecs; j++) {
00360 constVecPtr = dynamic_cast<const NOX::LAPACK::Vector*>(&(input[j]));
00361 for (int i=0; i<m; i++)
00362 B(i,j) = (*constVecPtr)(i);
00363 }
00364
00365 bool res = shiftedSolver.solve(false, nVecs, &B(0,0));
00366
00367 if (!res)
00368 return NOX::Abstract::Group::Failed;
00369
00370
00371 NOX::LAPACK::Vector* vecPtr;
00372 for (int j=0; j<nVecs; j++) {
00373 vecPtr = dynamic_cast<NOX::LAPACK::Vector*>(&(result[j]));
00374 for (int i=0; i<m; i++)
00375 (*vecPtr)(i) = B(i,j);
00376 }
00377
00378 return NOX::Abstract::Group::Ok;
00379 }
00380
00381 bool
00382 LOCA::LAPACK::Group::isComplex() const
00383 {
00384 return isValidComplex;
00385 }
00386
00387 NOX::Abstract::Group::ReturnType
00388 LOCA::LAPACK::Group::computeComplex(double frequency)
00389 {
00390 string callingFunction =
00391 "LOCA::LAPACK::computeComplex()";
00392 NOX::Abstract::Group::ReturnType finalStatus;
00393
00394 freq = frequency;
00395
00396
00397 finalStatus = computeJacobian();
00398 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00399
00400
00401 bool res =
00402 locaProblemInterface.computeShiftedMatrix(0.0, 1.0, xVector,
00403 shiftedSolver.getMatrix());
00404
00405
00406 NOX::LAPACK::Matrix<double>& jacobianMatrix = jacSolver.getMatrix();
00407 NOX::LAPACK::Matrix<double>& massMatrix = shiftedSolver.getMatrix();
00408 NOX::LAPACK::Matrix< std::complex<double> >& complexMatrix =
00409 complexSolver.getMatrix();
00410 int n = jacobianMatrix.numRows();
00411 for (int j=0; j<n; j++) {
00412 for (int i=0; i<n; i++) {
00413 complexMatrix(i,j) =
00414 std::complex<double>(jacobianMatrix(i,j), frequency*massMatrix(i,j));
00415 }
00416 }
00417
00418 if (finalStatus == NOX::Abstract::Group::Ok && res)
00419 isValidComplex = true;
00420
00421 if (res)
00422 return finalStatus;
00423 else
00424 return NOX::Abstract::Group::Failed;
00425 }
00426
00427 NOX::Abstract::Group::ReturnType
00428 LOCA::LAPACK::Group::applyComplex(const NOX::Abstract::Vector& input_real,
00429 const NOX::Abstract::Vector& input_imag,
00430 NOX::Abstract::Vector& result_real,
00431 NOX::Abstract::Vector& result_imag) const
00432 {
00433
00434 if (!isComplex())
00435 return NOX::Abstract::Group::BadDependency;
00436
00437 int n = complexSolver.getMatrix().numCols();
00438
00439
00440 std::vector< std::complex<double> > input(n);
00441 std::vector< std::complex<double> > result(n);
00442 const NOX::LAPACK::Vector& lapack_input_real =
00443 dynamic_cast<const NOX::LAPACK::Vector&>(input_real);
00444 const NOX::LAPACK::Vector& lapack_input_imag =
00445 dynamic_cast<const NOX::LAPACK::Vector&>(input_imag);
00446 for (int i=0; i<n; i++)
00447 input[i] = std::complex<double>(lapack_input_real(i),
00448 lapack_input_imag(i));
00449
00450
00451 complexSolver.apply(false, 1, &input[0], &result[0]);
00452
00453
00454 NOX::LAPACK::Vector& lapack_result_real =
00455 dynamic_cast<NOX::LAPACK::Vector&>(result_real);
00456 NOX::LAPACK::Vector& lapack_result_imag =
00457 dynamic_cast<NOX::LAPACK::Vector&>(result_imag);
00458 for (int i=0; i<n; i++) {
00459 lapack_result_real(i) = result[i].real();
00460 lapack_result_imag(i) = result[i].imag();
00461 }
00462
00463 return NOX::Abstract::Group::Ok;
00464 }
00465
00466 NOX::Abstract::Group::ReturnType
00467 LOCA::LAPACK::Group::applyComplexMultiVector(
00468 const NOX::Abstract::MultiVector& input_real,
00469 const NOX::Abstract::MultiVector& input_imag,
00470 NOX::Abstract::MultiVector& result_real,
00471 NOX::Abstract::MultiVector& result_imag) const
00472 {
00473
00474 if (!isComplex())
00475 return NOX::Abstract::Group::BadDependency;
00476
00477 int n = complexSolver.getMatrix().numRows();
00478 int p = input_real.numVectors();
00479
00480
00481 std::vector< std::complex<double> > input(n*p);
00482 std::vector< std::complex<double> > result(n*p);
00483 const NOX::LAPACK::Vector* lapack_input_real;
00484 const NOX::LAPACK::Vector* lapack_input_imag;
00485 for (int j=0; j<p; j++) {
00486 lapack_input_real =
00487 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_real[j]));
00488 lapack_input_imag =
00489 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_imag[j]));
00490 for (int i=0; i<n; i++)
00491 input[i+n*j] = std::complex<double>((*lapack_input_real)(i),
00492 (*lapack_input_imag)(i));
00493 }
00494
00495
00496 complexSolver.apply(false, p, &input[0], &result[0]);
00497
00498
00499 NOX::LAPACK::Vector* lapack_result_real;
00500 NOX::LAPACK::Vector* lapack_result_imag;
00501 for (int j=0; j<p; j++) {
00502 lapack_result_real =
00503 dynamic_cast<NOX::LAPACK::Vector*>(&(result_real[j]));
00504 lapack_result_imag =
00505 dynamic_cast<NOX::LAPACK::Vector*>(&(result_imag[j]));
00506 for (int i=0; i<n; i++) {
00507 (*lapack_result_real)(i) = result[i+n*j].real();
00508 (*lapack_result_imag)(i) = result[i+n*j].imag();
00509 }
00510 }
00511
00512 return NOX::Abstract::Group::Ok;
00513 }
00514
00515 NOX::Abstract::Group::ReturnType
00516 LOCA::LAPACK::Group::applyComplexInverseMultiVector(
00517 Teuchos::ParameterList& params,
00518 const NOX::Abstract::MultiVector& input_real,
00519 const NOX::Abstract::MultiVector& input_imag,
00520 NOX::Abstract::MultiVector& result_real,
00521 NOX::Abstract::MultiVector& result_imag) const
00522 {
00523
00524 if (!isComplex())
00525 return NOX::Abstract::Group::BadDependency;
00526
00527 int n = complexSolver.getMatrix().numRows();
00528 int p = input_real.numVectors();
00529
00530
00531 std::vector< std::complex<double> > input(n*p);
00532 const NOX::LAPACK::Vector* lapack_input_real;
00533 const NOX::LAPACK::Vector* lapack_input_imag;
00534 for (int j=0; j<p; j++) {
00535 lapack_input_real =
00536 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_real[j]));
00537 lapack_input_imag =
00538 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_imag[j]));
00539 for (int i=0; i<n; i++)
00540 input[i+n*j] = std::complex<double>((*lapack_input_real)(i),
00541 (*lapack_input_imag)(i));
00542 }
00543
00544
00545 bool res = complexSolver.solve(false, p, &input[0]);
00546
00547
00548 NOX::LAPACK::Vector* lapack_result_real;
00549 NOX::LAPACK::Vector* lapack_result_imag;
00550 for (int j=0; j<p; j++) {
00551 lapack_result_real =
00552 dynamic_cast<NOX::LAPACK::Vector*>(&(result_real[j]));
00553 lapack_result_imag =
00554 dynamic_cast<NOX::LAPACK::Vector*>(&(result_imag[j]));
00555 for (int i=0; i<n; i++) {
00556 (*lapack_result_real)(i) = input[i+n*j].real();
00557 (*lapack_result_imag)(i) = input[i+n*j].imag();
00558 }
00559 }
00560
00561 if (res)
00562 return NOX::Abstract::Group::Ok;
00563 else
00564 return NOX::Abstract::Group::Failed;
00565 }
00566
00567 NOX::Abstract::Group::ReturnType
00568 LOCA::LAPACK::Group::applyComplexTranspose(
00569 const NOX::Abstract::Vector& input_real,
00570 const NOX::Abstract::Vector& input_imag,
00571 NOX::Abstract::Vector& result_real,
00572 NOX::Abstract::Vector& result_imag) const
00573 {
00574
00575 if (!isComplex())
00576 return NOX::Abstract::Group::BadDependency;
00577
00578 int n = complexSolver.getMatrix().numCols();
00579
00580
00581 std::vector< std::complex<double> > input(n);
00582 std::vector< std::complex<double> > result(n);
00583 const NOX::LAPACK::Vector& lapack_input_real =
00584 dynamic_cast<const NOX::LAPACK::Vector&>(input_real);
00585 const NOX::LAPACK::Vector& lapack_input_imag =
00586 dynamic_cast<const NOX::LAPACK::Vector&>(input_imag);
00587 for (int i=0; i<n; i++)
00588 input[i] = std::complex<double>(lapack_input_real(i),
00589 lapack_input_imag(i));
00590
00591
00592 complexSolver.apply(true, 1, &input[0], &result[0]);
00593
00594
00595 NOX::LAPACK::Vector& lapack_result_real =
00596 dynamic_cast<NOX::LAPACK::Vector&>(result_real);
00597 NOX::LAPACK::Vector& lapack_result_imag =
00598 dynamic_cast<NOX::LAPACK::Vector&>(result_imag);
00599 for (int i=0; i<n; i++) {
00600 lapack_result_real(i) = result[i].real();
00601 lapack_result_imag(i) = result[i].imag();
00602 }
00603
00604 return NOX::Abstract::Group::Ok;
00605 }
00606
00607 NOX::Abstract::Group::ReturnType
00608 LOCA::LAPACK::Group::applyComplexTransposeMultiVector(
00609 const NOX::Abstract::MultiVector& input_real,
00610 const NOX::Abstract::MultiVector& input_imag,
00611 NOX::Abstract::MultiVector& result_real,
00612 NOX::Abstract::MultiVector& result_imag) const
00613 {
00614
00615 if (!isComplex())
00616 return NOX::Abstract::Group::BadDependency;
00617
00618 int n = complexSolver.getMatrix().numRows();
00619 int p = input_real.numVectors();
00620
00621
00622 std::vector< std::complex<double> > input(n*p);
00623 std::vector< std::complex<double> > result(n*p);
00624 const NOX::LAPACK::Vector* lapack_input_real;
00625 const NOX::LAPACK::Vector* lapack_input_imag;
00626 for (int j=0; j<p; j++) {
00627 lapack_input_real =
00628 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_real[j]));
00629 lapack_input_imag =
00630 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_imag[j]));
00631 for (int i=0; i<n; i++)
00632 input[i+n*j] = std::complex<double>((*lapack_input_real)(i),
00633 (*lapack_input_imag)(i));
00634 }
00635
00636
00637 complexSolver.apply(true, p, &input[0], &result[0]);
00638
00639
00640 NOX::LAPACK::Vector* lapack_result_real;
00641 NOX::LAPACK::Vector* lapack_result_imag;
00642 for (int j=0; j<p; j++) {
00643 lapack_result_real =
00644 dynamic_cast<NOX::LAPACK::Vector*>(&(result_real[j]));
00645 lapack_result_imag =
00646 dynamic_cast<NOX::LAPACK::Vector*>(&(result_imag[j]));
00647 for (int i=0; i<n; i++) {
00648 (*lapack_result_real)(i) = result[i+n*j].real();
00649 (*lapack_result_imag)(i) = result[i+n*j].imag();
00650 }
00651 }
00652
00653 return NOX::Abstract::Group::Ok;
00654 }
00655
00656 NOX::Abstract::Group::ReturnType
00657 LOCA::LAPACK::Group::applyComplexTransposeInverseMultiVector(
00658 Teuchos::ParameterList& params,
00659 const NOX::Abstract::MultiVector& input_real,
00660 const NOX::Abstract::MultiVector& input_imag,
00661 NOX::Abstract::MultiVector& result_real,
00662 NOX::Abstract::MultiVector& result_imag) const
00663 {
00664
00665 if (!isComplex())
00666 return NOX::Abstract::Group::BadDependency;
00667
00668 int n = complexSolver.getMatrix().numRows();
00669 int p = input_real.numVectors();
00670
00671
00672 std::vector< std::complex<double> > input(n*p);
00673 const NOX::LAPACK::Vector* lapack_input_real;
00674 const NOX::LAPACK::Vector* lapack_input_imag;
00675 for (int j=0; j<p; j++) {
00676 lapack_input_real =
00677 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_real[j]));
00678 lapack_input_imag =
00679 dynamic_cast<const NOX::LAPACK::Vector*>(&(input_imag[j]));
00680 for (int i=0; i<n; i++)
00681 input[i+n*j] = std::complex<double>((*lapack_input_real)(i),
00682 (*lapack_input_imag)(i));
00683 }
00684
00685
00686 bool res = complexSolver.solve(true, p, &input[0]);
00687
00688
00689 NOX::LAPACK::Vector* lapack_result_real;
00690 NOX::LAPACK::Vector* lapack_result_imag;
00691 for (int j=0; j<p; j++) {
00692 lapack_result_real =
00693 dynamic_cast<NOX::LAPACK::Vector*>(&(result_real[j]));
00694 lapack_result_imag =
00695 dynamic_cast<NOX::LAPACK::Vector*>(&(result_imag[j]));
00696 for (int i=0; i<n; i++) {
00697 (*lapack_result_real)(i) = input[i+n*j].real();
00698 (*lapack_result_imag)(i) = input[i+n*j].imag();
00699 }
00700 }
00701
00702 if (res)
00703 return NOX::Abstract::Group::Ok;
00704 else
00705 return NOX::Abstract::Group::Failed;
00706 }
00707
00708 NOX::Abstract::Group::ReturnType
00709 LOCA::LAPACK::Group::augmentJacobianForHomotopy(double a, double b)
00710 {
00711 NOX::LAPACK::Matrix<double>& jacobianMatrix = jacSolver.getMatrix();
00712 int size = jacobianMatrix.numRows();
00713
00714
00715 jacobianMatrix.scale(a);
00716
00717
00718 for (int i = 0; i < size; i++)
00719 jacobianMatrix(i,i) += b;
00720
00721 return NOX::Abstract::Group::Ok;
00722 }
00723
00724 void
00725 LOCA::LAPACK::Group::resetIsValid()
00726 {
00727 NOX::LAPACK::Group::resetIsValid();
00728 shiftedSolver.reset();
00729 isValidComplex = false;
00730 complexSolver.reset();
00731 }