AFEPack
|
00001 00006 #include "SparseMatrixTool.h" 00007 00008 namespace SparseMatrixTool { 00009 00021 template <typename number> 00022 void hCatSparseMatrix(const SparseMatrix<number>& m0, 00023 const SparseMatrix<number>& m1, 00024 SparsityPattern& sp, 00025 SparseMatrix<number>& m, 00026 bool is_pattern_ok) 00027 { 00028 const SparsityPattern& sp0 = m0.get_sparsity_pattern(); 00029 const SparsityPattern& sp1 = m1.get_sparsity_pattern(); 00030 if (! is_pattern_ok) 00031 hCatSparsityPattern(sp0, sp1, sp); 00032 m.reinit(sp); 00033 int i, j, n = sp0.n_cols(); 00034 const std::size_t * p_rowstart0 = sp0.get_rowstart_indices(); 00035 const u_int * p_column0 = sp0.get_column_numbers(); 00036 const std::size_t * p_rowstart1 = sp1.get_rowstart_indices(); 00037 const u_int * p_column1 = sp1.get_column_numbers(); 00038 for (i = 0;i < sp.n_rows();i ++) { 00039 for (j = p_rowstart0[i];j < p_rowstart0[i + 1];j ++) 00040 m.add(i, p_column0[j], m0.global_entry(j)); 00041 for (j = p_rowstart1[i];j < p_rowstart1[i + 1];j ++) 00042 m.add(i, n + p_column1[j], m1.global_entry(j)); 00043 } 00044 }; 00045 00057 template <typename number> 00058 void vCatSparseMatrix(const SparseMatrix<number>& m0, 00059 const SparseMatrix<number>& m1, 00060 SparsityPattern& sp, 00061 SparseMatrix<number>& m, 00062 bool is_pattern_ok) 00063 { 00064 const SparsityPattern& sp0 = m0.get_sparsity_pattern(); 00065 const SparsityPattern& sp1 = m1.get_sparsity_pattern(); 00066 if (! is_pattern_ok) 00067 vCatSparsityPattern(sp0, sp1, sp); 00068 m.reinit(sp); 00069 int i, j, n = sp0.n_rows(); 00070 const std::size_t * p_rowstart0 = sp0.get_rowstart_indices(); 00071 const u_int * p_column0 = sp0.get_column_numbers(); 00072 const std::size_t * p_rowstart1 = sp1.get_rowstart_indices(); 00073 const u_int * p_column1 = sp1.get_column_numbers(); 00074 for (i = 0;i < n;i ++) 00075 for (j = p_rowstart0[i];j < p_rowstart0[i + 1];j ++) 00076 m.add(i, p_column0[j], m0.global_entry(j)); 00077 for (i = 0;i < sp1.n_rows();i ++) 00078 for (j = p_rowstart1[i];j < p_rowstart1[i + 1];j ++) 00079 m.add(i + n, p_column1[j], m1.global_entry(j)); 00080 }; 00081 00102 template <typename number> 00103 void gammaCatSparseMatrix(const SparseMatrix<number>& m0, 00104 const SparseMatrix<number>& m1, 00105 SparsityPattern& sp, 00106 SparseMatrix<number>& m, 00107 bool is_pattern_ok) 00108 { 00109 const SparsityPattern& sp0 = m0.get_sparsity_pattern(); 00110 const SparsityPattern& sp1 = m1.get_sparsity_pattern(); 00111 if (! is_pattern_ok) 00112 gammaCatSparsityPattern(sp0, sp1, sp); 00113 m.reinit(sp); 00114 00115 int i, j, n = sp0.n_rows(); 00116 const std::size_t * p_rowstart0 = sp0.get_rowstart_indices(); 00117 const u_int * p_column0 = sp0.get_column_numbers(); 00118 const std::size_t * p_rowstart1 = sp1.get_rowstart_indices(); 00119 const u_int * p_column1 = sp1.get_column_numbers(); 00120 for (i = 0;i < n;i ++) { 00121 for (j = p_rowstart0[i];j < p_rowstart0[i + 1];j ++) 00122 m.add(i, p_column0[j], m0.global_entry(j)); 00123 for (j = p_rowstart1[i];j < p_rowstart1[i + 1];j ++) { 00124 m.add(i, n + p_column1[j], m1.global_entry(j)); 00125 m.add(n + p_column1[j], i, m1.global_entry(j)); 00126 } 00127 } 00128 }; 00129 00141 template <typename number> 00142 void dCatSparseMatrix(const SparseMatrix<number>& m0, 00143 const SparseMatrix<number>& m1, 00144 SparsityPattern& sp, 00145 SparseMatrix<number>& m, 00146 bool is_pattern_ok) 00147 { 00148 const SparsityPattern& sp0 = m0.get_sparsity_pattern(); 00149 const SparsityPattern& sp1 = m1.get_sparsity_pattern(); 00150 if (! is_pattern_ok) 00151 dCatSparsityPattern(sp0, sp1, sp); 00152 m.reinit(sp); 00153 00154 int i, j; 00155 int r0 = sp0.n_rows(), r1 = sp1.n_rows(); 00156 int c0 = sp0.n_cols(); 00157 const std::size_t * p_rowstart0 = sp0.get_rowstart_indices(); 00158 const u_int * p_column0 = sp0.get_column_numbers(); 00159 const std::size_t * p_rowstart1 = sp1.get_rowstart_indices(); 00160 const u_int * p_column1 = sp1.get_column_numbers(); 00161 for (i = 0;i < r0;i ++) 00162 for (j = p_rowstart0[i];j < p_rowstart0[i + 1];j ++) 00163 m.add(i, p_column0[j], m0.global_entry(j)); 00164 for (i = 0;i < r1;i ++) 00165 for (j = p_rowstart1[i];j < p_rowstart1[i + 1];j ++) 00166 m.add(r0 + i, c0 + p_column1[j], m1.global_entry(j)); 00167 }; 00168 00189 template <typename number> 00190 void fullCatSparseMatrix(const SparseMatrix<number>& m00, 00191 const SparseMatrix<number>& m01, 00192 const SparseMatrix<number>& m10, 00193 const SparseMatrix<number>& m11, 00194 SparsityPattern& sp, 00195 SparseMatrix<number>& m, 00196 bool is_pattern_ok) 00197 { 00198 SparsityPattern sp0, sp1; 00199 SparseMatrix<number> m0, m1; 00200 hCatSparseMatrix(m00, m01, sp0, m0, false); 00201 hCatSparseMatrix(m10, m11, sp1, m1, false); 00202 vCatSparseMatrix(m0, m1, sp, m, is_pattern_ok); 00203 }; 00204 00205 }; 00206