ColPack
|
00001 // An example of Bidirectional compression and recovery for Jacobian using Star Bicoloring 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 #include "ColPackHeaders.h" 00014 00015 using namespace ColPack; 00016 using namespace std; 00017 00018 #ifndef TOP_DIR 00019 #define TOP_DIR "." 00020 #endif 00021 00022 // baseDir should point to the directory (folder) containing the input file 00023 string baseDir=TOP_DIR; 00024 00025 #include "extra.h" //This .h file contains functions that are used in the below examples: 00026 //ReadMM(), MatrixMultiplication...(), Times2Plus1point5(), displayMatrix() and displayCompressedRowMatrix() 00027 00028 int main() 00029 { 00030 // s_InputFile = baseDir + <name of the input file> 00031 string s_InputFile; //path of the input file 00032 s_InputFile = baseDir; 00033 s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "column-compress.mtx"; 00034 00035 // Step 1: Determine sparsity structure of the Jacobian. 00036 // This step is done by an AD tool. For the purpose of illustration here, we read the structure from a file, 00037 // and store the structure in a Compressed Row Format. 00038 unsigned int *** uip3_SparsityPattern = new unsigned int **; //uip3_ means triple pointers of type unsigned int 00039 double*** dp3_Value = new double**; //dp3_ means triple pointers of type double. Other prefixes follow the same notation 00040 int rowCount, columnCount; 00041 ConvertMatrixMarketFormat2RowCompressedFormat(s_InputFile, uip3_SparsityPattern, dp3_Value,rowCount, columnCount); 00042 00043 cout<<"just for debugging purpose, display the 2 matrices: the matrix with SparsityPattern only and the matrix with Value"<<endl; 00044 cout<<fixed<<showpoint<<setprecision(2); //formatting output 00045 cout<<"(*uip3_SparsityPattern)"<<endl; 00046 displayCompressedRowMatrix((*uip3_SparsityPattern),rowCount); 00047 cout<<"(*dp3_Value)"<<endl; 00048 displayCompressedRowMatrix((*dp3_Value),rowCount); 00049 cout<<"Finish ConvertMatrixMarketFormat2RowCompressedFormat()"<<endl; 00050 Pause(); 00051 00052 //Step 2: Obtain the seed matrices via Star Bicoloring. 00053 double*** dp3_LeftSeed = new double**; 00054 int *ip1_LeftSeedRowCount = new int; 00055 int *ip1_LeftSeedColumnCount = new int; 00056 double*** dp3_RightSeed = new double**; 00057 int *ip1_RightSeedRowCount = new int; 00058 int *ip1_RightSeedColumnCount = new int; 00059 00060 //Step 2.1: Read the sparsity pattern of the given Jacobian matrix (compressed sparse rows format) 00061 //and create the corresponding bipartite graph 00062 BipartiteGraphBicoloringInterface *g = new BipartiteGraphBicoloringInterface(SRC_MEM_ADOLC, *uip3_SparsityPattern, rowCount, columnCount); 00063 00064 //Step 2.2: Color the graph based on the specified ordering and (Star) Bicoloring 00065 g->Bicoloring("SMALLEST_LAST", "IMPLICIT_COVERING__STAR_BICOLORING"); 00066 00067 //Step 2.3 (Option 1): From the coloring information, create and return the Left and Right seed matrices 00068 (*dp3_LeftSeed) = g->GetLeftSeedMatrix(ip1_LeftSeedRowCount, ip1_LeftSeedColumnCount ); 00069 (*dp3_RightSeed) = g->GetRightSeedMatrix(ip1_RightSeedRowCount, ip1_RightSeedColumnCount ); 00070 /* Notes: 00071 Step 2.3 (Option 2): From the coloring information, you can also get the vector of colorIDs of left and right vertices 00072 vector<int> vi_LeftVertexColors; 00073 g->GetLeftVertexColors(vi_LeftVertexColors); 00074 vector<int> RightVertexColors; 00075 g->GetRightVertexColors_Transformed(RightVertexColors); 00076 */ 00077 cout<<"Finish GenerateSeed()"<<endl; 00078 00079 //Display results of step 2 00080 printf(" dp3_LeftSeed %d x %d", *ip1_LeftSeedRowCount, *ip1_LeftSeedColumnCount); 00081 displayMatrix(*dp3_LeftSeed, *ip1_LeftSeedRowCount, *ip1_LeftSeedColumnCount); 00082 printf(" dp3_RightSeed %d x %d", *ip1_RightSeedRowCount, *ip1_RightSeedColumnCount); 00083 displayMatrix(*dp3_RightSeed, *ip1_RightSeedRowCount, *ip1_RightSeedColumnCount); 00084 Pause(); 00085 00086 // Step 3: Obtain the Jacobian-RightSeed and LeftSeed-Jacobian matrix products. 00087 // This step will also be done by an AD tool. For the purpose of illustration here: 00088 // - The left seed matrix LS is multiplied with the orginial matrix V (for Values). 00089 // The resulting matrix is stored in dp3_LeftCompressedMatrix. 00090 // - The orginial matrix V (for Values) is multiplied with the right seed matrix RS. 00091 // The resulting matrix is stored in dp3_RightCompressedMatrix. 00092 double*** dp3_LeftCompressedMatrix = new double**; 00093 double*** dp3_RightCompressedMatrix = new double**; 00094 cout<<"Start MatrixMultiplication() for both direction (left and right)"<<endl; 00095 MatrixMultiplication_SxV(*uip3_SparsityPattern, *dp3_Value, rowCount, columnCount, *dp3_LeftSeed, *ip1_LeftSeedRowCount, dp3_LeftCompressedMatrix); 00096 MatrixMultiplication_VxS(*uip3_SparsityPattern, *dp3_Value, rowCount, columnCount, *dp3_RightSeed, *ip1_RightSeedColumnCount, dp3_RightCompressedMatrix); 00097 cout<<"Finish MatrixMultiplication()"<<endl; 00098 00099 displayMatrix(*dp3_RightCompressedMatrix,rowCount,*ip1_RightSeedColumnCount); 00100 displayMatrix(*dp3_LeftCompressedMatrix,*ip1_LeftSeedRowCount, columnCount); 00101 Pause(); 00102 00103 //Step 4: Recover the numerical values of the original matrix from the compressed representation. 00104 // The new values are store in "dp3_NewValue" 00105 double*** dp3_NewValue = new double**; 00106 JacobianRecovery2D* jr2d = new JacobianRecovery2D; 00107 jr2d->DirectRecover_RowCompressedFormat(g, *dp3_LeftCompressedMatrix, *dp3_RightCompressedMatrix, *uip3_SparsityPattern, dp3_NewValue); 00108 cout<<"Finish Recover()"<<endl; 00109 00110 displayCompressedRowMatrix(*dp3_NewValue,rowCount); 00111 Pause(); 00112 00113 //Check for consistency, make sure the values in the 2 matrices are the same. 00114 if (CompressedRowMatricesAreEqual(*dp3_Value, *dp3_NewValue, rowCount,0)) cout<< "*dp3_Value == dp3_NewValue"<<endl; 00115 else cout<< "*dp3_Value != dp3_NewValue"<<endl; 00116 00117 Pause(); 00118 00119 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00120 00121 delete jr2d; 00122 jr2d = NULL; 00123 00124 delete dp3_NewValue; 00125 dp3_NewValue = NULL; 00126 00127 free_2DMatrix(dp3_RightCompressedMatrix, rowCount); 00128 dp3_RightCompressedMatrix = NULL; 00129 00130 free_2DMatrix(dp3_LeftCompressedMatrix, *ip1_LeftSeedRowCount); 00131 dp3_LeftCompressedMatrix = NULL; 00132 00133 delete dp3_RightSeed; 00134 dp3_RightSeed = NULL; 00135 00136 delete ip1_RightSeedColumnCount; 00137 ip1_RightSeedColumnCount = NULL; 00138 00139 delete ip1_RightSeedRowCount; 00140 ip1_RightSeedRowCount = NULL; 00141 00142 delete dp3_LeftSeed; 00143 dp3_LeftSeed = NULL; 00144 00145 delete ip1_LeftSeedColumnCount; 00146 ip1_LeftSeedColumnCount = NULL; 00147 00148 delete ip1_LeftSeedRowCount; 00149 ip1_LeftSeedRowCount = NULL; 00150 00151 free_2DMatrix(dp3_Value, rowCount); 00152 dp3_Value = NULL; 00153 00154 free_2DMatrix(uip3_SparsityPattern, rowCount); 00155 dp3_Value = NULL; 00156 00157 delete g; 00158 g=NULL; 00159 00160 return 0; 00161 }