ColPack
|
00001 // An example of Column compression and recovery for Jacobian 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 00035 int main() 00036 { 00037 // s_InputFile = baseDir + <name of the input file> 00038 string s_InputFile; //path of the input file 00039 s_InputFile = baseDir; 00040 s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "column-compress.mtx"; 00041 00042 // Step 1: Determine sparsity structure of the Jacobian. 00043 // This step is done by an AD tool. For the purpose of illustration here, we read the structure from a file, 00044 // and store the structure in a Compressed Row Format. 00045 unsigned int *** uip3_SparsityPattern = new unsigned int **; //uip3_ means triple pointers of type unsigned int 00046 double*** dp3_Value = new double**; //dp3_ means triple pointers of type double. Other prefixes follow the same notation 00047 int rowCount, columnCount; 00048 ConvertMatrixMarketFormat2RowCompressedFormat(s_InputFile, uip3_SparsityPattern, dp3_Value,rowCount, columnCount); 00049 00050 cout<<"just for debugging purpose, display the 2 matrices: the matrix with SparsityPattern only and the matrix with Value"<<endl; 00051 cout<<fixed<<showpoint<<setprecision(2); //formatting output 00052 cout<<"(*uip3_SparsityPattern)"<<endl; 00053 displayCompressedRowMatrix((*uip3_SparsityPattern),rowCount); 00054 cout<<"(*dp3_Value)"<<endl; 00055 displayCompressedRowMatrix((*dp3_Value),rowCount); 00056 cout<<"Finish ConvertMatrixMarketFormat2RowCompressedFormat()"<<endl; 00057 Pause(); 00058 00059 //Step 2: Obtain the seed matrix via coloring. 00060 double*** dp3_Seed = new double**; 00061 int *ip1_SeedRowCount = new int; 00062 int *ip1_SeedColumnCount = new int; 00063 int *ip1_ColorCount = new int; //The number of distinct colors used to color the graph 00064 00065 //Step 2.1: Read the sparsity pattern of the given Jacobian matrix (compressed sparse rows format) 00066 //and create the corresponding bipartite graph 00067 BipartiteGraphPartialColoringInterface *g = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, *uip3_SparsityPattern, rowCount, columnCount); 00068 00069 //Step 2.2: Do Partial-Distance-Two-Coloring the bipartite graph with the specified ordering 00070 g->PartialDistanceTwoColoring("SMALLEST_LAST", "COLUMN_PARTIAL_DISTANCE_TWO"); 00071 00072 //Step 2.3 (Option 1): From the coloring information, create and return the seed matrix 00073 (*dp3_Seed) = g->GetSeedMatrix(ip1_SeedRowCount, ip1_SeedColumnCount); 00074 /* Notes: 00075 Step 2.3 (Option 2): From the coloring information, you can also get the vector of colorIDs of left or right vertices (depend on the s_ColoringVariant that you choose) 00076 vector<int> vi_VertexPartialColors; 00077 g->GetVertexPartialColors(vi_VertexPartialColors); 00078 */ 00079 cout<<"Finish GenerateSeed()"<<endl; 00080 *ip1_ColorCount = *ip1_SeedColumnCount; 00081 00082 //Display results of step 2 00083 printf(" dp3_Seed %d x %d \n", *ip1_SeedRowCount, *ip1_SeedColumnCount); 00084 displayMatrix(*dp3_Seed, *ip1_SeedRowCount, *ip1_SeedColumnCount); 00085 Pause(); 00086 00087 // Step 3: Obtain the Jacobian-seed matrix product. 00088 // This step will also be done by an AD tool. For the purpose of illustration here, the orginial matrix V 00089 // (for Values) is multiplied with the seed matrix S. The resulting matrix is stored in dp3_CompressedMatrix. 00090 double*** dp3_CompressedMatrix = new double**; 00091 cout<<"Start MatrixMultiplication()"<<endl; 00092 MatrixMultiplication_VxS(*uip3_SparsityPattern, *dp3_Value, rowCount, columnCount, *dp3_Seed, *ip1_ColorCount, dp3_CompressedMatrix); 00093 cout<<"Finish MatrixMultiplication()"<<endl; 00094 00095 displayMatrix(*dp3_CompressedMatrix,rowCount,*ip1_ColorCount); 00096 Pause(); 00097 00098 //Step 4: Recover the numerical values of the original matrix from the compressed representation. 00099 // The new values are store in "dp2_JacobianValue" 00100 unsigned int** ip2_RowIndex = new unsigned int*; 00101 unsigned int** ip2_ColumnIndex = new unsigned int*; 00102 double** dp2_JacobianValue = new double*; 00103 JacobianRecovery1D* jr1d = new JacobianRecovery1D; 00104 jr1d->RecoverD2Cln_CoordinateFormat(g, *dp3_CompressedMatrix, *uip3_SparsityPattern, ip2_RowIndex, ip2_ColumnIndex, dp2_JacobianValue); 00105 cout<<"Finish Recover()"<<endl; 00106 00107 cout<<endl<<"Display result, the structure and values should be similar to the original one"<<endl; 00108 cout<<"Display *ip2_RowIndex"<<endl; 00109 displayVector(*ip2_RowIndex,g->GetEdgeCount()); 00110 cout<<"Display *ip2_ColumnIndex"<<endl; 00111 displayVector(*ip2_ColumnIndex, g->GetEdgeCount()); 00112 cout<<"Display *dp2_JacobianValue"<<endl; 00113 displayVector(*dp2_JacobianValue, g->GetEdgeCount()); 00114 Pause(); 00115 00116 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00117 00118 free_2DMatrix(uip3_SparsityPattern, rowCount); 00119 uip3_SparsityPattern=NULL; 00120 00121 free_2DMatrix(dp3_Value, rowCount); 00122 dp3_Value=NULL; 00123 00124 delete dp3_Seed; 00125 dp3_Seed = NULL; 00126 00127 delete ip1_SeedRowCount; 00128 ip1_SeedRowCount=NULL; 00129 00130 delete ip1_SeedColumnCount; 00131 ip1_SeedColumnCount = NULL; 00132 00133 free_2DMatrix(dp3_CompressedMatrix, rowCount); 00134 dp3_CompressedMatrix = NULL; 00135 00136 delete ip1_ColorCount; 00137 ip1_ColorCount = NULL; 00138 00139 delete jr1d; 00140 jr1d = NULL; 00141 00142 delete ip2_RowIndex; 00143 delete ip2_ColumnIndex; 00144 delete dp2_JacobianValue; 00145 ip2_RowIndex=NULL; 00146 ip2_ColumnIndex=NULL; 00147 dp2_JacobianValue=NULL; 00148 00149 delete g; 00150 g=NULL; 00151 00152 return 0; 00153 }