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