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