The R Project SVN R-packages

Rev

Rev 1282 | Rev 2254 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 1282 Rev 1283
Line 34... Line 34...
34 
{
 34 
{
35 
    char *nms[] = {"L", "U", "P", ""};
 35 
    char *nms[] = {"L", "U", "P", ""};
36 
    SEXP L, U, P, val = PROTECT(Matrix_make_named(VECSXP, nms)),
 36 
    SEXP L, U, P, val = PROTECT(Matrix_make_named(VECSXP, nms)),
37 
	lux = GET_SLOT(x, Matrix_xSym),
 37 
	lux = GET_SLOT(x, Matrix_xSym),
38 
	dd = GET_SLOT(x, Matrix_DimSym);
 38 
	dd = GET_SLOT(x, Matrix_DimSym);
39 
    int *perm, *pivot = INTEGER(GET_SLOT(x, Matrix_permSym)),
 39 
    int *iperm, *perm, *pivot = INTEGER(GET_SLOT(x, Matrix_permSym)),
40 
	i, n = INTEGER(dd)[0];
 40 
	i, n = INTEGER(dd)[0];
41 
 
 41 
 
42 
    SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("dtrMatrix")));
 42 
    SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("dtrMatrix")));
43 
    L = VECTOR_ELT(val, 0);
 43 
    L = VECTOR_ELT(val, 0);
44 
    SET_VECTOR_ELT(val, 1, NEW_OBJECT(MAKE_CLASS("dtrMatrix")));
 44 
    SET_VECTOR_ELT(val, 1, NEW_OBJECT(MAKE_CLASS("dtrMatrix")));
Line 54... Line 54...
54 
    SET_SLOT(U, Matrix_DimSym, duplicate(dd));
 54 
    SET_SLOT(U, Matrix_DimSym, duplicate(dd));
55 
    SET_SLOT(U, Matrix_uploSym, mkString("U"));
 55 
    SET_SLOT(U, Matrix_uploSym, mkString("U"));
56 
    SET_SLOT(U, Matrix_diagSym, mkString("N"));
 56 
    SET_SLOT(U, Matrix_diagSym, mkString("N"));
57 
    make_array_triangular(REAL(GET_SLOT(U, Matrix_xSym)), U);
 57 
    make_array_triangular(REAL(GET_SLOT(U, Matrix_xSym)), U);
58 
    SET_SLOT(P, Matrix_DimSym, duplicate(dd));
 58 
    SET_SLOT(P, Matrix_DimSym, duplicate(dd));
- 
 
 59 
    iperm = Calloc(n, int);
59 
    perm = INTEGER(ALLOC_SLOT(P, Matrix_permSym, INTSXP, n));
 60 
    perm = INTEGER(ALLOC_SLOT(P, Matrix_permSym, INTSXP, n));
- 
 
 61
60 
    for (i = 0; i < n; i++) perm[i] = i + 1;
 62 
    for (i = 0; i < n; i++) iperm[i] = i + 1; /* initialize permutation*/
61 
    for (i = 0; i < n; i++) {
 63 
    for (i = 0; i < n; i++) {	/* generate inverse permutation */
62 
	int newpos = pivot[i] - 1;
 64 
	int newpos = pivot[i] - 1;
63 
	if (newpos != i) {
 65 
	if (newpos != i) {
64 
	    int tmp = perm[i];
 66 
	    int tmp = iperm[i];
65 
 
 67 
 
66 
	    perm[i] = newpos + 1;
 68 
	    iperm[i] = iperm[newpos];
67 
	    perm[newpos] = tmp;
 69 
	    iperm[newpos] = tmp;
68 
	}
 70 
	}
69 
    }
 71 
    }
- 
 
 72 
				/* invert the inverse */
- 
 
 73 
    for (i = 0; i < n; i++) perm[iperm[i] - 1] = i + 1;
70 
    UNPROTECT(1);
 74 
    UNPROTECT(1);
71 
    return val;
 75 
    return val;
72 
}
 76 
}