00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "extra.h"
00022 #include "Pause.h"
00023 #include "mmio.h"
00024
00025 bool isValidOrdering(vector<int> & ordering, int offset) {
00026 vector<bool> isExist;
00027 int orderingNum = 0;
00028 isExist.resize(ordering.size(), false);
00029 for(int i=0; i<ordering.size(); i++) {
00030 orderingNum = ordering[i] - offset;
00031 if(orderingNum<0 || orderingNum>= ordering.size()) {
00032 cerr<<" This vertex # is not in the valid range [0, ordering.size()]. ordering[i]: "<<ordering[i]<<endl;
00033 return false;
00034 }
00035
00036 if(isExist[ orderingNum ]) {
00037 cerr<<"This vertex has been seen before. We have duplication!"<<endl;
00038 return false;
00039 }
00040
00041 isExist[ orderingNum ] = true;
00042 }
00043
00044 return true;
00045 }
00046
00047 int ReadRowCompressedFormat(string s_InputFile, unsigned int *** uip3_SparsityPattern, int& rowCount, int& columnCount) {
00048 string line;
00049 int lineCounter = 0,nz_counter = 0, nonzeros = 0, nnz_per_row = 0;
00050 unsigned int num = 0;
00051 istringstream in2;
00052 ifstream in (s_InputFile.c_str());
00053
00054 if(!in) {
00055 cout<<s_InputFile<<" not Found!"<<endl;
00056 exit(1);
00057 }
00058
00059 getline(in,line);
00060 lineCounter++;
00061 in2.str(line);
00062 in2 >> rowCount >> columnCount >> nonzeros;
00063
00064 (*uip3_SparsityPattern) = new unsigned int*[rowCount];
00065
00066 for(int i=0;i < rowCount; i++) {
00067 getline(in, line);
00068 lineCounter++;
00069 if(line!="")
00070 {
00071 in2.clear();
00072 in2.str(line);
00073 in2>>nnz_per_row;
00074 (*uip3_SparsityPattern)[i] = new unsigned int[nnz_per_row + 1];
00075 (*uip3_SparsityPattern)[i][0] = nnz_per_row;
00076
00077 for(int j=1; j<nnz_per_row+1; j++) {
00078 in2>>num;
00079 (*uip3_SparsityPattern)[i][j] = num;
00080 nz_counter++;
00081 }
00082 }
00083 else
00084 {
00085 cerr<<"* WARNING: ReadRowCompressedFormat()"<<endl;
00086 cerr<<"*\t line == \"\" at row "<<lineCounter<<". Empty line. Wrong input format. Can't process."<<endl;
00087 cerr<<"\t total non-zeros so far: "<<nz_counter<<endl;
00088 exit( -1);
00089 }
00090 }
00091
00092 if(nz_counter<nonzeros) {
00093 cerr<<"* WARNING: ReadRowCompressedFormat()"<<endl;
00094 cerr<<"*\t nz_counter<nonzeros+1. Wrong input format. Can't process."<<endl;
00095 cerr<<"\t total non-zeros so far: "<<nz_counter<<endl;
00096 exit( -1);
00097 }
00098
00099
00100
00101 return 0;
00102
00103 }
00104
00105 int RowCompressedFormat_2_SparseSolversFormat_StructureOnly (unsigned int ** uip2_HessianSparsityPattern, unsigned int ui_rowCount, unsigned int** ip2_RowIndex, unsigned int** ip2_ColumnIndex) {
00106
00107
00108 unsigned int nnz = 0;
00109 unsigned int nnz_in1Row = 0;
00110 (*ip2_RowIndex) = new unsigned int[ui_rowCount + 1];
00111 for (unsigned int i=0; i < ui_rowCount; i++) {
00112 nnz_in1Row = uip2_HessianSparsityPattern[i][0];
00113 (*ip2_RowIndex)[i] = nnz;
00114 for (unsigned int j = 1; j <= nnz_in1Row ; j++) {
00115 if (i <= uip2_HessianSparsityPattern[i][j]) nnz++;
00116 }
00117 }
00118 (*ip2_RowIndex)[ui_rowCount] = nnz;
00119
00120
00121 displayVector(*ip2_RowIndex,ui_rowCount+1);
00122
00123
00124 (*ip2_ColumnIndex) = new unsigned int[nnz];
00125 unsigned int count = 0;
00126 for (unsigned int i=0; i < ui_rowCount; i++) {
00127 nnz_in1Row = uip2_HessianSparsityPattern[i][0];
00128 for (unsigned int j = 1; j <= nnz_in1Row ; j++) {
00129 if (i <= uip2_HessianSparsityPattern[i][j]) {
00130 (*ip2_ColumnIndex)[count] = uip2_HessianSparsityPattern[i][j];
00131 count++;
00132 }
00133 }
00134 }
00135 if(count != nnz) {
00136 cerr<<"!!! count != nnz. count = "<<count<<endl;
00137 Pause();
00138 }
00139
00140 return nnz;
00141 }
00142
00143
00144 void ConvertDIMACSFormat2MatrixMarketFormat(string fileNameNoExt) {
00145 string inFileName = fileNameNoExt + ".gr";
00146 string outFileName = fileNameNoExt + ".mtx";
00147 string line, temp;
00148 ifstream in(inFileName.c_str());
00149 ofstream out(outFileName.c_str());
00150 istringstream iin;
00151
00152 while(in) {
00153 getline(in, line);
00154 if(line=="") break;
00155 switch(line[0]) {
00156 case 'a':
00157
00158 out<<line.substr(2,line.size()-2)<<endl;
00159 break;
00160 case 'c':
00161 break;
00162 default:
00163
00164 iin.str(line);
00165 iin>>temp>>temp>>temp;out<<temp<<" "<<temp<<" ";
00166 iin>>temp;out<<temp<<endl;
00167 break;
00168 }
00169 }
00170
00171 in.close();
00172 out.close();
00173 }
00174
00175 void randomOrdering(vector<int>& ordering) {
00176 srand(time(NULL));
00177 int size = ordering.size();
00178 int ran_num = 0;
00179 for(int i=0; i < size-1; i++) {
00180
00181 ran_num = (int)(((float) rand() / RAND_MAX) * (size -1 - i)) + i;
00182 swap(ordering[i],ordering[ran_num]);
00183 }
00184 }
00185
00186 string toUpper(string input) {
00187 string output = input;
00188
00189 for(int i = input.size() - 1; i>=0; i--) {
00190 if(input[i]==' ' || input[i]=='\t' || input[i]=='\n') {
00191 output[i] = '_';
00192 }
00193 else {
00194 output[i] = toupper(input[i]);
00195 }
00196 }
00197
00198 return output;
00199 }
00200
00201
00202 int Times2Plus1point5(double** dp2_Values, int i_RowCount, int i_ColumnCount) {
00203 for(int i=0; i < i_RowCount; i++) {
00204 for(int j=0; j < i_ColumnCount; j++) {
00205 if(dp2_Values[i][j] != 0.) dp2_Values[i][j] = dp2_Values[i][j]*2 + 1.5;
00206 }
00207
00208 }
00209 return 0;
00210 }
00211 int Times2(double** dp2_Values, int i_RowCount, int i_ColumnCount) {
00212 for(int i=0; i < i_RowCount; i++) {
00213 for(int j=0; j < i_ColumnCount; j++) {
00214 if(dp2_Values[i][j] != 0.) dp2_Values[i][j] = dp2_Values[i][j]*2;
00215 }
00216
00217 }
00218 return 0;
00219 }
00220
00221 int GenerateValues(unsigned int ** uip2_SparsityPattern, int rowCount, double*** dp3_Value) {
00222
00223 srand(0);
00224
00225 (*dp3_Value) = new double*[rowCount];
00226 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00227 unsigned int numOfNonZeros = uip2_SparsityPattern[i][0];
00228 (*dp3_Value)[i] = new double[numOfNonZeros + 1];
00229 (*dp3_Value)[i][0] = (double)numOfNonZeros;
00230 for(unsigned int j=1; j <= numOfNonZeros; j++) {
00231 (*dp3_Value)[i][j] = (rand()%2001 - 1000)/1000.0;
00232
00233 }
00234 }
00235
00236 return 0;
00237 }
00238
00239 int GenerateValuesForSymmetricMatrix(unsigned int ** uip2_SparsityPattern, int rowCount, double*** dp3_Value) {
00240
00241 srand(0);
00242
00243 int * nnzCount = new int[rowCount];
00244 for(unsigned int i=0; i < (unsigned int)rowCount; i++) nnzCount[i] = 0;
00245
00246 (*dp3_Value) = new double*[rowCount];
00247 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00248 unsigned int numOfNonZeros = uip2_SparsityPattern[i][0];
00249 (*dp3_Value)[i] = new double[numOfNonZeros + 1];
00250 (*dp3_Value)[i][0] = (double)numOfNonZeros;
00251 for(unsigned int j=1; j <= numOfNonZeros; j++) {
00252 if (uip2_SparsityPattern[i][j] >i) break;
00253 (*dp3_Value)[i][j] = (rand()%2001 - 1000)/1000.0; nnzCount[i]++;
00254 if (uip2_SparsityPattern[i][j] <i) {
00255 (*dp3_Value)[uip2_SparsityPattern[i][j]][nnzCount[uip2_SparsityPattern[i][j]]+1] = (*dp3_Value)[i][j]; nnzCount[uip2_SparsityPattern[i][j]]++;
00256 }
00257
00258 }
00259 }
00260
00261 delete[] nnzCount;
00262
00263 return 0;
00264 }
00265
00266 int ConvertMatrixMarketFormatToRowCompressedFormat(string s_InputFile, unsigned int *** uip3_SparsityPattern, double*** dp3_Value, int &rowCount, int &columnCount) {
00267
00268 string m_s_InputFile=s_InputFile;
00269
00270
00271 int rowCounter=0, nonzeros=0, rowIndex=0, colIndex=0, nz_counter=0;
00272
00273 float value;
00274 bool b_getValue, b_symmetric;
00275 istringstream in2;
00276 string line="";
00277 map<int,vector<int> > nodeList;
00278 map<int,vector<float> > valueList;
00279
00280
00281 MM_typecode matcode;
00282 FILE *f;
00283 if ((f = fopen(m_s_InputFile.c_str(), "r")) == NULL) {
00284 cout<<m_s_InputFile<<" not Found!"<<endl;
00285 exit(1);
00286 }
00287 else cout<<"Found file "<<m_s_InputFile<<endl;
00288
00289 if (mm_read_banner(f, &matcode) != 0)
00290 {
00291 printf("Could not process Matrix Market banner.\n");
00292 exit(1);
00293 }
00294
00295 if( mm_is_pattern(matcode) ) {
00296 b_getValue = false;
00297 }
00298 else b_getValue = true;
00299
00300 if(mm_is_symmetric(matcode)) {
00301 b_symmetric = true;
00302 }
00303 else b_symmetric = false;
00304
00305
00306 char * result = mm_typecode_to_str(matcode);
00307 printf("Graph of Market Market type: [%s]\n", result);
00308 free(result);
00309 if (b_getValue) printf("\t Graph structure and VALUES will be read\n");
00310 else printf("\t Read graph struture only. Values will NOT be read\n");
00311 if( !( mm_is_coordinate(matcode) && (mm_is_symmetric(matcode) || mm_is_general(matcode) ) && ( mm_is_real(matcode) || mm_is_pattern(matcode) || mm_is_integer(matcode) ) ) ) {
00312 printf("Sorry, this application does not support this type.");
00313 exit(1);
00314 }
00315
00316 fclose(f);
00317
00318
00319
00320 ifstream in (m_s_InputFile.c_str());
00321 if(!in) {
00322 cout<<m_s_InputFile<<" not Found!"<<endl;
00323 exit(1);
00324 }
00325 else {
00326
00327 }
00328
00329 getline(in,line);
00330 rowCounter++;
00331 while(line.size()>0&&line[0]=='%') {
00332 getline(in,line);
00333 }
00334 in2.str(line);
00335 in2 >> rowCount >> columnCount >> nonzeros;
00336
00337
00338
00339 while(!in.eof() && nz_counter<nonzeros)
00340 {
00341 getline(in,line);
00342 rowCounter++;
00343
00344 if(line!="")
00345 {
00346 in2.clear();
00347 in2.str(line);
00348 in2>>rowIndex>>colIndex;
00349 rowIndex--;
00350 colIndex--;
00351
00352 if(b_symmetric) {
00353 if(rowIndex > colIndex) {
00354
00355
00356 nodeList[rowIndex].push_back(colIndex);
00357 nodeList[colIndex].push_back(rowIndex);
00358 nz_counter += 2;
00359
00360 if(b_getValue) {
00361 in2>>value;
00362
00363 valueList[rowIndex].push_back(value);
00364 valueList[colIndex].push_back(value);
00365 }
00366 }
00367 else if (rowIndex == colIndex) {
00368
00369 nodeList[rowIndex].push_back(rowIndex);
00370 nz_counter++;
00371 if(b_getValue) {
00372 in2>>value;
00373 valueList[rowIndex].push_back(value);
00374 }
00375 }
00376 else {
00377 cerr<<"* WARNING: ConvertMatrixMarketFormatToRowCompressedFormat()"<<endl;
00378 cerr<<"\t Found a nonzero in the upper triangular. A symmetric Matrix Market file format should only specify the nonzeros in the lower triangular."<<endl;
00379 exit( -1);
00380 }
00381 }
00382 else {
00383 cout<<"\t"<<setw(4)<<rowIndex<<setw(4)<<colIndex<<setw(4)<<nz_counter<<endl;
00384 nodeList[rowIndex].push_back(colIndex);
00385 nz_counter++;
00386 if(b_getValue) {
00387 in2>>value;
00388 cout<<"Value = "<<value<<endl;
00389 valueList[rowIndex].push_back(value);
00390 }
00391 }
00392
00393 }
00394 else
00395 {
00396 cerr<<"* WARNING: ConvertMatrixMarketFormatToRowCompressedFormat()"<<endl;
00397 cerr<<"*\t line == \"\" at row "<<rowCounter<<". Empty line. Wrong input format. Can't process."<<endl;
00398 cerr<<"\t total non-zeros so far: "<<nz_counter<<endl;
00399 exit( -1);
00400 }
00401 }
00402 if(nz_counter<nonzeros) {
00403 cerr<<"* WARNING: ConvertMatrixMarketFormatToRowCompressedFormat()"<<endl;
00404 cerr<<"*\t nz_counter<nonzeros+1. Wrong input format. Can't process."<<endl;
00405 cerr<<"\t total non-zeros so far: "<<nz_counter<<endl;
00406 exit( -1);
00407 }
00408
00409
00410 (*uip3_SparsityPattern) = new unsigned int*[rowCount];
00411 if(b_getValue) (*dp3_Value) = new double*[rowCount];
00412 for(int i=0;i<rowCount; i++) {
00413 unsigned int numOfNonZeros = nodeList[i].size();
00414
00415
00416
00417 (*uip3_SparsityPattern)[i] = new unsigned int[numOfNonZeros+1];
00418 (*uip3_SparsityPattern)[i][0] = numOfNonZeros;
00419
00420 if(b_getValue) {
00421 (*dp3_Value)[i] = new double[numOfNonZeros+1];
00422 (*dp3_Value)[i][0] = (double)numOfNonZeros;
00423 }
00424
00425 for(unsigned int j=0; j < numOfNonZeros; j++) {
00426 (*uip3_SparsityPattern)[i][j+1] = nodeList[i][j];
00427
00428 }
00429
00430 if(b_getValue) for(unsigned int j=0; j < numOfNonZeros; j++) {
00431 (*dp3_Value)[i][j+1] = valueList[i][j];
00432 }
00433
00434 }
00435
00436
00437 return(0);
00438 }
00439
00440 int MatrixMultiplication_VxS(unsigned int ** uip3_SparsityPattern, double** dp3_Value, int rowCount, int columnCount, double** dp2_seed, int colorCount, double*** dp3_CompressedMatrix) {
00441
00442
00443
00444 (*dp3_CompressedMatrix) = new double*[rowCount];
00445 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00446 (*dp3_CompressedMatrix)[i] = new double[colorCount];
00447 for(unsigned int j=0; j < (unsigned int)colorCount; j++) {
00448 (*dp3_CompressedMatrix)[i][j] = 0.;
00449 }
00450 }
00451
00452
00453
00454 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00455 unsigned int numOfNonZeros = uip3_SparsityPattern[i][0];
00456 for(unsigned int j=1; j <= numOfNonZeros; j++) {
00457 for(unsigned int k=0; k < (unsigned int)colorCount; k++) {
00458
00459 (*dp3_CompressedMatrix)[i][k] += dp3_Value[i][j]*dp2_seed[uip3_SparsityPattern[i][j]][k];
00460 }
00461 }
00462 }
00463
00464 return 0;
00465 }
00466
00467 int MatrixMultiplication_SxV(unsigned int ** uip3_SparsityPattern, double** dp3_Value, int rowCount, int columnCount, double** dp2_seed, int colorCount, double*** dp3_CompressedMatrix) {
00468
00469
00470
00471 (*dp3_CompressedMatrix) = new double*[colorCount];
00472 for(unsigned int i=0; i < (unsigned int)colorCount; i++) {
00473 (*dp3_CompressedMatrix)[i] = new double[columnCount];
00474 for(unsigned int j=0; j < (unsigned int)columnCount; j++) {
00475 (*dp3_CompressedMatrix)[i][j] = 0.;
00476 }
00477 }
00478
00479
00480
00481 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00482 unsigned int numOfNonZeros = uip3_SparsityPattern[i][0];
00483 for(unsigned int j=1; j <= numOfNonZeros; j++) {
00484 for(unsigned int k=0; k < (unsigned int)colorCount; k++) {
00485
00486 (*dp3_CompressedMatrix)[k][uip3_SparsityPattern[i][j]] += dp2_seed[k][i]*dp3_Value[i][j];
00487 }
00488 }
00489 }
00490
00491 return 0;
00492 }
00493
00494 bool CompressedRowMatricesREqual(double** dp3_Value, double** dp3_NewValue, int rowCount, bool compare_exact, bool print_all) {
00495 double ratio = 1.;
00496 int none_equal_count = 0;
00497
00498 for(unsigned int i=0; i < (unsigned int)rowCount; i++) {
00499 unsigned int numOfNonZeros = (unsigned int)dp3_Value[i][0];
00500 if (numOfNonZeros != (unsigned int)dp3_NewValue[i][0]) {
00501 printf("Number of non-zeros in row %d are not equal. dp3_Value[i][0] = %d; dp3_NewValue[i][0] = %d; \n",i,(unsigned int)dp3_Value[i][0],(unsigned int)dp3_NewValue[i][0]);
00502 if (print_all) {
00503 none_equal_count++;
00504 continue;
00505 }
00506 else return false;
00507 }
00508 for(unsigned int j=0; j <= numOfNonZeros; j++) {
00509 if (compare_exact) {
00510 if (dp3_Value[i][j] != dp3_NewValue[i][j]) {
00511 printf("At row %d, column %d, dp3_Value[i][j](%f) != dp3_NewValue[i][j](%f) \n",i,j,dp3_Value[i][j],dp3_NewValue[i][j]);
00512 if (print_all) {
00513 none_equal_count++;
00514 }
00515 else {
00516 printf("You may want to set the flag \"compare_exact\" to 0 to compare the values approximately\n");
00517 return false;
00518 }
00519 }
00520 }
00521 else {
00522 if(dp3_NewValue[i][j] == 0.) {
00523 if(dp3_Value[i][j] != 0.) {
00524 printf("At row %d, column %d, dp3_Value[i][j](%f) != dp3_NewValue[i][j](0) \n",i,j,dp3_Value[i][j]);
00525 if (print_all) {
00526 none_equal_count++;
00527 }
00528 else return false;
00529 }
00530 }
00531 else {
00532 ratio = dp3_Value[i][j] / dp3_NewValue[i][j];
00533 if( ratio < .99 || ratio > 1.02) {
00534 printf("At row %d, column %d, dp3_Value[i][j](%f) != dp3_NewValue[i][j](%f) ; dp3_Value[i][j] / dp3_NewValue[i][j]=%f\n",i,j,dp3_Value[i][j],dp3_NewValue[i][j], ratio);
00535 if (print_all) {
00536 none_equal_count++;
00537 }
00538 else return false;
00539 }
00540 }
00541 }
00542 }
00543 }
00544
00545 if(none_equal_count!=0) {
00546 printf("Total: %d lines. (The total # of non-equals can be greater)\n",none_equal_count);
00547 if (compare_exact) printf("You may want to set the flag \"compare_exact\" to 0 to compare the values approximately\n");
00548 return false;
00549 }
00550 else return true;
00551 }
00552
00553