The R Project SVN – R-packages – Diff – /branches/Matrix-APIchange/src/dppMatrix.c


#include "dppMatrix.h"
 
SEXP dppMatrix_validate(SEXP obj)
{
/*     int i, n = INTEGER(GET_SLOT(obj, Matrix_DimSym))[0]; */
/*     double *x = REAL(GET_SLOT(obj, Matrix_xSym)); */
 
    /* quick but nondefinitive check on positive definiteness */
/*     for (i = 0; i < n; i++) */
/* 	if (x[i * np1] < 0) */
/* 	    return mkString(_("dppMatrix is not positive definite")); */
    return dspMatrix_validate(obj);
}
 
SEXP dppMatrix_chol(SEXP x)
{
    SEXP val = get_factors(x, "pCholesky"),
	dimP = GET_SLOT(x, Matrix_DimSym),
	uploP = GET_SLOT(x, Matrix_uploSym);
    const char *uplo = CHAR(STRING_ELT(uploP, 0));
    int *dims = INTEGER(dimP), info;
 
    if (val != R_NilValue) return val;
    dims = INTEGER(dimP);
    val = PROTECT(NEW_OBJECT(MAKE_CLASS("pCholesky")));
    SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
    SET_SLOT(val, Matrix_diagSym, mkString("N"));
    SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
    SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
    F77_CALL(dpptrf)(uplo, dims, REAL(GET_SLOT(val, Matrix_xSym)), &info);
    if (info) {
	if(info > 0) /* e.g. x singular */
	    error(_("the leading minor of order %d is not positive definite"),
		    info);
	else /* should never happen! */
	    error(_("Lapack routine %s returned error code %d"), "dpptrf", info);
    }
    UNPROTECT(1);
    return set_factors(x, val, "pCholesky");
}
 
SEXP dppMatrix_rcond(SEXP obj, SEXP type)
{
    SEXP Chol = dppMatrix_chol(obj);
    char typnm[] = {'O', '\0'};	/* always use the one norm */
    int *dims = INTEGER(GET_SLOT(Chol, Matrix_DimSym)), info;
    double anorm = get_norm_sp(obj, typnm), rcond;
 
    F77_CALL(dppcon)(uplo_P(Chol), dims,
		     REAL(GET_SLOT(Chol, Matrix_xSym)), &anorm, &rcond,
		     (double *) R_alloc(3*dims[0], sizeof(double)),
		     (int *) R_alloc(dims[0], sizeof(int)), &info);
    return ScalarReal(rcond);
}
 
SEXP dppMatrix_solve(SEXP x)
{
    SEXP Chol = dppMatrix_chol(x);
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dppMatrix")));
    int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), info;
 
    SET_SLOT(val, Matrix_uploSym, duplicate(GET_SLOT(Chol, Matrix_uploSym)));
    SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(Chol, Matrix_xSym)));
    SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(Chol, Matrix_DimSym)));
    F77_CALL(dpptri)(uplo_P(val), dims,
		     REAL(GET_SLOT(val, Matrix_xSym)), &info);
    UNPROTECT(1);
    return val;
}
 
SEXP dppMatrix_matrix_solve(SEXP a, SEXP b)
{
    SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
    SEXP Chol = dppMatrix_chol(a);
    int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
	*bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
    int n = bdims[0], nrhs = bdims[1], info;
 
    if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
	error(_("Dimensions of system to be solved are inconsistent"));
    F77_CALL(dpptrs)(uplo_P(Chol), &n, &nrhs,
		     REAL(GET_SLOT(Chol, Matrix_xSym)),
		     REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);
    UNPROTECT(1);
    return val;
}

The R Project SVN R-packages

(root)/branches/Matrix-APIchange/src/dppMatrix.c – Rev 4418 → 4529

Rev 4418		Rev 4529
1	`#include "dppMatrix.h"`	1	`#include "dppMatrix.h"`
2		2
3	`SEXP dppMatrix_validate(SEXP obj)`	3	`SEXP dppMatrix_validate(SEXP obj)`
4	`{`	4	`{`
5	`/* int i, n = INTEGER(GET_SLOT(obj, Matrix_DimSym))[0]; */`	5	`/* int i, n = INTEGER(GET_SLOT(obj, Matrix_DimSym))[0]; */`
6	`/* double x = REAL(GET_SLOT(obj, Matrix_xSym)); /`	6	`/* double x = REAL(GET_SLOT(obj, Matrix_xSym)); /`
7		7
8	`/* quick but nondefinitive check on positive definiteness */`	8	`/* quick but nondefinitive check on positive definiteness */`
9	`/* for (i = 0; i < n; i++) */`	9	`/* for (i = 0; i < n; i++) */`
10	`/* if (x[i * np1] < 0) */`	10	`/* if (x[i * np1] < 0) */`
11	`/* return mkString(_("dppMatrix is not positive definite")); */`	11	`/* return mkString(_("dppMatrix is not positive definite")); */`
12	`return dspMatrix_validate(obj);`	12	`return dspMatrix_validate(obj);`
13	`}`	13	`}`
14		14
15	`SEXP dppMatrix_chol(SEXP x)`	15	`SEXP dppMatrix_chol(SEXP x)`
16	`{`	16	`{`
17	`SEXP val = get_factors(x, "pCholesky"),`	17	`SEXP val = get_factors(x, "pCholesky"),`
18	`dimP = GET_SLOT(x, Matrix_DimSym),`	18	`dimP = GET_SLOT(x, Matrix_DimSym),`
19	`uploP = GET_SLOT(x, Matrix_uploSym);`	19	`uploP = GET_SLOT(x, Matrix_uploSym);`
20	`const char *uplo = CHAR(STRING_ELT(uploP, 0));`	20	`const char *uplo = CHAR(STRING_ELT(uploP, 0));`
21	`int *dims = INTEGER(dimP), info;`	21	`int *dims = INTEGER(dimP), info;`
22		22
23	`if (val != R_NilValue) return val;`	23	`if (val != R_NilValue) return val;`
24	`dims = INTEGER(dimP);`	24	`dims = INTEGER(dimP);`
25	`val = PROTECT(NEW_OBJECT(MAKE_CLASS("pCholesky")));`	25	`val = PROTECT(NEW_OBJECT(MAKE_CLASS("pCholesky")));`
26	`SET_SLOT(val, Matrix_uploSym, duplicate(uploP));`	26	`SET_SLOT(val, Matrix_uploSym, duplicate(uploP));`
27	`SET_SLOT(val, Matrix_diagSym, mkString("N"));`	27	`SET_SLOT(val, Matrix_diagSym, mkString("N"));`
28	`SET_SLOT(val, Matrix_DimSym, duplicate(dimP));`	28	`SET_SLOT(val, Matrix_DimSym, duplicate(dimP));`
29	`SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));`	29	`SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));`
30	`F77_CALL(dpptrf)(uplo, dims, REAL(GET_SLOT(val, Matrix_xSym)), &info);`	30	`F77_CALL(dpptrf)(uplo, dims, REAL(GET_SLOT(val, Matrix_xSym)), &info);`
31	`if (info) {`	31	`if (info) {`
32	`if(info > 0) /* e.g. x singular */`	32	`if(info > 0) /* e.g. x singular */`
33	`error(_("the leading minor of order %d is not positive definite"),`	33	`error(_("the leading minor of order %d is not positive definite"),`
34	`info);`	34	`info);`
35	`else /* should never happen! */`	35	`else /* should never happen! */`
36	`error(_("Lapack routine %s returned error code %d"), "dpptrf", info);`	36	`error(_("Lapack routine %s returned error code %d"), "dpptrf", info);`
37	`}`	37	`}`
38	`UNPROTECT(1);`	38	`UNPROTECT(1);`
39	`return set_factors(x, val, "pCholesky");`	39	`return set_factors(x, val, "pCholesky");`
40	`}`	40	`}`
41		41
42	`SEXP dppMatrix_rcond(SEXP obj, SEXP type)`	42	`SEXP dppMatrix_rcond(SEXP obj, SEXP type)`
43	`{`	43	`{`
44	`SEXP Chol = dppMatrix_chol(obj);`	44	`SEXP Chol = dppMatrix_chol(obj);`
45	`char typnm[] = {'O', '\0'}; /* always use the one norm */`	45	`char typnm[] = {'O', '\0'}; /* always use the one norm */`
46	`int *dims = INTEGER(GET_SLOT(Chol, Matrix_DimSym)), info;`	46	`int *dims = INTEGER(GET_SLOT(Chol, Matrix_DimSym)), info;`
47	`double anorm = get_norm_sp(obj, typnm), rcond;`	47	`double anorm = get_norm_sp(obj, typnm), rcond;`
48		48
49	`F77_CALL(dppcon)(uplo_P(Chol), dims,`	49	`F77_CALL(dppcon)(uplo_P(Chol), dims,`
50	`REAL(GET_SLOT(Chol, Matrix_xSym)), &anorm, &rcond,`	50	`REAL(GET_SLOT(Chol, Matrix_xSym)), &anorm, &rcond,`
51	`(double ) R_alloc(3dims[0], sizeof(double)),`	51	`(double ) R_alloc(3dims[0], sizeof(double)),`
52	`(int *) R_alloc(dims[0], sizeof(int)), &info);`	52	`(int *) R_alloc(dims[0], sizeof(int)), &info);`
53	`return ScalarReal(rcond);`	53	`return ScalarReal(rcond);`
54	`}`	54	`}`
55		55
56	`SEXP dppMatrix_solve(SEXP x)`	56	`SEXP dppMatrix_solve(SEXP x)`
57	`{`	57	`{`
58	`SEXP Chol = dppMatrix_chol(x);`	58	`SEXP Chol = dppMatrix_chol(x);`
59	`SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dppMatrix")));`	59	`SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dppMatrix")));`
60	`int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), info;`	60	`int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), info;`
61		61
62	`SET_SLOT(val, Matrix_uploSym, duplicate(GET_SLOT(Chol, Matrix_uploSym)));`	62	`SET_SLOT(val, Matrix_uploSym, duplicate(GET_SLOT(Chol, Matrix_uploSym)));`
63	`SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(Chol, Matrix_xSym)));`	63	`SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(Chol, Matrix_xSym)));`
64	`SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(Chol, Matrix_DimSym)));`	64	`SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(Chol, Matrix_DimSym)));`
65	`F77_CALL(dpptri)(uplo_P(val), dims,`	65	`F77_CALL(dpptri)(uplo_P(val), dims,`
66	`REAL(GET_SLOT(val, Matrix_xSym)), &info);`	66	`REAL(GET_SLOT(val, Matrix_xSym)), &info);`
67	`UNPROTECT(1);`	67	`UNPROTECT(1);`
68	`return val;`	68	`return val;`
69	`}`	69	`}`
70		70
71	`SEXP dppMatrix_matrix_solve(SEXP a, SEXP b)`	71	`SEXP dppMatrix_matrix_solve(SEXP a, SEXP b)`
72	`{`	72	`{`
73	`SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));`	73	`SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));`
74	`SEXP Chol = dppMatrix_chol(a);`	74	`SEXP Chol = dppMatrix_chol(a);`
75	`int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),`	75	`int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),`
76	`*bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));`	76	`*bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));`
77	`int n = bdims[0], nrhs = bdims[1], info;`	77	`int n = bdims[0], nrhs = bdims[1], info;`
78		78
79	`if (adims != bdims \|\| bdims[1] < 1 \|\| *adims < 1)`	79	`if (adims != bdims \|\| bdims[1] < 1 \|\| *adims < 1)`
80	`error(_("Dimensions of system to be solved are inconsistent"));`	80	`error(_("Dimensions of system to be solved are inconsistent"));`
81	`F77_CALL(dpptrs)(uplo_P(Chol), &n, &nrhs,`	81	`F77_CALL(dpptrs)(uplo_P(Chol), &n, &nrhs,`
82	`REAL(GET_SLOT(Chol, Matrix_xSym)),`	82	`REAL(GET_SLOT(Chol, Matrix_xSym)),`
83	`REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);`	83	`REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);`
84	`UNPROTECT(1);`	84	`UNPROTECT(1);`
85	`return val;`	85	`return val;`
86	`}`	86	`}`