The R Project SVN R-packages

Rev

Rev 4909 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 4909 Rev 4972
Line 57... Line 57...
57 
			      /* uploT = */ (*CHAR(asChar(uplo)) == 'U')? 1: -1,
 57 
			      /* uploT = */ (*CHAR(asChar(uplo)) == 'U')? 1: -1,
58 
			      Rkind,
 58 
			      Rkind,
59 
			      /* diag = */ CHAR(STRING_ELT(diag, 0)),
 59 
			      /* diag = */ CHAR(STRING_ELT(diag, 0)),
60 
			      GET_SLOT(x, Matrix_DimNamesSym));
 60 
			      GET_SLOT(x, Matrix_DimNamesSym));
61 
}
 61 
}
- 
 
 62 
 
- 
 
 63 
SEXP Tsparse_diagU2N(SEXP x)
- 
 
 64 
{
- 
 
 65 
    char *valid[] = {"dtTMatrix", /* 0 */
- 
 
 66 
		     "ltTMatrix", /* 1 */
- 
 
 67 
		     "ntTMatrix", /* 2 : no x slot */
- 
 
 68 
		     "ztTMatrix", /* 3 */
- 
 
 69 
		     ""};
- 
 
 70 
/* #define xSXP(iTyp) ((iTyp == 0) ? REALSXP : ((iTyp == 1) ? LGLSXP : /\* else *\/ CPLXSXP)); */
- 
 
 71 
/* #define xTYPE(iTyp) ((iTyp == 0) ? double : ((iTyp == 1) ? int : /\* else *\/ Rcomplex)); */
- 
 
 72 
    int ctype = Matrix_check_class(class_P(x), valid);
- 
 
 73 
 
- 
 
 74 
    if (ctype < 0 || *diag_P(x) != 'U') {
- 
 
 75 
	/* "trivially fast" when not triangular (<==> no 'diag' slot),
- 
 
 76 
	   or not *unit* triangular */
- 
 
 77 
	return (x);
- 
 
 78 
    }
- 
 
 79 
    else { /* instead of going to Csparse -> Cholmod -> Csparse -> Tsparse, work directly: */
- 
 
 80 
	int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
- 
 
 81 
	    nnz = length(GET_SLOT(x, Matrix_iSym)),
- 
 
 82 
	    new_n = nnz + n;
- 
 
 83 
	SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS(class_P(x))));
- 
 
 84 
	int *islot = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, new_n)),
- 
 
 85 
	    *jslot = INTEGER(ALLOC_SLOT(ans, Matrix_jSym, INTSXP, new_n));
- 
 
 86 
 
- 
 
 87 
	slot_dup(ans, x, Matrix_DimSym);
- 
 
 88 
	SET_DimNames(ans, x);
- 
 
 89 
	slot_dup(ans, x, Matrix_uploSym);
- 
 
 90 
	SET_SLOT(ans, Matrix_diagSym, mkString("N"));
- 
 
 91 
 
- 
 
 92 
	/* Build the new i- and j- slots : first copy the current : */
- 
 
 93 
	Memcpy(islot, INTEGER(GET_SLOT(x, Matrix_iSym)), nnz);
- 
 
 94 
	Memcpy(jslot, INTEGER(GET_SLOT(x, Matrix_jSym)), nnz);
- 
 
 95 
	/* then, add the new (i,j) slot entries: */
- 
 
 96 
	for(i = 0; i < n; i++) {
- 
 
 97 
	    islot[i + nnz] = i;
- 
 
 98 
	    jslot[i + nnz] = i;
- 
 
 99 
	}
- 
 
 100 
 
- 
 
 101 
	/* build the new x-slot : */
- 
 
 102 
	switch(ctype) {
- 
 
 103 
	case 0: { /* "d" */
- 
 
 104 
	    double *x_new = REAL(ALLOC_SLOT(ans, Matrix_xSym,
- 
 
 105 
					    REALSXP, new_n));
- 
 
 106 
	    Memcpy(x_new, REAL(GET_SLOT(x, Matrix_xSym)), nnz);
- 
 
 107 
	    for(i = 0; i < n; i++) /* add  x[i,i] = 1. */
- 
 
 108 
		x_new[i + nnz] = 1.;
- 
 
 109 
	    break;
- 
 
 110 
	}
- 
 
 111 
	case 1: { /* "l" */
- 
 
 112 
	    int *x_new = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym,
- 
 
 113 
					    LGLSXP, new_n));
- 
 
 114 
	    Memcpy(x_new, LOGICAL(GET_SLOT(x, Matrix_xSym)), nnz);
- 
 
 115 
	    for(i = 0; i < n; i++) /* add  x[i,i] = 1 (= TRUE) */
- 
 
 116 
		x_new[i + nnz] = 1;
- 
 
 117 
	    break;
- 
 
 118 
	}
- 
 
 119 
	case 2: /* "n" */
- 
 
 120 
		/* nothing to do here */
- 
 
 121 
	    break;
- 
 
 122 
 
- 
 
 123 
	case 3: { /* "z" */
- 
 
 124 
	    Rcomplex *x_new = COMPLEX(ALLOC_SLOT(ans, Matrix_xSym,
- 
 
 125 
						 CPLXSXP, new_n));
- 
 
 126 
	    Memcpy(x_new, COMPLEX(GET_SLOT(x, Matrix_xSym)), nnz);
- 
 
 127 
	    for(i = 0; i < n; i++) /* add  x[i,i] = 1 (= TRUE) */
- 
 
 128 
		x_new[i + nnz] = (Rcomplex) {1., 0.};
- 
 
 129 
	    break;
- 
 
 130 
	}
- 
 
 131 
 
- 
 
 132 
	}/* switch() */
- 
 
 133 
 
- 
 
 134 
	UNPROTECT(1);
- 
 
 135 
	return ans;
- 
 
 136 
    }
- 
 
 137 
}