ColPack
|
00001 // An example of (Column) compression (by using Star coloring) and direct recovery for Hessian 00002 /* How to compile this driver manually: 00003 Please make sure that "baseDir" point to the directory (folder) containing the input matrix file, and 00004 s_InputFile should point to the input file that you want to use 00005 To compile the code, replace the Main.cpp file in Main directory with this file 00006 and run "make" in ColPack installation directory. Make will generate "ColPack.exe" executable 00007 Run "ColPack.exe" 00008 00009 Note: If you got "symbol lookup error ... undefined symbol " 00010 Please make sure that your LD_LIBRARY_PATH contains libColPack.so 00011 00012 //*/ 00013 00014 #include "ColPackHeaders.h" 00015 00016 using namespace ColPack; 00017 using namespace std; 00018 00019 #ifndef TOP_DIR 00020 #define TOP_DIR "." 00021 #endif 00022 00023 // baseDir should point to the directory (folder) containing the input file 00024 string baseDir=TOP_DIR; 00025 00026 #include "extra.h" //This .h file contains functions that are used in the below examples: 00027 //ReadMM(), MatrixMultiplication...(), Times2Plus1point5(), displayMatrix() and displayCompressedRowMatrix() 00028 00029 int main() 00030 { 00031 // s_InputFile = baseDir + <name of the input file> 00032 string s_InputFile; //path of the input file 00033 s_InputFile = baseDir; 00034 s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "mtx-spear-head.mtx"; 00035 00036 // Step 1: Determine sparsity structure of the Jacobian. 00037 // This step is done by an AD tool. For the purpose of illustration here, we read the structure from a file, 00038 // and store the structure in a Compressed Row Format. 00039 unsigned int *** uip3_SparsityPattern = new unsigned int **; 00040 double*** dp3_Value = new double**; 00041 int rowCount, columnCount; 00042 ConvertMatrixMarketFormat2RowCompressedFormat(s_InputFile, uip3_SparsityPattern, dp3_Value,rowCount, columnCount); 00043 00044 cout<<"just for debugging purpose, display the 2 matrices: the matrix with SparsityPattern only and the matrix with Value"<<endl; 00045 cout<<fixed<<showpoint<<setprecision(2); //formatting output 00046 cout<<"(*uip3_SparsityPattern)"<<endl; 00047 displayCompressedRowMatrix((*uip3_SparsityPattern),rowCount); 00048 cout<<"(*dp3_Value)"<<endl; 00049 displayCompressedRowMatrix((*dp3_Value),rowCount); 00050 cout<<"Finish ConvertMatrixMarketFormat2RowCompressedFormat()"<<endl; 00051 Pause(); 00052 00053 //Step 2: Obtain the seed matrix via coloring. 00054 double*** dp3_Seed = new double**; 00055 int *ip1_SeedRowCount = new int; 00056 int *ip1_SeedColumnCount = new int; 00057 int *ip1_ColorCount = new int; //The number of distinct colors used to color the graph 00058 00059 //Step 2.1: Read the sparsity pattern of the given Hessian matrix (compressed sparse rows format) 00060 //and create the corresponding graph 00061 GraphColoringInterface *g = new GraphColoringInterface(SRC_MEM_ADOLC,*uip3_SparsityPattern, rowCount); 00062 00063 //Step 2.2: Color the bipartite graph with the specified ordering 00064 g->Coloring("SMALLEST_LAST", "STAR"); 00065 00066 //Step 2.3 (Option 1): From the coloring information, create and return the seed matrix 00067 (*dp3_Seed) = g->GetSeedMatrix(ip1_SeedRowCount, ip1_SeedColumnCount); 00068 /* Notes: 00069 Step 2.3 (Option 2): From the coloring information, you can also get the vector of colorIDs of vertices 00070 vector<int> vi_VertexColors; 00071 g->GetVertexColors(vi_VertexColors); 00072 */ 00073 cout<<"Finish GenerateSeed()"<<endl; 00074 *ip1_ColorCount = *ip1_SeedColumnCount; 00075 00076 displayMatrix(*dp3_Seed,columnCount,*ip1_SeedColumnCount); 00077 Pause(); 00078 00079 // Step 3: Obtain the Hessian-seed matrix product. 00080 // This step will also be done by an AD tool. For the purpose of illustration here, the orginial matrix V 00081 // (for Values) is multiplied with the seed matrix S. The resulting matrix is stored in dp3_CompressedMatrix. 00082 double*** dp3_CompressedMatrix = new double**; 00083 cout<<"Start MatrixMultiplication()"<<endl; 00084 MatrixMultiplication_VxS(*uip3_SparsityPattern, *dp3_Value, rowCount, columnCount, *dp3_Seed, *ip1_ColorCount, dp3_CompressedMatrix); 00085 cout<<"Finish MatrixMultiplication()"<<endl; 00086 00087 displayMatrix(*dp3_CompressedMatrix,rowCount,*ip1_ColorCount); 00088 Pause(); 00089 00090 //Step 4: Recover the numerical values of the original matrix from the compressed representation. 00091 // The new values are store in "dp3_NewValue" 00092 double*** dp3_NewValue = new double**; 00093 HessianRecovery* hr = new HessianRecovery; 00094 hr->DirectRecover_RowCompressedFormat(g, *dp3_CompressedMatrix, *uip3_SparsityPattern, dp3_NewValue); 00095 cout<<"Finish Recover()"<<endl; 00096 00097 displayCompressedRowMatrix(*dp3_NewValue,rowCount); 00098 Pause(); 00099 00100 //Check for consistency, make sure the values in the 2 matrices are the same. 00101 if (CompressedRowMatricesAreEqual(*dp3_Value, *dp3_NewValue, rowCount,0,1)) cout<< "*dp3_Value == dp3_NewValue"<<endl; 00102 else cout<< "*dp3_Value != dp3_NewValue"<<endl; 00103 00104 Pause(); 00105 00106 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00107 00108 free_2DMatrix(uip3_SparsityPattern, rowCount); 00109 uip3_SparsityPattern=NULL; 00110 00111 free_2DMatrix(dp3_Value, rowCount); 00112 dp3_Value=NULL; 00113 00114 delete dp3_Seed; 00115 dp3_Seed = NULL; 00116 00117 delete ip1_SeedRowCount; 00118 ip1_SeedRowCount=NULL; 00119 00120 delete ip1_SeedColumnCount; 00121 ip1_SeedColumnCount = NULL; 00122 00123 free_2DMatrix(dp3_CompressedMatrix, rowCount); 00124 dp3_CompressedMatrix = NULL; 00125 00126 delete ip1_ColorCount; 00127 ip1_ColorCount = NULL; 00128 00129 delete hr; 00130 hr = NULL; 00131 00132 delete dp3_NewValue; 00133 dp3_NewValue=NULL; 00134 00135 delete g; 00136 g=NULL; 00137 00138 return 0; 00139 }