SCM

SCM Repository

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

Annotation of /pkg/src/pdMat.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 10 #include "Mutils.h"
2 :     #include <R_ext/Lapack.h>
3 :    
4 :     SEXP pdCompSymm_pdFactor(SEXP pd)
5 :     {
6 :     int i, j, ncol = asInteger(GET_SLOT(pd, install("ncol")));
7 :     double *param = REAL(GET_SLOT(pd, install("param"))), sd,
8 :     lcorr, corr, corr2, *L, offdiag;
9 :     SEXP retval;
10 :    
11 :     retval = PROTECT(allocMatrix(REALSXP, ncol, ncol));
12 :     L = REAL(retval);
13 :     sd = exp(param[0]);
14 :     lcorr = exp(param[1]); /* logistic of the correlation */
15 :     corr = (lcorr - 1.)/((double) ncol - 1.)/(lcorr + 1.);
16 :     corr2 = sd * sqrt(1. - corr);
17 :     corr = sd * sqrt((1. + (ncol - 1.) * corr) / ((double) ncol));
18 :    
19 :     for (i = 0; i < ncol; i++) {
20 :     L[i * ncol] = corr;
21 :     }
22 :    
23 :     for (i = 1; i < ncol; i++) {
24 :     offdiag = -corr2/sqrt((double) (i * (i + 1)));
25 :     for (j = 0; j < i; j++) {
26 :     L[i + (j * ncol)] = offdiag;
27 :     }
28 :     L[i * (ncol + 1)] = -sd * i;
29 :     }
30 :     UNPROTECT(1);
31 :     return retval;
32 :     }
33 :    
34 :     SEXP nlme_Chol(SEXP A)
35 :     {
36 :     SEXP ans = PROTECT((TYPEOF(A) == REALSXP)?duplicate(A):
37 :     coerceVector(A, REALSXP));
38 :     SEXP adims = getAttrib(A, R_DimSymbol);
39 :     int m = INTEGER(adims)[0];
40 :     int n = INTEGER(adims)[1];
41 :     int i, j;
42 :    
43 :     if (!(isMatrix(A)))
44 :     error("A must be a numeric (real) matrix");
45 :     if (m != n) error("A must be a square matrix");
46 :     for (j = 0; j < n; j++) { /* zero the lower triangle */
47 :     for (i = j+1; i < n; i++) {
48 :     REAL(ans)[i + j * n] = 0.;
49 :     }
50 :     }
51 :    
52 :     F77_CALL(dpotrf)("Upper", &m, REAL(ans), &m, &i);
53 :     nlme_check_Lapack_error(i, "dpotrf");
54 :     UNPROTECT(1);
55 :     return ans;
56 :     }
57 :    
58 :    

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