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 #include "Epetra_config.h"
00040 #include "Epetra_Map.h"
00041 #include "Epetra_Vector.h"
00042 #include "Epetra_Comm.h"
00043 #include "Epetra_Import.h"
00044
00045 #include "LOCA_Epetra_AugmentedOp.H"
00046 #include "LOCA_GlobalData.H"
00047 #include "LOCA_ErrorCheck.H"
00048
00049 LOCA::Epetra::AugmentedOp::AugmentedOp(
00050 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00051 const Teuchos::RCP<Epetra_Operator>& jac,
00052 const Teuchos::RCP<const Epetra_MultiVector>& a_,
00053 const Teuchos::RCP<const Epetra_MultiVector>& b_,
00054 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> c_)
00055 : globalData(global_data),
00056 label("LOCA::Epetra::AugmentedOp"),
00057 jacOperator(jac),
00058 underlyingMap(jacOperator->OperatorDomainMap()),
00059 underlyingComm(underlyingMap.Comm()),
00060 localMap(a_->NumVectors(), 0, jacOperator->Comm()),
00061 a(a_),
00062 b(b_),
00063 c(View, localMap, c_->values(), c_->stride(), c_->numCols()),
00064 underlyingLength(a->MyLength()),
00065 numConstraints(a->NumVectors()),
00066 extendedMapPtr(NULL),
00067 extendedImportMapPtr(NULL),
00068 extendedImporter(NULL),
00069 importedInput(NULL),
00070 result_y(NULL),
00071 tmp(NULL),
00072 haveParamComponent(false),
00073 useTranspose(false),
00074 dlapack()
00075 {
00076
00077
00078 haveParamComponent = (underlyingComm.MyPID() == 0);
00079
00080
00081 buildExtendedMap(underlyingMap, extendedMapPtr, false,
00082 haveParamComponent);
00083
00084
00085 buildExtendedMap(underlyingMap, extendedImportMapPtr, true,
00086 haveParamComponent);
00087
00088
00089 extendedImporter = new Epetra_Import(*extendedImportMapPtr,
00090 *extendedMapPtr);
00091 }
00092
00093 LOCA::Epetra::AugmentedOp::~AugmentedOp()
00094 {
00095 delete extendedMapPtr;
00096 delete extendedImportMapPtr;
00097 delete extendedImporter;
00098 delete importedInput;
00099 delete result_y;
00100 delete tmp;
00101 }
00102
00103 int
00104 LOCA::Epetra::AugmentedOp::SetUseTranspose(bool UseTranspose)
00105 {
00106 useTranspose = UseTranspose;
00107 jacOperator->SetUseTranspose(UseTranspose);
00108 return 0;
00109 }
00110
00111 int
00112 LOCA::Epetra::AugmentedOp::Apply(const Epetra_MultiVector& Input,
00113 Epetra_MultiVector& Result) const
00114 {
00115
00116
00117 int n = Input.NumVectors();
00118
00119
00120 if (importedInput == NULL || n != importedInput->NumVectors())
00121 globalData->locaErrorCheck->throwError("LOCA::Epetra::AugmentedOp::Apply()"
00122 "Must call init() before Apply()!");
00123
00124
00125 importedInput->Import(Input, *extendedImporter, Insert);
00126
00127
00128 double **input_view;
00129 double **result_view;
00130 importedInput->ExtractView(&input_view);
00131 Result.ExtractView(&result_view);
00132
00133
00134
00135
00136 double **input_view_y = new double*[n];
00137 for (int i=0; i<n; i++)
00138 input_view_y[i] = input_view[i]+underlyingLength;
00139
00140
00141 Epetra_MultiVector input_x(View, underlyingMap, input_view, n);
00142 Epetra_MultiVector input_y(View, localMap, input_view_y, n);
00143 Epetra_MultiVector result_x(View, underlyingMap, result_view, n);
00144
00145
00146 jacOperator->Apply(input_x, result_x);
00147
00148 if (useTranspose) {
00149
00150
00151 result_x.Multiply('N', 'N', 1.0, *b, input_y, 1.0);
00152
00153
00154 result_y->Multiply('T', 'N', 1.0, *a, input_x, 0.0);
00155 result_y->Multiply('T', 'N', 1.0, c, input_y, 1.0);
00156
00157 }
00158 else {
00159
00160
00161 result_x.Multiply('N', 'N', 1.0, *a, input_y, 1.0);
00162
00163
00164 result_y->Multiply('T', 'N', 1.0, *b, input_x, 0.0);
00165 result_y->Multiply('N', 'N', 1.0, c, input_y, 1.0);
00166
00167 }
00168
00169
00170 if (haveParamComponent)
00171 for (int j=0; j<n; j++)
00172 for (int i=0; i<numConstraints; i++)
00173 result_view[j][underlyingLength+i] = (*result_y)[j][i];
00174
00175 delete [] input_view_y;
00176
00177 return 0;
00178 }
00179
00180 int
00181 LOCA::Epetra::AugmentedOp::ApplyInverse(const Epetra_MultiVector& Input,
00182 Epetra_MultiVector& Result) const
00183 {
00184
00185
00186 int n = Input.NumVectors();
00187
00188
00189 if (importedInput == NULL || n != importedInput->NumVectors())
00190 globalData->locaErrorCheck->throwError(
00191 "LOCA::Epetra::AugmentedOp::ApplyInverse()"
00192 "Must call init() before ApplyInverse()!");
00193
00194
00195 importedInput->Import(Input, *extendedImporter, Insert);
00196
00197
00198 double **input_view;
00199 double **result_view;
00200 importedInput->ExtractView(&input_view);
00201 Result.ExtractView(&result_view);
00202
00203
00204
00205
00206 double **input_view_y = new double*[n];
00207 for (int i=0; i<n; i++)
00208 input_view_y[i] = input_view[i]+underlyingLength;
00209
00210
00211 Epetra_MultiVector input_x(View, underlyingMap, input_view, n);
00212 Epetra_MultiVector input_y(View, localMap, input_view_y, n);
00213 Epetra_MultiVector result_x(View, underlyingMap, result_view, n);
00214
00215
00216 result_y->Scale(1.0, input_y);
00217
00218
00219 Epetra_MultiVector tmp2(c);
00220 tmp->PutScalar(0.0);
00221
00222
00223 jacOperator->ApplyInverse(input_x, result_x);
00224
00225 if (useTranspose) {
00226
00227
00228 jacOperator->ApplyInverse(*b, *tmp);
00229
00230
00231 result_y->Multiply('T', 'N', -1.0, *a, result_x, 1.0);
00232
00233
00234 tmp2.Multiply('T', 'N', -1.0, *a, *tmp, 1.0);
00235
00236 }
00237 else {
00238
00239
00240 jacOperator->ApplyInverse(*a, *tmp);
00241
00242
00243 result_y->Multiply('T', 'N', -1.0, *b, result_x, 1.0);
00244
00245
00246 tmp2.Multiply('T', 'N', -1.0, *b, *tmp, 1.0);
00247
00248 }
00249
00250
00251 int *ipiv = new int[numConstraints];
00252 int info;
00253 double *result_y_view;
00254 int result_y_lda;
00255 result_y->ExtractView(&result_y_view, &result_y_lda);
00256 double *tmp2_view;
00257 int tmp2_lda;
00258 tmp2.ExtractView(&tmp2_view, &tmp2_lda);
00259 dlapack.GESV(numConstraints, n, tmp2_view, tmp2_lda, ipiv, result_y_view,
00260 result_y_lda, &info);
00261 delete [] ipiv;
00262 if (info != 0) {
00263 globalData->locaErrorCheck->throwError(
00264 "LOCA::Epetra::AugmentedOp::ApplyInverse()"
00265 "Solve of dense matrix failed!");
00266 }
00267
00268
00269 result_x.Multiply('N', 'N', -1.0, *tmp, *result_y, 1.0);
00270
00271
00272 if (haveParamComponent)
00273 for (int j=0; j<n; j++)
00274 for (int i=0; i<numConstraints; i++)
00275 result_view[j][underlyingLength+i] = (*result_y)[j][i];
00276
00277 delete [] input_view_y;
00278
00279 return 0;
00280 }
00281
00282 double
00283 LOCA::Epetra::AugmentedOp::NormInf() const
00284 {
00285 double Jn, an, bn;
00286 double *t = new double[numConstraints];
00287
00288 Jn = jacOperator->NormInf();
00289
00290 a->NormInf(t);
00291 an = t[0];
00292 for (int i=1; i<numConstraints; i++)
00293 if (t[i] > an)
00294 an = t[i];
00295
00296 b->NormInf(t);
00297 bn = t[0];
00298 for (int i=1; i<numConstraints; i++)
00299 if (t[i] > bn)
00300 bn = t[i];
00301
00302 delete [] t;
00303
00304 return Jn + an + bn;
00305 }
00306
00307
00308 const char*
00309 LOCA::Epetra::AugmentedOp::Label () const
00310 {
00311 return const_cast<char*>(label.c_str());
00312 }
00313
00314 bool
00315 LOCA::Epetra::AugmentedOp::UseTranspose() const
00316 {
00317 return useTranspose;
00318 }
00319
00320 bool
00321 LOCA::Epetra::AugmentedOp::HasNormInf() const
00322 {
00323 return jacOperator->HasNormInf();
00324 }
00325
00326 const Epetra_Comm&
00327 LOCA::Epetra::AugmentedOp::Comm() const
00328 {
00329 return underlyingComm;
00330 }
00331 const Epetra_Map&
00332 LOCA::Epetra::AugmentedOp::OperatorDomainMap() const
00333 {
00334 return *extendedMapPtr;
00335 }
00336
00337 const Epetra_Map&
00338 LOCA::Epetra::AugmentedOp::OperatorRangeMap() const
00339 {
00340 return *extendedMapPtr;
00341 }
00342
00343 void
00344 LOCA::Epetra::AugmentedOp::init(const Epetra_MultiVector& x)
00345 {
00346 if (importedInput == NULL || importedInput->NumVectors() != x.NumVectors()) {
00347 if (importedInput != NULL) {
00348 delete importedInput;
00349 delete result_y;
00350 delete tmp;
00351 }
00352 importedInput = new Epetra_MultiVector(*extendedImportMapPtr,
00353 x.NumVectors());
00354 result_y = new Epetra_MultiVector(localMap, x.NumVectors());
00355 tmp = new Epetra_MultiVector(underlyingMap, x.NumVectors());
00356 }
00357 }
00358
00359 Teuchos::RCP<Epetra_MultiVector>
00360 LOCA::Epetra::AugmentedOp::buildEpetraAugmentedMultiVec(
00361 const Epetra_MultiVector& x,
00362 const NOX::Abstract::MultiVector::DenseMatrix *y,
00363 bool doCopy) const
00364 {
00365 Teuchos::RCP<Epetra_MultiVector> extVec =
00366 Teuchos::rcp(new Epetra_MultiVector(*extendedMapPtr, x.NumVectors()));
00367
00368 if (doCopy) {
00369 for (int j=0; j<x.NumVectors(); j++) {
00370 for (int i=0; i<underlyingLength; i++)
00371 (*extVec)[j][i] = x[j][i];
00372
00373 if (haveParamComponent && y != NULL)
00374 for (int i=0; i<numConstraints; i++)
00375 (*extVec)[j][underlyingLength+i] = (*y)(i,j);
00376 }
00377 }
00378 return extVec;
00379 }
00380
00381 void
00382 LOCA::Epetra::AugmentedOp::setEpetraAugmentedMultiVec(
00383 Epetra_MultiVector& x,
00384 NOX::Abstract::MultiVector::DenseMatrix& y,
00385 const Epetra_MultiVector& augMultiVec) const
00386 {
00387
00388 if (importedInput == NULL || x.NumVectors() != importedInput->NumVectors())
00389 globalData->locaErrorCheck->throwError(
00390 "LOCA::Epetra::AugmentedOp::setEpetraAugmentedMultiVec()"
00391 "Must call init() before setEpetraAugmentedMultiVec()!");
00392
00393
00394 importedInput->Import(augMultiVec, *extendedImporter, Insert);
00395
00396 for (int j=0; j<x.NumVectors(); j++) {
00397 for (int i=0; i<underlyingLength; i++)
00398 x[j][i] = (*importedInput)[j][i];
00399 for (int i=0; i<numConstraints; i++)
00400 y(i,j) = (*importedInput)[j][underlyingLength+i];
00401 }
00402 }
00403
00404 void
00405 LOCA::Epetra::AugmentedOp::buildExtendedMap(const Epetra_BlockMap& uMap,
00406 Epetra_Map*& eMapPtr,
00407 bool buildImporter,
00408 bool haveParam)
00409 {
00410 Epetra_BlockMap& nonconstUnderlyingMap = const_cast<Epetra_BlockMap&>(uMap);
00411
00412
00413 Epetra_Map* uPointMapPtr =
00414 dynamic_cast<Epetra_Map*>(&nonconstUnderlyingMap);
00415 bool allocatedPointMap = false;
00416 if (uPointMapPtr == NULL) {
00417 allocatedPointMap = true;
00418 blockMap2PointMap(uMap, uPointMapPtr);
00419 }
00420
00421 int max_gid = uPointMapPtr->MaxAllGID();
00422 int num_global_elements = uPointMapPtr->NumGlobalElements();
00423 int num_my_elements = uPointMapPtr->NumMyElements();
00424 int *global_elements = uPointMapPtr->MyGlobalElements();
00425 const Epetra_Comm& comm = uPointMapPtr->Comm();
00426 int index_base = uPointMapPtr->IndexBase();
00427
00428 int ext_num_global_elements;
00429 int ext_num_my_elements;
00430 int *ext_global_elements;
00431
00432
00433 if (buildImporter)
00434 ext_num_global_elements =
00435 num_global_elements + numConstraints*comm.NumProc();
00436 else
00437 ext_num_global_elements = num_global_elements + numConstraints;
00438
00439
00440 if (buildImporter || haveParam)
00441 ext_num_my_elements = num_my_elements + numConstraints;
00442 else
00443 ext_num_my_elements = num_my_elements;
00444
00445
00446 ext_global_elements = new int[ext_num_my_elements];
00447
00448
00449 for (int i=0; i<num_my_elements; i++) {
00450 ext_global_elements[i] = global_elements[i];
00451 }
00452 if (buildImporter || haveParam)
00453 for (int i=0; i<numConstraints; i++)
00454 ext_global_elements[num_my_elements+i] = max_gid + 1 + i;
00455
00456
00457 eMapPtr = new Epetra_Map(ext_num_global_elements, ext_num_my_elements,
00458 ext_global_elements, index_base, comm);
00459
00460
00461 delete [] ext_global_elements;
00462 if (allocatedPointMap)
00463 delete uPointMapPtr;
00464 }
00465
00466 int
00467 LOCA::Epetra::AugmentedOp::blockMap2PointMap(const Epetra_BlockMap& BlockMap,
00468 Epetra_Map*& PointMap) const
00469 {
00470
00471
00472
00473
00474
00475 int MaxElementSize = BlockMap.MaxElementSize();
00476 int PtNumMyElements = BlockMap.NumMyPoints();
00477 int * PtMyGlobalElements = 0;
00478 if (PtNumMyElements>0) PtMyGlobalElements = new int[PtNumMyElements];
00479
00480 int NumMyElements = BlockMap.NumMyElements();
00481
00482 int curID = 0;
00483 for (int i=0; i<NumMyElements; i++) {
00484 int StartID = BlockMap.GID(i)*MaxElementSize;
00485 int ElementSize = BlockMap.ElementSize(i);
00486 for (int j=0; j<ElementSize; j++) PtMyGlobalElements[curID++] = StartID+j;
00487 }
00488 assert(curID==PtNumMyElements);
00489
00490 PointMap = new Epetra_Map(-1, PtNumMyElements, PtMyGlobalElements, BlockMap.IndexBase(), BlockMap.Comm());
00491
00492 if (PtNumMyElements>0) delete [] PtMyGlobalElements;
00493
00494 if (!BlockMap.PointSameAs(*PointMap)) {EPETRA_CHK_ERR(-1);}
00495 return(0);
00496 }