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_Comm.h"
00040 #include "Epetra_Map.h"
00041 #include "Epetra_MapColoring.h"
00042 #include "Epetra_Import.h"
00043 #include "Epetra_Vector.h"
00044 #include "Epetra_IntVector.h"
00045 #include "Epetra_CrsGraph.h"
00046 #include "Epetra_CrsMatrix.h"
00047
00048 #include "Epetra_Time.h"
00049
00050 #include "NOX_Epetra_Interface_Required.H"
00051 #include "NOX_Utils.H"
00052
00053 #include "NOX_Epetra_FiniteDifferenceColoring.H"
00054
00055 using namespace NOX;
00056 using namespace NOX::Epetra;
00057
00058
00059
00060 FiniteDifferenceColoring::FiniteDifferenceColoring(
00061 Teuchos::ParameterList& printingParams,
00062 const Teuchos::RCP<Interface::Required>& i,
00063 const NOX::Epetra::Vector& x,
00064 const Teuchos::RCP<Epetra_MapColoring>& colorMap_,
00065 const Teuchos::RCP< vector<Epetra_IntVector> >& columns_,
00066 bool parallelColoring,
00067 bool distance1_,
00068 double beta_, double alpha_) :
00069 FiniteDifference(printingParams, i, x, beta_, alpha_),
00070 coloringType(NOX_SERIAL),
00071 distance1(distance1_),
00072 colorMap(colorMap_),
00073 columns(columns_),
00074 numColors(colorMap->NumColors()),
00075 maxNumColors(colorMap->MaxNumColors()),
00076 colorList(colorMap->ListOfColors()),
00077 cMap(0),
00078 Importer(0),
00079 colorVect(0),
00080 betaColorVect(0),
00081 mappedColorVect(0),
00082 xCol_perturb(0),
00083 columnMap(0),
00084 rowColImporter(0)
00085 {
00086 label = "NOX::FiniteDifferenceColoring Jacobian";
00087
00088 if( parallelColoring )
00089 coloringType = NOX_PARALLEL;
00090
00091 createColorContainers();
00092 }
00093
00094 FiniteDifferenceColoring::FiniteDifferenceColoring(
00095 Teuchos::ParameterList& printingParams,
00096 const Teuchos::RCP<Interface::Required>& i,
00097 const NOX::Epetra::Vector& x,
00098 const Teuchos::RCP<Epetra_CrsGraph>& rawGraph_,
00099 const Teuchos::RCP<Epetra_MapColoring>& colorMap_,
00100 const Teuchos::RCP< vector<Epetra_IntVector> >& columns_,
00101 bool parallelColoring,
00102 bool distance1_,
00103 double beta_, double alpha_) :
00104 FiniteDifference(printingParams, i, x, rawGraph_, beta_, alpha_),
00105 coloringType(NOX_SERIAL),
00106 distance1(distance1_),
00107 colorMap(colorMap_),
00108 columns(columns_),
00109 numColors(colorMap->NumColors()),
00110 maxNumColors(colorMap->MaxNumColors()),
00111 colorList(colorMap->ListOfColors()),
00112 cMap(0),
00113 Importer(0),
00114 colorVect(0),
00115 betaColorVect(0),
00116 mappedColorVect(new Epetra_Vector(rawGraph_->ColMap())),
00117 xCol_perturb(new Epetra_Vector(rawGraph_->ColMap())),
00118 columnMap(&(rawGraph_->ColMap())),
00119 rowColImporter(new Epetra_Import(*columnMap, fo.Map()))
00120 {
00121 label = "NOX::FiniteDifferenceColoring Jacobian";
00122
00123 if( parallelColoring )
00124 coloringType = NOX_PARALLEL;
00125
00126 createColorContainers();
00127 }
00128
00129 FiniteDifferenceColoring::~FiniteDifferenceColoring()
00130 {
00131 delete cMap; cMap = 0;
00132 delete Importer; Importer = 0;
00133 delete colorVect; colorVect = 0;
00134 delete betaColorVect; betaColorVect = 0;
00135 delete rowColImporter; rowColImporter = 0;
00136 delete xCol_perturb; xCol_perturb = 0;
00137 delete mappedColorVect; mappedColorVect = 0;
00138 }
00139
00140 bool FiniteDifferenceColoring::computeJacobian(const Epetra_Vector& x)
00141 {
00142 return( computeJacobian(x, *this));
00143 }
00144
00145 bool FiniteDifferenceColoring::computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
00146 {
00147
00148
00149 FiniteDifferenceColoring* testMatrix =
00150 dynamic_cast<FiniteDifferenceColoring*>(&Jac);
00151 if (testMatrix == 0) {
00152 cout << "ERROR: NOX::Epetra::FiniteDifferenceColoring::computeJacobian() - "
00153 << "Jacobian to evaluate is not a FiniteDifferenceColoring object!"
00154 << endl;
00155 throw "NOX Error";
00156 }
00157
00158 const Epetra_BlockMap& map = fo.Map();
00159
00160
00161 Epetra_Time fillTimer(x.Comm());
00162
00163
00164
00165 Epetra_CrsMatrix& jac = *testMatrix->jacobian;
00166
00167
00168 jac.PutScalar(0.0);
00169
00170
00171 if ( diffType == Centered )
00172 if ( Teuchos::is_null(fmPtr) )
00173 fmPtr = Teuchos::rcp(new Epetra_Vector(x));
00174
00175 double scaleFactor = 1.0;
00176 if ( diffType == Backward )
00177 scaleFactor = -1.0;
00178
00179 int myMin = map.MinMyGID();
00180 int myMax = map.MaxMyGID();
00181
00182
00183
00184
00185 xCol_perturb->Import(x, *rowColImporter, Insert);
00186
00187
00188 if( coloringType == NOX_SERIAL )
00189 computeF(*xCol_perturb, fo, NOX::Epetra::Interface::Required::FD_Res);
00190 else {
00191
00192 for (int i=0; i<x_perturb.MyLength(); i++)
00193 x_perturb[i] = (*xCol_perturb)[columnMap->LID(map.GID(i))];
00194 computeF(x_perturb, fo, NOX::Epetra::Interface::Required::FD_Res);
00195 }
00196
00197
00198 list<int>::iterator allBegin = listOfAllColors.begin(),
00199 allEnd = listOfAllColors.end(),
00200 allIter;
00201
00202 std::map<int, int>::iterator mapEnd = colorToNumMap.end(),
00203 myMapIter;
00204 for ( allIter = allBegin; allIter != allEnd; ++allIter ) {
00205 bool skipIt = true;
00206 int k = -1;
00207 myMapIter = colorToNumMap.find( *allIter );
00208 if( myMapIter != mapEnd ) {
00209 skipIt = false;
00210 k = (*myMapIter).second;
00211 }
00212
00213
00214
00215
00216 int color = -1;
00217 if( !skipIt ) color = colorList[k];
00218 cMap = colorMap->GenerateMap(color);
00219 colorVect = new Epetra_Vector(*cMap);
00220 betaColorVect = new Epetra_Vector(*cMap);
00221 betaColorVect->PutScalar(beta);
00222
00223
00224
00225
00226
00227
00228 for (int i=0; i<colorVect->MyLength(); i++)
00229 (*colorVect)[i] = (*xCol_perturb)[columnMap->LID(cMap->GID(i))];
00230 colorVect->Abs(*colorVect);
00231 colorVect->Update(1.0, *betaColorVect, alpha);
00232
00233
00234
00235 mappedColorVect->PutScalar(0.0);
00236
00237
00238 for (int i=0; i<colorVect->MyLength(); i++)
00239 (*mappedColorVect)[columnMap->LID(cMap->GID(i))] = (*colorVect)[i];
00240 xCol_perturb->Update(scaleFactor, *mappedColorVect, 1.0);
00241
00242
00243 if( coloringType == NOX_SERIAL )
00244 computeF(*xCol_perturb, fp, NOX::Epetra::Interface::Required::FD_Res);
00245 else {
00246
00247 for (int i=0; i<x_perturb.MyLength(); i++)
00248 x_perturb[i] = (*xCol_perturb)[columnMap->LID(map.GID(i))];
00249 computeF(x_perturb, fp, NOX::Epetra::Interface::Required::FD_Res);
00250 }
00251
00252 if ( diffType == Centered ) {
00253 xCol_perturb->Update(-2.0, *mappedColorVect, 1.0);
00254 if( coloringType == NOX_SERIAL )
00255 computeF(*xCol_perturb, *fmPtr,
00256 NOX::Epetra::Interface::Required::FD_Res);
00257 else {
00258
00259 for (int i=0; i<x_perturb.MyLength(); i++)
00260 x_perturb[i] = (*xCol_perturb)[columnMap->LID(map.GID(i))];
00261 computeF(x_perturb, *fmPtr, NOX::Epetra::Interface::Required::FD_Res);
00262 }
00263 }
00264
00265
00266
00267 if ( diffType != Centered ) {
00268 Jc.Update(1.0, fp, -1.0, fo, 0.0);
00269 }
00270 else {
00271 Jc.Update(1.0, fp, -1.0, *fmPtr, 0.0);
00272 }
00273
00274
00275 if( !skipIt ) {
00276 for (int j = myMin; j < myMax+1; j++) {
00277
00278 if (!map.MyGID(j))
00279 continue;
00280 int globalColumnID = (*columns)[k][map.LID(j)];
00281
00282
00283 if( distance1 && (j != globalColumnID) )
00284 continue;
00285
00286 if( globalColumnID >= 0) {
00287
00288
00289
00290 if ( diffType != Centered )
00291 Jc[map.LID(j)] /=
00292 (scaleFactor * (*mappedColorVect)[columnMap->LID(globalColumnID)]);
00293 else
00294 Jc[map.LID(j)] /=
00295 (2.0 * (*mappedColorVect)[columnMap->LID(globalColumnID)]);
00296
00297 int err = jac.ReplaceGlobalValues(j,1,&Jc[map.LID(j)],&globalColumnID);
00298 if(err) {
00299 cout << "ERROR (" << map.Comm().MyPID() << ") : "
00300 << "Inserting global value with indices (" << j << ","
00301 << globalColumnID << ") = " << Jc[map.LID(j)] << endl;
00302 }
00303 }
00304 }
00305 }
00306
00307
00308 delete Importer; Importer = 0;
00309 delete betaColorVect; betaColorVect = 0;
00310 delete colorVect; colorVect = 0;
00311 delete cMap; cMap = 0;
00312
00313
00314 xCol_perturb->Import(x, *rowColImporter, Insert);
00315 }
00316
00317 double fillTime = fillTimer.ElapsedTime();
00318 x.Comm().Barrier();
00319
00320 if (utils.isPrintType(Utils::Details)) {
00321 for(int n = 0; n < map.Comm().NumProc(); n++) {
00322 if(map.Comm().MyPID() == n)
00323 cout << "\tTime to fill Jacobian [" << n << "] --> "
00324 << fillTime << " sec." << endl;
00325 x.Comm().Barrier();
00326 }
00327 cout << endl;
00328 }
00329
00330 jac.FillComplete();
00331
00332 return true;
00333 }
00334
00335 void FiniteDifferenceColoring::createColorContainers()
00336 {
00337
00338 int sumNumColors = numColors;
00339 colorMap->Comm().SumAll(&numColors, &sumNumColors, 1);
00340 int * allColorList = new int[maxNumColors*colorMap->Comm().NumProc()];
00341 int * myColorList = new int[maxNumColors];
00342 for( int i=0; i< maxNumColors; i++)
00343 if( i<numColors )
00344 myColorList[i] = colorList[i];
00345 else
00346 myColorList[i] = -1;
00347
00348 colorMap->Comm().GatherAll( myColorList, allColorList, maxNumColors );
00349
00350
00351 for( int i = 0; i < maxNumColors*colorMap->Comm().NumProc(); i++ )
00352 listOfAllColors.push_back( allColorList[i] );
00353
00354 listOfAllColors.remove(-1);
00355 listOfAllColors.sort();
00356 listOfAllColors.unique();
00357
00358
00359 for( int i = 0; i < numColors; i++ )
00360 colorToNumMap[ colorList[i] ] = i;
00361
00362
00363 delete [] myColorList;
00364 delete [] allColorList;
00365 }