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