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 int rowCount_for_dp3_NewValue = jr2d->DirectRecover_RowCompressedFormat_unmanaged(g, *dp3_LeftCompressedMatrix, *dp3_RightCompressedMatrix, *uip3_SparsityPattern, dp3_NewValue); 00108 /* RecoverD2Cln_RowCompressedFormat_unmanaged is called instead of RecoverD2Cln_RowCompressedFormat so that 00109 we could manage the memory deallocation for dp3_NewValue by ourselves. 00110 This way, we can reuse (*dp3_NewValue) to store new values if RecoverD2Cln_RowCompressedFormat...() need to be called again. 00111 //*/ 00112 00113 cout<<"Finish Recover()"<<endl; 00114 00115 displayCompressedRowMatrix(*dp3_NewValue,rowCount); 00116 Pause(); 00117 00118 //Check for consistency, make sure the values in the 2 matrices are the same. 00119 if (CompressedRowMatricesAreEqual(*dp3_Value, *dp3_NewValue, rowCount,0)) cout<< "*dp3_Value == dp3_NewValue"<<endl; 00120 else cout<< "*dp3_Value != dp3_NewValue"<<endl; 00121 00122 Pause(); 00123 00124 /* Let say that we have new matrices with the same sparsity structure (only the values changed), 00125 We can take advantage of the memory already allocated to (*dp3_NewValue) and use (*dp3_NewValue) to stored the new values 00126 by calling DirectRecover_RowCompressedFormat_usermem recovery function. 00127 This function works in the same way as DirectRecover_RowCompressedFormat_unmanaged except that the memory for (*dp3_NewValue) 00128 is reused and therefore, takes less time than the _unmanaged counterpart. 00129 //*/ 00130 for(int i=0; i<3;i++) { 00131 jr2d->DirectRecover_RowCompressedFormat_usermem(g, *dp3_LeftCompressedMatrix, *dp3_RightCompressedMatrix, *uip3_SparsityPattern, dp3_NewValue); 00132 } 00133 00134 //Deallocate memory for 2-dimensional array (*dp3_NewValue) 00135 for(int i=0; i<rowCount;i++) { 00136 free((*dp3_NewValue)[i]); 00137 } 00138 free(*dp3_NewValue); 00139 00140 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00141 00142 delete jr2d; 00143 jr2d = NULL; 00144 00145 delete dp3_NewValue; 00146 dp3_NewValue = NULL; 00147 00148 free_2DMatrix(dp3_RightCompressedMatrix, rowCount); 00149 dp3_RightCompressedMatrix = NULL; 00150 00151 free_2DMatrix(dp3_LeftCompressedMatrix, *ip1_LeftSeedRowCount); 00152 dp3_LeftCompressedMatrix = NULL; 00153 00154 delete dp3_RightSeed; 00155 dp3_RightSeed = NULL; 00156 00157 delete ip1_RightSeedColumnCount; 00158 ip1_RightSeedColumnCount = NULL; 00159 00160 delete ip1_RightSeedRowCount; 00161 ip1_RightSeedRowCount = NULL; 00162 00163 delete dp3_LeftSeed; 00164 dp3_LeftSeed = NULL; 00165 00166 delete ip1_LeftSeedColumnCount; 00167 ip1_LeftSeedColumnCount = NULL; 00168 00169 delete ip1_LeftSeedRowCount; 00170 ip1_LeftSeedRowCount = NULL; 00171 00172 free_2DMatrix(dp3_Value, rowCount); 00173 dp3_Value = NULL; 00174 00175 free_2DMatrix(uip3_SparsityPattern, rowCount); 00176 dp3_Value = NULL; 00177 00178 delete g; 00179 g=NULL; 00180 00181 return 0; 00182 }