SCM

SCM Repository

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

View of /pkg/src/pdMat.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (download) (as text) (annotate)
Mon Mar 22 20:20:05 2004 UTC (15 years, 3 months ago) by bates
File size: 1513 byte(s)
Initial import
#include "Mutils.h"
#include <R_ext/Lapack.h>

SEXP pdCompSymm_pdFactor(SEXP pd)
{
    int i, j, ncol = asInteger(GET_SLOT(pd, install("ncol")));
    double *param = REAL(GET_SLOT(pd, install("param"))), sd,
        lcorr, corr, corr2, *L, offdiag;
    SEXP retval;
    
    retval = PROTECT(allocMatrix(REALSXP, ncol, ncol));
    L = REAL(retval);
    sd = exp(param[0]);
    lcorr = exp(param[1]);	/* logistic of the correlation */
    corr = (lcorr - 1.)/((double) ncol - 1.)/(lcorr + 1.);
    corr2 = sd * sqrt(1. - corr);
    corr = sd * sqrt((1. + (ncol - 1.) * corr) / ((double) ncol));
    
    for (i = 0; i < ncol; i++) {
	L[i * ncol] = corr;
    }
    
    for (i = 1; i < ncol; i++) {
	offdiag = -corr2/sqrt((double) (i * (i + 1)));
	for (j = 0; j < i; j++) {
	    L[i + (j * ncol)] = offdiag;
	}
	L[i * (ncol + 1)] = -sd * i;
    }
    UNPROTECT(1);
    return retval;
}

SEXP nlme_Chol(SEXP A)
{
    SEXP ans = PROTECT((TYPEOF(A) == REALSXP)?duplicate(A):
		       coerceVector(A, REALSXP));
    SEXP adims = getAttrib(A, R_DimSymbol);
    int m = INTEGER(adims)[0];
    int n = INTEGER(adims)[1];
    int i, j;

    if (!(isMatrix(A)))
	error("A must be a numeric (real) matrix");
    if (m != n) error("A must be a square matrix");
    for (j = 0; j < n; j++) {	/* zero the lower triangle */
	for (i = j+1; i < n; i++) {
	    REAL(ans)[i + j * n] = 0.;
	}
    }

    F77_CALL(dpotrf)("Upper", &m, REAL(ans), &m, &i);
    nlme_check_Lapack_error(i, "dpotrf");
    UNPROTECT(1);
    return ans;
}



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