ColPack: Recovery/HessianRecovery.cpp Source File

Go to the documentation of this file.
00001 /************************************************************************************
00002     Copyright (C) 2005-2008 Assefaw H. Gebremedhin, Arijit Tarafdar, Duc Nguyen,
00003     Alex Pothen
00004 
00005     This file is part of ColPack.
00006 
00007     ColPack is free software: you can redistribute it and/or modify
00008     it under the terms of the GNU Lesser General Public License as published
00009     by the Free Software Foundation, either version 3 of the License, or
00010     (at your option) any later version.
00011 
00012     ColPack is distributed in the hope that it will be useful,
00013     but WITHOUT ANY WARRANTY; without even the implied warranty of
00014     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015     GNU Lesser General Public License for more details.
00016 
00017     You should have received a copy of the GNU Lesser General Public License
00018     along with ColPack.  If not, see <http://www.gnu.org/licenses/>.
00019 ************************************************************************************/
00020 
00021 #include "ColPackHeaders.h"
00022 
00023 using namespace std;
00024 
00025 namespace ColPack
00026 {
00027 
00028         int HessianRecovery::DirectRecover_RowCompressedFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, double*** dp3_HessianValue) {
00029                 if(g==NULL) {
00030                         cerr<<"g==NULL"<<endl;
00031                         return _FALSE;
00032                 }
00033 
00034                 if(AF_available) reset();
00035 
00036                 int rowCount = g->GetVertexCount();
00037                 int colorCount = g->GetVertexColorCount();
00038                 //cout<<"colorCount="<<colorCount<<endl;
00039                 vector<int> vi_VertexColors;
00040                 g->GetVertexColors(vi_VertexColors);
00041 
00042                 //Do (column-)color statistic for each row, i.e., see how many elements in that row has color 0, color 1 ...
00043                 int** colorStatistic = new int*[rowCount];      //color statistic for each row. For example, colorStatistic[0] is color statistic for row 0
00044                                                                                                         //If row 0 has 5 columns with color 3 => colorStatistic[0][3] = 5;
00045                 //Allocate memory for colorStatistic[rowCount][colorCount] and initilize the matrix
00046                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00047                         colorStatistic[i] = new int[colorCount];
00048                         for(unsigned int j=0; j < (unsigned int)colorCount; j++) colorStatistic[i][j] = 0;
00049                 }
00050 
00051                 //populate colorStatistic
00052                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00053                         int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00054                         for(unsigned int j=1; j <= (unsigned int)numOfNonZeros; j++) {
00055                                 //non-zero in the Hessian: [i][uip2_HessianSparsityPattern[i][j]]
00056                                 //color of that column: vi_VertexColors[uip2_HessianSparsityPattern[i][j]]
00057                                 colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]++;
00058                         }
00059                 }
00060 
00061                 //allocate memory for *dp3_HessianValue. The dp3_HessianValue and uip2_HessianSparsityPattern matrices should have the same size
00062                 *dp3_HessianValue = new double*[rowCount];
00063                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00064                         unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00065                         (*dp3_HessianValue)[i] = new double[numOfNonZeros+1];
00066                         (*dp3_HessianValue)[i][0] = numOfNonZeros; //initialize value of the 1st entry
00067                         for(unsigned int j=1; j <= numOfNonZeros; j++) (*dp3_HessianValue)[i][j] = 0.; //initialize value of other entries
00068                 }
00069 
00070                 //Now, go to the main part, recover the values of non-zero entries in the Hessian
00071                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00072                         unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00073                         for(unsigned int j=1; j <= numOfNonZeros; j++) {
00074                                 if(i == uip2_HessianSparsityPattern[i][j]) { // the non-zero is in the diagonal of the matrix
00075                                         (*dp3_HessianValue)[i][j] = dp2_CompressedMatrix[i][vi_VertexColors[i]];
00076                                         //printf("Recover diagonal (*dp3_HessianValue)[%d][%d] = %f from dp2_CompressedMatrix[%d][%d] \n", i, j, dp2_CompressedMatrix[i][vi_VertexColors[i]], i, vi_VertexColors[i]);         
00077                                 }
00078                                 else {// i != uip2_HessianSparsityPattern[i][j] // the non-zero is NOT in the diagonal of the matrix
00079                                         if(colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]==1) {
00080                                                 (*dp3_HessianValue)[i][j] = dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]];
00081                                                 //printf("Recover (*dp3_HessianValue)[%d][%d] = %f from dp2_CompressedMatrix[%d][%d] \n", i, j, dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]], i, vi_VertexColors[uip2_HessianSparsityPattern[i][j]]);         
00082                                         }
00083                                         else {
00084                                                 (*dp3_HessianValue)[i][j] = dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]];
00085                                                 //printf("Recover (*dp3_HessianValue)[%d][%d] = %f from dp2_CompressedMatrix[%d][%d] \n", i, j, dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]], uip2_HessianSparsityPattern[i][j], vi_VertexColors[i]);         
00086                                         }
00087                                 }
00088                         }
00089                 }
00090 
00091                 free_2DMatrix(colorStatistic, rowCount);
00092                 colorStatistic = NULL;
00093 
00094                 AF_available = true;
00095                 i_AF_rowCount = rowCount;
00096                 dp2_AF_Value = *dp3_HessianValue;
00097 
00098                 return (_TRUE);
00099         }
00100 
00101 
00102         int HessianRecovery::DirectRecover_SparseSolversFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) {
00103                 if(g==NULL) {
00104                         cerr<<"g==NULL"<<endl;
00105                         return _FALSE;
00106                 }
00107 
00108                 if(SSF_available) {
00109 cout<<"SSF_available="<<SSF_available<<endl; Pause();
00110                         reset();
00111                 }
00112 
00113                 int rowCount = g->GetVertexCount();
00114                 int colorCount = g->GetVertexColorCount();
00115                 //cout<<"colorCount="<<colorCount<<endl;
00116                 vector<int> vi_VertexColors;
00117                 g->GetVertexColors(vi_VertexColors);
00118 
00119                 //g->PrintGraph();
00120 
00121                 unsigned int numOfNonZerosInHessianValue = RowCompressedFormat_2_SparseSolversFormat_StructureOnly(uip2_HessianSparsityPattern, rowCount, ip2_RowIndex, ip2_ColumnIndex);
00122 
00123                 //Do (column-)color statistic for each row, i.e., see how many elements in that row has color 0, color 1 ...
00124                 int** colorStatistic = new int*[rowCount];      //color statistic for each row. For example, colorStatistic[0] is color statistic for row 0
00125                                                                                                         //If row 0 has 5 columns with color 3 => colorStatistic[0][3] = 5;
00126                 //Allocate memory for colorStatistic[rowCount][colorCount] and initilize the matrix
00127                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00128                         colorStatistic[i] = new int[colorCount];
00129                         for(unsigned int j=0; j < (unsigned int)colorCount; j++) colorStatistic[i][j] = 0;
00130                 }
00131 
00132                 //populate colorStatistic
00133                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00134                         int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00135                         for(unsigned int j=1; j <= (unsigned int)numOfNonZeros; j++) {
00136                                 //non-zero in the Hessian: [i][uip2_HessianSparsityPattern[i][j]]
00137                                 //color of that column: vi_VertexColors[uip2_HessianSparsityPattern[i][j]]
00138                                 colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]++;
00139                         }
00140                 }
00141 
00142                 //cout<<"allocate memory for *dp2_HessianValue rowCount="<<rowCount<<endl;
00143                 //printf("i=%d\t numOfNonZerosInHessianValue=%d \n", i, numOfNonZerosInHessianValue);
00144                 (*dp2_HessianValue) = new double[numOfNonZerosInHessianValue]; //allocate memory for *dp2_JacobianValue.
00145                 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) (*dp2_HessianValue)[i] = 0.; //initialize value of other entries
00146 
00147 
00148                 //Now, go to the main part, recover the values of non-zero entries in the Hessian
00149                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00150                         unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00151                         unsigned int offset = 0;
00152                         //printf("\ti=%d, \t NumOfNonzeros=%d,\t  \n",i,numOfNonZeros);
00153                         for(unsigned int j=1; j <= numOfNonZeros; j++) {
00154                                 if (i > uip2_HessianSparsityPattern[i][j]) {
00155                                   offset++;
00156                                   continue;
00157                                 }
00158                                 else if(i == uip2_HessianSparsityPattern[i][j]) { // the non-zero is in the diagonal of the matrix
00159                                         (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = dp2_CompressedMatrix[i][vi_VertexColors[i]];
00160                                         //printf("DIAGONAL Recover (*dp2_HessianValue)[%d] = %f from dp2_CompressedMatrix[%d][%d] \n", (*ip2_RowIndex)[i] + j - offset - 1, (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - 1],i,vi_VertexColors[i]);
00161                                         //printf("\t (*ip2_RowIndex)[i = %d] = %d, j = %d, offset = %d \n", i, (*ip2_RowIndex)[i], j, offset);
00162                                 }
00163                                 else {// i != uip2_HessianSparsityPattern[i][j] // the non-zero is NOT in the diagonal of the matrix
00164                                         if(colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]==1) {
00165                                                 (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]];
00166                                           //printf("Recover (*dp2_HessianValue)[%d] = %f from dp2_CompressedMatrix[%d][%d] \n", (*ip2_RowIndex)[i] + j - offset - 1, (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - 1], i , vi_VertexColors[uip2_HessianSparsityPattern[i][j]]);
00167                                         }
00168                                         else {
00169                                                 (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]];
00170                                           //printf("Recover (*dp2_HessianValue)[%d] = %f from dp2_CompressedMatrix[%d][%d] \n", (*ip2_RowIndex)[i] + j - offset - 1, (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - 1], uip2_HessianSparsityPattern[i][j], vi_VertexColors[i]);
00171                                         }
00172                                 }
00173                         }
00174                 }
00175 
00176                 free_2DMatrix(colorStatistic, rowCount);
00177                 colorStatistic = NULL;
00178 
00179 
00180                 //Making the array indices to start at 1 instead of 0 to conform with theIntel MKL sparse storage scheme for the direct sparse solvers
00181                 for(unsigned int i=0; i <= (unsigned int) rowCount ; i++) {
00182                   (*ip2_RowIndex)[i]++;
00183                 }
00184                 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) {
00185                   (*ip2_ColumnIndex)[i]++;
00186                 }
00187 
00188 
00189                 SSF_available = true;
00190                 i_SSF_rowCount = rowCount;
00191                 ip_SSF_RowIndex = *ip2_RowIndex;
00192                 ip_SSF_ColumnIndex = *ip2_ColumnIndex;
00193                 dp_SSF_Value = *dp2_HessianValue;
00194 
00195                 return (_TRUE);
00196         }
00197 
00198 
00199         int HessianRecovery::DirectRecover_CoordinateFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) {
00200                 if(g==NULL) {
00201                         cerr<<"g==NULL"<<endl;
00202                         return _FALSE;
00203                 }
00204 
00205                 if(CF_available) reset();
00206 
00207                 int rowCount = g->GetVertexCount();
00208                 int colorCount = g->GetVertexColorCount();
00209                 //cout<<"colorCount="<<colorCount<<endl;
00210                 vector<int> vi_VertexColors;
00211                 g->GetVertexColors(vi_VertexColors);
00212 
00213                 //Do (column-)color statistic for each row, i.e., see how many elements in that row has color 0, color 1 ...
00214                 int** colorStatistic = new int*[rowCount];      //color statistic for each row. For example, colorStatistic[0] is color statistic for row 0
00215                                                                                                         //If row 0 has 5 columns with color 3 => colorStatistic[0][3] = 5;
00216                 //Allocate memory for colorStatistic[rowCount][colorCount] and initilize the matrix
00217                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00218                         colorStatistic[i] = new int[colorCount];
00219                         for(unsigned int j=0; j < (unsigned int)colorCount; j++) colorStatistic[i][j] = 0;
00220                 }
00221 
00222                 //populate colorStatistic
00223                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00224                         int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00225                         for(unsigned int j=1; j <= (unsigned int)numOfNonZeros; j++) {
00226                                 //non-zero in the Hessian: [i][uip2_HessianSparsityPattern[i][j]]
00227                                 //color of that column: vi_VertexColors[uip2_HessianSparsityPattern[i][j]]
00228                                 colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]++;
00229                         }
00230                 }
00231 
00232                 vector<unsigned int> RowIndex;
00233                 vector<unsigned int> ColumnIndex;
00234                 vector<double> HessianValue;
00235 
00236                 //Now, go to the main part, recover the values of non-zero entries in the Hessian
00237                 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00238                         unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00239                         for(unsigned int j=1; j <= numOfNonZeros; j++) {
00240                                 if(uip2_HessianSparsityPattern[i][j]<i) continue;
00241 
00242                                 if(i == uip2_HessianSparsityPattern[i][j]) { // the non-zero is in the diagonal of the matrix
00243                                         HessianValue.push_back(dp2_CompressedMatrix[i][vi_VertexColors[i]]);
00244                                 }
00245                                 else {// i != uip2_HessianSparsityPattern[i][j] // the non-zero is NOT in the diagonal of the matrix
00246                                         if(colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]==1) {
00247                                                 HessianValue.push_back(dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]);
00248                                         }
00249                                         else {
00250                                                 HessianValue.push_back(dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]]);
00251                                         }
00252                                 }
00253                                 RowIndex.push_back(i);
00254                                 ColumnIndex.push_back(uip2_HessianSparsityPattern[i][j]);
00255                         }
00256                 }
00257 
00258                 unsigned int numOfNonZeros = RowIndex.size();
00259                 (*ip2_RowIndex) = new unsigned int[numOfNonZeros];
00260                 (*ip2_ColumnIndex) = new unsigned int[numOfNonZeros];
00261                 (*dp2_HessianValue) = new double[numOfNonZeros]; //allocate memory for *dp2_HessianValue.
00262 
00263                 for(int i=0; i < numOfNonZeros; i++) {
00264                         (*ip2_RowIndex)[i] = RowIndex[i];
00265                         (*ip2_ColumnIndex)[i] = ColumnIndex[i];
00266                         (*dp2_HessianValue)[i] = HessianValue[i];
00267                 }
00268 
00269                 free_2DMatrix(colorStatistic, rowCount);
00270                 colorStatistic = NULL;
00271 
00272 
00273                 CF_available = true;
00274                 i_CF_rowCount = rowCount;
00275                 ip_CF_RowIndex = *ip2_RowIndex;
00276                 ip_CF_ColumnIndex = *ip2_ColumnIndex;
00277                 dp_CF_Value = *dp2_HessianValue;
00278 
00279                 return (numOfNonZeros);
00280         }
00281 
00282         int HessianRecovery::IndirectRecover_RowCompressedFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, double*** dp3_HessianValue) {
00283                 if(g==NULL) {
00284                         cerr<<"g==NULL"<<endl;
00285                         return _FALSE;
00286                 }
00287 
00288                 if(AF_available) reset();
00289 
00290                 int i=0,j=0;
00291                 int i_VertexCount = _UNKNOWN;
00292                 int i_EdgeID, i_SetID;
00293                 vector<int> vi_Sets;
00294                 map< int, vector<int> > mivi_VertexSets;
00295 
00296                 vector<int> vi_Vertices;
00297                 g->GetVertices(vi_Vertices);
00298 
00299                 vector<int> vi_Edges;
00300                 g->GetEdges(vi_Edges);
00301 
00302                 vector<int> vi_VertexColors;
00303                 g->GetVertexColors(vi_VertexColors);
00304 
00305                 map< int, map< int, int> > mimi2_VertexEdgeMap;
00306                 g->GetVertexEdgeMap(mimi2_VertexEdgeMap);
00307 
00308                 DisjointSets ds_DisjointSets;
00309                 g->GetDisjointSets(ds_DisjointSets);
00310 
00311                 //populate vi_Sets & mivi_VertexSets
00312                 vi_Sets.clear();
00313                 mivi_VertexSets.clear();
00314 
00315                 i_VertexCount = g->GetVertexCount();
00316 
00317                 for(i=0; i<i_VertexCount; i++) // for each vertex A (indexed by i)
00318                 {
00319                 for(j=vi_Vertices[i]; j<vi_Vertices[STEP_UP(i)]; j++) // for each of the vertex B that connect to A
00320                         {
00321                                 if(i < vi_Edges[j]) // if the index of A (i) is less than the index of B (vi_Edges[j])
00322                                                                                 //basicly each edge is represented by (vertex with smaller ID, vertex with larger ID). This way, we don't insert a specific edge twice
00323                                 {
00324                                         i_EdgeID = mimi2_VertexEdgeMap[i][vi_Edges[j]];
00325 
00326                                         i_SetID = ds_DisjointSets.FindAndCompress(i_EdgeID);
00327 
00328                                         if(i_SetID == i_EdgeID) // that edge is the root of the set => create new set
00329                                         {
00330                                                 vi_Sets.push_back(i_SetID);
00331                                         }
00332 
00333                                         mivi_VertexSets[i_SetID].push_back(i);
00334                                         mivi_VertexSets[i_SetID].push_back(vi_Edges[j]);
00335                                 }
00336                         }
00337                 }
00338 
00339                 int i_MaximumVertexDegree;
00340 
00341                 int i_HighestInducedVertexDegree;
00342 
00343                 int i_LeafVertex, i_ParentVertex, i_PresentVertex;
00344 
00345                 int i_VertexDegree;
00346 
00347                 int i_SetCount, i_SetSize;
00348                 //i_SetCount = vi_Sets.size();
00349                 //i_SetSize: size (number of edges?) in a bicolored tree
00350 
00351                 double d_Value;
00352 
00353                 vector<int> vi_EvaluatedDiagonals;
00354 
00355                 vector<int> vi_InducedVertexDegrees;
00356 
00357                 vector<double> vd_IncludedVertices;
00358 
00359                 vector< vector<int> > v2i_VertexAdjacency;
00360 
00361                 vector< vector<double> > v2d_NonzeroAdjacency;
00362 
00363                 vector< list<int> > vli_GroupedInducedVertexDegrees;
00364 
00365                 vector< list<int>::iterator > vlit_VertexLocations;
00366 
00367                 i_MaximumVertexDegree = g->GetMaximumVertexDegree();
00368 
00369         #if DEBUG == 5103
00370 
00371                 cout<<endl;
00372                 cout<<"DEBUG 5103 | Hessian Evaluation | Bicolored Sets"<<endl;
00373                 cout<<endl;
00374 
00375                 i_SetCount = (signed) vi_Sets.size();
00376 
00377                 for(i=0; i<i_SetCount; i++)
00378                 {
00379                         cout<<STEP_UP(vi_Sets[i])<<"\t"<<" : ";
00380 
00381                         i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size();
00382 
00383                         for(j=0; j<i_SetSize; j++)
00384                         {
00385                                 if(j == STEP_DOWN(i_SetSize))
00386                                 {
00387                                 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<" ("<<i_SetSize<<")"<<endl;
00388                                 }
00389                                 else
00390                                 {
00391                                 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<", ";
00392                                 }
00393                         }
00394                 }
00395 
00396                 cout<<endl;
00397                 cout<<"[Set Count = "<<i_SetCount<<"]"<<endl;
00398                 cout<<endl;
00399 
00400         #endif
00401 
00402                 //Step 5: from here on
00403                 i_VertexCount = g->GetVertexCount();
00404 
00405                 v2i_VertexAdjacency.clear();
00406                 v2i_VertexAdjacency.resize((unsigned) i_VertexCount);
00407 
00408                 v2d_NonzeroAdjacency.clear();
00409                 v2d_NonzeroAdjacency.resize((unsigned) i_VertexCount);
00410 
00411                 vi_EvaluatedDiagonals.clear();
00412                 vi_EvaluatedDiagonals.resize((unsigned) i_VertexCount, _FALSE);
00413 
00414                 vi_InducedVertexDegrees.clear();
00415                 vi_InducedVertexDegrees.resize((unsigned) i_VertexCount, _FALSE);
00416 
00417                 vd_IncludedVertices.clear();
00418                 vd_IncludedVertices.resize((unsigned) i_VertexCount, _UNKNOWN);
00419 
00420                 i_ParentVertex = _UNKNOWN;
00421 
00422                 i_SetCount = (signed) vi_Sets.size();
00423 
00424                 for(i=0; i<i_SetCount; i++)
00425                 {
00426                         vli_GroupedInducedVertexDegrees.clear();
00427                         vli_GroupedInducedVertexDegrees.resize((unsigned) STEP_UP(i_MaximumVertexDegree));
00428 
00429                         vlit_VertexLocations.clear();
00430                         vlit_VertexLocations.resize((unsigned) i_VertexCount);
00431 
00432                         i_HighestInducedVertexDegree = _UNKNOWN;
00433 
00434                         i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size();
00435 
00436                         for(j=0; j<i_SetSize; j++)
00437                         {
00438                                 i_PresentVertex = mivi_VertexSets[vi_Sets[i]][j];
00439 
00440                                 vd_IncludedVertices[i_PresentVertex] = _FALSE;
00441 
00442                                 if(vi_InducedVertexDegrees[i_PresentVertex] != _FALSE)
00443                                 {
00444                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].erase(vlit_VertexLocations[i_PresentVertex]);
00445                                 }
00446 
00447                                 vi_InducedVertexDegrees[i_PresentVertex]++;
00448 
00449                                 if(i_HighestInducedVertexDegree < vi_InducedVertexDegrees[i_PresentVertex])
00450                                 {
00451                                         i_HighestInducedVertexDegree = vi_InducedVertexDegrees[i_PresentVertex];
00452                                 }
00453 
00454                                 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].push_front(i_PresentVertex);
00455 
00456                                 vlit_VertexLocations[i_PresentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].begin();
00457                         }
00458 
00459         #if DEBUG == 5103
00460 
00461                         int k;
00462 
00463                         list<int>::iterator lit_ListIterator;
00464 
00465                         cout<<endl;
00466                         cout<<"DEBUG 5103 | Hessian Evaluation | Induced Vertex Degrees | Set "<<STEP_UP(i)<<endl;
00467                         cout<<endl;
00468 
00469                         for(j=0; j<STEP_UP(i_HighestInducedVertexDegree); j++)
00470                         {
00471                                 i_SetSize = (signed) vli_GroupedInducedVertexDegrees[j].size();
00472 
00473                                 if(i_SetSize == _FALSE)
00474                                 {
00475                                         continue;
00476                                 }
00477 
00478                                 k = _FALSE;
00479 
00480                                 cout<<"Degree "<<j<<"\t"<<" : ";
00481 
00482                                 for(lit_ListIterator=vli_GroupedInducedVertexDegrees[j].begin(); lit_ListIterator!=vli_GroupedInducedVertexDegrees[j].end(); lit_ListIterator++)
00483                                 {
00484                                         if(k == STEP_DOWN(i_SetSize))
00485                                         {
00486                                                 cout<<STEP_UP(*lit_ListIterator)<<" ("<<i_SetSize<<")"<<endl;
00487                                         }
00488                                         else
00489                                         {
00490                                                 cout<<STEP_UP(*lit_ListIterator)<<", ";
00491                                         }
00492 
00493                                         k++;
00494                                 }
00495                         }
00496 
00497         #endif
00498 
00499         #if DEBUG == 5103
00500 
00501                         cout<<endl;
00502                         cout<<"DEBUG 5103 | Hessian Evaluation | Retrieved Elements"<<"| Set "<<STEP_UP(i)<<endl;
00503                         cout<<endl;
00504 
00505         #endif
00506 //#define DEBUG 5103
00507                         //get the diagonal values
00508                         for (int index = 0; index < i_VertexCount; index++) {
00509                                 if(vi_EvaluatedDiagonals[index] == _FALSE)
00510                                 {
00511                                         d_Value = dp2_CompressedMatrix[index][vi_VertexColors[index]];
00512 
00513         #if DEBUG == 5103
00514 
00515                                         cout<<"Element["<<STEP_UP(index)<<"]["<<STEP_UP(index)<<"] = "<<d_Value<<endl;
00516 
00517         #endif
00518                                         v2i_VertexAdjacency[index].push_back(index);
00519                                         v2d_NonzeroAdjacency[index].push_back(d_Value);
00520 
00521                                         vi_EvaluatedDiagonals[index] = _TRUE;
00522 
00523                                 }
00524                         }
00525 
00526                         for ( ; ; )
00527                         {
00528                                 if(vli_GroupedInducedVertexDegrees[_TRUE].empty()) // If there is no leaf left on the color tree
00529                                 {
00530                                         i_LeafVertex = vli_GroupedInducedVertexDegrees[_FALSE].front();
00531 
00532                                         vi_InducedVertexDegrees[i_LeafVertex] = _FALSE;
00533 
00534                                         vd_IncludedVertices[i_LeafVertex] = _UNKNOWN;
00535 
00536                                         break;
00537                                 }
00538 
00539                                 i_LeafVertex = vli_GroupedInducedVertexDegrees[_TRUE].front();
00540 
00541                                 vli_GroupedInducedVertexDegrees[_TRUE].pop_front();
00542 
00543 
00544                                 //Find i_ParentVertex
00545                                 for(j=vi_Vertices[i_LeafVertex]; j<vi_Vertices[STEP_UP(i_LeafVertex)]; j++)
00546                                 {
00547                                         if(vd_IncludedVertices[vi_Edges[j]] != _UNKNOWN)
00548                                         {
00549                                                 i_ParentVertex = vi_Edges[j];
00550 
00551                                                 break;
00552                                         }
00553                                 }
00554 
00555                                 d_Value = dp2_CompressedMatrix[i_LeafVertex][vi_VertexColors[i_ParentVertex]] - vd_IncludedVertices[i_LeafVertex];
00556 
00557                                 vd_IncludedVertices[i_ParentVertex] += d_Value;
00558 
00559                                 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE;
00560                                 vd_IncludedVertices[i_LeafVertex] = _UNKNOWN;
00561                                 if(vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].size()>1) {
00562                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].erase(vlit_VertexLocations[i_ParentVertex]);
00563                                 }
00564                                 else {
00565                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].pop_back();
00566                                 }
00567 
00568                                 vi_InducedVertexDegrees[i_ParentVertex]--;
00569                                 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].push_back(i_ParentVertex);
00570 
00571                                 //Update position of the iterator pointing to i_ParentVertex in "InducedVertexDegrees" structure
00572                                 vlit_VertexLocations[i_ParentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].end();
00573                                 --vlit_VertexLocations[i_ParentVertex];
00574 
00575                                 v2i_VertexAdjacency[i_LeafVertex].push_back(i_ParentVertex);
00576                                 v2d_NonzeroAdjacency[i_LeafVertex].push_back(d_Value);
00577 
00578                                 v2i_VertexAdjacency[i_ParentVertex].push_back(i_LeafVertex);
00579                                 v2d_NonzeroAdjacency[i_ParentVertex].push_back(d_Value);
00580 
00581 
00582         #if DEBUG == 5103
00583 
00584                                 cout<<"Element["<<STEP_UP(i_LeafVertex)<<"]["<<STEP_UP(i_ParentVertex)<<"] = "<<d_Value<<endl;
00585         #endif
00586 
00587                         }
00588                 }
00589 
00590 
00591                 //allocate memory for *dp3_HessianValue. The dp3_HessianValue and uip2_HessianSparsityPattern matrices should have the same size
00592                 *dp3_HessianValue = new double*[i_VertexCount];
00593                 for(unsigned int i=0; i < (unsigned int)i_VertexCount; i++) {
00594                         unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0];
00595                         (*dp3_HessianValue)[i] = new double[numOfNonZeros+1];
00596                         (*dp3_HessianValue)[i][0] = numOfNonZeros; //initialize value of the 1st entry
00597                         for(unsigned int j=1; j <= numOfNonZeros; j++) (*dp3_HessianValue)[i][j] = 0.; //initialize value of other entries
00598                 }
00599 
00600                 //populate dp3_HessianValue row by row, column by column
00601                 for(i=0; i<i_VertexCount; i++) {
00602                         int NumOfNonzeros = uip2_HessianSparsityPattern[i][0];
00603                         i_VertexDegree = (signed) v2i_VertexAdjacency[i].size();
00604                         for(j=1; j<=NumOfNonzeros; j++) {
00605                                 int targetColumnID = uip2_HessianSparsityPattern[i][j];
00606                                 for (int k=0; k<i_VertexDegree; k++) {// search through the v2i_VertexAdjacency matrix to find the correct column
00607                                         if(targetColumnID == v2i_VertexAdjacency[i][k]) { //found it
00608                                                 (*dp3_HessianValue)[i][j] = v2d_NonzeroAdjacency[i][k];
00609                                                 break;
00610                                         }
00611                                 }
00612                         }
00613                 }
00614 
00615         #undef DEBUG
00616 
00617                 AF_available = true;
00618                 i_AF_rowCount = i_VertexCount;
00619                 dp2_AF_Value = *dp3_HessianValue;
00620 
00621                 return (_TRUE);
00622         }
00623 
00624 
00625         int HessianRecovery::IndirectRecover_SparseSolversFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) {
00626                 if(g==NULL) {
00627                         cerr<<"g==NULL"<<endl;
00628                         return _FALSE;
00629                 }
00630 
00631 
00632                 if(SSF_available) {
00633 cout<<"SSF_available="<<SSF_available<<endl; Pause();
00634                         reset();
00635                 }
00636 
00637                 unsigned int numOfNonZerosInHessianValue = RowCompressedFormat_2_SparseSolversFormat_StructureOnly(uip2_HessianSparsityPattern, g->GetVertexCount(), ip2_RowIndex, ip2_ColumnIndex);
00638 
00639                 int i=0,j=0;
00640                 int i_VertexCount = _UNKNOWN;
00641                 int i_EdgeID, i_SetID;
00642                 vector<int> vi_Sets;
00643                 map< int, vector<int> > mivi_VertexSets;
00644 
00645                 vector<int> vi_Vertices;
00646                 g->GetVertices(vi_Vertices);
00647 
00648                 vector<int> vi_Edges;
00649                 g->GetEdges(vi_Edges);
00650 
00651                 vector<int> vi_VertexColors;
00652                 g->GetVertexColors(vi_VertexColors);
00653 
00654                 map< int, map< int, int> > mimi2_VertexEdgeMap;
00655                 g->GetVertexEdgeMap(mimi2_VertexEdgeMap);
00656 
00657                 DisjointSets ds_DisjointSets;
00658                 g->GetDisjointSets(ds_DisjointSets);
00659 
00660                 //populate vi_Sets & mivi_VertexSets
00661                 vi_Sets.clear();
00662                 mivi_VertexSets.clear();
00663 
00664                 i_VertexCount = g->GetVertexCount();
00665 
00666                 for(i=0; i<i_VertexCount; i++) // for each vertex A (indexed by i)
00667                 {
00668                         for(j=vi_Vertices[i]; j<vi_Vertices[STEP_UP(i)]; j++) // for each of the vertex B that connect to A
00669                         {
00670                                 if(i < vi_Edges[j]) // if the index of A (i) is less than the index of B (vi_Edges[j])
00671                                                                                 //basic each edge is represented by (vertex with smaller ID, vertex with larger ID). This way, we don't insert a specific edge twice
00672                                 {
00673                                         i_EdgeID = mimi2_VertexEdgeMap[i][vi_Edges[j]];
00674 
00675                                         i_SetID = ds_DisjointSets.FindAndCompress(i_EdgeID);
00676 
00677                                         if(i_SetID == i_EdgeID) // that edge is the root of the set => create new set
00678                                         {
00679                                                 vi_Sets.push_back(i_SetID);
00680                                         }
00681 
00682                                         mivi_VertexSets[i_SetID].push_back(i);
00683                                         mivi_VertexSets[i_SetID].push_back(vi_Edges[j]);
00684                                 }
00685                         }
00686                 }
00687 
00688                 int i_MaximumVertexDegree;
00689 
00690                 int i_HighestInducedVertexDegree;
00691 
00692                 int i_LeafVertex, i_ParentVertex, i_PresentVertex;
00693 
00694                 int i_VertexDegree;
00695 
00696                 int i_SetCount, i_SetSize;
00697 
00698                 double d_Value;
00699 
00700                 vector<int> vi_EvaluatedDiagonals;
00701 
00702                 vector<int> vi_InducedVertexDegrees;
00703 
00704                 vector<double> vd_IncludedVertices;
00705 
00706                 vector< vector<int> > v2i_VertexAdjacency;
00707 
00708                 vector< vector<double> > v2d_NonzeroAdjacency;
00709 
00710                 vector< list<int> > vli_GroupedInducedVertexDegrees;
00711 
00712                 vector< list<int>::iterator > vlit_VertexLocations;
00713 
00714                 i_MaximumVertexDegree = g->GetMaximumVertexDegree();
00715 
00716         #if DEBUG == 5103
00717 
00718                 cout<<endl;
00719                 cout<<"DEBUG 5103 | Hessian Evaluation | Bicolored Sets"<<endl;
00720                 cout<<endl;
00721 
00722                 i_SetCount = (signed) vi_Sets.size();
00723 
00724                 for(i=0; i<i_SetCount; i++)
00725                 {
00726                         cout<<STEP_UP(vi_Sets[i])<<"\t"<<" : ";
00727 
00728                         i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size();
00729 
00730                         for(j=0; j<i_SetSize; j++)
00731                         {
00732                                 if(j == STEP_DOWN(i_SetSize))
00733                                 {
00734                                 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<" ("<<i_SetSize<<")"<<endl;
00735                                 }
00736                                 else
00737                                 {
00738                                 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<", ";
00739                                 }
00740                         }
00741                 }
00742 
00743                 cout<<endl;
00744                 cout<<"[Set Count = "<<i_SetCount<<"]"<<endl;
00745                 cout<<endl;
00746 
00747         #endif
00748 
00749                 //Step 5: from here on
00750                 i_VertexCount = g->GetVertexCount();
00751 
00752                 v2i_VertexAdjacency.clear();
00753                 v2i_VertexAdjacency.resize((unsigned) i_VertexCount);
00754 
00755                 v2d_NonzeroAdjacency.clear();
00756                 v2d_NonzeroAdjacency.resize((unsigned) i_VertexCount);
00757 
00758                 vi_EvaluatedDiagonals.clear();
00759                 vi_EvaluatedDiagonals.resize((unsigned) i_VertexCount, _FALSE);
00760 
00761                 vi_InducedVertexDegrees.clear();
00762                 vi_InducedVertexDegrees.resize((unsigned) i_VertexCount, _FALSE);
00763 
00764                 vd_IncludedVertices.clear();
00765                 vd_IncludedVertices.resize((unsigned) i_VertexCount, _UNKNOWN);
00766 
00767                 i_ParentVertex = _UNKNOWN;
00768 
00769                 i_SetCount = (signed) vi_Sets.size();
00770 
00771                 for(i=0; i<i_SetCount; i++)
00772                 {
00773                         vli_GroupedInducedVertexDegrees.clear();
00774                         vli_GroupedInducedVertexDegrees.resize((unsigned) STEP_UP(i_MaximumVertexDegree));
00775 
00776                         vlit_VertexLocations.clear();
00777                         vlit_VertexLocations.resize((unsigned) i_VertexCount);
00778 
00779                         i_HighestInducedVertexDegree = _UNKNOWN;
00780 
00781                         i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size();
00782 
00783                         for(j=0; j<i_SetSize; j++)
00784                         {
00785                                 i_PresentVertex = mivi_VertexSets[vi_Sets[i]][j];
00786 
00787                                 vd_IncludedVertices[i_PresentVertex] = _FALSE;
00788 
00789                                 if(vi_InducedVertexDegrees[i_PresentVertex] != _FALSE)
00790                                 {
00791                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].erase(vlit_VertexLocations[i_PresentVertex]);
00792                                 }
00793 
00794                                 vi_InducedVertexDegrees[i_PresentVertex]++;
00795 
00796                                 if(i_HighestInducedVertexDegree < vi_InducedVertexDegrees[i_PresentVertex])
00797                                 {
00798                                         i_HighestInducedVertexDegree = vi_InducedVertexDegrees[i_PresentVertex];
00799                                 }
00800 
00801                                 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].push_front(i_PresentVertex);
00802 
00803                                 vlit_VertexLocations[i_PresentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].begin();
00804                         }
00805 
00806         #if DEBUG == 5103
00807 
00808                         int k;
00809 
00810                         list<int>::iterator lit_ListIterator;
00811 
00812                         cout<<endl;
00813                         cout<<"DEBUG 5103 | Hessian Evaluation | Induced Vertex Degrees | Set "<<STEP_UP(i)<<endl;
00814                         cout<<endl;
00815 
00816                         for(j=0; j<STEP_UP(i_HighestInducedVertexDegree); j++)
00817                         {
00818                                 i_SetSize = (signed) vli_GroupedInducedVertexDegrees[j].size();
00819 
00820                                 if(i_SetSize == _FALSE)
00821                                 {
00822                                         continue;
00823                                 }
00824 
00825                                 k = _FALSE;
00826 
00827                                 cout<<"Degree "<<j<<"\t"<<" : ";
00828 
00829                                 for(lit_ListIterator=vli_GroupedInducedVertexDegrees[j].begin(); lit_ListIterator!=vli_GroupedInducedVertexDegrees[j].end(); lit_ListIterator++)
00830                                 {
00831                                         if(k == STEP_DOWN(i_SetSize))
00832                                         {
00833                                                 cout<<STEP_UP(*lit_ListIterator)<<" ("<<i_SetSize<<")"<<endl;
00834                                         }
00835                                         else
00836                                         {
00837                                                 cout<<STEP_UP(*lit_ListIterator)<<", ";
00838                                         }
00839 
00840                                         k++;
00841                                 }
00842                         }
00843 
00844         #endif
00845 
00846         #if DEBUG == 5103
00847 
00848                         cout<<endl;
00849                         cout<<"DEBUG 5103 | Hessian Evaluation | Retrieved Elements"<<"| Set "<<STEP_UP(i)<<endl;
00850                         cout<<endl;
00851 
00852         #endif
00853 //#define DEBUG 5103
00854                         //get the diagonal values
00855                         for (int index = 0; index < i_VertexCount; index++) {
00856                                 if(vi_EvaluatedDiagonals[index] == _FALSE)
00857                                 {
00858                                         d_Value = dp2_CompressedMatrix[index][vi_VertexColors[index]];
00859 
00860         #if DEBUG == 5103
00861 
00862                                         cout<<"Element["<<STEP_UP(index)<<"]["<<STEP_UP(index)<<"] = "<<d_Value<<endl;
00863 
00864         #endif
00865                                         v2i_VertexAdjacency[index].push_back(index);
00866                                         v2d_NonzeroAdjacency[index].push_back(d_Value);
00867 
00868                                         vi_EvaluatedDiagonals[index] = _TRUE;
00869 
00870                                 }
00871                         }
00872 
00873                         for ( ; ; )
00874                         {
00875                                 if(vli_GroupedInducedVertexDegrees[_TRUE].empty()) // If there is no leaf left on the color tree
00876                                 {
00877                                         i_LeafVertex = vli_GroupedInducedVertexDegrees[_FALSE].front();
00878 
00879                                         vi_InducedVertexDegrees[i_LeafVertex] = _FALSE;
00880 
00881                                         vd_IncludedVertices[i_LeafVertex] = _UNKNOWN;
00882 
00883                                         break;
00884                                 }
00885 
00886                                 i_LeafVertex = vli_GroupedInducedVertexDegrees[_TRUE].front();
00887 
00888                                 vli_GroupedInducedVertexDegrees[_TRUE].pop_front();
00889 
00890                                 //Find i_ParentVertex
00891                                 for(j=vi_Vertices[i_LeafVertex]; j<vi_Vertices[STEP_UP(i_LeafVertex)]; j++)
00892                                 {
00893                                         if(vd_IncludedVertices[vi_Edges[j]] != _UNKNOWN)
00894                                         {
00895                                                 i_ParentVertex = vi_Edges[j];
00896 
00897                                                 break;
00898                                         }
00899                                 }
00900 
00901                                 d_Value = dp2_CompressedMatrix[i_LeafVertex][vi_VertexColors[i_ParentVertex]] - vd_IncludedVertices[i_LeafVertex];
00902 
00903                                 vd_IncludedVertices[i_ParentVertex] += d_Value;
00904 
00905                                 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE;
00906                                 vd_IncludedVertices[i_LeafVertex] = _UNKNOWN;
00907                                 if(vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].size()>1) {
00908                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].erase(vlit_VertexLocations[i_ParentVertex]);
00909                                 }
00910                                 else {
00911                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].pop_back();
00912                                 }
00913 
00914                                 vi_InducedVertexDegrees[i_ParentVertex]--;
00915                                 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].push_back(i_ParentVertex);
00916 
00917                                 //Update position of the iterator pointing to i_ParentVertex in "InducedVertexDegrees" structure
00918                                 vlit_VertexLocations[i_ParentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].end();
00919                                 --vlit_VertexLocations[i_ParentVertex];
00920 
00921                                 v2i_VertexAdjacency[i_LeafVertex].push_back(i_ParentVertex);
00922                                 v2d_NonzeroAdjacency[i_LeafVertex].push_back(d_Value);
00923 
00924                                 v2i_VertexAdjacency[i_ParentVertex].push_back(i_LeafVertex);
00925                                 v2d_NonzeroAdjacency[i_ParentVertex].push_back(d_Value);
00926 
00927         #if DEBUG == 5103
00928 
00929                                 cout<<"Element["<<STEP_UP(i_LeafVertex)<<"]["<<STEP_UP(i_ParentVertex)<<"] = "<<d_Value<<endl;
00930         #endif
00931 
00932                         }
00933                 }
00934 
00935 
00936 
00937                 //cout<<"allocate memory for *dp2_HessianValue rowCount="<<rowCount<<endl;
00938                 //printf("i=%d\t numOfNonZerosInHessianValue=%d \n", i, numOfNonZerosInHessianValue);
00939                 (*dp2_HessianValue) = new double[numOfNonZerosInHessianValue]; //allocate memory for *dp2_JacobianValue.
00940                 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) (*dp2_HessianValue)[i] = 0.; //initialize value of other entries
00941 
00942                 //populate dp2_HessianValue row by row, column by column
00943                 for(i=0; i<i_VertexCount; i++) {
00944                         int NumOfNonzeros = uip2_HessianSparsityPattern[i][0];
00945                         int offset = 0;
00946                         i_VertexDegree = (signed) v2i_VertexAdjacency[i].size();
00947                         for(j=1; j<=NumOfNonzeros; j++) {
00948                                 if( i > uip2_HessianSparsityPattern[i][j] ) {
00949                                   offset++;
00950                                   continue;
00951                                 }
00952                                 int targetColumnID = uip2_HessianSparsityPattern[i][j];
00953                                 for (int k=0; k<i_VertexDegree; k++) {// search through the v2i_VertexAdjacency matrix to find the correct column
00954                                         if(targetColumnID == v2i_VertexAdjacency[i][k]) { //found it
00955                                                 (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = v2d_NonzeroAdjacency[i][k];
00956                                                 break;
00957                                         }
00958                                 }
00959                         }
00960                 }
00961 
00962         #undef DEBUG
00963 
00964 
00965                 //Making the array indices to start at 1 instead of 0 to conform with theIntel MKL sparse storage scheme for the direct sparse solvers
00966                 for(unsigned int i=0; i <= (unsigned int) i_VertexCount ; i++) {
00967                   (*ip2_RowIndex)[i]++;
00968                 }
00969                 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) {
00970                   (*ip2_ColumnIndex)[i]++;
00971                 }
00972 
00973                 SSF_available = true;
00974                 i_SSF_rowCount = i_VertexCount;
00975                 ip_SSF_RowIndex = *ip2_RowIndex;
00976                 ip_SSF_ColumnIndex = *ip2_ColumnIndex;
00977                 dp_SSF_Value = *dp2_HessianValue;
00978 
00979                 return (_TRUE);
00980         }
00981 
00982         int HessianRecovery::IndirectRecover_CoordinateFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) {
00983 //#define DEBUG 5103
00984                 if(g==NULL) {
00985                         cerr<<"g==NULL"<<endl;
00986                         return _FALSE;
00987                 }
00988 
00989                 if(CF_available) reset();
00990 
00991                 int i=0,j=0;
00992                 int i_VertexCount = _UNKNOWN;
00993                 int i_EdgeID, i_SetID;
00994                 vector<int> vi_Sets;
00995                 map< int, vector<int> > mivi_VertexSets;
00996                 vector<int> vi_Vertices;
00997                 g->GetVertices(vi_Vertices);
00998                 vector<int> vi_Edges;
00999                 g->GetEdges(vi_Edges);
01000                 vector<int> vi_VertexColors;
01001                 g->GetVertexColors(vi_VertexColors);
01002                 map< int, map< int, int> > mimi2_VertexEdgeMap;
01003                 g->GetVertexEdgeMap(mimi2_VertexEdgeMap);
01004                 DisjointSets ds_DisjointSets;
01005                 g->GetDisjointSets(ds_DisjointSets);
01006 
01007                 //populate vi_Sets & mivi_VertexSets
01008                 vi_Sets.clear();
01009                 mivi_VertexSets.clear();
01010 
01011                 i_VertexCount = g->GetVertexCount();
01012 
01013                 for(i=0; i<i_VertexCount; i++) // for each vertex A (indexed by i)
01014                 {
01015                         for(j=vi_Vertices[i]; j<vi_Vertices[STEP_UP(i)]; j++) // for each of the vertex B that connect to A
01016                         {
01017                                 if(i < vi_Edges[j]) // if the index of A (i) is less than the index of B (vi_Edges[j])
01018                                                                                 //basic each edge is represented by (vertex with smaller ID, vertex with larger ID). This way, we don't insert a specific edge twice
01019                                 {
01020                                         i_EdgeID = mimi2_VertexEdgeMap[i][vi_Edges[j]];
01021 
01022                                         i_SetID = ds_DisjointSets.FindAndCompress(i_EdgeID);
01023 
01024                                         if(i_SetID == i_EdgeID) // that edge is the root of the set => create new set
01025                                         {
01026                                                 vi_Sets.push_back(i_SetID);
01027                                         }
01028 
01029                                         mivi_VertexSets[i_SetID].push_back(i);
01030                                         mivi_VertexSets[i_SetID].push_back(vi_Edges[j]);
01031                                 }
01032                         }
01033                 }
01034 
01035                 int i_MaximumVertexDegree;
01036 
01037                 int i_HighestInducedVertexDegree;
01038 
01039                 int i_LeafVertex, i_ParentVertex, i_PresentVertex;
01040 
01041                 int i_VertexDegree;
01042 
01043                 int i_SetCount, i_SetSize;
01044 
01045                 double d_Value;
01046 
01047                 vector<int> vi_EvaluatedDiagonals;
01048 
01049                 vector<int> vi_InducedVertexDegrees;
01050 
01051                 vector<double> vd_IncludedVertices;
01052 
01053                 vector< vector<int> > v2i_VertexAdjacency;
01054 
01055                 vector< vector<double> > v2d_NonzeroAdjacency;
01056 
01057                 vector< list<int> > vli_GroupedInducedVertexDegrees;
01058 
01059                 vector< list<int>::iterator > vlit_VertexLocations;
01060 
01061                 i_MaximumVertexDegree = g->GetMaximumVertexDegree();
01062 
01063         #if DEBUG == 5103
01064 
01065                 cout<<endl;
01066                 cout<<"DEBUG 5103 | Hessian Evaluation | Bicolored Sets"<<endl;
01067                 cout<<endl;
01068 
01069                 i_SetCount = (signed) vi_Sets.size();
01070 
01071                 for(i=0; i<i_SetCount; i++)
01072                 {
01073                         cout<<STEP_UP(vi_Sets[i])<<"\t"<<" : ";
01074 
01075                         i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size();
01076 
01077                         for(j=0; j<i_SetSize; j++)
01078                         {
01079                                 if(j == STEP_DOWN(i_SetSize))
01080                                 {
01081                                 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<" ("<<i_SetSize<<")"<<endl;
01082                                 }
01083                                 else
01084                                 {
01085                                 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<", ";
01086                                 }
01087                         }
01088                 }
01089 
01090                 cout<<endl;
01091                 cout<<"[Set Count = "<<i_SetCount<<"]"<<endl;
01092                 cout<<endl;
01093 
01094         #endif
01095 
01096                 //Step 5: from here on
01097                 i_VertexCount = g->GetVertexCount();
01098 
01099                 v2i_VertexAdjacency.clear();
01100                 v2i_VertexAdjacency.resize((unsigned) i_VertexCount);
01101 
01102                 v2d_NonzeroAdjacency.clear();
01103                 v2d_NonzeroAdjacency.resize((unsigned) i_VertexCount);
01104 
01105                 vi_EvaluatedDiagonals.clear();
01106                 vi_EvaluatedDiagonals.resize((unsigned) i_VertexCount, _FALSE);
01107 
01108                 vi_InducedVertexDegrees.clear();
01109                 vi_InducedVertexDegrees.resize((unsigned) i_VertexCount, _FALSE);
01110 
01111                 vd_IncludedVertices.clear();
01112                 vd_IncludedVertices.resize((unsigned) i_VertexCount, _UNKNOWN);
01113 
01114                 i_ParentVertex = _UNKNOWN;
01115 
01116                 i_SetCount = (signed) vi_Sets.size();
01117 
01118                 for(i=0; i<i_SetCount; i++)
01119                 {
01120                         vli_GroupedInducedVertexDegrees.clear();
01121                         vli_GroupedInducedVertexDegrees.resize((unsigned) STEP_UP(i_MaximumVertexDegree));
01122 
01123                         vlit_VertexLocations.clear();
01124                         vlit_VertexLocations.resize((unsigned) i_VertexCount);
01125 
01126                         i_HighestInducedVertexDegree = _UNKNOWN;
01127 
01128                         i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size();
01129 
01130                         for(j=0; j<i_SetSize; j++)
01131                         {
01132                                 i_PresentVertex = mivi_VertexSets[vi_Sets[i]][j];
01133 
01134                                 vd_IncludedVertices[i_PresentVertex] = _FALSE;
01135 
01136                                 if(vi_InducedVertexDegrees[i_PresentVertex] != _FALSE)
01137                                 {
01138                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].erase(vlit_VertexLocations[i_PresentVertex]);
01139                                 }
01140 
01141                                 vi_InducedVertexDegrees[i_PresentVertex]++;
01142 
01143                                 if(i_HighestInducedVertexDegree < vi_InducedVertexDegrees[i_PresentVertex])
01144                                 {
01145                                         i_HighestInducedVertexDegree = vi_InducedVertexDegrees[i_PresentVertex];
01146                                 }
01147                                 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].push_front(i_PresentVertex);
01148 
01149                                 vlit_VertexLocations[i_PresentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].begin();
01150                         }
01151 
01152         #if DEBUG == 5103
01153 
01154                         int k;
01155 
01156                         list<int>::iterator lit_ListIterator;
01157 
01158                         cout<<endl;
01159                         cout<<"DEBUG 5103 | Hessian Evaluation | Induced Vertex Degrees | Set "<<STEP_UP(i)<<endl;
01160                         cout<<endl;
01161 
01162                         for(j=0; j<STEP_UP(i_HighestInducedVertexDegree); j++)
01163                         {
01164                                 i_SetSize = (signed) vli_GroupedInducedVertexDegrees[j].size();
01165 
01166                                 if(i_SetSize == _FALSE)
01167                                 {
01168                                         continue;
01169                                 }
01170 
01171                                 k = _FALSE;
01172 
01173                                 cout<<"Degree "<<j<<"\t"<<" : ";
01174 
01175                                 for(lit_ListIterator=vli_GroupedInducedVertexDegrees[j].begin(); lit_ListIterator!=vli_GroupedInducedVertexDegrees[j].end(); lit_ListIterator++)
01176                                 {
01177                                         if(k == STEP_DOWN(i_SetSize))
01178                                         {
01179                                                 cout<<STEP_UP(*lit_ListIterator)<<" ("<<i_SetSize<<")"<<endl;
01180                                         }
01181                                         else
01182                                         {
01183                                                 cout<<STEP_UP(*lit_ListIterator)<<", ";
01184                                         }
01185 
01186                                         k++;
01187                                 }
01188                         }
01189 
01190         #endif
01191 
01192         #if DEBUG == 5103
01193 
01194                         cout<<endl;
01195                         cout<<"DEBUG 5103 | Hessian Evaluation | Retrieved Elements"<<"| Set "<<STEP_UP(i)<<endl;
01196                         cout<<endl;
01197 
01198         #endif
01199                         //get the diagonal values
01200                         for (int index = 0; index < i_VertexCount; index++) {
01201                                 if(vi_EvaluatedDiagonals[index] == _FALSE)
01202                                 {
01203                                         d_Value = dp2_CompressedMatrix[index][vi_VertexColors[index]];
01204 
01205         #if DEBUG == 5103
01206 
01207                                         cout<<"Element["<<STEP_UP(index)<<"]["<<STEP_UP(index)<<"] = "<<d_Value<<endl;
01208 
01209         #endif
01210                                         v2i_VertexAdjacency[index].push_back(index);
01211                                         v2d_NonzeroAdjacency[index].push_back(d_Value);
01212 
01213                                         vi_EvaluatedDiagonals[index] = _TRUE;
01214 
01215                                 }
01216                         }
01217 
01218                         for ( ; ; )
01219                         {
01220                                 if(vli_GroupedInducedVertexDegrees[_TRUE].empty()) // If there is no leaf left on the color tree
01221                                 {
01222                                         i_LeafVertex = vli_GroupedInducedVertexDegrees[_FALSE].front();
01223 
01224                                         vi_InducedVertexDegrees[i_LeafVertex] = _FALSE;
01225 
01226                                         vd_IncludedVertices[i_LeafVertex] = _UNKNOWN;
01227 
01228                                         break;
01229                                 }
01230 
01231                                 i_LeafVertex = vli_GroupedInducedVertexDegrees[_TRUE].front();
01232 
01233                                 vli_GroupedInducedVertexDegrees[_TRUE].pop_front();
01234 
01235                                 //Find i_ParentVertex
01236                                 for(j=vi_Vertices[i_LeafVertex]; j<vi_Vertices[STEP_UP(i_LeafVertex)]; j++)
01237                                 {
01238                                         if(vd_IncludedVertices[vi_Edges[j]] != _UNKNOWN)
01239                                         {
01240                                                 i_ParentVertex = vi_Edges[j];
01241 
01242                                                 break;
01243                                         }
01244                                 }
01245 
01246                                 d_Value = dp2_CompressedMatrix[i_LeafVertex][vi_VertexColors[i_ParentVertex]] - vd_IncludedVertices[i_LeafVertex];
01247 
01248                                 vd_IncludedVertices[i_ParentVertex] += d_Value;
01249 
01250                                 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE;
01251                                 vd_IncludedVertices[i_LeafVertex] = _UNKNOWN;
01252                                 if(vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].size()>1) {
01253                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].erase(vlit_VertexLocations[i_ParentVertex]);
01254                                 }
01255                                 else {
01256                                         vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].pop_back();
01257                                 }
01258 
01259                                 vi_InducedVertexDegrees[i_ParentVertex]--;
01260                                 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].push_back(i_ParentVertex);
01261 
01262                                 //Update position of the iterator pointing to i_ParentVertex in "InducedVertexDegrees" structure
01263                                 vlit_VertexLocations[i_ParentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].end();
01264                                 --vlit_VertexLocations[i_ParentVertex];
01265 
01266                                 v2i_VertexAdjacency[i_LeafVertex].push_back(i_ParentVertex);
01267                                 v2d_NonzeroAdjacency[i_LeafVertex].push_back(d_Value);
01268 
01269                                 v2i_VertexAdjacency[i_ParentVertex].push_back(i_LeafVertex);
01270                                 v2d_NonzeroAdjacency[i_ParentVertex].push_back(d_Value);
01271 
01272 
01273         #if DEBUG == 5103
01274 
01275                                 cout<<"Element["<<STEP_UP(i_LeafVertex)<<"]["<<STEP_UP(i_ParentVertex)<<"] = "<<d_Value<<endl;
01276         #endif
01277 
01278                         }
01279                 }
01280 
01281                 vector<unsigned int> RowIndex;
01282                 vector<unsigned int> ColumnIndex;
01283                 vector<double> HessianValue;
01284 
01285                 //populate dp3_HessianValue row by row, column by column
01286                 for(i=0; i<i_VertexCount; i++) {
01287                         int NumOfNonzeros = uip2_HessianSparsityPattern[i][0];
01288                         i_VertexDegree = (signed) v2i_VertexAdjacency[i].size();
01289                         for(j=1; j<=NumOfNonzeros; j++) {
01290                                 int targetColumnID = uip2_HessianSparsityPattern[i][j];
01291                                 if(targetColumnID<i) continue;
01292                                 for (int k=0; k<i_VertexDegree; k++) {// search through the v2i_VertexAdjacency matrix to find the correct column
01293                                         if(targetColumnID == v2i_VertexAdjacency[i][k]) { //found it
01294                                                 HessianValue.push_back(v2d_NonzeroAdjacency[i][k]);
01295                                                 break;
01296                                         }
01297                                 }
01298                                 RowIndex.push_back(i);
01299                                 ColumnIndex.push_back(uip2_HessianSparsityPattern[i][j]);
01300                         }
01301                 }
01302 
01303                 unsigned int numOfNonZeros = RowIndex.size();
01304                 (*ip2_RowIndex) = new unsigned int[numOfNonZeros];
01305                 (*ip2_ColumnIndex) = new unsigned int[numOfNonZeros];
01306                 (*dp2_HessianValue) = new double[numOfNonZeros]; //allocate memory for *dp2_HessianValue.
01307 
01308                 for(int i=0; i < numOfNonZeros; i++) {
01309                         (*ip2_RowIndex)[i] = RowIndex[i];
01310                         (*ip2_ColumnIndex)[i] = ColumnIndex[i];
01311                         (*dp2_HessianValue)[i] = HessianValue[i];
01312                 }
01313 
01314 
01315                 CF_available = true;
01316                 i_CF_rowCount = numOfNonZeros;
01317                 ip_CF_RowIndex = *ip2_RowIndex;
01318                 ip_CF_ColumnIndex = *ip2_ColumnIndex;
01319                 dp_CF_Value = *dp2_HessianValue;
01320 
01321                 return (numOfNonZeros);
01322         }
01323 }