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: a matrix 00013 double*** dp3_NewValue; 00014 //*/ 00015 00016 #include "ColPackHeaders.h" 00017 00018 using namespace ColPack; 00019 using namespace std; 00020 00021 #ifndef TOP_DIR 00022 #define TOP_DIR "." 00023 #endif 00024 00025 // baseDir should point to the directory (folder) containing the input file 00026 string baseDir=TOP_DIR; 00027 00028 #include "extra.h" //This .h file contains functions that are used in the below examples 00029 00030 int main() 00031 { 00032 // s_InputFile = baseDir + <name of the input file> 00033 string s_InputFile; //path of the input file 00034 s_InputFile = baseDir; 00035 s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "column-compress.mtx"; 00036 //s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "hess_pat.mtx"; 00037 00038 // Step 1: Determine sparsity structure of the Jacobian. 00039 // This step is done by an AD tool. For the purpose of illustration here, we read the structure from a file, 00040 // and store the structure in a Compressed Row Format and then ADIC format. 00041 unsigned int *** uip3_SparsityPattern = new unsigned int **; //uip3_ means triple pointers of type unsigned int 00042 double*** dp3_Value = new double**; //dp3_ means triple pointers of type double. Other prefixes follow the same notation 00043 int rowCount, columnCount; 00044 ConvertMatrixMarketFormat2RowCompressedFormat(s_InputFile, uip3_SparsityPattern, dp3_Value,rowCount, columnCount); 00045 00046 cout<<"just for debugging purpose, display the 2 matrices: the matrix with SparsityPattern only and the matrix with Value"<<endl; 00047 cout<<"Matrix rowCount = "<<rowCount<<"; columnCount = "<<columnCount<<endl; 00048 cout<<fixed<<showpoint<<setprecision(2); //formatting output 00049 cout<<"(*uip3_SparsityPattern)"<<endl; 00050 displayCompressedRowMatrix((*uip3_SparsityPattern),rowCount, true); 00051 cout<<"(*dp3_Value)"<<endl; 00052 displayCompressedRowMatrix((*dp3_Value),rowCount); 00053 cout<<"Finish ConvertMatrixMarketFormat2RowCompressedFormat()"<<endl; 00054 Pause(); 00055 00056 std::list<std::set<int> > lsi_SparsityPattern; 00057 std::list<std::vector<double> > lvd_Value; 00058 ConvertRowCompressedFormat2ADIC( (*uip3_SparsityPattern) , rowCount, (*dp3_Value), lsi_SparsityPattern, lvd_Value); 00059 00060 cout<<"just for debugging purpose, display the matrix in ADIC format rowCount = "<<rowCount<<endl; 00061 cout<<"Display lsi_SparsityPattern"<<endl; 00062 DisplayADICFormat_Sparsity(lsi_SparsityPattern); 00063 cout<<"Display lvd_Value"<<endl; 00064 DisplayADICFormat_Value(lvd_Value); 00065 cout<<"Finish ConvertRowCompressedFormat2CSR()"<<endl; 00066 Pause(); 00067 00068 //Step 2: Coloring. 00069 int *ip1_ColorCount = new int; //The number of distinct colors used to color the graph 00070 00071 //Step 2.1: Read the sparsity pattern of the given Jacobian matrix (ADIC format) 00072 //and create the corresponding bipartite graph 00073 BipartiteGraphPartialColoringInterface *g = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADIC, &lsi_SparsityPattern, columnCount); 00074 00075 //Step 2.2: Do Partial-Distance-Two-Coloring the bipartite graph with the specified ordering 00076 g->PartialDistanceTwoColoring("SMALLEST_LAST", "COLUMN_PARTIAL_DISTANCE_TWO"); 00077 00078 //Step 2.3: From the coloring information, you can get the vector of colorIDs of left or right vertices (depend on the s_ColoringVariant that you choose) 00079 vector<int> vi_VertexPartialColors; 00080 g->GetVertexPartialColors(vi_VertexPartialColors); 00081 *ip1_ColorCount = g->GetRightVertexColorCount(); 00082 cout<<"Finish GetVertexPartialColors()"<<endl; 00083 00084 //Display results of step 2 00085 printf(" Display vi_VertexPartialColors *ip1_ColorCount=%d \n",*ip1_ColorCount); 00086 displayVector(vi_VertexPartialColors); 00087 Pause(); 00088 00089 // Step 3: Obtain the Jacobian-seed matrix product. 00090 // This step will also be done by an AD tool. For the purpose of illustration here, the orginial matrix V 00091 // (for Values) is multiplied with the seed matrix S (represented as a vector of colors vi_VertexPartialColors). 00092 // The resulting matrix is stored in dp3_CompressedMatrix. 00093 double*** dp3_CompressedMatrix = new double**; 00094 cout<<"Start MatrixMultiplication()"<<endl; 00095 MatrixMultiplication_VxS__usingVertexPartialColors(lsi_SparsityPattern, lvd_Value, columnCount, vi_VertexPartialColors, *ip1_ColorCount, dp3_CompressedMatrix); 00096 cout<<"Finish MatrixMultiplication()"<<endl; 00097 00098 displayMatrix(*dp3_CompressedMatrix,rowCount,*ip1_ColorCount); 00099 Pause(); 00100 00101 //Step 4: Recover the numerical values of the original matrix from the compressed representation. 00102 // The new values are store in "lvd_NewValue" 00103 std::list<std::vector<double> > lvd_NewValue; 00104 JacobianRecovery1D* jr1d = new JacobianRecovery1D; 00105 jr1d->RecoverD2Cln_ADICFormat(g, *dp3_CompressedMatrix, lsi_SparsityPattern, lvd_NewValue); 00106 cout<<"Finish Recover()"<<endl; 00107 00108 DisplayADICFormat_Value(lvd_NewValue); 00109 Pause(); 00110 00111 //Check for consistency, make sure the values in the 2 matrices are the same. 00112 if (ADICMatricesAreEqual(lvd_Value, lvd_NewValue,0)) cout<< "lvd_Value == lvd_NewValue"<<endl; 00113 else cout<< "lvd_Value != lvd_NewValue"<<endl; 00114 00115 Pause(); 00116 00117 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00118 00119 free_2DMatrix(uip3_SparsityPattern, rowCount); 00120 uip3_SparsityPattern=NULL; 00121 00122 free_2DMatrix(dp3_Value, rowCount); 00123 dp3_Value=NULL; 00124 00125 free_2DMatrix(dp3_CompressedMatrix, rowCount); 00126 dp3_CompressedMatrix = NULL; 00127 00128 delete ip1_ColorCount; 00129 ip1_ColorCount = NULL; 00130 00131 delete jr1d; 00132 jr1d = NULL; 00133 00134 delete g; 00135 g=NULL; 00136 00137 return 0; 00138 }