SCM

SCM Repository

[matrix] Annotation of /pkg/src/Mutils.c
ViewVC logotype

Annotation of /pkg/src/Mutils.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (view) (download) (as text)

1 : bates 10 #include "Mutils.h"
2 : bates 47 #include "triplet_to_col.h"
3 : bates 10 #include <R_ext/Lapack.h>
4 :    
5 :     SEXP
6 :     Matrix_DSym,
7 :     Matrix_DIsqrtSym,
8 :     Matrix_DimSym,
9 :     Matrix_GpSym,
10 :     Matrix_LIiSym,
11 :     Matrix_LIpSym,
12 :     Matrix_LIxSym,
13 :     Matrix_LiSym,
14 :     Matrix_LpSym,
15 :     Matrix_LxSym,
16 :     Matrix_OmegaSym,
17 :     Matrix_ParentSym,
18 :     Matrix_RXXSym,
19 :     Matrix_RZXSym,
20 :     Matrix_XtXSym,
21 :     Matrix_ZtXSym,
22 :     Matrix_bVarSym,
23 :     Matrix_devianceSym,
24 :     Matrix_devCompSym,
25 :     Matrix_diagSym,
26 :     Matrix_iSym,
27 :     Matrix_ipermSym,
28 :     Matrix_jSym,
29 :     Matrix_matSym,
30 :     Matrix_ncSym,
31 :     Matrix_pSym,
32 :     Matrix_permSym,
33 :     Matrix_statusSym,
34 :     Matrix_uploSym,
35 :     Matrix_xSym,
36 :     Matrix_zSym;
37 :    
38 :     SEXP Matrix_init(void)
39 :     {
40 :     Matrix_DSym = install("D");
41 :     Matrix_DIsqrtSym = install("DIsqrt");
42 :     Matrix_DimSym = install("Dim");
43 :     Matrix_GpSym = install("Gp");
44 :     Matrix_LIiSym = install("LIi");
45 :     Matrix_LIpSym = install("LIp");
46 :     Matrix_LIxSym = install("LIx");
47 :     Matrix_LiSym = install("Li");
48 :     Matrix_LpSym = install("Lp");
49 :     Matrix_LxSym = install("Lx");
50 :     Matrix_OmegaSym = install("Omega");
51 :     Matrix_ParentSym = install("Parent");
52 :     Matrix_RXXSym = install("RXX");
53 :     Matrix_RZXSym = install("RZX");
54 :     Matrix_XtXSym = install("XtX");
55 :     Matrix_ZtXSym = install("ZtX");
56 :     Matrix_bVarSym = install("bVar");
57 :     Matrix_devianceSym = install("deviance");
58 :     Matrix_devCompSym = install("devComp");
59 :     Matrix_diagSym = install("diag");
60 :     Matrix_iSym = install("i");
61 :     Matrix_ipermSym = install("iperm");
62 :     Matrix_jSym = install("j");
63 :     Matrix_matSym = install("mat");
64 :     Matrix_ncSym = install("nc");
65 :     Matrix_pSym = install("p");
66 :     Matrix_permSym = install("perm");
67 :     Matrix_statusSym = install("status");
68 :     Matrix_uploSym = install("uplo");
69 :     Matrix_xSym = install("x");
70 :     Matrix_zSym = install("z");
71 :     return R_NilValue;
72 :     }
73 :    
74 :     char norm_type(char *typstr)
75 :     {
76 :     char typup;
77 :    
78 :     if (strlen(typstr) != 1)
79 :     error("argument type[1]='%s' must be a character string of string length 1",
80 :     typstr);
81 :     typup = toupper(*typstr);
82 :     if (typup == '1') typup = 'O'; /* aliases */
83 :     if (typup == 'E') typup = 'F';
84 :     if (typup != 'M' && typup != 'O' && typup != 'I' && typup != 'F')
85 :     error("argument type[1]='%s' must be one of 'M','1','O','I','F' or 'E'",
86 :     typstr);
87 :     return typup;
88 :     }
89 :    
90 :     char rcond_type(char *typstr)
91 :     {
92 :     char typup;
93 :    
94 :     if (strlen(typstr) != 1)
95 :     error("argument type[1]='%s' must be a character string of string length 1",
96 :     typstr);
97 :     typup = toupper(*typstr);
98 :     if (typup == '1') typup = 'O'; /* alias */
99 :     if (typup != 'O' && typup != 'I')
100 :     error("argument type[1]='%s' must be one of '1','O', or 'I'",
101 :     typstr);
102 :     return typup;
103 :     }
104 :    
105 :     double get_double_by_name(SEXP obj, char *nm)
106 :     {
107 :     SEXP nms = getAttrib(obj, R_NamesSymbol);
108 :     int i, len = length(obj);
109 :    
110 :     if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue))
111 :     error("object must be a named, numeric vector");
112 :     for (i = 0; i < len; i++) {
113 :     if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
114 :     return REAL(obj)[i];
115 :     }
116 :     }
117 :     return R_NaReal;
118 :     }
119 :    
120 :     SEXP
121 :     set_double_by_name(SEXP obj, double val, char *nm)
122 :     {
123 :     SEXP nms = getAttrib(obj, R_NamesSymbol);
124 :     int i, len = length(obj);
125 :    
126 :     if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue))
127 :     error("object must be a named, numeric vector");
128 :     for (i = 0; i < len; i++) {
129 :     if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
130 :     REAL(obj)[i] = val;
131 :     return obj;
132 :     }
133 :     }
134 :     {SEXP nx = PROTECT(allocVector(REALSXP, len + 1)),
135 :     nnms = allocVector(STRSXP, len + 1);
136 :    
137 :     setAttrib(nx, R_NamesSymbol, nnms);
138 :     for (i = 0; i < len; i++) {
139 :     REAL(nx)[i] = REAL(obj)[i];
140 :     SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i)));
141 :     }
142 :     REAL(nx)[len] = val;
143 :     SET_STRING_ELT(nnms, len, mkChar(nm));
144 :     UNPROTECT(1);
145 :     return nx;
146 :     }
147 :     }
148 :    
149 :     SEXP as_det_obj(double val, int log, int sign)
150 :     {
151 :     SEXP det = PROTECT(allocVector(VECSXP, 2)),
152 :     nms = allocVector(STRSXP, 2),
153 :     vv = ScalarReal(val);
154 :    
155 :     setAttrib(det, R_NamesSymbol, nms);
156 :     SET_STRING_ELT(nms, 0, mkChar("modulus"));
157 :     SET_STRING_ELT(nms, 1, mkChar("sign"));
158 :     setAttrib(vv, install("logarithm"), ScalarLogical(log));
159 :     SET_VECTOR_ELT(det, 0, vv);
160 :     SET_VECTOR_ELT(det, 1, ScalarInteger(sign));
161 :     setAttrib(det, R_ClassSymbol, ScalarString(mkChar("det")));
162 :     UNPROTECT(1);
163 :     return det;
164 :     }
165 :    
166 :     SEXP get_factorization(SEXP obj, char *nm)
167 :     {
168 :     SEXP fac = GET_SLOT(obj, install("factorization")),
169 :     nms = getAttrib(fac, R_NamesSymbol);
170 :     int i, len = length(fac);
171 :    
172 :     if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue))
173 :     error("factorization slot must be a named list");
174 :     for (i = 0; i < len; i++) {
175 :     if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
176 :     return VECTOR_ELT(fac, i);
177 :     }
178 :     }
179 :     return R_NilValue;
180 :     }
181 :    
182 :     SEXP set_factorization(SEXP obj, SEXP val, char *nm)
183 :     {
184 :     SEXP fac = GET_SLOT(obj, install("factorization")),
185 :     nms = getAttrib(fac, R_NamesSymbol), nfac, nnms;
186 :     int i, len = length(fac);
187 :    
188 :     if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue))
189 :     error("factorization slot must be a named list");
190 :     for (i = 0; i < len; i++) {
191 :     if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
192 :     SET_VECTOR_ELT(fac, i, val);
193 :     return val;
194 :     }
195 :     }
196 :     nfac = allocVector(VECSXP, len + 1);
197 :     nnms = allocVector(STRSXP, len + 1);
198 :     setAttrib(nfac, R_NamesSymbol, nnms);
199 :     for (i = 0; i < len; i++) {
200 :     SET_VECTOR_ELT(nfac, i, VECTOR_ELT(fac, i));
201 :     SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i)));
202 :     }
203 :     SET_VECTOR_ELT(nfac, len, val);
204 :     SET_STRING_ELT(nnms, len, mkChar(nm));
205 :     SET_SLOT(obj, install("factorization"), nfac);
206 :     return val;
207 :     }
208 :    
209 :     SEXP cscMatrix_set_Dim(SEXP x, int nrow)
210 :     {
211 :     int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
212 :    
213 :     dims[0] = nrow;
214 :     dims[1] = length(GET_SLOT(x, Matrix_pSym)) - 1;
215 :     return x;
216 :     }
217 :    
218 :     int csc_unsorted_columns(int ncol, const int p[], const int i[])
219 :     {
220 :     int j;
221 :     for (j = 0; j < ncol; j++) {
222 :     int ind, lst = p[j+1] - 1;
223 :     for (ind = p[j]; ind < lst; ind++) {
224 :     if (i[ind] > i[ind+1]) return 1;
225 :     }
226 :     }
227 :     return 0;
228 :     }
229 :    
230 :     void csc_sort_columns(int ncol, const int p[], int i[], double x[])
231 :     {
232 :     int j, maxdiff, *ord;
233 :     double *dd;
234 :    
235 :     maxdiff = -1;
236 :     for (j = 0; j < ncol; j++) {
237 :     int diff = p[j+1] - p[j];
238 :     if (diff > maxdiff) maxdiff = diff;
239 :     }
240 :     ord = Calloc(maxdiff, int);
241 :     dd = Calloc(maxdiff, double);
242 :     for (j = 0; j < ncol; j++) {
243 :     int cLen = p[j+1] - p[j];
244 :     if (cLen > 1) {
245 :     int k, offset = p[j];
246 :     for (k = 0; k < cLen; k++) ord[k] = k;
247 :     R_qsort_int_I(i + offset, ord, 1, cLen);
248 :     for (k = 0; k < cLen; k++) dd[k] = x[ord[k] + offset];
249 :     Memcpy(x + offset, dd, cLen);
250 :     }
251 :     }
252 :     Free(ord); Free(dd);
253 :     }
254 :    
255 :     SEXP csc_check_column_sorting(SEXP m)
256 :     {
257 :     int *mp = INTEGER(GET_SLOT(m, Matrix_pSym)),
258 :     *mi = INTEGER(GET_SLOT(m, Matrix_iSym)),
259 :     ncol = INTEGER(GET_SLOT(m, Matrix_DimSym))[1];
260 :    
261 :     if (csc_unsorted_columns(ncol, mp, mi))
262 :     csc_sort_columns(ncol, mp, mi, REAL(GET_SLOT(m, Matrix_xSym)));
263 :     return m;
264 :     }
265 :    
266 :     SEXP triple_as_SEXP(int nrow, int ncol, int nz,
267 :     const int Ti [], const int Tj [], const double Tx [],
268 :     char *Rclass)
269 :     {
270 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS(Rclass)));
271 :     int *Ai, *Ap;
272 :     double *Ax;
273 :    
274 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
275 :     Ap = INTEGER(GET_SLOT(val, Matrix_pSym));
276 :     Ai = Calloc(nz, int); Ax = Calloc(nz, double);
277 :     triplet_to_col(nrow, ncol, nz, Ti, Tj, Tx, Ap, Ai, Ax);
278 :     nz = Ap[ncol];
279 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nz));
280 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), Ai, nz); Free(Ai);
281 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nz));
282 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), Ax, nz); Free(Ax);
283 :     UNPROTECT(1);
284 :     return cscMatrix_set_Dim(val, nrow);
285 :     }
286 :    
287 :     void csc_components_transpose(int m, int n, int nnz,
288 :     const int xp[], const int xi[],
289 :     const double xx[],
290 :     int ap[], int ai[], double ax[])
291 :     {
292 :     int k, kk,
293 :     *ind = (int *) R_alloc(nnz, sizeof(int)),
294 :     *aj = (int *) R_alloc(nnz, sizeof(int));
295 :    
296 :     Memcpy(aj, xi, nnz); /* copy xi into aj and sort */
297 :     for (k = 0; k < nnz; k++) ind[k] = k;
298 :     R_qsort_int_I(aj, ind, 1, nnz);
299 :    
300 :     ap[0] = 0; kk = 0; /* generate ap from aj */
301 :     for (k = 1; k < m; k++) {
302 :     while (aj[kk] < k) kk++;
303 :     ap[k] = kk;
304 :     }
305 :     ap[m] = nnz;
306 :    
307 :     for (k = 0; k < n; k++) { /* overwrite aj with (implicit) xj */
308 :     for (kk = xp[k]; kk < xp[k+1]; kk++) aj[kk] = k;
309 :     }
310 :     for (k = 0; k < nnz; k++) { /* write ax and ai from xx and xj */
311 :     kk = ind[k];
312 :     ax[k] = xx[kk];
313 :     ai[k] = aj[kk];
314 :     }
315 :     if (csc_unsorted_columns(m, ap, ai)) csc_sort_columns(m, ap, ai, ax);
316 :     }
317 :    
318 :     void ssc_symbolic_permute(int n, int upper, const int perm[],
319 :     int Ap[], int Ai[])
320 :     {
321 :     int
322 :     j, k,
323 :     nnz = Ap[n],
324 :     *Aj = Calloc(nnz, int),
325 :     *ord = Calloc(nnz, int),
326 :     *ii = Calloc(nnz, int);
327 :    
328 :     for (j = 0; j < n; j++) {
329 :     int pj = perm[j];
330 :     for (k = Ap[j]; k < Ap[j+1]; k++) {
331 :     Aj[k] = pj;
332 :     }
333 :     }
334 :     for (k = 0; k < nnz; k++) {
335 :     Ai[k] = perm[Ai[k]];
336 :     ord[k] = k;
337 :     if ((upper && Ai[k] > Aj[k]) || (!upper && Ai[k] < Aj[k])) {
338 :     int tmp = Ai[k]; Ai[k] = Aj[k]; Aj[k] = tmp;
339 :     warning("in part of code that should not be used");
340 :     }
341 :     }
342 :     R_qsort_int_I(Aj, ord, 1, nnz); /* sort Aj carrying along ind */
343 :    
344 :     k = nnz - 1;
345 :     for (j = n - 1; j >= 0; j--) { /* generate new Ap */
346 :     for(; Aj[k] >= j && k >= 0; k--) Ap[j] = k;
347 :     }
348 :     for (k = 0; k < nnz; k++) ii[k] = Ai[ord[k]];
349 :     Memcpy(Ai, ii, nnz);
350 :     for (j = 0; j < n; j++) R_isort(Ai + Ap[j], Ap[j+1] - Ap[j]);
351 :     Free(Aj); Free(ord); Free(ii);
352 :     }
353 :    
354 :    
355 :     /**
356 :     * Symmetrize a matrix by copying the strict upper triangle into the
357 :     * lower triangle.
358 :     *
359 :     * @param a pointer to a matrix in Fortran storage mode
360 :     * @param nc number of columns (and rows and leading dimension) in the matrix
361 :     *
362 :     * @return a, symmetrized
363 :     */
364 :     double *
365 :     nlme_symmetrize(double *a, const int nc)
366 :     {
367 :     int i, j;
368 :    
369 :     for (i = 1; i < nc; i++)
370 :     for (j = 0; j < i; j++)
371 :     a[i + j*nc] = a[j + i*nc];
372 :     return a;
373 :     }
374 :    
375 :     /**
376 :     * Check the error code returned by an Lapack routine and create an
377 :     * appropriate error message.
378 :     *
379 :     * @param info Error code as returned from the Lapack routine
380 :     * @param laName Character string containing the name of the Lapack routine
381 :     */
382 :     void
383 :     nlme_check_Lapack_error(int info, const char *laName)
384 :     {
385 :     if (info != 0) {
386 :     if (info > 0)
387 :     error("error code %d from Lapack routine %s", info, laName);
388 :     error("argument no. %d to Lapack routine %s is illegal",
389 :     -info, laName);
390 :     }
391 :     }
392 :    
393 :     /**
394 :     * Calculate the inner product of vec(nlev*D^{-1} - A'A)/2 and the
395 :     * pdgradient array regarded as a nc*nc by plen matrix. This
396 :     * calculation is used in several of the LMEgradient methods.
397 :     *
398 :     * @param factor The nc by nc factor of the pdMat object
399 :     * @param A The nc by nc matrix A from the LME decomposition.
400 :     * @param nlev The number of groups associated with the random effect
401 :     * @param nc The number of columns in the matrix
402 :     * @param pdgradient A pdgradient object of dimension nc by nc by plen
403 :     * @param value An array of length plen in which the gradient will be returned
404 :     *
405 :     * @return value, with the LME gradient
406 :     */
407 :     double *
408 :     LMEgradient(const double* factor, const double* A, const int nlev,
409 :     const int nc, const double* pdgradient, const int plen,
410 :     double* value)
411 :     {
412 :     int info, ncsq = nc*nc, one_i = 1;
413 :     double nlev2_d = ((double) nlev)/2., mhalf_d = -0.5, one_d = 1.0,
414 :     zero_d = 0.0;
415 :     double *fact = Calloc(nc * nc, double);
416 :    
417 :     F77_CALL(dlacpy)("U", &nc, &nc, factor, &nc, fact, &nc);
418 :     F77_CALL(dpotri)("U", &nc, fact, &nc, &info);
419 :     nlme_check_Lapack_error(info, "dpotri");
420 :     F77_CALL(dsyrk)("U", "T", &nc, &nc, &mhalf_d, A, &nc, &nlev2_d,
421 :     fact, &nc);
422 :     nlme_symmetrize(fact, nc);
423 :     F77_CALL(dgemv)("T", &ncsq, &plen, &one_d, pdgradient, &ncsq,
424 :     fact, &one_i, &zero_d, value, &one_i);
425 :     Free(fact);
426 :     return value;
427 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge