The R Project SVN R-packages

Rev

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

Rev 5010 Rev 5012
Line 19... Line 19...
19 
 * update the numeric values before returning.
 19 
 * update the numeric values before returning.
20 
 *
 20 
 *
21 
 * If no cached copy is available then evaluate one, cache it (for
 21 
 * If no cached copy is available then evaluate one, cache it (for
22 
 * zero Imult), and return a copy.
 22 
 * zero Imult), and return a copy.
23 
 *
 23 
 *
24 
 * @param A      dsCMatrix object
 24 
 * @param Ap     dsCMatrix object
25 
 * @param perm   integer indicating if permutation is required (>0),
 - 
 
26 
 *               forbidden (0) or optional (<0)
 - 
 
27 
 * @param perm   integer indicating if permutation is required (>0),
 25 
 * @param perm   integer indicating if permutation is required (>0),
28 
 *               forbidden (0) or optional (<0)
 26 
 *               forbidden (0) or optional (<0)
29 
 * @param LDL    integer indicating if the LDL' form is required (>0),
 27 
 * @param LDL    integer indicating if the LDL' form is required (>0),
30 
 *               forbidden (0) or optional (<0)
 28 
 *               forbidden (0) or optional (<0)
31 
 * @param super  integer indicating if the supernodal form is required (>0),
 29 
 * @param super  integer indicating if the supernodal form is required (>0),
32 
 *               forbidden (0) or optional (<0)
 30 
 *               forbidden (0) or optional (<0)
- 
 
 31 
 * @param Imult  numeric multiplier of I in  |A + Imult * I|
33 
 */
 32 
 */
34 
static CHM_FR
 33 
static CHM_FR
35 
internal_chm_factor(SEXP Ap, int perm, int LDL, int super, double Imult)
 34 
internal_chm_factor(SEXP Ap, int perm, int LDL, int super, double Imult)
36 
{
 35 
{
37 
    SEXP facs = GET_SLOT(Ap, Matrix_factorSym);
 36 
    SEXP facs = GET_SLOT(Ap, Matrix_factorSym);
Line 127... Line 126...
127 
 * Fast version of getting at the diagonal matrix D of the
 126 
 * Fast version of getting at the diagonal matrix D of the
128 
 * (generalized) simplicial Cholesky LDL' decomposition of a
 127 
 * (generalized) simplicial Cholesky LDL' decomposition of a
129 
 * (sparse symmetric) dsCMatrix.
 128 
 * (sparse symmetric) dsCMatrix.
130 
 *
 129 
 *
131 
 * @param Ap  symmetric CsparseMatrix
 130 
 * @param Ap  symmetric CsparseMatrix
132 
 * @param permp  logical indicating if permutation is allowed
 131 
 * @param permpP  logical indicating if permutation is allowed
133 
 *
 132 
 *
134 
 * @return SEXP containing either the vector diagonal entries of D,
 133 
 * @return SEXP containing either the vector diagonal entries of D,
135 
 *         or just  sum_i D[i], prod_i D[i] or  sum_i log(D[i]).
 134 
 *         or just  sum_i D[i], prod_i D[i] or  sum_i log(D[i]).
136 
 */
 135 
 */
137 
/* Started as copy + modification of dsCMatrix_Cholesky();
 - 
 
138 
 * builds strongly on diag_tC(...., SEXP resultKind) in Csparse.c
 - 
 
139 
*/
 - 
 
140 
SEXP dsCMatrix_LDL_D(SEXP Ap, SEXP permP, SEXP resultKind)
 136 
SEXP dsCMatrix_LDL_D(SEXP Ap, SEXP permP, SEXP resultKind)
141 
{
 137 
{
142 
    char fname[12] = "spDCholesky"; /* template for factorization name */
 138 
    CHM_FR L = internal_chm_factor(Ap, asLogical(permP),
143 
    /* S|s : super or not
 139 
			   /*LDL*/ 1, /*super*/0, /*Imult*/0.);
144 
     * P|p : permuted or not
 140 
    return diag_tC_ptr(L->n,
145 
     * D|d :  LDL' or not (= LL')
 - 
 
146 
     */
 - 
 
147 
    const int perm = asLogical(permP), LDL = 1, super = 0;
 - 
 
148 
    SEXP Chol;
 141 
		       L->p,
149 
    CHM_SP A;
 - 
 
150 
    CHM_FR L;
 142 
		       L->x,
151 
    int sup, ll;
 143 
		       L->Perm,
152 
 
 - 
 
153 
    /* if (super) fname[0] = 'S'; */
 - 
 
154 
    if (perm) fname[1] = 'P';
 - 
 
155 
    /* if (LDL) fname[2] = 'D'; */
 - 
 
156 
    Chol = get_factors(Ap, fname);
 - 
 
157 
    if (Chol != R_NilValue) { /* use the *cached* Cholesky() factor */
 - 
 
158 
	return diag_tC(GET_SLOT(Chol, Matrix_pSym),
 - 
 
159 
		       GET_SLOT(Chol, Matrix_xSym),
 - 
 
160 
		       GET_SLOT(Chol, Matrix_permSym),
 - 
 
161 
		       resultKind);
 144 
		       resultKind);
162 
    }
 - 
 
163 
    else {
 - 
 
164 
	A = AS_CHM_SP(Ap);
 - 
 
165 
	R_CheckStack();
 - 
 
166 
	if (!A->stype)
 - 
 
167 
	    error("Non-symmetric matrix passed to dsCMatrix_LDL_D");
 - 
 
168 
 
 - 
 
169 
	sup = c.supernodal;
 - 
 
170 
	ll = c.final_ll;
 - 
 
171 
 
 - 
 
172 
	c.final_ll = !LDL;	/* leave as LL' or form LDL' */
 - 
 
173 
	c.supernodal = super ? CHOLMOD_SUPERNODAL : CHOLMOD_SIMPLICIAL;
 - 
 
174 
 
 - 
 
175 
	if (perm) {
 - 
 
176 
	    L = cholmod_analyze(A, &c); /* get fill-reducing permutation */
 - 
 
177 
	} else {			/* require identity permutation */
 - 
 
178 
	    int nmethods = c.nmethods,
 - 
 
179 
		ord0 = c.method[0].ordering, postorder = c.postorder;
 - 
 
180 
	    c.nmethods = 1;
 - 
 
181 
	    c.method[0].ordering = CHOLMOD_NATURAL; c.postorder = FALSE;
 - 
 
182 
	    L = cholmod_analyze(A, &c);
 - 
 
183 
	    c.nmethods = nmethods;
 - 
 
184 
	    c.method[0].ordering = ord0;            c.postorder = postorder;
 - 
 
185 
	}
 - 
 
186 
	if (!cholmod_factorize(A, L, &c))
 - 
 
187 
	    error(_("Cholesky factorization failed"));
 - 
 
188 
	/* restore previous setting */
 - 
 
189 
	c.supernodal = sup;
 - 
 
190 
	c.final_ll = ll;
 - 
 
191 
 
 - 
 
192 
	/* Now, instead of
 - 
 
193 
	 *   Chol = set_factors(Ap, chm_factor_to_SEXP(L, 1), fname);
 - 
 
194 
	 *   return Chol;
 - 
 
195 
	 *
 - 
 
196 
	 * get the correct entries from the CHOLMOD 'L' struct directly : */
 - 
 
197 
 
 - 
 
198 
	return diag_tC_ptr(L->n,
 - 
 
199 
			   L->p,
 - 
 
200 
			   L->x,
 - 
 
201 
			   L->Perm,
 - 
 
202 
			   resultKind);
 - 
 
203 
    }
 - 
 
204 
}
 145 
}
205 
 
 146 
 
206 
SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b)
 147 
SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b)
207 
{
 148 
{
208 
    CHM_FR L = internal_chm_factor(a, -1, -1, -1, 0.);
 149 
    CHM_FR L = internal_chm_factor(a, -1, -1, -1, 0.);