1
2 /* @cond INNERDOC */
3 /**
4 * @file
5 * @brief
6 * Auxiliary functions.
7 */
8
9 /*
10
11 Copyright (C) 2008-2019 Michele Martone
12
13 This file is part of librsb.
14
15 librsb is free software; you can redistribute it and/or modify it
16 under the terms of the GNU Lesser General Public License as published
17 by the Free Software Foundation; either version 3 of the License, or
18 (at your option) any later version.
19
20 librsb is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23 License for more details.
24
25 You should have received a copy of the GNU Lesser General Public
26 License along with librsb; see the file COPYING.
27 If not, see <http://www.gnu.org/licenses/>.
28
29 */
30 /*
31 The code in this file was generated automatically by an M4 script.
32 It is not meant to be used as an API (Application Programming Interface).
33 p.s.: right now, only row major matrix access is considered.
34
35 */
36
37
38 #ifdef __cplusplus
39 extern "C" {
40 #endif /* __cplusplus */
41
42 #define RSB_WANT_OMP 1
43 #define RSB_MAX_OMP_THREADS 4
44 #include <omp.h> /* OpenMP parallelism (EXPERIMENTAL) */
45
46
47 #include "rsb_common.h"
rsb_do_csr_ilu0_DOUBLE(struct rsb_coo_matrix_t * coop)48 rsb_err_t rsb_do_csr_ilu0_DOUBLE(struct rsb_coo_matrix_t * coop){
49 /**
50 * \ingroup gr_internals
51 FIXME: INCOMPLETE, EXPERIMENTAL, TEMPORARILY HERE
52 On exit, the matrix will contain the L and U factors of a pattern preserving incomplete LU factorization (ILU 0).
53 */
54 rsb_err_t errval = RSB_ERR_NO_ERROR;
55 rsb_coo_idx_t i;
56
57 {
58 double *VA = coop->VA;
59 const rsb_coo_idx_t *PA = coop->IA;
60 const rsb_coo_idx_t *JA = coop->JA;
61 for(i=1;i<coop->nr;++i)
62 {
63 const rsb_nnz_idx_t ifp = PA[i],ilp = PA[i+1],irnz = ilp-ifp;
64 rsb_nnz_idx_t idp = RSB_MARKER_NNZ_VALUE,ikp = RSB_MARKER_NNZ_VALUE;
65 if(irnz)
66 {
67
68 idp = rsb__nnz_split_coo_bsearch(JA+ifp,i,irnz)+ifp;
69 assert(idp<=ilp);
70 assert(idp>=ifp);
71 for(ikp=ifp;ikp<idp;++ikp)// k = 1...i-1
72 {
73 /* FIXME: write a sparse vectors dot product macro and apply it here */
74 const rsb_nnz_idx_t k = JA[ikp],kfp = PA[k],klp = PA[k+1],krnz = klp-kfp;
75 const int kdp = rsb__nnz_split_coo_bsearch(JA+kfp,k,krnz)+kfp;
76 rsb_nnz_idx_t kjp = kfp,ijp = ikp+1;
77 VA[ikp]/=VA[kdp];
78 /* FIXME: to optimize this phase, we should loop on the shorter row */
79 for(;ijp<ilp;++ijp)// j = k+1...n
80 {
81 for(;JA[kjp]<JA[ijp] && kjp<klp;++kjp)
82 ;
83 if(kjp==klp)
84 goto out;
85 /* JA[kjp]>=JA[ijp] */
86 for(;JA[kjp]>JA[ijp] && ijp<ilp;++ijp)
87 ;
88 if(ijp==ilp)
89 goto out;
90 /* JA[kjp]==JA[ijp] */
91 VA[ijp]-=VA[ikp]*VA[kjp];
92 }
93 out:
94 RSB_NULL_STATEMENT_FOR_COMPILER_HAPPINESS
95
96 }
97 }
98 }
99 }
100 RSB_DO_ERR_RETURN(errval)
101 }
rsb_do_csr_ilu0_FLOAT(struct rsb_coo_matrix_t * coop)102 rsb_err_t rsb_do_csr_ilu0_FLOAT(struct rsb_coo_matrix_t * coop){
103 /**
104 * \ingroup gr_internals
105 FIXME: INCOMPLETE, EXPERIMENTAL, TEMPORARILY HERE
106 On exit, the matrix will contain the L and U factors of a pattern preserving incomplete LU factorization (ILU 0).
107 */
108 rsb_err_t errval = RSB_ERR_NO_ERROR;
109 rsb_coo_idx_t i;
110
111 {
112 float *VA = coop->VA;
113 const rsb_coo_idx_t *PA = coop->IA;
114 const rsb_coo_idx_t *JA = coop->JA;
115 for(i=1;i<coop->nr;++i)
116 {
117 const rsb_nnz_idx_t ifp = PA[i],ilp = PA[i+1],irnz = ilp-ifp;
118 rsb_nnz_idx_t idp = RSB_MARKER_NNZ_VALUE,ikp = RSB_MARKER_NNZ_VALUE;
119 if(irnz)
120 {
121
122 idp = rsb__nnz_split_coo_bsearch(JA+ifp,i,irnz)+ifp;
123 assert(idp<=ilp);
124 assert(idp>=ifp);
125 for(ikp=ifp;ikp<idp;++ikp)// k = 1...i-1
126 {
127 /* FIXME: write a sparse vectors dot product macro and apply it here */
128 const rsb_nnz_idx_t k = JA[ikp],kfp = PA[k],klp = PA[k+1],krnz = klp-kfp;
129 const int kdp = rsb__nnz_split_coo_bsearch(JA+kfp,k,krnz)+kfp;
130 rsb_nnz_idx_t kjp = kfp,ijp = ikp+1;
131 VA[ikp]/=VA[kdp];
132 /* FIXME: to optimize this phase, we should loop on the shorter row */
133 for(;ijp<ilp;++ijp)// j = k+1...n
134 {
135 for(;JA[kjp]<JA[ijp] && kjp<klp;++kjp)
136 ;
137 if(kjp==klp)
138 goto out;
139 /* JA[kjp]>=JA[ijp] */
140 for(;JA[kjp]>JA[ijp] && ijp<ilp;++ijp)
141 ;
142 if(ijp==ilp)
143 goto out;
144 /* JA[kjp]==JA[ijp] */
145 VA[ijp]-=VA[ikp]*VA[kjp];
146 }
147 out:
148 RSB_NULL_STATEMENT_FOR_COMPILER_HAPPINESS
149
150 }
151 }
152 }
153 }
154 RSB_DO_ERR_RETURN(errval)
155 }
rsb_do_csr_ilu0_FLOAT_COMPLEX(struct rsb_coo_matrix_t * coop)156 rsb_err_t rsb_do_csr_ilu0_FLOAT_COMPLEX(struct rsb_coo_matrix_t * coop){
157 /**
158 * \ingroup gr_internals
159 FIXME: INCOMPLETE, EXPERIMENTAL, TEMPORARILY HERE
160 On exit, the matrix will contain the L and U factors of a pattern preserving incomplete LU factorization (ILU 0).
161 */
162 rsb_err_t errval = RSB_ERR_NO_ERROR;
163 rsb_coo_idx_t i;
164
165 {
166 float complex *VA = coop->VA;
167 const rsb_coo_idx_t *PA = coop->IA;
168 const rsb_coo_idx_t *JA = coop->JA;
169 for(i=1;i<coop->nr;++i)
170 {
171 const rsb_nnz_idx_t ifp = PA[i],ilp = PA[i+1],irnz = ilp-ifp;
172 rsb_nnz_idx_t idp = RSB_MARKER_NNZ_VALUE,ikp = RSB_MARKER_NNZ_VALUE;
173 if(irnz)
174 {
175
176 idp = rsb__nnz_split_coo_bsearch(JA+ifp,i,irnz)+ifp;
177 assert(idp<=ilp);
178 assert(idp>=ifp);
179 for(ikp=ifp;ikp<idp;++ikp)// k = 1...i-1
180 {
181 /* FIXME: write a sparse vectors dot product macro and apply it here */
182 const rsb_nnz_idx_t k = JA[ikp],kfp = PA[k],klp = PA[k+1],krnz = klp-kfp;
183 const int kdp = rsb__nnz_split_coo_bsearch(JA+kfp,k,krnz)+kfp;
184 rsb_nnz_idx_t kjp = kfp,ijp = ikp+1;
185 VA[ikp]/=VA[kdp];
186 /* FIXME: to optimize this phase, we should loop on the shorter row */
187 for(;ijp<ilp;++ijp)// j = k+1...n
188 {
189 for(;JA[kjp]<JA[ijp] && kjp<klp;++kjp)
190 ;
191 if(kjp==klp)
192 goto out;
193 /* JA[kjp]>=JA[ijp] */
194 for(;JA[kjp]>JA[ijp] && ijp<ilp;++ijp)
195 ;
196 if(ijp==ilp)
197 goto out;
198 /* JA[kjp]==JA[ijp] */
199 VA[ijp]-=VA[ikp]*VA[kjp];
200 }
201 out:
202 RSB_NULL_STATEMENT_FOR_COMPILER_HAPPINESS
203
204 }
205 }
206 }
207 }
208 RSB_DO_ERR_RETURN(errval)
209 }
rsb_do_csr_ilu0_DOUBLE_COMPLEX(struct rsb_coo_matrix_t * coop)210 rsb_err_t rsb_do_csr_ilu0_DOUBLE_COMPLEX(struct rsb_coo_matrix_t * coop){
211 /**
212 * \ingroup gr_internals
213 FIXME: INCOMPLETE, EXPERIMENTAL, TEMPORARILY HERE
214 On exit, the matrix will contain the L and U factors of a pattern preserving incomplete LU factorization (ILU 0).
215 */
216 rsb_err_t errval = RSB_ERR_NO_ERROR;
217 rsb_coo_idx_t i;
218
219 {
220 double complex *VA = coop->VA;
221 const rsb_coo_idx_t *PA = coop->IA;
222 const rsb_coo_idx_t *JA = coop->JA;
223 for(i=1;i<coop->nr;++i)
224 {
225 const rsb_nnz_idx_t ifp = PA[i],ilp = PA[i+1],irnz = ilp-ifp;
226 rsb_nnz_idx_t idp = RSB_MARKER_NNZ_VALUE,ikp = RSB_MARKER_NNZ_VALUE;
227 if(irnz)
228 {
229
230 idp = rsb__nnz_split_coo_bsearch(JA+ifp,i,irnz)+ifp;
231 assert(idp<=ilp);
232 assert(idp>=ifp);
233 for(ikp=ifp;ikp<idp;++ikp)// k = 1...i-1
234 {
235 /* FIXME: write a sparse vectors dot product macro and apply it here */
236 const rsb_nnz_idx_t k = JA[ikp],kfp = PA[k],klp = PA[k+1],krnz = klp-kfp;
237 const int kdp = rsb__nnz_split_coo_bsearch(JA+kfp,k,krnz)+kfp;
238 rsb_nnz_idx_t kjp = kfp,ijp = ikp+1;
239 VA[ikp]/=VA[kdp];
240 /* FIXME: to optimize this phase, we should loop on the shorter row */
241 for(;ijp<ilp;++ijp)// j = k+1...n
242 {
243 for(;JA[kjp]<JA[ijp] && kjp<klp;++kjp)
244 ;
245 if(kjp==klp)
246 goto out;
247 /* JA[kjp]>=JA[ijp] */
248 for(;JA[kjp]>JA[ijp] && ijp<ilp;++ijp)
249 ;
250 if(ijp==ilp)
251 goto out;
252 /* JA[kjp]==JA[ijp] */
253 VA[ijp]-=VA[ikp]*VA[kjp];
254 }
255 out:
256 RSB_NULL_STATEMENT_FOR_COMPILER_HAPPINESS
257
258 }
259 }
260 }
261 }
262 RSB_DO_ERR_RETURN(errval)
263 }
264
rsb__prec_ilu0(struct rsb_mtx_t * mtxAp)265 rsb_err_t rsb__prec_ilu0(struct rsb_mtx_t * mtxAp){
266 rsb_err_t errval = RSB_ERR_NO_ERROR;
267 struct rsb_coo_matrix_t coo;
268
269 if(!mtxAp || !rsb__is_terminal_recursive_matrix(mtxAp) ||
270 !rsb__is_css_matrix(mtxAp) || (mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES) ||
271 /*mtxAp->typecode != RSB_NUMERICAL_TYPE_DOUBLE || */!rsb__is_square(mtxAp) || rsb__is_symmetric(mtxAp) ||
272 RSB_DO_FLAG_HAS(mtxAp->flags,RSB_FLAG_UNIT_DIAG_IMPLICIT)
273 )
274 {
275 RSB_ERROR(RSB_ERRM_ES);
276 errval = RSB_ERR_BADARGS;
277 goto err;
278 }
279 if(mtxAp->nr==1)
280 goto err;
281 if((errval = rsb__project_rsb_to_coo(mtxAp,&coo))!=RSB_ERR_NO_ERROR)
282 goto err;
283 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
284 if( mtxAp->typecode == RSB_NUMERICAL_TYPE_DOUBLE )
285 return rsb_do_csr_ilu0_DOUBLE(&coo);
286 else
287 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
288 #ifdef RSB_NUMERICAL_TYPE_FLOAT
289 if( mtxAp->typecode == RSB_NUMERICAL_TYPE_FLOAT )
290 return rsb_do_csr_ilu0_FLOAT(&coo);
291 else
292 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
293 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
294 if( mtxAp->typecode == RSB_NUMERICAL_TYPE_FLOAT_COMPLEX )
295 return rsb_do_csr_ilu0_FLOAT_COMPLEX(&coo);
296 else
297 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
298 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
299 if( mtxAp->typecode == RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX )
300 return rsb_do_csr_ilu0_DOUBLE_COMPLEX(&coo);
301 else
302 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
303 errval = RSB_ERR_INTERNAL_ERROR;
304 err:
305 RSB_DO_ERR_RETURN(errval)
306 }
307
rsb__prec_csr_ilu0(struct rsb_coo_matrix_t * coop)308 rsb_err_t rsb__prec_csr_ilu0(struct rsb_coo_matrix_t * coop){
309 // FIXME: termporary
310 if(coop->nr==1)
311 goto err;
312 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
313 if( coop->typecode == RSB_NUMERICAL_TYPE_DOUBLE )
314 return rsb_do_csr_ilu0_DOUBLE(coop);
315 else
316 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
317 #ifdef RSB_NUMERICAL_TYPE_FLOAT
318 if( coop->typecode == RSB_NUMERICAL_TYPE_FLOAT )
319 return rsb_do_csr_ilu0_FLOAT(coop);
320 else
321 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
322 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
323 if( coop->typecode == RSB_NUMERICAL_TYPE_FLOAT_COMPLEX )
324 return rsb_do_csr_ilu0_FLOAT_COMPLEX(coop);
325 else
326 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
327 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
328 if( coop->typecode == RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX )
329 return rsb_do_csr_ilu0_DOUBLE_COMPLEX(coop);
330 else
331 #endif /* RSB_M4_NUMERICAL_TYPE_PREPROCESSOR_SYMBOL(mtype) */
332 err:
333 return RSB_ERR_INTERNAL_ERROR;
334 }
335
336 #ifdef __cplusplus
337 }
338 #endif /* __cplusplus */
339
340
341 /* @endcond */
342