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 //ReadMM(), MatrixMultiplication...(), Times2Plus1point5(), displayMatrix() and displayCompressedRowMatrix() 00030 00031 int main() 00032 { 00033 // s_InputFile = baseDir + <name of the input file> 00034 string s_InputFile; //path of the input file 00035 s_InputFile = baseDir; 00036 s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "column-compress.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. 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<<fixed<<showpoint<<setprecision(2); //formatting output 00048 cout<<"(*uip3_SparsityPattern)"<<endl; 00049 displayCompressedRowMatrix((*uip3_SparsityPattern),rowCount); 00050 cout<<"(*dp3_Value)"<<endl; 00051 displayCompressedRowMatrix((*dp3_Value),rowCount); 00052 cout<<"Finish ConvertMatrixMarketFormat2RowCompressedFormat()"<<endl; 00053 Pause(); 00054 00055 //Step 2: Obtain the seed matrix via coloring. 00056 double*** dp3_Seed = new double**; 00057 int *ip1_SeedRowCount = new int; 00058 int *ip1_SeedColumnCount = new int; 00059 int *ip1_ColorCount = new int; //The number of distinct colors used to color the graph 00060 00061 //Step 2.1: Read the sparsity pattern of the given Jacobian matrix (compressed sparse rows format) 00062 //and create the corresponding bipartite graph 00063 BipartiteGraphPartialColoringInterface *g = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, *uip3_SparsityPattern, rowCount, columnCount); 00064 00065 //Step 2.2: Do Partial-Distance-Two-Coloring the bipartite graph with the specified ordering 00066 g->PartialDistanceTwoColoring("SMALLEST_LAST", "COLUMN_PARTIAL_DISTANCE_TWO"); 00067 00068 //Step 2.3 (Option 1): From the coloring information, create and return the seed matrix 00069 (*dp3_Seed) = g->GetSeedMatrix(ip1_SeedRowCount, ip1_SeedColumnCount); 00070 /* Notes: 00071 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) 00072 vector<int> vi_VertexPartialColors; 00073 g->GetVertexPartialColors(vi_VertexPartialColors); 00074 */ 00075 cout<<"Finish GenerateSeed()"<<endl; 00076 *ip1_ColorCount = *ip1_SeedColumnCount; 00077 00078 //Display results of step 2 00079 printf(" dp3_Seed %d x %d \n", *ip1_SeedRowCount, *ip1_SeedColumnCount); 00080 displayMatrix(*dp3_Seed, *ip1_SeedRowCount, *ip1_SeedColumnCount); 00081 Pause(); 00082 00083 // Step 3: Obtain the Jacobian-seed matrix product. 00084 // This step will also be done by an AD tool. For the purpose of illustration here, the orginial matrix V 00085 // (for Values) is multiplied with the seed matrix S. The resulting matrix is stored in dp3_CompressedMatrix. 00086 double*** dp3_CompressedMatrix = new double**; 00087 cout<<"Start MatrixMultiplication()"<<endl; 00088 MatrixMultiplication_VxS(*uip3_SparsityPattern, *dp3_Value, rowCount, columnCount, *dp3_Seed, *ip1_ColorCount, dp3_CompressedMatrix); 00089 cout<<"Finish MatrixMultiplication()"<<endl; 00090 00091 displayMatrix(*dp3_CompressedMatrix,rowCount,*ip1_ColorCount); 00092 Pause(); 00093 00094 //Step 4: Recover the numerical values of the original matrix from the compressed representation. 00095 // The new values are store in "dp3_NewValue" 00096 double*** dp3_NewValue = new double**; 00097 JacobianRecovery1D* jr1d = new JacobianRecovery1D; 00098 int rowCount_for_dp3_NewValue = jr1d->RecoverD2Cln_RowCompressedFormat_unmanaged(g, *dp3_CompressedMatrix, *uip3_SparsityPattern, dp3_NewValue); 00099 /* RecoverD2Cln_RowCompressedFormat_unmanaged is called instead of RecoverD2Cln_RowCompressedFormat so that 00100 we could manage the memory deallocation for dp3_NewValue by ourselves. 00101 This way, we can reuse (*dp3_NewValue) to store new values if RecoverD2Cln_RowCompressedFormat...() need to be called again. 00102 //*/ 00103 00104 cout<<"Finish Recover()"<<endl; 00105 00106 displayCompressedRowMatrix(*dp3_NewValue,rowCount); 00107 Pause(); 00108 00109 //Check for consistency, make sure the values in the 2 matrices are the same. 00110 if (CompressedRowMatricesAreEqual(*dp3_Value, *dp3_NewValue, rowCount,0)) cout<< "*dp3_Value == dp3_NewValue"<<endl; 00111 else cout<< "*dp3_Value != dp3_NewValue"<<endl; 00112 00113 Pause(); 00114 00115 /* Let say that we have new matrices with the same sparsity structure (only the values changed), 00116 We can take advantage of the memory already allocated to (*dp3_NewValue) and use (*dp3_NewValue) to stored the new values 00117 by calling RecoverD2Cln_RowCompressedFormat_usermem recovery function. 00118 This function works in the same way as RecoverD2Cln_RowCompressedFormat_unmanaged except that the memory for (*dp3_NewValue) 00119 is reused and therefore, takes less time than the _unmanaged counterpart. 00120 //*/ 00121 for(int i=0; i<3;i++) { 00122 jr1d->RecoverD2Cln_RowCompressedFormat_usermem(g, *dp3_CompressedMatrix, *uip3_SparsityPattern, dp3_NewValue); 00123 } 00124 00125 //Deallocate memory for 2-dimensional array (*dp3_NewValue) 00126 for(int i=0; i<rowCount_for_dp3_NewValue;i++) { 00127 free((*dp3_NewValue)[i]); 00128 } 00129 free(*dp3_NewValue); 00130 00131 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00132 00133 free_2DMatrix(uip3_SparsityPattern, rowCount); 00134 uip3_SparsityPattern=NULL; 00135 00136 free_2DMatrix(dp3_Value, rowCount); 00137 dp3_Value=NULL; 00138 00139 delete dp3_Seed; 00140 dp3_Seed = NULL; 00141 00142 delete ip1_SeedRowCount; 00143 ip1_SeedRowCount=NULL; 00144 00145 delete ip1_SeedColumnCount; 00146 ip1_SeedColumnCount = NULL; 00147 00148 free_2DMatrix(dp3_CompressedMatrix, rowCount); 00149 dp3_CompressedMatrix = NULL; 00150 00151 delete ip1_ColorCount; 00152 ip1_ColorCount = NULL; 00153 00154 delete jr1d; 00155 jr1d = NULL; 00156 00157 delete dp3_NewValue; 00158 dp3_NewValue=NULL; 00159 00160 delete g; 00161 g=NULL; 00162 00163 return 0; 00164 }