SCM Repository
Annotation of /pkg/Matrix/src/Mutils.h
Parent Directory
|
Revision Log
Revision 1200 -
(view)
(download)
(as text)
Original Path: pkg/src/Mutils.h
1 : | bates | 10 | #ifndef MATRIX_MUTILS_H |
2 : | #define MATRIX_MUTILS_H | ||
3 : | |||
4 : | bates | 582 | #ifdef __cplusplus |
5 : | extern "C" { | ||
6 : | #endif | ||
7 : | maechler | 890 | |
8 : | bates | 10 | #include <Rdefines.h> |
9 : | bates | 268 | #include <Rconfig.h> |
10 : | bates | 582 | #include <R.h> /* to include Rconfig.h */ |
11 : | maechler | 890 | |
12 : | bates | 582 | #ifdef ENABLE_NLS |
13 : | #include <libintl.h> | ||
14 : | #define _(String) dgettext ("Matrix", String) | ||
15 : | #else | ||
16 : | #define _(String) (String) | ||
17 : | #endif | ||
18 : | maechler | 890 | |
19 : | SEXP triangularMatrix_validate(SEXP obj); | ||
20 : | SEXP symmetricMatrix_validate(SEXP obj); | ||
21 : | maechler | 1164 | SEXP dense_nonpacked_validate(SEXP obj); |
22 : | maechler | 890 | |
23 : | bates | 582 | /* enum constants from cblas.h and some short forms */ |
24 : | enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102}; | ||
25 : | enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113}; | ||
26 : | enum CBLAS_UPLO {CblasUpper=121, CblasLower=122}; | ||
27 : | enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132}; | ||
28 : | enum CBLAS_SIDE {CblasLeft=141, CblasRight=142}; | ||
29 : | bates | 447 | #define RMJ CblasRowMajor |
30 : | #define CMJ CblasColMajor | ||
31 : | #define NTR CblasNoTrans | ||
32 : | #define TRN CblasTrans | ||
33 : | #define CTR CblasConjTrans | ||
34 : | #define UPP CblasUpper | ||
35 : | #define LOW CblasLower | ||
36 : | #define NUN CblasNonUnit | ||
37 : | #define UNT CblasUnit | ||
38 : | #define LFT CblasLeft | ||
39 : | #define RGT CblasRight | ||
40 : | |||
41 : | bates | 10 | char norm_type(char *typstr); |
42 : | char rcond_type(char *typstr); | ||
43 : | double get_double_by_name(SEXP obj, char *nm); | ||
44 : | SEXP set_double_by_name(SEXP obj, double val, char *nm); | ||
45 : | SEXP as_det_obj(double val, int log, int sign); | ||
46 : | bates | 476 | SEXP get_factors(SEXP obj, char *nm); |
47 : | SEXP set_factors(SEXP obj, SEXP val, char *nm); | ||
48 : | bates | 478 | SEXP dgCMatrix_set_Dim(SEXP x, int nrow); |
49 : | maechler | 943 | char uplo_value(SEXP x); |
50 : | char diag_value(SEXP x); | ||
51 : | |||
52 : | bates | 10 | int csc_unsorted_columns(int ncol, const int p[], const int i[]); |
53 : | void csc_sort_columns(int ncol, const int p[], int i[], double x[]); | ||
54 : | SEXP triple_as_SEXP(int nrow, int ncol, int nz, | ||
55 : | const int Ti [], const int Tj [], const double Tx [], | ||
56 : | char *Rclass); | ||
57 : | SEXP csc_check_column_sorting(SEXP A); | ||
58 : | bates | 588 | void csc_compTr(int m, int n, int nnz, |
59 : | const int xp[], const int xi[], const double xx[], | ||
60 : | int ap[], int ai[], double ax[]); | ||
61 : | bates | 10 | void ssc_symbolic_permute(int n, int upper, const int perm[], |
62 : | int Ap[], int Ai[]); | ||
63 : | bates | 493 | SEXP Matrix_make_named(int TYP, char **names); |
64 : | bates | 592 | SEXP check_scalar_string(SEXP sP, char *vals, char *nm); |
65 : | bates | 597 | double *packed_getDiag(double *dest, SEXP x); |
66 : | bates | 862 | SEXP Matrix_getElement(SEXP list, char *nm); |
67 : | bates | 592 | |
68 : | maechler | 952 | #define PACKED_TO_FULL(TYPE) \ |
69 : | TYPE *packed_to_full_ ## TYPE(TYPE *dest, const TYPE *src, \ | ||
70 : | int n, enum CBLAS_UPLO uplo) | ||
71 : | PACKED_TO_FULL(double); | ||
72 : | PACKED_TO_FULL(int); | ||
73 : | #undef PACKED_TO_FULL | ||
74 : | bates | 738 | |
75 : | maechler | 952 | #define FULL_TO_PACKED(TYPE) \ |
76 : | TYPE *full_to_packed_ ## TYPE(TYPE *dest, const TYPE *src, int n, \ | ||
77 : | enum CBLAS_UPLO uplo, enum CBLAS_DIAG diag) | ||
78 : | FULL_TO_PACKED(double); | ||
79 : | FULL_TO_PACKED(int); | ||
80 : | #undef FULL_TO_PACKED | ||
81 : | |||
82 : | |||
83 : | bates | 592 | extern /* stored pointers to symbols initialized in R_init_Matrix */ |
84 : | bates | 329 | #include "Syms.h" |
85 : | bates | 10 | |
86 : | bates | 432 | /* zero an array */ |
87 : | bates | 441 | #define AZERO(x, n) {int _I_, _SZ_ = (n); for(_I_ = 0; _I_ < _SZ_; _I_++) (x)[_I_] = 0;} |
88 : | bates | 432 | |
89 : | bates | 597 | /* number of elements in one triangle of a square matrix of order n */ |
90 : | #define PACKED_LENGTH(n) ((n) * ((n) + 1))/2 | ||
91 : | |||
92 : | bates | 738 | /* duplicate the slot with name given by sym from src to dest */ |
93 : | #define slot_dup(dest, src, sym) SET_SLOT(dest, sym, duplicate(GET_SLOT(src, sym))) | ||
94 : | |||
95 : | maechler | 951 | #define uplo_P(_x_) CHAR(STRING_ELT(GET_SLOT(_x_, Matrix_uploSym), 0)) |
96 : | #define diag_P(_x_) CHAR(STRING_ELT(GET_SLOT(_x_, Matrix_diagSym), 0)) | ||
97 : | |||
98 : | |||
99 : | maechler | 890 | /** |
100 : | bates | 597 | * Check for valid length of a packed triangular array and return the |
101 : | * corresponding number of columns | ||
102 : | maechler | 890 | * |
103 : | bates | 597 | * @param len length of a packed triangular array |
104 : | maechler | 890 | * |
105 : | bates | 597 | * @return number of columns |
106 : | */ | ||
107 : | static R_INLINE | ||
108 : | maechler | 890 | int packed_ncol(int len) |
109 : | bates | 597 | { |
110 : | int disc = 8 * len + 1; /* discriminant */ | ||
111 : | int sqrtd = (int) sqrt((double) disc); | ||
112 : | |||
113 : | if (len < 0 || disc != sqrtd * sqrtd) | ||
114 : | error(_("invalid 'len' = %d in packed_ncol")); | ||
115 : | return (sqrtd - 1)/2; | ||
116 : | } | ||
117 : | |||
118 : | maechler | 890 | /** |
119 : | bates | 536 | * Allocate an SEXP of given type and length, assign it as slot nm in |
120 : | * the object, and return the SEXP. The validity of this function | ||
121 : | * depends on SET_SLOT not duplicating val when NAMED(val) == 0. If | ||
122 : | * this behavior changes then ALLOC_SLOT must use SET_SLOT followed by | ||
123 : | * GET_SLOT to ensure that the value returned is indeed the SEXP in | ||
124 : | * the slot. | ||
125 : | maechler | 890 | * |
126 : | bates | 536 | * @param obj object in which to assign the slot |
127 : | * @param nm name of the slot, as an R name object | ||
128 : | * @param type type of SEXP to allocate | ||
129 : | * @param length length of SEXP to allocate | ||
130 : | maechler | 890 | * |
131 : | bates | 536 | * @return SEXP of given type and length assigned as slot nm in obj |
132 : | */ | ||
133 : | static R_INLINE | ||
134 : | SEXP ALLOC_SLOT(SEXP obj, SEXP nm, SEXPTYPE type, int length) | ||
135 : | { | ||
136 : | SEXP val = allocVector(type, length); | ||
137 : | bates | 441 | |
138 : | bates | 536 | SET_SLOT(obj, nm, val); |
139 : | return val; | ||
140 : | } | ||
141 : | |||
142 : | maechler | 890 | /** |
143 : | bates | 679 | * Expand compressed pointers in the array mp into a full set of indices |
144 : | bates | 536 | * in the array mj. |
145 : | maechler | 890 | * |
146 : | bates | 679 | * @param ncol number of columns (or rows) |
147 : | bates | 536 | * @param mp column pointer vector of length ncol + 1 |
148 : | bates | 679 | * @param mj vector of length mp[ncol] to hold the result |
149 : | maechler | 890 | * |
150 : | bates | 536 | * @return mj |
151 : | */ | ||
152 : | static R_INLINE | ||
153 : | bates | 679 | int* expand_cmprPt(int ncol, const int mp[], int mj[]) |
154 : | bates | 536 | { |
155 : | int j; | ||
156 : | for (j = 0; j < ncol; j++) { | ||
157 : | int j2 = mp[j+1], jj; | ||
158 : | for (jj = mp[j]; jj < j2; jj++) mj[jj] = j; | ||
159 : | } | ||
160 : | return mj; | ||
161 : | } | ||
162 : | |||
163 : | |||
164 : | maechler | 890 | /** |
165 : | bates | 536 | * Return the linear index of the (row, col) entry in a csc structure. |
166 : | * If the entry is not found and missing is 0 an error is signaled; | ||
167 : | * otherwise the missing value is returned. | ||
168 : | maechler | 890 | * |
169 : | bates | 536 | * @param p vector of column pointers |
170 : | * @param i vector of row indices | ||
171 : | * @param row row index | ||
172 : | * @param col column index | ||
173 : | * @param missing the value to return is the row, col entry does not | ||
174 : | * exist. If this is zero and the row, col entry does not exist an | ||
175 : | * error is signaled. | ||
176 : | maechler | 890 | * |
177 : | bates | 536 | * @return index of element at (row, col) if it exists, otherwise missing |
178 : | */ | ||
179 : | static R_INLINE int | ||
180 : | check_csc_index(const int p[], const int i[], int row, int col, int missing) | ||
181 : | { | ||
182 : | int k, k2 = p[col + 1]; | ||
183 : | /* linear search - perhaps replace by bsearch */ | ||
184 : | for (k = p[col]; k < k2; k++) if (i[k] == row) return k; | ||
185 : | if (!missing) | ||
186 : | error("row %d and column %d not defined in rowind and colptr", | ||
187 : | row, col); | ||
188 : | return missing; | ||
189 : | } | ||
190 : | |||
191 : | SEXP alloc3Darray(SEXPTYPE mode, int nrow, int ncol, int nface); | ||
192 : | |||
193 : | /** | ||
194 : | * Calculate the zero-based index in a row-wise packed lower triangular matrix. | ||
195 : | * This is used for the arrays of blocked sparse matrices. | ||
196 : | * | ||
197 : | * @param i column number (zero-based) | ||
198 : | * @param k row number (zero-based) | ||
199 : | * | ||
200 : | * @return The index of the (k,i) element of a packed lower triangular matrix | ||
201 : | */ | ||
202 : | static R_INLINE | ||
203 : | int Lind(int k, int i) | ||
204 : | { | ||
205 : | if (k < i) error("Lind(k = %d, i = %d) must have k >= i", k, i); | ||
206 : | return (k * (k + 1))/2 + i; | ||
207 : | } | ||
208 : | |||
209 : | maechler | 890 | /** |
210 : | bates | 536 | * Check for a complete match on matrix dimensions |
211 : | maechler | 890 | * |
212 : | bates | 536 | * @param xd dimensions of first matrix |
213 : | * @param yd dimensions of second matrix | ||
214 : | maechler | 890 | * |
215 : | bates | 536 | * @return 1 if dimensions match, otherwise 0 |
216 : | */ | ||
217 : | static R_INLINE | ||
218 : | int match_mat_dims(const int xd[], const int yd[]) | ||
219 : | { | ||
220 : | return xd[0] == yd[0] && xd[1] == yd[1]; | ||
221 : | } | ||
222 : | |||
223 : | double *expand_csc_column(double *dest, int m, int j, | ||
224 : | const int Ap[], const int Ai[], const double Ax[]); | ||
225 : | |||
226 : | maechler | 890 | /** |
227 : | bates | 565 | * Apply a permutation to an integer vector |
228 : | maechler | 890 | * |
229 : | bates | 565 | * @param i vector of 0-based indices |
230 : | * @param n length of vector i | ||
231 : | * @param perm 0-based permutation vector of length max(i) + 1 | ||
232 : | */ | ||
233 : | static R_INLINE void | ||
234 : | int_permute(int i[], int n, const int perm[]) | ||
235 : | { | ||
236 : | int j; | ||
237 : | for (j = 0; j < n; j++) i[j] = perm[i[j]]; | ||
238 : | } | ||
239 : | |||
240 : | maechler | 890 | /** |
241 : | bates | 565 | * Force index pairs to be in the upper triangle of a matrix |
242 : | maechler | 890 | * |
243 : | bates | 565 | * @param i vector of 0-based row indices |
244 : | * @param j vector of 0-based column indices | ||
245 : | * @param nnz length of index vectors | ||
246 : | */ | ||
247 : | static R_INLINE void | ||
248 : | make_upper_triangular(int i[], int j[], int nnz) | ||
249 : | { | ||
250 : | int k; | ||
251 : | for (k = 0; k < nnz; k++) { | ||
252 : | if (i[k] > j[k]) { | ||
253 : | int tmp = i[k]; | ||
254 : | i[k] = j[k]; | ||
255 : | j[k] = tmp; | ||
256 : | } | ||
257 : | } | ||
258 : | } | ||
259 : | bates | 572 | |
260 : | maechler | 1200 | void make_d_matrix_triangular(double *x, SEXP from); |
261 : | void make_i_matrix_triangular( int *x, SEXP from); | ||
262 : | bates | 582 | |
263 : | maechler | 1200 | void make_d_matrix_symmetric(double *to, SEXP from); |
264 : | void make_i_matrix_symmetric( int *to, SEXP from); | ||
265 : | |||
266 : | bates | 738 | SEXP Matrix_expand_pointers(SEXP pP); |
267 : | |||
268 : | bates | 796 | |
269 : | maechler | 890 | /** |
270 : | bates | 796 | * Elementwise increment dest by src |
271 : | maechler | 890 | * |
272 : | bates | 796 | * @param dest vector to be incremented |
273 : | * @param src vector to be added to dest | ||
274 : | * @param n length of vectors | ||
275 : | maechler | 890 | * |
276 : | bates | 796 | * @return dest |
277 : | */ | ||
278 : | static R_INLINE double* | ||
279 : | vecIncrement(double dest[], const double src[], int n) { | ||
280 : | int i; | ||
281 : | for (i = 0; i < n; i++) dest[i] += src[i]; | ||
282 : | return dest; | ||
283 : | } | ||
284 : | |||
285 : | maechler | 890 | /** |
286 : | bates | 796 | * Elementwise sum of src1 and src2 into dest |
287 : | maechler | 890 | * |
288 : | bates | 796 | * @param dest vector to be incremented |
289 : | * @param src1 vector to be added | ||
290 : | * @param src1 second vector to be added | ||
291 : | * @param n length of vectors | ||
292 : | maechler | 890 | * |
293 : | bates | 796 | * @return dest |
294 : | */ | ||
295 : | static R_INLINE double* | ||
296 : | vecSum(double dest[], const double src1[], const double src2[], | ||
297 : | int n) { | ||
298 : | int i; | ||
299 : | for (i = 0; i < n; i++) dest[i] = src1[i] + src2[i]; | ||
300 : | return dest; | ||
301 : | } | ||
302 : | |||
303 : | bates | 963 | SEXP alloc_real_classed_matrix(char *class, int nrow, int ncol); |
304 : | bates | 1162 | SEXP alloc_dgeMatrix(int m, int n, SEXP rownms, SEXP colnms); |
305 : | SEXP alloc_dpoMatrix(int n, char *uplo, SEXP rownms, SEXP colnms); | ||
306 : | SEXP alloc_dtrMatrix(int n, char *uplo, char *diag, SEXP rownms, SEXP colnms); | ||
307 : | SEXP alloc_dsCMatrix(int n, int nz, char *uplo, SEXP rownms, SEXP colnms); | ||
308 : | bates | 963 | |
309 : | bates | 582 | #ifdef __cplusplus |
310 : | } | ||
311 : | bates | 572 | #endif |
312 : | bates | 582 | |
313 : | #endif /* MATRIX_MUTILS_H_ */ |
root@r-forge.r-project.org | ViewVC Help |
Powered by ViewVC 1.0.0 |