|
ColPack
|
|
ColPack
|
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 //#define DEBUG 5103 00025 namespace ColPack 00026 { 00027 00028 int HessianRecovery::DirectRecover_RowCompressedFormat_usermem(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 int rowCount = g->GetVertexCount(); 00035 int colorCount = g->GetVertexColorCount(); 00036 //cout<<"colorCount="<<colorCount<<endl; 00037 vector<int> vi_VertexColors; 00038 g->GetVertexColors(vi_VertexColors); 00039 00040 //Do (column-)color statistic for each row, i.e., see how many elements in that row has color 0, color 1 ... 00041 int** colorStatistic = new int*[rowCount]; //color statistic for each row. For example, colorStatistic[0] is color statistic for row 0 00042 //If row 0 has 5 columns with color 3 => colorStatistic[0][3] = 5; 00043 //Allocate memory for colorStatistic[rowCount][colorCount] and initilize the matrix 00044 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00045 colorStatistic[i] = new int[colorCount]; 00046 for(unsigned int j=0; j < (unsigned int)colorCount; j++) colorStatistic[i][j] = 0; 00047 } 00048 00049 //populate colorStatistic 00050 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00051 int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00052 for(unsigned int j=1; j <= (unsigned int)numOfNonZeros; j++) { 00053 //non-zero in the Hessian: [i][uip2_HessianSparsityPattern[i][j]] 00054 //color of that column: vi_VertexColors[uip2_HessianSparsityPattern[i][j]] 00055 colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]++; 00056 } 00057 } 00058 00059 //Now, go to the main part, recover the values of non-zero entries in the Hessian 00060 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00061 unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00062 for(unsigned int j=1; j <= numOfNonZeros; j++) { 00063 if(i == uip2_HessianSparsityPattern[i][j]) { // the non-zero is in the diagonal of the matrix 00064 (*dp3_HessianValue)[i][j] = dp2_CompressedMatrix[i][vi_VertexColors[i]]; 00065 //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]); 00066 } 00067 else {// i != uip2_HessianSparsityPattern[i][j] // the non-zero is NOT in the diagonal of the matrix 00068 if(colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]==1) { 00069 (*dp3_HessianValue)[i][j] = dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]; 00070 //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]]); 00071 } 00072 else { 00073 (*dp3_HessianValue)[i][j] = dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]]; 00074 //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]); 00075 } 00076 } 00077 } 00078 } 00079 00080 free_2DMatrix(colorStatistic, rowCount); 00081 colorStatistic = NULL; 00082 00083 return (rowCount); 00084 } 00085 00086 int HessianRecovery::DirectRecover_RowCompressedFormat_unmanaged(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, double*** dp3_HessianValue) { 00087 if(g==NULL) { 00088 cerr<<"g==NULL"<<endl; 00089 return _FALSE; 00090 } 00091 00092 int rowCount = g->GetVertexCount(); 00093 00094 //allocate memory for *dp3_HessianValue. The dp3_HessianValue and uip2_HessianSparsityPattern matrices should have the same size 00095 *dp3_HessianValue = (double**) malloc(rowCount * sizeof(double*)); 00096 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00097 unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00098 (*dp3_HessianValue)[i] = (double*) malloc( (numOfNonZeros+1) * sizeof(double) ); 00099 (*dp3_HessianValue)[i][0] = numOfNonZeros; //initialize value of the 1st entry 00100 for(unsigned int j=1; j <= numOfNonZeros; j++) (*dp3_HessianValue)[i][j] = 0.; //initialize value of other entries 00101 } 00102 00103 return DirectRecover_RowCompressedFormat_usermem(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, dp3_HessianValue); 00104 } 00105 00106 int HessianRecovery::DirectRecover_RowCompressedFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, double*** dp3_HessianValue) { 00107 00108 int returnValue = DirectRecover_RowCompressedFormat_unmanaged(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, dp3_HessianValue); 00109 00110 if(AF_available) reset(); 00111 00112 AF_available = true; 00113 i_AF_rowCount = g->GetVertexCount(); 00114 dp2_AF_Value = *dp3_HessianValue; 00115 00116 return (returnValue); 00117 } 00118 00119 00120 int HessianRecovery::DirectRecover_SparseSolversFormat_usermem(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue, unsigned int numOfNonZerosInHessianValue) { 00121 00122 if(g==NULL) { 00123 cerr<<"g==NULL"<<endl; 00124 return _FALSE; 00125 } 00126 00127 int rowCount = g->GetVertexCount(); 00128 int colorCount = g->GetVertexColorCount(); 00129 //cout<<"colorCount="<<colorCount<<endl; 00130 vector<int> vi_VertexColors; 00131 g->GetVertexColors(vi_VertexColors); 00132 00133 //g->PrintGraph(); 00134 00135 //Making the array indices to start at 0 instead of 1 00136 for(unsigned int i=0; i <= (unsigned int) rowCount ; i++) { 00137 (*ip2_RowIndex)[i]--; 00138 } 00139 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) { 00140 (*ip2_ColumnIndex)[i]--; 00141 } 00142 00143 //Do (column-)color statistic for each row, i.e., see how many elements in that row has color 0, color 1 ... 00144 int** colorStatistic = new int*[rowCount]; //color statistic for each row. For example, colorStatistic[0] is color statistic for row 0 00145 //If row 0 has 5 columns with color 3 => colorStatistic[0][3] = 5; 00146 //Allocate memory for colorStatistic[rowCount][colorCount] and initilize the matrix 00147 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00148 colorStatistic[i] = new int[colorCount]; 00149 for(unsigned int j=0; j < (unsigned int)colorCount; j++) colorStatistic[i][j] = 0; 00150 } 00151 00152 //populate colorStatistic 00153 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00154 int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00155 for(unsigned int j=1; j <= (unsigned int)numOfNonZeros; j++) { 00156 //non-zero in the Hessian: [i][uip2_HessianSparsityPattern[i][j]] 00157 //color of that column: vi_VertexColors[uip2_HessianSparsityPattern[i][j]] 00158 colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]++; 00159 } 00160 } 00161 00162 00163 //Now, go to the main part, recover the values of non-zero entries in the Hessian 00164 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00165 unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00166 unsigned int offset = 0; 00167 //printf("\ti=%d, \t NumOfNonzeros=%d,\t \n",i,numOfNonZeros); 00168 for(unsigned int j=1; j <= numOfNonZeros; j++) { 00169 if (i > uip2_HessianSparsityPattern[i][j]) { 00170 offset++; 00171 continue; 00172 } 00173 else if(i == uip2_HessianSparsityPattern[i][j]) { // the non-zero is in the diagonal of the matrix 00174 (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = dp2_CompressedMatrix[i][vi_VertexColors[i]]; 00175 //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]); 00176 //printf("\t (*ip2_RowIndex)[i = %d] = %d, j = %d, offset = %d \n", i, (*ip2_RowIndex)[i], j, offset); 00177 } 00178 else {// i != uip2_HessianSparsityPattern[i][j] // the non-zero is NOT in the diagonal of the matrix 00179 if(colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]==1) { 00180 (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]; 00181 //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]]); 00182 } 00183 else { 00184 (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]]; 00185 //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]); 00186 } 00187 } 00188 } 00189 } 00190 00191 free_2DMatrix(colorStatistic, rowCount); 00192 colorStatistic = NULL; 00193 00194 //Making the array indices to start at 1 instead of 0 to conform with the Intel MKL sparse storage scheme for the direct sparse solvers 00195 for(unsigned int i=0; i <= (unsigned int) rowCount ; i++) { 00196 (*ip2_RowIndex)[i]++; 00197 } 00198 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) { 00199 (*ip2_ColumnIndex)[i]++; 00200 } 00201 00202 return (rowCount); 00203 } 00204 00205 int HessianRecovery::DirectRecover_SparseSolversFormat_unmanaged(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue, unsigned int numOfNonZerosInHessianValue) { 00206 00207 if(g==NULL) { 00208 cerr<<"g==NULL"<<endl; 00209 return _FALSE; 00210 } 00211 00212 int rowCount = g->GetVertexCount(); 00213 00214 if (numOfNonZerosInHessianValue < 1) { 00215 numOfNonZerosInHessianValue = ConvertRowCompressedFormat2SparseSolversFormat_StructureOnly(uip2_HessianSparsityPattern, rowCount, ip2_RowIndex, ip2_ColumnIndex); 00216 00217 //Making the array indices to start at 1 instead of 0 to conform with the Intel MKL sparse storage scheme for the direct sparse solvers 00218 for(unsigned int i=0; i <= (unsigned int) rowCount ; i++) { 00219 (*ip2_RowIndex)[i]++; 00220 } 00221 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) { 00222 (*ip2_ColumnIndex)[i]++; 00223 } 00224 } 00225 00226 //cout<<"allocate memory for *dp2_HessianValue rowCount="<<rowCount<<endl; 00227 //printf("i=%d\t numOfNonZerosInHessianValue=%d \n", i, numOfNonZerosInHessianValue); 00228 (*dp2_HessianValue) = (double*) malloc(numOfNonZerosInHessianValue * sizeof(double)); //allocate memory for *dp2_JacobianValue. 00229 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) (*dp2_HessianValue)[i] = 0.; //initialize value of other entries 00230 00231 int returnValue = DirectRecover_SparseSolversFormat_usermem(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue, numOfNonZerosInHessianValue); 00232 00233 return returnValue; 00234 } 00235 00236 int HessianRecovery::DirectRecover_SparseSolversFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 00237 int returnValue = DirectRecover_SparseSolversFormat_unmanaged(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue); 00238 00239 if(SSF_available) { 00240 //cout<<"SSF_available="<<SSF_available<<endl; Pause(); 00241 reset(); 00242 } 00243 00244 SSF_available = true; 00245 i_SSF_rowCount = g->GetVertexCount(); 00246 ip_SSF_RowIndex = *ip2_RowIndex; 00247 ip_SSF_ColumnIndex = *ip2_ColumnIndex; 00248 dp_SSF_Value = *dp2_HessianValue; 00249 00250 return (returnValue); 00251 } 00252 00253 int HessianRecovery::DirectRecover_CoordinateFormat_vectors_OMP(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, vector<unsigned int> &RowIndex, vector<unsigned int> &ColumnIndex, vector<double> &HessianValue) { 00254 int rowCount = g->GetVertexCount(); 00255 int colorCount = g->GetVertexColorCount(); 00256 //cout<<"colorCount="<<colorCount<<endl; 00257 vector<int> vi_VertexColors; 00258 g->GetVertexColors(vi_VertexColors); 00259 00260 //Do (column-)color statistic for each row, i.e., see how many elements in that row has color 0, color 1 ... 00261 int** colorStatistic = new int*[rowCount]; //color statistic for each row. For example, colorStatistic[0] is color statistic for row 0 00262 //If row 0 has 5 columns with color 3 => colorStatistic[0][3] = 5; 00263 //Allocate memory for colorStatistic[rowCount][colorCount] and initilize the matrix 00264 #pragma omp parallel for default(none) schedule(static) shared(colorStatistic, colorCount, rowCount) 00265 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00266 colorStatistic[i] = new int[colorCount]; 00267 for(unsigned int j=0; j < (unsigned int)colorCount; j++) colorStatistic[i][j] = 0; 00268 } 00269 00270 //populate colorStatistic 00271 #pragma omp parallel for default(none) schedule(static) shared(rowCount, uip2_HessianSparsityPattern, colorStatistic, vi_VertexColors) 00272 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00273 int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00274 for(unsigned int j=1; j <= (unsigned int)numOfNonZeros; j++) { 00275 //non-zero in the Hessian: [i][uip2_HessianSparsityPattern[i][j]] 00276 //color of that column: vi_VertexColors[uip2_HessianSparsityPattern[i][j]] 00277 colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]++; 00278 } 00279 } 00280 00281 int* row_size = new int[rowCount]; 00282 int total_length; 00283 row_size[0]=0; 00284 00285 #pragma omp parallel for default(none) schedule(static) shared(rowCount,row_size, uip2_HessianSparsityPattern) reduction(+:total_length) 00286 // Find out the number of elements in each rows. Note: row_size[1] will has the number of elements in the first row (!!!NOT the 2nd row) 00287 for(unsigned int i=0; i < (unsigned int)rowCount-1; i++) { 00288 unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00289 int rsize=0; 00290 for(unsigned int j=1; j <= numOfNonZeros; j++) { 00291 if(uip2_HessianSparsityPattern[i][j]<i) continue; 00292 rsize++; 00293 } 00294 total_length = total_length + rsize; 00295 row_size[i+1]=rsize; 00296 } 00297 for(unsigned int j=1; j <= uip2_HessianSparsityPattern[rowCount-1][0]; j++) { 00298 if(uip2_HessianSparsityPattern[rowCount-1][j]<rowCount-1) continue; 00299 total_length++; 00300 } 00301 00302 // Allocate memory for RowIndex, ColumnIndex and HessianValue. Also modify row_size 00303 #pragma omp sections 00304 { 00305 #pragma omp section 00306 { 00307 RowIndex.resize(total_length); 00308 } 00309 #pragma omp section 00310 { 00311 ColumnIndex.resize(total_length); 00312 } 00313 #pragma omp section 00314 { 00315 HessianValue.resize(total_length); 00316 } 00317 #pragma omp section 00318 { 00319 //Modify row_size so that each element will keep the starting index of RowIndex, ColumnIndex and HessianValue 00320 for(unsigned int i=1; i < (unsigned int)rowCount; i++) { 00321 row_size[i] += row_size[i-1]; 00322 } 00323 } 00324 } 00325 00326 //Now, go to the main part, recover the values of non-zero entries in the Hessian 00328 #pragma omp parallel for default(none) schedule(static) shared(rowCount, uip2_HessianSparsityPattern, HessianValue, RowIndex, ColumnIndex, dp2_CompressedMatrix, vi_VertexColors, colorStatistic, row_size) 00329 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00330 unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00331 int index = 0; 00332 for(unsigned int j=1; j <= numOfNonZeros; j++) { 00333 if(uip2_HessianSparsityPattern[i][j]<i) continue; 00334 00335 if(i == uip2_HessianSparsityPattern[i][j]) { // the non-zero is in the diagonal of the matrix 00336 HessianValue.push_back(dp2_CompressedMatrix[i][vi_VertexColors[i]]); 00337 } 00338 else {// i != uip2_HessianSparsityPattern[i][j] // the non-zero is NOT in the diagonal of the matrix 00339 if(colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]==1) { 00340 HessianValue[row_size[i]+index] =dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]; 00341 } 00342 else { 00343 HessianValue[row_size[i]+index] =dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]]; 00344 } 00345 } 00346 RowIndex[row_size[i]+index] = i; 00347 ColumnIndex[row_size[i]+index] = uip2_HessianSparsityPattern[i][j]; 00348 index++; 00349 } 00350 } 00351 00352 free_2DMatrix(colorStatistic, rowCount); 00353 colorStatistic = NULL; 00354 00355 return (RowIndex.size()); 00356 } 00357 00358 int HessianRecovery::DirectRecover_CoordinateFormat_vectors(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, vector<unsigned int> &RowIndex, vector<unsigned int> &ColumnIndex, vector<double> &HessianValue) { 00359 int rowCount = g->GetVertexCount(); 00360 int colorCount = g->GetVertexColorCount(); 00361 //cout<<"colorCount="<<colorCount<<endl; 00362 vector<int> vi_VertexColors; 00363 g->GetVertexColors(vi_VertexColors); 00364 00365 //Do (column-)color statistic for each row, i.e., see how many elements in that row has color 0, color 1 ... 00366 int** colorStatistic = new int*[rowCount]; //color statistic for each row. For example, colorStatistic[0] is color statistic for row 0 00367 //If row 0 has 5 columns with color 3 => colorStatistic[0][3] = 5; 00368 //Allocate memory for colorStatistic[rowCount][colorCount] and initilize the matrix 00369 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00370 colorStatistic[i] = new int[colorCount]; 00371 for(unsigned int j=0; j < (unsigned int)colorCount; j++) colorStatistic[i][j] = 0; 00372 } 00373 00374 //populate colorStatistic 00375 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00376 int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00377 for(unsigned int j=1; j <= (unsigned int)numOfNonZeros; j++) { 00378 //non-zero in the Hessian: [i][uip2_HessianSparsityPattern[i][j]] 00379 //color of that column: vi_VertexColors[uip2_HessianSparsityPattern[i][j]] 00380 colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]++; 00381 } 00382 } 00383 00384 //Now, go to the main part, recover the values of non-zero entries in the Hessian 00385 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00386 unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00387 for(unsigned int j=1; j <= numOfNonZeros; j++) { 00388 if(uip2_HessianSparsityPattern[i][j]<i) continue; 00389 00390 if(i == uip2_HessianSparsityPattern[i][j]) { // the non-zero is in the diagonal of the matrix 00391 HessianValue.push_back(dp2_CompressedMatrix[i][vi_VertexColors[i]]); 00392 } 00393 else {// i != uip2_HessianSparsityPattern[i][j] // the non-zero is NOT in the diagonal of the matrix 00394 if(colorStatistic[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]==1) { 00395 HessianValue.push_back(dp2_CompressedMatrix[i][vi_VertexColors[uip2_HessianSparsityPattern[i][j]]]); 00396 } 00397 else { 00398 HessianValue.push_back(dp2_CompressedMatrix[uip2_HessianSparsityPattern[i][j]][vi_VertexColors[i]]); 00399 } 00400 } 00401 RowIndex.push_back(i); 00402 ColumnIndex.push_back(uip2_HessianSparsityPattern[i][j]); 00403 } 00404 } 00405 00406 free_2DMatrix(colorStatistic, rowCount); 00407 colorStatistic = NULL; 00408 00409 return (RowIndex.size()); 00410 } 00411 00412 int HessianRecovery::DirectRecover_CoordinateFormat_usermem(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 00413 if(g==NULL) { 00414 cerr<<"g==NULL"<<endl; 00415 return _FALSE; 00416 } 00417 00418 vector<unsigned int> RowIndex; 00419 vector<unsigned int> ColumnIndex; 00420 vector<double> HessianValue; 00421 00422 int returnValue = DirectRecover_CoordinateFormat_vectors_OMP(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, RowIndex, ColumnIndex, HessianValue); 00423 00424 unsigned int numOfNonZeros = RowIndex.size(); 00425 00426 #pragma omp parallel for default(none) schedule(static) shared(numOfNonZeros, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue, RowIndex, ColumnIndex, HessianValue) 00427 for(int i=0; i < numOfNonZeros; i++) { 00428 (*ip2_RowIndex)[i] = RowIndex[i]; 00429 (*ip2_ColumnIndex)[i] = ColumnIndex[i]; 00430 (*dp2_HessianValue)[i] = HessianValue[i]; 00431 } 00432 00433 return (returnValue); 00434 } 00435 00436 int HessianRecovery::DirectRecover_CoordinateFormat_usermem_serial(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 00437 if(g==NULL) { 00438 cerr<<"g==NULL"<<endl; 00439 return _FALSE; 00440 } 00441 00442 vector<unsigned int> RowIndex; 00443 vector<unsigned int> ColumnIndex; 00444 vector<double> HessianValue; 00445 00446 int returnValue = DirectRecover_CoordinateFormat_vectors(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, RowIndex, ColumnIndex, HessianValue); 00447 00448 unsigned int numOfNonZeros = RowIndex.size(); 00449 00450 for(int i=0; i < numOfNonZeros; i++) { 00451 (*ip2_RowIndex)[i] = RowIndex[i]; 00452 (*ip2_ColumnIndex)[i] = ColumnIndex[i]; 00453 (*dp2_HessianValue)[i] = HessianValue[i]; 00454 } 00455 00456 return (returnValue); 00457 } 00458 00459 int HessianRecovery::DirectRecover_CoordinateFormat_unmanaged_OMP(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 00460 if(g==NULL) { 00461 cerr<<"g==NULL"<<endl; 00462 return _FALSE; 00463 } 00464 00465 vector<unsigned int> RowIndex; 00466 vector<unsigned int> ColumnIndex; 00467 vector<double> HessianValue; 00468 00469 int returnValue = DirectRecover_CoordinateFormat_vectors_OMP(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, RowIndex, ColumnIndex, HessianValue); 00470 00471 unsigned int numOfNonZeros = RowIndex.size(); 00472 #pragma omp sections 00473 { 00474 #pragma omp section 00475 { 00476 (*ip2_RowIndex) = (unsigned int*) malloc(numOfNonZeros * sizeof(unsigned int)); 00477 } 00478 #pragma omp section 00479 { 00480 (*ip2_ColumnIndex) = (unsigned int*) malloc(numOfNonZeros * sizeof(unsigned int)); 00481 } 00482 #pragma omp section 00483 { 00484 (*dp2_HessianValue) = (double*) malloc(numOfNonZeros * sizeof(double)); //allocate memory for *dp2_HessianValue. 00485 } 00486 } 00487 00488 #pragma omp parallel for default(none) schedule(static) shared(numOfNonZeros, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue, RowIndex, ColumnIndex, HessianValue) 00489 for(int i=0; i < numOfNonZeros; i++) { 00490 (*ip2_RowIndex)[i] = RowIndex[i]; 00491 (*ip2_ColumnIndex)[i] = ColumnIndex[i]; 00492 (*dp2_HessianValue)[i] = HessianValue[i]; 00493 } 00494 00495 return (returnValue); 00496 } 00497 00498 int HessianRecovery::DirectRecover_CoordinateFormat_unmanaged(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 00499 if(g==NULL) { 00500 cerr<<"g==NULL"<<endl; 00501 return _FALSE; 00502 } 00503 00504 vector<unsigned int> RowIndex; 00505 vector<unsigned int> ColumnIndex; 00506 vector<double> HessianValue; 00507 00508 int returnValue = DirectRecover_CoordinateFormat_vectors(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, RowIndex, ColumnIndex, HessianValue); 00509 00510 unsigned int numOfNonZeros = RowIndex.size(); 00511 (*ip2_RowIndex) = (unsigned int*) malloc(numOfNonZeros * sizeof(unsigned int)); 00512 (*ip2_ColumnIndex) = (unsigned int*) malloc(numOfNonZeros * sizeof(unsigned int)); 00513 (*dp2_HessianValue) = (double*) malloc(numOfNonZeros * sizeof(double)); //allocate memory for *dp2_HessianValue. 00514 00515 for(int i=0; i < numOfNonZeros; i++) { 00516 (*ip2_RowIndex)[i] = RowIndex[i]; 00517 (*ip2_ColumnIndex)[i] = ColumnIndex[i]; 00518 (*dp2_HessianValue)[i] = HessianValue[i]; 00519 } 00520 00521 return (returnValue); 00522 } 00523 00524 int HessianRecovery::DirectRecover_CoordinateFormat_OMP(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 00525 int numOfNonZeros = DirectRecover_CoordinateFormat_unmanaged_OMP(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, ip2_RowIndex,ip2_ColumnIndex, dp2_HessianValue); 00526 00527 if(CF_available) reset(); 00528 00529 CF_available = true; 00530 i_CF_rowCount = g->GetVertexCount(); 00531 ip_CF_RowIndex = *ip2_RowIndex; 00532 ip_CF_ColumnIndex = *ip2_ColumnIndex; 00533 dp_CF_Value = *dp2_HessianValue; 00534 00535 return (numOfNonZeros); 00536 } 00537 00538 int HessianRecovery::DirectRecover_CoordinateFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 00539 int numOfNonZeros = DirectRecover_CoordinateFormat_unmanaged(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, ip2_RowIndex,ip2_ColumnIndex, dp2_HessianValue); 00540 00541 if(CF_available) reset(); 00542 00543 CF_available = true; 00544 i_CF_rowCount = g->GetVertexCount(); 00545 ip_CF_RowIndex = *ip2_RowIndex; 00546 ip_CF_ColumnIndex = *ip2_ColumnIndex; 00547 dp_CF_Value = *dp2_HessianValue; 00548 00549 return (numOfNonZeros); 00550 } 00551 00552 int HessianRecovery::IndirectRecover_RowCompressedFormat_usermem(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, double*** dp3_HessianValue) { 00553 if(g==NULL) { 00554 cerr<<"g==NULL"<<endl; 00555 return _FALSE; 00556 } 00557 00558 int i=0,j=0; 00559 int i_VertexCount = _UNKNOWN; 00560 int i_EdgeID, i_SetID; 00561 vector<int> vi_Sets; 00562 map< int, vector<int> > mivi_VertexSets; 00563 00564 vector<int> vi_Vertices; 00565 g->GetVertices(vi_Vertices); 00566 00567 vector<int> vi_Edges; 00568 g->GetEdges(vi_Edges); 00569 00570 vector<int> vi_VertexColors; 00571 g->GetVertexColors(vi_VertexColors); 00572 00573 map< int, map< int, int> > mimi2_VertexEdgeMap; 00574 g->GetVertexEdgeMap(mimi2_VertexEdgeMap); 00575 00576 DisjointSets ds_DisjointSets; 00577 g->GetDisjointSets(ds_DisjointSets); 00578 00579 //populate vi_Sets & mivi_VertexSets 00580 vi_Sets.clear(); 00581 mivi_VertexSets.clear(); 00582 00583 i_VertexCount = g->GetVertexCount(); 00584 00585 for(i=0; i<i_VertexCount; i++) // for each vertex A (indexed by i) 00586 { 00587 for(j=vi_Vertices[i]; j<vi_Vertices[STEP_UP(i)]; j++) // for each of the vertex B that connect to A 00588 { 00589 if(i < vi_Edges[j]) // if the index of A (i) is less than the index of B (vi_Edges[j]) 00590 //basicly each edge is represented by (vertex with smaller ID, vertex with larger ID). This way, we don't insert a specific edge twice 00591 { 00592 i_EdgeID = mimi2_VertexEdgeMap[i][vi_Edges[j]]; 00593 00594 i_SetID = ds_DisjointSets.FindAndCompress(i_EdgeID); 00595 00596 if(i_SetID == i_EdgeID) // that edge is the root of the set => create new set 00597 { 00598 vi_Sets.push_back(i_SetID); 00599 } 00600 00601 mivi_VertexSets[i_SetID].push_back(i); 00602 mivi_VertexSets[i_SetID].push_back(vi_Edges[j]); 00603 } 00604 } 00605 } 00606 00607 int i_MaximumVertexDegree; 00608 00609 int i_HighestInducedVertexDegree; 00610 00611 int i_LeafVertex, i_ParentVertex, i_PresentVertex; 00612 00613 int i_VertexDegree; 00614 00615 int i_SetCount, i_SetSize; 00616 //i_SetCount = vi_Sets.size(); 00617 //i_SetSize: size (number of edges?) in a bicolored tree 00618 00619 double d_Value; 00620 00621 vector<int> vi_EvaluatedDiagonals; 00622 00623 vector<int> vi_InducedVertexDegrees; 00624 00625 vector<double> vd_IncludedVertices; 00626 00627 vector< vector<int> > v2i_VertexAdjacency; 00628 00629 vector< vector<double> > v2d_NonzeroAdjacency; 00630 00631 vector< list<int> > vli_GroupedInducedVertexDegrees; 00632 00633 vector< list<int>::iterator > vlit_VertexLocations; 00634 00635 i_MaximumVertexDegree = g->GetMaximumVertexDegree(); 00636 00637 #if DEBUG == 5103 00638 00639 cout<<endl; 00640 cout<<"DEBUG 5103 | Hessian Evaluation | Bicolored Sets"<<endl; 00641 cout<<endl; 00642 00643 i_SetCount = (signed) vi_Sets.size(); 00644 00645 for(i=0; i<i_SetCount; i++) 00646 { 00647 cout<<STEP_UP(vi_Sets[i])<<"\t"<<" : "; 00648 00649 i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size(); 00650 00651 for(j=0; j<i_SetSize; j++) 00652 { 00653 if(j == STEP_DOWN(i_SetSize)) 00654 { 00655 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<" ("<<i_SetSize<<")"<<endl; 00656 } 00657 else 00658 { 00659 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<", "; 00660 } 00661 } 00662 } 00663 00664 cout<<endl; 00665 cout<<"[Set Count = "<<i_SetCount<<"]"<<endl; 00666 cout<<endl; 00667 00668 #endif 00669 00670 //Step 5: from here on 00671 i_VertexCount = g->GetVertexCount(); 00672 00673 v2i_VertexAdjacency.clear(); 00674 v2i_VertexAdjacency.resize((unsigned) i_VertexCount); 00675 00676 v2d_NonzeroAdjacency.clear(); 00677 v2d_NonzeroAdjacency.resize((unsigned) i_VertexCount); 00678 00679 vi_EvaluatedDiagonals.clear(); 00680 vi_EvaluatedDiagonals.resize((unsigned) i_VertexCount, _FALSE); 00681 00682 vi_InducedVertexDegrees.clear(); 00683 vi_InducedVertexDegrees.resize((unsigned) i_VertexCount, _FALSE); 00684 00685 vd_IncludedVertices.clear(); 00686 vd_IncludedVertices.resize((unsigned) i_VertexCount, _UNKNOWN); 00687 00688 i_ParentVertex = _UNKNOWN; 00689 00690 i_SetCount = (signed) vi_Sets.size(); 00691 00692 for(i=0; i<i_SetCount; i++) 00693 { 00694 vli_GroupedInducedVertexDegrees.clear(); 00695 vli_GroupedInducedVertexDegrees.resize((unsigned) STEP_UP(i_MaximumVertexDegree)); 00696 00697 vlit_VertexLocations.clear(); 00698 vlit_VertexLocations.resize((unsigned) i_VertexCount); 00699 00700 i_HighestInducedVertexDegree = _UNKNOWN; 00701 00702 i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size(); 00703 00704 for(j=0; j<i_SetSize; j++) 00705 { 00706 i_PresentVertex = mivi_VertexSets[vi_Sets[i]][j]; 00707 00708 vd_IncludedVertices[i_PresentVertex] = _FALSE; 00709 00710 if(vi_InducedVertexDegrees[i_PresentVertex] != _FALSE) 00711 { 00712 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].erase(vlit_VertexLocations[i_PresentVertex]); 00713 } 00714 00715 vi_InducedVertexDegrees[i_PresentVertex]++; 00716 00717 if(i_HighestInducedVertexDegree < vi_InducedVertexDegrees[i_PresentVertex]) 00718 { 00719 i_HighestInducedVertexDegree = vi_InducedVertexDegrees[i_PresentVertex]; 00720 } 00721 00722 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].push_front(i_PresentVertex); 00723 00724 vlit_VertexLocations[i_PresentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].begin(); 00725 } 00726 00727 #if DEBUG == 5103 00728 00729 int k; 00730 00731 list<int>::iterator lit_ListIterator; 00732 00733 cout<<endl; 00734 cout<<"DEBUG 5103 | Hessian Evaluation | Induced Vertex Degrees | Set "<<STEP_UP(i)<<endl; 00735 cout<<endl; 00736 00737 for(j=0; j<STEP_UP(i_HighestInducedVertexDegree); j++) 00738 { 00739 i_SetSize = (signed) vli_GroupedInducedVertexDegrees[j].size(); 00740 00741 if(i_SetSize == _FALSE) 00742 { 00743 continue; 00744 } 00745 00746 k = _FALSE; 00747 00748 cout<<"Degree "<<j<<"\t"<<" : "; 00749 00750 for(lit_ListIterator=vli_GroupedInducedVertexDegrees[j].begin(); lit_ListIterator!=vli_GroupedInducedVertexDegrees[j].end(); lit_ListIterator++) 00751 { 00752 if(k == STEP_DOWN(i_SetSize)) 00753 { 00754 cout<<STEP_UP(*lit_ListIterator)<<"{"<<vd_IncludedVertices[*lit_ListIterator]<<"} ("<<i_SetSize<<")"<<endl; 00755 } 00756 else 00757 { 00758 cout<<STEP_UP(*lit_ListIterator)<<"{"<<vd_IncludedVertices[*lit_ListIterator]<<"}, "; 00759 } 00760 00761 k++; 00762 } 00763 } 00764 00765 #endif 00766 00767 #if DEBUG == 5103 00768 00769 cout<<endl; 00770 cout<<"DEBUG 5103 | Hessian Evaluation | Retrieved Elements"<<"| Set "<<STEP_UP(i)<<endl; 00771 cout<<endl; 00772 00773 #endif 00774 //#define DEBUG 5103 00775 //get the diagonal values 00776 for (int index = 0; index < i_VertexCount; index++) { 00777 if(vi_EvaluatedDiagonals[index] == _FALSE) 00778 { 00779 d_Value = dp2_CompressedMatrix[index][vi_VertexColors[index]]; 00780 00781 #if DEBUG == 5103 00782 00783 cout<<"Element["<<STEP_UP(index)<<"]["<<STEP_UP(index)<<"] = "<<d_Value<<endl; 00784 00785 #endif 00786 v2i_VertexAdjacency[index].push_back(index); 00787 v2d_NonzeroAdjacency[index].push_back(d_Value); 00788 00789 vi_EvaluatedDiagonals[index] = _TRUE; 00790 00791 } 00792 } 00793 00794 for ( ; ; ) 00795 { 00796 if(vli_GroupedInducedVertexDegrees[_TRUE].empty()) // If there is no leaf left on the color tree 00797 { 00798 i_LeafVertex = vli_GroupedInducedVertexDegrees[_FALSE].front(); 00799 00800 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE; 00801 00802 vd_IncludedVertices[i_LeafVertex] = _UNKNOWN; 00803 00804 break; 00805 } 00806 00807 i_LeafVertex = vli_GroupedInducedVertexDegrees[_TRUE].front(); 00808 00809 vli_GroupedInducedVertexDegrees[_TRUE].pop_front(); 00810 00811 00812 //Find i_ParentVertex 00813 for(j=vi_Vertices[i_LeafVertex]; j<vi_Vertices[STEP_UP(i_LeafVertex)]; j++) 00814 { 00815 if(vd_IncludedVertices[vi_Edges[j]] != _UNKNOWN) 00816 { 00817 i_ParentVertex = vi_Edges[j]; 00818 00819 break; 00820 } 00821 } 00822 00823 d_Value = dp2_CompressedMatrix[i_LeafVertex][vi_VertexColors[i_ParentVertex]] - vd_IncludedVertices[i_LeafVertex]; 00824 00825 vd_IncludedVertices[i_ParentVertex] += d_Value; 00826 00827 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE; 00828 vd_IncludedVertices[i_LeafVertex] = _UNKNOWN; 00829 if(vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].size()>1) { 00830 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].erase(vlit_VertexLocations[i_ParentVertex]); 00831 } 00832 else { 00833 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].pop_back(); 00834 } 00835 00836 vi_InducedVertexDegrees[i_ParentVertex]--; 00837 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].push_back(i_ParentVertex); 00838 00839 //Update position of the iterator pointing to i_ParentVertex in "InducedVertexDegrees" structure 00840 vlit_VertexLocations[i_ParentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].end(); 00841 --vlit_VertexLocations[i_ParentVertex]; 00842 00843 v2i_VertexAdjacency[i_LeafVertex].push_back(i_ParentVertex); 00844 v2d_NonzeroAdjacency[i_LeafVertex].push_back(d_Value); 00845 00846 v2i_VertexAdjacency[i_ParentVertex].push_back(i_LeafVertex); 00847 v2d_NonzeroAdjacency[i_ParentVertex].push_back(d_Value); 00848 00849 00850 #if DEBUG == 5103 00851 00852 cout<<"Element["<<STEP_UP(i_LeafVertex)<<"]["<<STEP_UP(i_ParentVertex)<<"] = "<<d_Value<<endl; 00853 #endif 00854 00855 } 00856 } 00857 00858 00859 //allocate memory for *dp3_HessianValue. The dp3_HessianValue and uip2_HessianSparsityPattern matrices should have the same size 00860 //*dp3_HessianValue = new double*[i_VertexCount]; 00861 //for(unsigned int i=0; i < (unsigned int)i_VertexCount; i++) { 00862 // unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00863 // (*dp3_HessianValue)[i] = new double[numOfNonZeros+1]; 00864 // (*dp3_HessianValue)[i][0] = numOfNonZeros; //initialize value of the 1st entry 00865 // for(unsigned int j=1; j <= numOfNonZeros; j++) (*dp3_HessianValue)[i][j] = 0.; //initialize value of other entries 00866 //} 00867 00868 //populate dp3_HessianValue row by row, column by column 00869 for(i=0; i<i_VertexCount; i++) { 00870 int NumOfNonzeros = uip2_HessianSparsityPattern[i][0]; 00871 i_VertexDegree = (signed) v2i_VertexAdjacency[i].size(); 00872 for(j=1; j<=NumOfNonzeros; j++) { 00873 int targetColumnID = uip2_HessianSparsityPattern[i][j]; 00874 for (int k=0; k<i_VertexDegree; k++) {// search through the v2i_VertexAdjacency matrix to find the correct column 00875 if(targetColumnID == v2i_VertexAdjacency[i][k]) { //found it 00876 (*dp3_HessianValue)[i][j] = v2d_NonzeroAdjacency[i][k]; 00877 break; 00878 } 00879 } 00880 } 00881 } 00882 00883 00884 return (i_VertexCount); 00885 } 00886 00887 int HessianRecovery::IndirectRecover_RowCompressedFormat_unmanaged(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, double*** dp3_HessianValue) { 00888 if(g==NULL) { 00889 cerr<<"g==NULL"<<endl; 00890 return _FALSE; 00891 } 00892 00893 int rowCount = g->GetVertexCount(); 00894 00895 //allocate memory for *dp3_HessianValue. The dp3_HessianValue and uip2_HessianSparsityPattern matrices should have the same size 00896 *dp3_HessianValue = (double**) malloc(rowCount * sizeof(double*)); 00897 for(unsigned int i=0; i < (unsigned int)rowCount; i++) { 00898 unsigned int numOfNonZeros = uip2_HessianSparsityPattern[i][0]; 00899 (*dp3_HessianValue)[i] = (double*) malloc( (numOfNonZeros+1) * sizeof(double) ); 00900 (*dp3_HessianValue)[i][0] = numOfNonZeros; //initialize value of the 1st entry 00901 for(unsigned int j=1; j <= numOfNonZeros; j++) (*dp3_HessianValue)[i][j] = 0.; //initialize value of other entries 00902 } 00903 00904 return IndirectRecover_RowCompressedFormat_usermem(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, dp3_HessianValue); 00905 } 00906 00907 int HessianRecovery::IndirectRecover_RowCompressedFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, double*** dp3_HessianValue) { 00908 00909 int returnValue = IndirectRecover_RowCompressedFormat_unmanaged( g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, dp3_HessianValue); 00910 00911 if(AF_available) reset(); 00912 00913 AF_available = true; 00914 i_AF_rowCount = g->GetVertexCount(); 00915 dp2_AF_Value = *dp3_HessianValue; 00916 00917 return (returnValue); 00918 } 00919 00920 int HessianRecovery::IndirectRecover_SparseSolversFormat_usermem(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue, unsigned int numOfNonZerosInHessianValue) { 00921 00922 if(g==NULL) { 00923 cerr<<"g==NULL"<<endl; 00924 return _FALSE; 00925 } 00926 00927 //unsigned int numOfNonZerosInHessianValue = ConvertRowCompressedFormat2SparseSolversFormat_StructureOnly(uip2_HessianSparsityPattern, g->GetVertexCount(), ip2_RowIndex, ip2_ColumnIndex); 00928 00929 int i_VertexCount = g->GetVertexCount(); 00930 00931 //Making the array indices to start at 0 instead of 1 00932 for(unsigned int i=0; i <= (unsigned int) i_VertexCount ; i++) { 00933 (*ip2_RowIndex)[i]--; 00934 } 00935 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) { 00936 (*ip2_ColumnIndex)[i]--; 00937 } 00938 00939 int i=0,j=0; 00940 int i_EdgeID, i_SetID; 00941 vector<int> vi_Sets; 00942 map< int, vector<int> > mivi_VertexSets; 00943 00944 vector<int> vi_Vertices; 00945 g->GetVertices(vi_Vertices); 00946 00947 vector<int> vi_Edges; 00948 g->GetEdges(vi_Edges); 00949 00950 vector<int> vi_VertexColors; 00951 g->GetVertexColors(vi_VertexColors); 00952 00953 map< int, map< int, int> > mimi2_VertexEdgeMap; 00954 g->GetVertexEdgeMap(mimi2_VertexEdgeMap); 00955 00956 DisjointSets ds_DisjointSets; 00957 g->GetDisjointSets(ds_DisjointSets); 00958 00959 //populate vi_Sets & mivi_VertexSets 00960 vi_Sets.clear(); 00961 mivi_VertexSets.clear(); 00962 00963 00964 for(i=0; i<i_VertexCount; i++) // for each vertex A (indexed by i) 00965 { 00966 for(j=vi_Vertices[i]; j<vi_Vertices[STEP_UP(i)]; j++) // for each of the vertex B that connect to A 00967 { 00968 if(i < vi_Edges[j]) // if the index of A (i) is less than the index of B (vi_Edges[j]) 00969 //basic each edge is represented by (vertex with smaller ID, vertex with larger ID). This way, we don't insert a specific edge twice 00970 { 00971 i_EdgeID = mimi2_VertexEdgeMap[i][vi_Edges[j]]; 00972 00973 i_SetID = ds_DisjointSets.FindAndCompress(i_EdgeID); 00974 00975 if(i_SetID == i_EdgeID) // that edge is the root of the set => create new set 00976 { 00977 vi_Sets.push_back(i_SetID); 00978 } 00979 00980 mivi_VertexSets[i_SetID].push_back(i); 00981 mivi_VertexSets[i_SetID].push_back(vi_Edges[j]); 00982 } 00983 } 00984 } 00985 00986 int i_MaximumVertexDegree; 00987 00988 int i_HighestInducedVertexDegree; 00989 00990 int i_LeafVertex, i_ParentVertex, i_PresentVertex; 00991 00992 int i_VertexDegree; 00993 00994 int i_SetCount, i_SetSize; 00995 00996 double d_Value; 00997 00998 vector<int> vi_EvaluatedDiagonals; 00999 01000 vector<int> vi_InducedVertexDegrees; 01001 01002 vector<double> vd_IncludedVertices; 01003 01004 vector< vector<int> > v2i_VertexAdjacency; 01005 01006 vector< vector<double> > v2d_NonzeroAdjacency; 01007 01008 vector< list<int> > vli_GroupedInducedVertexDegrees; 01009 01010 vector< list<int>::iterator > vlit_VertexLocations; 01011 01012 i_MaximumVertexDegree = g->GetMaximumVertexDegree(); 01013 01014 #if DEBUG == 5103 01015 01016 cout<<endl; 01017 cout<<"DEBUG 5103 | Hessian Evaluation | Bicolored Sets"<<endl; 01018 cout<<endl; 01019 01020 i_SetCount = (signed) vi_Sets.size(); 01021 01022 for(i=0; i<i_SetCount; i++) 01023 { 01024 cout<<STEP_UP(vi_Sets[i])<<"\t"<<" : "; 01025 01026 i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size(); 01027 01028 for(j=0; j<i_SetSize; j++) 01029 { 01030 if(j == STEP_DOWN(i_SetSize)) 01031 { 01032 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<" ("<<i_SetSize<<")"<<endl; 01033 } 01034 else 01035 { 01036 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<", "; 01037 } 01038 } 01039 } 01040 01041 cout<<endl; 01042 cout<<"[Set Count = "<<i_SetCount<<"]"<<endl; 01043 cout<<endl; 01044 01045 #endif 01046 01047 //Step 5: from here on 01048 i_VertexCount = g->GetVertexCount(); 01049 01050 v2i_VertexAdjacency.clear(); 01051 v2i_VertexAdjacency.resize((unsigned) i_VertexCount); 01052 01053 v2d_NonzeroAdjacency.clear(); 01054 v2d_NonzeroAdjacency.resize((unsigned) i_VertexCount); 01055 01056 vi_EvaluatedDiagonals.clear(); 01057 vi_EvaluatedDiagonals.resize((unsigned) i_VertexCount, _FALSE); 01058 01059 vi_InducedVertexDegrees.clear(); 01060 vi_InducedVertexDegrees.resize((unsigned) i_VertexCount, _FALSE); 01061 01062 vd_IncludedVertices.clear(); 01063 vd_IncludedVertices.resize((unsigned) i_VertexCount, _UNKNOWN); 01064 01065 i_ParentVertex = _UNKNOWN; 01066 01067 i_SetCount = (signed) vi_Sets.size(); 01068 01069 for(i=0; i<i_SetCount; i++) 01070 { 01071 vli_GroupedInducedVertexDegrees.clear(); 01072 vli_GroupedInducedVertexDegrees.resize((unsigned) STEP_UP(i_MaximumVertexDegree)); 01073 01074 vlit_VertexLocations.clear(); 01075 vlit_VertexLocations.resize((unsigned) i_VertexCount); 01076 01077 i_HighestInducedVertexDegree = _UNKNOWN; 01078 01079 i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size(); 01080 01081 for(j=0; j<i_SetSize; j++) 01082 { 01083 i_PresentVertex = mivi_VertexSets[vi_Sets[i]][j]; 01084 01085 vd_IncludedVertices[i_PresentVertex] = _FALSE; 01086 01087 if(vi_InducedVertexDegrees[i_PresentVertex] != _FALSE) 01088 { 01089 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].erase(vlit_VertexLocations[i_PresentVertex]); 01090 } 01091 01092 vi_InducedVertexDegrees[i_PresentVertex]++; 01093 01094 if(i_HighestInducedVertexDegree < vi_InducedVertexDegrees[i_PresentVertex]) 01095 { 01096 i_HighestInducedVertexDegree = vi_InducedVertexDegrees[i_PresentVertex]; 01097 } 01098 01099 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].push_front(i_PresentVertex); 01100 01101 vlit_VertexLocations[i_PresentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].begin(); 01102 } 01103 01104 #if DEBUG == 5103 01105 01106 int k; 01107 01108 list<int>::iterator lit_ListIterator; 01109 01110 cout<<endl; 01111 cout<<"DEBUG 5103 | Hessian Evaluation | Induced Vertex Degrees | Set "<<STEP_UP(i)<<endl; 01112 cout<<endl; 01113 01114 for(j=0; j<STEP_UP(i_HighestInducedVertexDegree); j++) 01115 { 01116 i_SetSize = (signed) vli_GroupedInducedVertexDegrees[j].size(); 01117 01118 if(i_SetSize == _FALSE) 01119 { 01120 continue; 01121 } 01122 01123 k = _FALSE; 01124 01125 cout<<"Degree "<<j<<"\t"<<" : "; 01126 01127 for(lit_ListIterator=vli_GroupedInducedVertexDegrees[j].begin(); lit_ListIterator!=vli_GroupedInducedVertexDegrees[j].end(); lit_ListIterator++) 01128 { 01129 if(k == STEP_DOWN(i_SetSize)) 01130 { 01131 cout<<STEP_UP(*lit_ListIterator)<<" ("<<i_SetSize<<")"<<endl; 01132 } 01133 else 01134 { 01135 cout<<STEP_UP(*lit_ListIterator)<<", "; 01136 } 01137 01138 k++; 01139 } 01140 } 01141 01142 #endif 01143 01144 #if DEBUG == 5103 01145 01146 cout<<endl; 01147 cout<<"DEBUG 5103 | Hessian Evaluation | Retrieved Elements"<<"| Set "<<STEP_UP(i)<<endl; 01148 cout<<endl; 01149 01150 #endif 01151 //#define DEBUG 5103 01152 //get the diagonal values 01153 for (int index = 0; index < i_VertexCount; index++) { 01154 if(vi_EvaluatedDiagonals[index] == _FALSE) 01155 { 01156 d_Value = dp2_CompressedMatrix[index][vi_VertexColors[index]]; 01157 01158 #if DEBUG == 5103 01159 01160 cout<<"Element["<<STEP_UP(index)<<"]["<<STEP_UP(index)<<"] = "<<d_Value<<endl; 01161 01162 #endif 01163 v2i_VertexAdjacency[index].push_back(index); 01164 v2d_NonzeroAdjacency[index].push_back(d_Value); 01165 01166 vi_EvaluatedDiagonals[index] = _TRUE; 01167 01168 } 01169 } 01170 01171 for ( ; ; ) 01172 { 01173 if(vli_GroupedInducedVertexDegrees[_TRUE].empty()) // If there is no leaf left on the color tree 01174 { 01175 i_LeafVertex = vli_GroupedInducedVertexDegrees[_FALSE].front(); 01176 01177 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE; 01178 01179 vd_IncludedVertices[i_LeafVertex] = _UNKNOWN; 01180 01181 break; 01182 } 01183 01184 i_LeafVertex = vli_GroupedInducedVertexDegrees[_TRUE].front(); 01185 01186 vli_GroupedInducedVertexDegrees[_TRUE].pop_front(); 01187 01188 //Find i_ParentVertex 01189 for(j=vi_Vertices[i_LeafVertex]; j<vi_Vertices[STEP_UP(i_LeafVertex)]; j++) 01190 { 01191 if(vd_IncludedVertices[vi_Edges[j]] != _UNKNOWN) 01192 { 01193 i_ParentVertex = vi_Edges[j]; 01194 01195 break; 01196 } 01197 } 01198 01199 d_Value = dp2_CompressedMatrix[i_LeafVertex][vi_VertexColors[i_ParentVertex]] - vd_IncludedVertices[i_LeafVertex]; 01200 01201 vd_IncludedVertices[i_ParentVertex] += d_Value; 01202 01203 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE; 01204 vd_IncludedVertices[i_LeafVertex] = _UNKNOWN; 01205 if(vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].size()>1) { 01206 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].erase(vlit_VertexLocations[i_ParentVertex]); 01207 } 01208 else { 01209 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].pop_back(); 01210 } 01211 01212 vi_InducedVertexDegrees[i_ParentVertex]--; 01213 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].push_back(i_ParentVertex); 01214 01215 //Update position of the iterator pointing to i_ParentVertex in "InducedVertexDegrees" structure 01216 vlit_VertexLocations[i_ParentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].end(); 01217 --vlit_VertexLocations[i_ParentVertex]; 01218 01219 v2i_VertexAdjacency[i_LeafVertex].push_back(i_ParentVertex); 01220 v2d_NonzeroAdjacency[i_LeafVertex].push_back(d_Value); 01221 01222 v2i_VertexAdjacency[i_ParentVertex].push_back(i_LeafVertex); 01223 v2d_NonzeroAdjacency[i_ParentVertex].push_back(d_Value); 01224 01225 #if DEBUG == 5103 01226 01227 cout<<"Element["<<STEP_UP(i_LeafVertex)<<"]["<<STEP_UP(i_ParentVertex)<<"] = "<<d_Value<<endl; 01228 #endif 01229 01230 } 01231 } 01232 01233 01234 01235 //cout<<"allocate memory for *dp2_HessianValue rowCount="<<rowCount<<endl; 01236 //printf("i=%d\t numOfNonZerosInHessianValue=%d \n", i, numOfNonZerosInHessianValue); 01237 //(*dp2_HessianValue) = new double[numOfNonZerosInHessianValue]; //allocate memory for *dp2_JacobianValue. 01238 //for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) (*dp2_HessianValue)[i] = 0.; //initialize value of other entries 01239 01240 //populate dp2_HessianValue row by row, column by column 01241 for(i=0; i<i_VertexCount; i++) { 01242 int NumOfNonzeros = uip2_HessianSparsityPattern[i][0]; 01243 int offset = 0; 01244 i_VertexDegree = (signed) v2i_VertexAdjacency[i].size(); 01245 for(j=1; j<=NumOfNonzeros; j++) { 01246 if( i > uip2_HessianSparsityPattern[i][j] ) { 01247 offset++; 01248 continue; 01249 } 01250 int targetColumnID = uip2_HessianSparsityPattern[i][j]; 01251 for (int k=0; k<i_VertexDegree; k++) {// search through the v2i_VertexAdjacency matrix to find the correct column 01252 if(targetColumnID == v2i_VertexAdjacency[i][k]) { //found it 01253 (*dp2_HessianValue)[(*ip2_RowIndex)[i] + j - offset - 1] = v2d_NonzeroAdjacency[i][k]; 01254 break; 01255 } 01256 } 01257 } 01258 } 01259 01260 01261 01262 //Making the array indices to start at 1 instead of 0 to conform with theIntel MKL sparse storage scheme for the direct sparse solvers 01263 for(unsigned int i=0; i <= (unsigned int) i_VertexCount ; i++) { 01264 (*ip2_RowIndex)[i]++; 01265 } 01266 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) { 01267 (*ip2_ColumnIndex)[i]++; 01268 } 01269 01270 return (i_VertexCount); 01271 } 01272 01273 int HessianRecovery::IndirectRecover_SparseSolversFormat_unmanaged(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue, unsigned int numOfNonZerosInHessianValue) { 01274 01275 if(g==NULL) { 01276 cerr<<"g==NULL"<<endl; 01277 return _FALSE; 01278 } 01279 01280 int rowCount = g->GetVertexCount(); 01281 01282 if (numOfNonZerosInHessianValue < 1) { 01283 numOfNonZerosInHessianValue = ConvertRowCompressedFormat2SparseSolversFormat_StructureOnly(uip2_HessianSparsityPattern, rowCount, ip2_RowIndex, ip2_ColumnIndex); 01284 01285 //Making the array indices to start at 1 instead of 0 to conform with the Intel MKL sparse storage scheme for the direct sparse solvers 01286 for(unsigned int i=0; i <= (unsigned int) rowCount ; i++) { 01287 (*ip2_RowIndex)[i]++; 01288 } 01289 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) { 01290 (*ip2_ColumnIndex)[i]++; 01291 } 01292 } 01293 01294 //cout<<"allocate memory for *dp2_HessianValue rowCount="<<rowCount<<endl; 01295 //printf("i=%d\t numOfNonZerosInHessianValue=%d \n", i, numOfNonZerosInHessianValue); 01296 (*dp2_HessianValue) = (double*) malloc(numOfNonZerosInHessianValue * sizeof(double)); //allocate memory for *dp2_JacobianValue. 01297 for(unsigned int i=0; i < numOfNonZerosInHessianValue; i++) (*dp2_HessianValue)[i] = 0.; //initialize value of other entries 01298 01299 int returnValue = IndirectRecover_SparseSolversFormat_usermem(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue, numOfNonZerosInHessianValue); 01300 01301 return returnValue; 01302 } 01303 01304 int HessianRecovery::IndirectRecover_SparseSolversFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 01305 01306 int returnValue = IndirectRecover_SparseSolversFormat_unmanaged( g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue); 01307 01308 if(SSF_available) { 01309 //cout<<"SSF_available="<<SSF_available<<endl; Pause(); 01310 reset(); 01311 } 01312 01313 SSF_available = true; 01314 i_SSF_rowCount = g->GetVertexCount(); 01315 ip_SSF_RowIndex = *ip2_RowIndex; 01316 ip_SSF_ColumnIndex = *ip2_ColumnIndex; 01317 dp_SSF_Value = *dp2_HessianValue; 01318 01319 return (returnValue); 01320 } 01321 01322 int HessianRecovery::IndirectRecover_CoordinateFormat_unmanaged(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** uip2_RowIndex, unsigned int** uip2_ColumnIndex, double** dp2_HessianValue) { 01323 #ifdef _COLPACK_CHECKPOINT_ 01324 cout<<"IN 1 HessianRecovery::IndirectRecover_CoordinateFormat_unmanaged()"<<endl; 01325 string s_postfix = "-IndirectRecover_CoordinateFormat_vectors"; 01326 cout<<"*WriteMatrixMarket_ADOLCInput("<<s_postfix<<", 1, uip2_HessianSparsityPattern, "<< g->GetVertexCount() <<", " << g->GetVertexCount() <<", dp2_CompressedMatrix, " << g->GetVertexCount() <<", " << g->GetVertexColorCount() <<endl; 01327 WriteMatrixMarket_ADOLCInput(s_postfix, 1, uip2_HessianSparsityPattern, g->GetVertexCount(), g->GetVertexCount() , dp2_CompressedMatrix, g->GetVertexCount(), g->GetVertexColorCount() ); 01328 #endif 01329 if(g==NULL) { 01330 cerr<<"g==NULL"<<endl; 01331 return _FALSE; 01332 } 01333 01334 vector<unsigned int> RowIndex; 01335 vector<unsigned int> ColumnIndex; 01336 vector<double> HessianValue; 01337 01338 int returnValue = IndirectRecover_CoordinateFormat_vectors(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, RowIndex, ColumnIndex, HessianValue); 01339 01340 unsigned int numOfNonZeros = returnValue; 01341 (*uip2_RowIndex) = (unsigned int*) malloc(numOfNonZeros * sizeof(unsigned int)); 01342 (*uip2_ColumnIndex) = (unsigned int*) malloc(numOfNonZeros * sizeof(unsigned int)); 01343 (*dp2_HessianValue) = (double*) malloc(numOfNonZeros * sizeof(double)); //allocate memory for *dp2_HessianValue. 01344 01345 for(int i=0; i < numOfNonZeros; i++) { 01346 (*uip2_RowIndex)[i] = RowIndex[i]; 01347 (*uip2_ColumnIndex)[i] = ColumnIndex[i]; 01348 (*dp2_HessianValue)[i] = HessianValue[i]; 01349 } 01350 01351 #ifdef _COLPACK_CHECKPOINT_ 01352 cout<<"IN 2 HessianRecovery::IndirectRecover_CoordinateFormat_unmanaged()"<<endl; 01353 unsigned int*** dp3_Pattern = (unsigned int***) malloc( sizeof(unsigned int**) ); 01354 double*** dp3_Values = (double***) malloc( sizeof(double**) ); 01355 ConvertCoordinateFormat2RowCompressedFormat((*uip2_RowIndex), (*uip2_ColumnIndex), (*dp2_HessianValue), g->GetVertexCount(), numOfNonZeros, dp3_Pattern, dp3_Values ); 01356 01357 s_postfix = "-IndirectRecover_CoordinateFormat_vectors"; 01358 cout<<"*WriteMatrixMarket_ADOLCInput("<<s_postfix<<", 2, uip2_HessianSparsityPattern, "<< g->GetVertexCount() <<", " << g->GetVertexCount() <<", dp2_CompressedMatrix, " << g->GetVertexCount() <<", " << g->GetVertexColorCount()<<", dp3_Values" <<endl; 01359 WriteMatrixMarket_ADOLCInput(s_postfix, 2, uip2_HessianSparsityPattern, g->GetVertexCount(), g->GetVertexCount() , dp2_CompressedMatrix, g->GetVertexCount(), g->GetVertexColorCount(), dp3_Values ); 01360 01361 //Deallocate dp3_Pattern & dp3_Values 01362 freeMatrix(dp3_Pattern, g->GetVertexCount()); 01363 freeMatrix(dp3_Values, g->GetVertexCount()); 01364 #endif 01365 return (returnValue); 01366 } 01367 01368 // !!! Note: the speed of this function could be improved by reserve the vectors size for (RowIndex, ColumnIndex, HessianValue) ahead of time 01369 int HessianRecovery::IndirectRecover_CoordinateFormat_vectors(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, vector<unsigned int> &RowIndex, vector<unsigned int> &ColumnIndex, vector<double> &HessianValue) { 01370 int i=0,j=0; 01371 int i_VertexCount = _UNKNOWN; 01372 int i_EdgeID, i_SetID; 01373 vector<int> vi_Sets; 01374 map< int, vector<int> > mivi_VertexSets; 01375 vector<int> vi_Vertices; 01376 g->GetVertices(vi_Vertices); 01377 vector<int> vi_Edges; 01378 g->GetEdges(vi_Edges); 01379 vector<int> vi_VertexColors; 01380 g->GetVertexColors(vi_VertexColors); 01381 map< int, map< int, int> > mimi2_VertexEdgeMap; 01382 g->GetVertexEdgeMap(mimi2_VertexEdgeMap); 01383 DisjointSets ds_DisjointSets; 01384 g->GetDisjointSets(ds_DisjointSets); 01385 //populate vi_Sets & mivi_VertexSets 01386 vi_Sets.clear(); 01387 mivi_VertexSets.clear(); 01388 01389 i_VertexCount = g->GetVertexCount(); 01390 01391 for(i=0; i<i_VertexCount; i++) // for each vertex A (indexed by i) 01392 { 01393 for(j=vi_Vertices[i]; j<vi_Vertices[STEP_UP(i)]; j++) // for each of the vertex B that connect to A 01394 { 01395 if(i < vi_Edges[j]) // if the index of A (i) is less than the index of B (vi_Edges[j]) 01396 //basically each edge is represented by (vertex with smaller ID, vertex with larger ID). This way, we don't insert a specific edge twice 01397 { 01398 i_EdgeID = mimi2_VertexEdgeMap[i][vi_Edges[j]]; 01399 01400 i_SetID = ds_DisjointSets.FindAndCompress(i_EdgeID); 01401 01402 if(i_SetID == i_EdgeID) // that edge is the root of the set => create new set 01403 { 01404 vi_Sets.push_back(i_SetID); 01405 } 01406 01407 mivi_VertexSets[i_SetID].push_back(i); 01408 mivi_VertexSets[i_SetID].push_back(vi_Edges[j]); 01409 } 01410 } 01411 } 01412 01413 int i_MaximumVertexDegree; 01414 01415 int i_HighestInducedVertexDegree; 01416 01417 int i_LeafVertex, i_ParentVertex, i_PresentVertex; 01418 01419 int i_VertexDegree; 01420 01421 int i_SetCount, i_SetSize; 01422 01423 double d_Value; 01424 01425 vector<int> vi_EvaluatedDiagonals; 01426 01427 vector<int> vi_InducedVertexDegrees; 01428 01429 vector<double> vd_IncludedVertices; 01430 vector<bool> vb_IncludedVertices; 01431 01432 vector< vector<int> > v2i_VertexAdjacency; 01433 01434 vector< vector<double> > v2d_NonzeroAdjacency; 01435 01436 vector< list<int> > vli_GroupedInducedVertexDegrees; 01437 01438 vector< list<int>::iterator > vlit_VertexLocations; 01439 01440 i_MaximumVertexDegree = g->GetMaximumVertexDegree(); 01441 #if DEBUG == 5103 01442 01443 cout<<endl; 01444 cout<<"DEBUG 5103 | Hessian Evaluation | Bicolored Sets"<<endl; 01445 cout<<endl; 01446 01447 i_SetCount = (signed) vi_Sets.size(); 01448 01449 for(i=0; i<i_SetCount; i++) 01450 { 01451 cout<<STEP_UP(vi_Sets[i])<<"\t"<<" : "; 01452 01453 i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size(); 01454 01455 for(j=0; j<i_SetSize; j++) 01456 { 01457 if(j == STEP_DOWN(i_SetSize)) 01458 { 01459 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<" ("<<i_SetSize<<")"<<endl; 01460 } 01461 else 01462 { 01463 cout<<STEP_UP(mivi_VertexSets[vi_Sets[i]][j])<<", "; 01464 } 01465 } 01466 } 01467 01468 cout<<endl; 01469 cout<<"[Set Count = "<<i_SetCount<<"]"<<endl; 01470 cout<<endl; 01471 01472 #endif 01473 01474 //Step 5: from here on 01475 i_VertexCount = g->GetVertexCount(); 01476 01477 v2i_VertexAdjacency.clear(); 01478 v2i_VertexAdjacency.resize((unsigned) i_VertexCount); 01479 01480 v2d_NonzeroAdjacency.clear(); 01481 v2d_NonzeroAdjacency.resize((unsigned) i_VertexCount); 01482 01483 vi_EvaluatedDiagonals.clear(); 01484 vi_EvaluatedDiagonals.resize((unsigned) i_VertexCount, _FALSE); 01485 01486 vi_InducedVertexDegrees.clear(); 01487 vi_InducedVertexDegrees.resize((unsigned) i_VertexCount, _FALSE); 01488 01489 vd_IncludedVertices.clear(); 01490 vd_IncludedVertices.resize((unsigned) i_VertexCount, 0.); 01491 vb_IncludedVertices.clear(); 01492 vb_IncludedVertices.resize((unsigned) i_VertexCount, false); 01493 01494 i_ParentVertex = _UNKNOWN; 01495 01496 i_SetCount = (signed) vi_Sets.size(); 01497 01498 for(i=0; i<i_SetCount; i++) 01499 { 01500 vli_GroupedInducedVertexDegrees.clear(); 01501 vli_GroupedInducedVertexDegrees.resize((unsigned) STEP_UP(i_MaximumVertexDegree)); 01502 01503 vlit_VertexLocations.clear(); 01504 vlit_VertexLocations.resize((unsigned) i_VertexCount); 01505 01506 i_HighestInducedVertexDegree = _UNKNOWN; 01507 01508 i_SetSize = (signed) mivi_VertexSets[vi_Sets[i]].size(); 01509 01510 for(j=0; j<i_SetSize; j++) 01511 { 01512 i_PresentVertex = mivi_VertexSets[vi_Sets[i]][j]; 01513 01514 vb_IncludedVertices[i_PresentVertex] = true; 01515 vd_IncludedVertices[i_PresentVertex] = 0; 01516 01517 if(vi_InducedVertexDegrees[i_PresentVertex] != _FALSE) 01518 { 01519 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].erase(vlit_VertexLocations[i_PresentVertex]); 01520 } 01521 01522 vi_InducedVertexDegrees[i_PresentVertex]++; 01523 01524 if(i_HighestInducedVertexDegree < vi_InducedVertexDegrees[i_PresentVertex]) 01525 { 01526 i_HighestInducedVertexDegree = vi_InducedVertexDegrees[i_PresentVertex]; 01527 } 01528 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].push_front(i_PresentVertex); 01529 01530 vlit_VertexLocations[i_PresentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_PresentVertex]].begin(); 01531 } 01532 01533 #if DEBUG == 5103 01534 01535 int k; 01536 01537 list<int>::iterator lit_ListIterator; 01538 01539 cout<<endl; 01540 cout<<"DEBUG 5103 | Hessian Evaluation | Induced Vertex Degrees | Set "<<STEP_UP(i)<<endl; 01541 cout<<endl; 01542 01543 for(j=0; j<STEP_UP(i_HighestInducedVertexDegree); j++) 01544 { 01545 i_SetSize = (signed) vli_GroupedInducedVertexDegrees[j].size(); 01546 01547 if(i_SetSize == _FALSE) 01548 { 01549 continue; 01550 } 01551 01552 k = _FALSE; 01553 01554 cout<<"Degree "<<j<<"\t"<<" : "; 01555 01556 for(lit_ListIterator=vli_GroupedInducedVertexDegrees[j].begin(); lit_ListIterator!=vli_GroupedInducedVertexDegrees[j].end(); lit_ListIterator++) 01557 { 01558 if(k == STEP_DOWN(i_SetSize)) 01559 { 01560 cout<<STEP_UP(*lit_ListIterator)<<"{"<<vd_IncludedVertices[*lit_ListIterator]<<"} ("<<i_SetSize<<")"<<endl; 01561 } 01562 else 01563 { 01564 cout<<STEP_UP(*lit_ListIterator)<<"{"<<vd_IncludedVertices[*lit_ListIterator]<<"}, "; 01565 } 01566 01567 k++; 01568 } 01569 } 01570 01571 #endif 01572 01573 #if DEBUG == 5103 01574 01575 cout<<endl; 01576 cout<<"DEBUG 5103 | Hessian Evaluation | Retrieved Elements"<<"| Set "<<STEP_UP(i)<<endl; 01577 cout<<endl; 01578 01579 #endif 01580 //get the diagonal values 01581 for (int index = 0; index < i_VertexCount; index++) { 01582 if(vi_EvaluatedDiagonals[index] == _FALSE) 01583 { 01584 d_Value = dp2_CompressedMatrix[index][vi_VertexColors[index]]; 01585 01586 #if DEBUG == 5103 01587 01588 cout<<"Element["<<STEP_UP(index)<<"]["<<STEP_UP(index)<<"] = "<<d_Value<<endl; 01589 01590 #endif 01591 v2i_VertexAdjacency[index].push_back(index); 01592 v2d_NonzeroAdjacency[index].push_back(d_Value); 01593 01594 vi_EvaluatedDiagonals[index] = _TRUE; 01595 01596 } 01597 } 01598 01599 for ( ; ; ) 01600 { 01601 if(vli_GroupedInducedVertexDegrees[_TRUE].empty()) // If there is no leaf left on the color tree 01602 { 01603 i_LeafVertex = vli_GroupedInducedVertexDegrees[_FALSE].front(); 01604 01605 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE; 01606 01607 vb_IncludedVertices[i_LeafVertex] = false; 01608 01609 break; 01610 } 01611 01612 i_LeafVertex = vli_GroupedInducedVertexDegrees[_TRUE].front(); 01613 01614 vli_GroupedInducedVertexDegrees[_TRUE].pop_front(); 01615 01616 //Find i_ParentVertex 01617 for(j=vi_Vertices[i_LeafVertex]; j<vi_Vertices[STEP_UP(i_LeafVertex)]; j++) 01618 { 01619 if(vb_IncludedVertices[vi_Edges[j]]) 01620 { 01621 i_ParentVertex = vi_Edges[j]; 01622 01623 break; 01624 } 01625 } 01626 01627 d_Value = dp2_CompressedMatrix[i_LeafVertex][vi_VertexColors[i_ParentVertex]] - vd_IncludedVertices[i_LeafVertex]; 01628 01629 vd_IncludedVertices[i_ParentVertex] += d_Value; 01630 01631 vi_InducedVertexDegrees[i_LeafVertex] = _FALSE; 01632 vb_IncludedVertices[i_LeafVertex] = false; 01633 if(vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].size()>1) { 01634 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].erase(vlit_VertexLocations[i_ParentVertex]); 01635 } 01636 else { 01637 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].pop_back(); 01638 } 01639 01640 vi_InducedVertexDegrees[i_ParentVertex]--; 01641 vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].push_back(i_ParentVertex); 01642 01643 //Update position of the iterator pointing to i_ParentVertex in "InducedVertexDegrees" structure 01644 vlit_VertexLocations[i_ParentVertex] = vli_GroupedInducedVertexDegrees[vi_InducedVertexDegrees[i_ParentVertex]].end(); 01645 --vlit_VertexLocations[i_ParentVertex]; 01646 01647 v2i_VertexAdjacency[i_LeafVertex].push_back(i_ParentVertex); 01648 v2d_NonzeroAdjacency[i_LeafVertex].push_back(d_Value); 01649 01650 v2i_VertexAdjacency[i_ParentVertex].push_back(i_LeafVertex); 01651 v2d_NonzeroAdjacency[i_ParentVertex].push_back(d_Value); 01652 01653 01654 #if DEBUG == 5103 01655 01656 cout<<"Element["<<STEP_UP(i_LeafVertex)<<"]["<<STEP_UP(i_ParentVertex)<<"] = "<<d_Value<<endl; 01657 #endif 01658 01659 } 01660 } 01661 01662 //populate dp3_HessianValue row by row, column by column 01663 for(i=0; i<i_VertexCount; i++) { 01664 int NumOfNonzeros = uip2_HessianSparsityPattern[i][0]; 01665 i_VertexDegree = (signed) v2i_VertexAdjacency[i].size(); 01666 for(j=1; j<=NumOfNonzeros; j++) { 01667 int targetColumnID = uip2_HessianSparsityPattern[i][j]; 01668 if(targetColumnID<i) continue; 01669 for (int k=0; k<i_VertexDegree; k++) {// search through the v2i_VertexAdjacency matrix to find the correct column 01670 if(targetColumnID == v2i_VertexAdjacency[i][k]) { //found it 01671 HessianValue.push_back(v2d_NonzeroAdjacency[i][k]); 01672 break; 01673 } 01674 } 01675 RowIndex.push_back(i); 01676 ColumnIndex.push_back(uip2_HessianSparsityPattern[i][j]); 01677 } 01678 } 01679 01680 return ( RowIndex.size() ); 01681 } 01682 01683 int HessianRecovery::IndirectRecover_CoordinateFormat_usermem(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 01684 if(g==NULL) { 01685 cerr<<"g==NULL"<<endl; 01686 return _FALSE; 01687 } 01688 01689 vector<unsigned int> RowIndex; 01690 vector<unsigned int> ColumnIndex; 01691 vector<double> HessianValue; 01692 01693 01694 int returnValue = IndirectRecover_CoordinateFormat_vectors(g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, RowIndex, ColumnIndex, HessianValue); 01695 01696 unsigned int numOfNonZeros = returnValue; 01697 01698 for(int i=0; i < numOfNonZeros; i++) { 01699 (*ip2_RowIndex)[i] = RowIndex[i]; 01700 (*ip2_ColumnIndex)[i] = ColumnIndex[i]; 01701 (*dp2_HessianValue)[i] = HessianValue[i]; 01702 } 01703 01704 return (returnValue); 01705 } 01706 01707 int HessianRecovery::IndirectRecover_CoordinateFormat(GraphColoringInterface* g, double** dp2_CompressedMatrix, unsigned int ** uip2_HessianSparsityPattern, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex, double** dp2_HessianValue) { 01708 //#define DEBUG 5103 01709 int returnValue = IndirectRecover_CoordinateFormat_unmanaged( g, dp2_CompressedMatrix, uip2_HessianSparsityPattern, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue); 01710 01711 if(CF_available) reset(); 01712 01713 CF_available = true; 01714 i_CF_rowCount = returnValue; 01715 ip_CF_RowIndex = *ip2_RowIndex; 01716 ip_CF_ColumnIndex = *ip2_ColumnIndex; 01717 dp_CF_Value = *dp2_HessianValue; 01718 01719 return (returnValue); 01720 } 01721 }
1.7.6.1 
