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 Return by recovery routine: three vectors in "Coordinate Format" (zero-based indexing) 00013 http://www.intel.com/software/products/mkl/docs/webhelp/appendices/mkl_appA_SMSF.html#mkl_appA_SMSF_5 00014 unsigned int** ip2_RowIndex 00015 unsigned int** ip2_ColumnIndex 00016 double** dp2_JacobianValue // corresponding non-zero values 00017 //*/ 00018 00019 #include "ColPackHeaders.h" 00020 00021 using namespace ColPack; 00022 using namespace std; 00023 00024 #ifndef TOP_DIR 00025 #define TOP_DIR "." 00026 #endif 00027 00028 // baseDir should point to the directory (folder) containing the input file 00029 string baseDir=TOP_DIR; 00030 00031 #include "extra.h" //This .h file contains functions that are used in the below examples: 00032 //ReadMM(), MatrixMultiplication...(), Times2Plus1point5(), displayMatrix() and displayCompressedRowMatrix() 00033 00034 int main() 00035 { 00036 // s_InputFile = baseDir + <name of the input file> 00037 string s_InputFile; //path of the input file 00038 s_InputFile = baseDir; 00039 s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "mtx-spear-head.mtx"; 00040 00041 // Step 1: Determine sparsity structure of the Jacobian. 00042 // This step is done by an AD tool. For the purpose of illustration here, we read the structure from a file, 00043 // and store the structure in a Compressed Row Format. 00044 unsigned int *** uip3_SparsityPattern = new unsigned int **; 00045 double*** dp3_Value = new double**; 00046 int rowCount, columnCount; 00047 ConvertMatrixMarketFormat2RowCompressedFormat(s_InputFile, uip3_SparsityPattern, dp3_Value,rowCount, columnCount); 00048 00049 cout<<"just for debugging purpose, display the 2 matrices: the matrix with SparsityPattern only and the matrix with Value"<<endl; 00050 cout<<fixed<<showpoint<<setprecision(2); //formatting output 00051 cout<<"(*uip3_SparsityPattern)"<<endl; 00052 displayCompressedRowMatrix((*uip3_SparsityPattern),rowCount); 00053 cout<<"(*dp3_Value)"<<endl; 00054 displayCompressedRowMatrix((*dp3_Value),rowCount); 00055 cout<<"Finish ConvertMatrixMarketFormat2RowCompressedFormat()"<<endl; 00056 Pause(); 00057 00058 //Step 2: Obtain the seed matrix via coloring. 00059 double*** dp3_Seed = new double**; 00060 int *ip1_SeedRowCount = new int; 00061 int *ip1_SeedColumnCount = new int; 00062 int *ip1_ColorCount = new int; //The number of distinct colors used to color the graph 00063 00064 //Step 2.1: Read the sparsity pattern of the given Hessian matrix (compressed sparse rows format) 00065 //and create the corresponding graph 00066 GraphColoringInterface *g = new GraphColoringInterface(SRC_MEM_ADOLC,*uip3_SparsityPattern, rowCount); 00067 00068 //Step 2.2: Color the bipartite graph with the specified ordering 00069 g->Coloring("SMALLEST_LAST", "STAR"); 00070 00071 //Step 2.3 (Option 1): From the coloring information, create and return the seed matrix 00072 (*dp3_Seed) = g->GetSeedMatrix(ip1_SeedRowCount, ip1_SeedColumnCount); 00073 /* Notes: 00074 Step 2.3 (Option 2): From the coloring information, you can also get the vector of colorIDs of vertices 00075 vector<int> vi_VertexColors; 00076 g->GetVertexColors(vi_VertexColors); 00077 */ 00078 cout<<"Finish GenerateSeed()"<<endl; 00079 *ip1_ColorCount = *ip1_SeedColumnCount; 00080 00081 displayMatrix(*dp3_Seed,columnCount,*ip1_SeedColumnCount); 00082 Pause(); 00083 00084 // Step 3: Obtain the Hessian-seed matrix product. 00085 // This step will also be done by an AD tool. For the purpose of illustration here, the orginial matrix V 00086 // (for Values) is multiplied with the seed matrix S. The resulting matrix is stored in dp3_CompressedMatrix. 00087 double*** dp3_CompressedMatrix = new double**; 00088 cout<<"Start MatrixMultiplication()"<<endl; 00089 MatrixMultiplication_VxS(*uip3_SparsityPattern, *dp3_Value, rowCount, columnCount, *dp3_Seed, *ip1_ColorCount, dp3_CompressedMatrix); 00090 cout<<"Finish MatrixMultiplication()"<<endl; 00091 00092 displayMatrix(*dp3_CompressedMatrix,rowCount,*ip1_ColorCount); 00093 Pause(); 00094 00095 //Step 4: Recover the numerical values of the original matrix from the compressed representation. 00096 // The new values are store in "dp2_HessianValue" 00097 unsigned int** ip2_RowIndex = new unsigned int*; 00098 unsigned int** ip2_ColumnIndex = new unsigned int*; 00099 double** dp2_HessianValue = new double*; 00100 HessianRecovery* hr = new HessianRecovery; 00101 int i_arraySize = hr->DirectRecover_CoordinateFormat(g, *dp3_CompressedMatrix, *uip3_SparsityPattern, ip2_RowIndex, ip2_ColumnIndex, dp2_HessianValue); 00102 cout<<"Finish Recover()"<<endl; 00103 00104 cout<<endl<<"Display result, the structure and values should be similar to the upper half of the original matrix"<<endl; 00105 cout<<"Display *ip2_RowIndex"<<endl; 00106 displayVector(*ip2_RowIndex,i_arraySize); 00107 cout<<"Display *ip2_ColumnIndex"<<endl; 00108 displayVector(*ip2_ColumnIndex, i_arraySize); 00109 cout<<"Display *dp2_HessianValue"<<endl; 00110 displayVector(*dp2_HessianValue, i_arraySize); 00111 00112 Pause(); 00113 00114 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00115 00116 free_2DMatrix(uip3_SparsityPattern, rowCount); 00117 uip3_SparsityPattern=NULL; 00118 00119 free_2DMatrix(dp3_Value, rowCount); 00120 dp3_Value=NULL; 00121 00122 delete dp3_Seed; 00123 dp3_Seed = NULL; 00124 00125 delete ip1_SeedRowCount; 00126 ip1_SeedRowCount=NULL; 00127 00128 delete ip1_SeedColumnCount; 00129 ip1_SeedColumnCount = NULL; 00130 00131 free_2DMatrix(dp3_CompressedMatrix, rowCount); 00132 dp3_CompressedMatrix = NULL; 00133 00134 delete ip1_ColorCount; 00135 ip1_ColorCount = NULL; 00136 00137 delete hr; 00138 hr = NULL; 00139 00140 delete ip2_RowIndex; 00141 delete ip2_ColumnIndex; 00142 delete dp2_HessianValue; 00143 ip2_RowIndex = NULL; 00144 ip2_ColumnIndex = NULL; 00145 dp2_HessianValue = NULL; 00146 00147 delete g; 00148 g=NULL; 00149 00150 return 0; 00151 }