ColPack
|
00001 // An example of Column compression and recovery for Jacobian 00002 /* 00003 This example simulates the situation where the matrix input for BipartiteGraphPartialColoringInterface() constructor in Step 2 is in CSR format 00004 (in ColPack, we also call this format SSF, Sparse Solvers Format). 00005 Step 1 in this example is used to build the sparse matrix in CSR format. In reality, the data structure in CSR format is the output of some AD tools. 00006 00007 How to compile this driver manually: 00008 Please make sure that "baseDir" point to the directory (folder) containing the input matrix file, and 00009 s_InputFile should point to the input file that you want to use 00010 To compile the code, replace the Main.cpp file in Main directory with this file 00011 and run "make" in ColPack installation directory. Make will generate "ColPack.exe" executable 00012 Run "ColPack.exe" 00013 00014 Note: If you got "symbol lookup error ... undefined symbol " 00015 Please make sure that your LD_LIBRARY_PATH contains libColPack.so 00016 00017 Return by recovery routine: a matrix 00018 double*** dp3_NewValue; 00019 //*/ 00020 00021 #include "ColPackHeaders.h" 00022 00023 using namespace ColPack; 00024 using namespace std; 00025 00026 #ifndef TOP_DIR 00027 #define TOP_DIR "." 00028 #endif 00029 00030 // baseDir should point to the directory (folder) containing the input file 00031 string baseDir=TOP_DIR; 00032 00033 #include "extra.h" //This .h file contains functions that are used in the below examples: 00034 //ReadMM(), MatrixMultiplication...(), Times2Plus1point5(), displayMatrix() and displayCompressedRowMatrix() 00035 00036 int main() 00037 { 00038 // s_InputFile = baseDir + <name of the input file> 00039 string s_InputFile; //path of the input file 00040 s_InputFile = baseDir; 00041 s_InputFile += DIR_SEPARATOR; s_InputFile += "Graphs"; s_InputFile += DIR_SEPARATOR; s_InputFile += "column-compress.mtx"; 00042 00043 // Step 1: Determine sparsity structure of the Jacobian. 00044 // This step is done by an AD tool. For the purpose of illustration here, we read the structure from a file, 00045 // and store the structure in a Compressed Row Format and then CSR format. 00046 unsigned int *** uip3_SparsityPattern = new unsigned int **; //uip3_ means triple pointers of type unsigned int 00047 double*** dp3_Value = new double**; //dp3_ means double pointers of type double. Other prefixes follow the same notation 00048 int rowCount, columnCount; 00049 ConvertMatrixMarketFormat2RowCompressedFormat(s_InputFile, uip3_SparsityPattern, dp3_Value,rowCount, columnCount); 00050 00051 cout<<"just for debugging purpose, display the 2 matrices: the matrix with SparsityPattern only and the matrix with Value"<<endl; 00052 cout<<fixed<<showpoint<<setprecision(2); //formatting output 00053 cout<<"(*uip3_SparsityPattern)"<<endl; 00054 displayCompressedRowMatrix((*uip3_SparsityPattern),rowCount); 00055 cout<<"(*dp3_Value)"<<endl; 00056 displayCompressedRowMatrix((*dp3_Value),rowCount); 00057 cout<<"Finish ConvertMatrixMarketFormat2RowCompressedFormat()"<<endl; 00058 Pause(); 00059 00060 int** ip_RowIndex = new int*; 00061 int** ip_ColumnIndex = new int*; 00062 ConvertRowCompressedFormat2CSR( (*uip3_SparsityPattern) , rowCount, ip_RowIndex, ip_ColumnIndex); 00063 00064 cout<<"just for debugging purpose, display the matrix in CSR format rowCount = "<<rowCount<<endl; 00065 cout<<"Display *ip_RowIndex"<<endl; 00066 displayVector(*ip_RowIndex,rowCount+1); 00067 cout<<"Display *ip_ColumnIndex"<<endl; 00068 displayVector(*ip_ColumnIndex, (*ip_RowIndex)[rowCount]); 00069 cout<<"Finish ConvertRowCompressedFormat2CSR()"<<endl; 00070 Pause(); 00071 00072 //Step 2: Obtain the seed matrix via coloring. 00073 double*** dp3_Seed = new double**; 00074 int *ip1_SeedRowCount = new int; 00075 int *ip1_SeedColumnCount = new int; 00076 int *ip1_ColorCount = new int; //The number of distinct colors used to color the graph 00077 00078 //Step 2.1: Read the sparsity pattern of the given Jacobian matrix (compressed sparse rows format) 00079 //and create the corresponding bipartite graph 00080 BipartiteGraphPartialColoringInterface *g = new BipartiteGraphPartialColoringInterface(SRC_MEM_CSR, *ip_RowIndex, rowCount, columnCount, *ip_ColumnIndex); 00081 00082 //Step 2.2: Do Partial-Distance-Two-Coloring the bipartite graph with the specified ordering 00083 g->PartialDistanceTwoColoring("SMALLEST_LAST", "COLUMN_PARTIAL_DISTANCE_TWO"); 00084 00085 //Step 2.3 (Option 1): From the coloring information, create and return the seed matrix 00086 (*dp3_Seed) = g->GetSeedMatrix(ip1_SeedRowCount, ip1_SeedColumnCount); 00087 /* Notes: 00088 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) 00089 vector<int> vi_VertexPartialColors; 00090 g->GetVertexPartialColors(vi_VertexPartialColors); 00091 */ 00092 cout<<"Finish GenerateSeed()"<<endl; 00093 *ip1_ColorCount = *ip1_SeedColumnCount; 00094 00095 //Display results of step 2 00096 printf(" dp3_Seed %d x %d \n", *ip1_SeedRowCount, *ip1_SeedColumnCount); 00097 displayMatrix(*dp3_Seed, *ip1_SeedRowCount, *ip1_SeedColumnCount); 00098 Pause(); 00099 00100 // Step 3: Obtain the Jacobian-seed matrix product. 00101 // This step will also be done by an AD tool. For the purpose of illustration here, the orginial matrix V 00102 // (for Values) is multiplied with the seed matrix S. The resulting matrix is stored in dp3_CompressedMatrix. 00103 double*** dp3_CompressedMatrix = new double**; 00104 cout<<"Start MatrixMultiplication()"<<endl; 00105 MatrixMultiplication_VxS(*uip3_SparsityPattern, *dp3_Value, rowCount, columnCount, *dp3_Seed, *ip1_ColorCount, dp3_CompressedMatrix); 00106 cout<<"Finish MatrixMultiplication()"<<endl; 00107 00108 displayMatrix(*dp3_CompressedMatrix,rowCount,*ip1_ColorCount); 00109 Pause(); 00110 00111 //Step 4: Recover the numerical values of the original matrix from the compressed representation. 00112 // The new values are store in "dp3_NewValue" 00113 double*** dp3_NewValue = new double**; 00114 JacobianRecovery1D* jr1d = new JacobianRecovery1D; 00115 jr1d->RecoverD2Cln_RowCompressedFormat(g, *dp3_CompressedMatrix, *uip3_SparsityPattern, dp3_NewValue); 00116 cout<<"Finish Recover()"<<endl; 00117 00118 displayCompressedRowMatrix(*dp3_NewValue,rowCount); 00119 Pause(); 00120 00121 //Check for consistency, make sure the values in the 2 matrices are the same. 00122 if (CompressedRowMatricesAreEqual(*dp3_Value, *dp3_NewValue, rowCount,0)) cout<< "*dp3_Value == dp3_NewValue"<<endl; 00123 else cout<< "*dp3_Value != dp3_NewValue"<<endl; 00124 00125 Pause(); 00126 00127 //Deallocate memory using functions in Utilities/MatrixDeallocation.h 00128 00129 free_2DMatrix(uip3_SparsityPattern, rowCount); 00130 uip3_SparsityPattern=NULL; 00131 00132 free_2DMatrix(dp3_Value, rowCount); 00133 dp3_Value=NULL; 00134 00135 delete dp3_Seed; 00136 dp3_Seed = NULL; 00137 00138 delete ip1_SeedRowCount; 00139 ip1_SeedRowCount=NULL; 00140 00141 delete ip1_SeedColumnCount; 00142 ip1_SeedColumnCount = NULL; 00143 00144 free_2DMatrix(dp3_CompressedMatrix, rowCount); 00145 dp3_CompressedMatrix = NULL; 00146 00147 delete ip1_ColorCount; 00148 ip1_ColorCount = NULL; 00149 00150 delete jr1d; 00151 jr1d = NULL; 00152 00153 delete dp3_NewValue; 00154 dp3_NewValue=NULL; 00155 00156 delete g; 00157 g=NULL; 00158 00159 delete[] (*ip_RowIndex); 00160 delete ip_RowIndex; 00161 delete[] (*ip_ColumnIndex); 00162 delete ip_ColumnIndex; 00163 00164 return 0; 00165 }