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



    CHM_FR N = AS_CHM_FR(dsCMatrix_Cholesky(x, pivot, ScalarLogical(FALSE),
					    ScalarLogical(FALSE)));
    /* Must use a copy; cholmod_factor_to_sparse modifies first arg. */
    CHM_FR Ncp = cholmod_copy_factor(N, &c);
    CHM_SP L = cholmod_factor_to_sparse(Ncp, &c);
    CHM_SP R = cholmod_transpose(L, /*values*/ 1, &c);
    SEXP ans = PROTECT(chm_sparse_to_SEXP(R, 1/*do_free*/, 1/*uploT*/, 0/*Rkind*/,
					  "N"/*diag*/, GET_SLOT(x, Matrix_DimNamesSym)));
    R_CheckStack();
 
    cholmod_free_factor(&Ncp, &c);
    cholmod_free_sparse(&L, &c);
    if (asLogical(pivot)) {
	SEXP piv = PROTECT(allocVector(INTSXP, N->n));
	int *dest = INTEGER(piv), *src = (int*)N->Perm, i;

    if (LDL) fname[2] = 'D';
    Chol = get_factors(Ap, fname);
    if (Chol != R_NilValue) return Chol; /* return a cached factor */
 
    A = AS_CHM_SP(Ap);
    R_CheckStack();
    if (!A->stype)
	error("Non-symmetric matrix passed to dsCMatrix_chol");
 
    sup = c.supernodal;
    ll = c.final_ll;

SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b)
{
    SEXP Chol = get_factor_pattern(a, "...Cholesky", 3);
    CHM_FR L;
    CHM_SP cx, cb = AS_CHM_SP(b);
    R_CheckStack();
 
    if (Chol == R_NilValue) /* compute (and cache) "sPDCholesky" */
	Chol = dsCMatrix_Cholesky(a,
				  ScalarLogical(1),  /* permuted  : "P" */
				  ScalarLogical(1),  /* LDL'	  : "D" */
				  ScalarLogical(0)); /* simplicial: "s" */
    L = AS_CHM_FR(Chol);
    R_CheckStack();
    cx = cholmod_spsolve(CHOLMOD_A, L, cb, &c);
    return chm_sparse_to_SEXP(cx, /*do_free*/ 1, /*uploT*/ 0,
			      /*Rkind*/ 0, /*diag*/ "N",
			      /*dimnames = */ R_NilValue);
}

SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b)
{
    SEXP Chol = get_factor_pattern(a, "...Cholesky", 3);
    CHM_FR L;
    CHM_DN cx, cb = AS_CHM_DN(PROTECT(mMatrix_as_dgeMatrix(b)));
    R_CheckStack();
 
    if (Chol == R_NilValue) /* compute (and cache) "sPDCholesky" */
	Chol = dsCMatrix_Cholesky(a,
				  ScalarLogical(1),  /* permuted  : "P" */
				  ScalarLogical(1),  /* LDL'      : "D" */
				  ScalarLogical(0)); /* simplicial: "s" */
    L = AS_CHM_FR(Chol);
    R_CheckStack();
    cx = cholmod_solve(CHOLMOD_A, L, cb, &c);
    UNPROTECT(1);
    return chm_dense_to_SEXP(cx, 1, 0, /*dimnames = */ R_NilValue);
}
 

SEXP dsCMatrix_to_dgTMatrix(SEXP x)
{
    CHM_SP A = AS_CHM_SP(x);
    CHM_SP Afull = cholmod_copy(A, /*stype*/ 0, /*mode*/ 1, &c);
    CHM_TR At = cholmod_sparse_to_triplet(Afull, &c);
    R_CheckStack();
 
    if (!A->stype)
	error("Non-symmetric matrix passed to dsCMatrix_to_dgTMatrix");
    cholmod_free_sparse(&Afull, &c);
    return chm_triplet_to_SEXP(At, 1, /*uploT*/ 0, /*Rkind*/ 0, "",

The R Project SVN R-packages

(root)/branches/Matrix-APIchange/src/dsCMatrix.c @ 4560 – Rev 4530 → 4558

Rev 4530		Rev 4558
Line 5...		Line 5...
5	`CHM_FR N = AS_CHM_FR(dsCMatrix_Cholesky(x, pivot, ScalarLogical(FALSE),`	5	`CHM_FR N = AS_CHM_FR(dsCMatrix_Cholesky(x, pivot, ScalarLogical(FALSE),`
6	`ScalarLogical(FALSE)));`	6	`ScalarLogical(FALSE)));`
7	`/* Must use a copy; cholmod_factor_to_sparse modifies first arg. */`	7	`/* Must use a copy; cholmod_factor_to_sparse modifies first arg. */`
8	`CHM_FR Ncp = cholmod_copy_factor(N, &c);`	8	`CHM_FR Ncp = cholmod_copy_factor(N, &c);`
9	`CHM_SP L = cholmod_factor_to_sparse(Ncp, &c);`	9	`CHM_SP L = cholmod_factor_to_sparse(Ncp, &c);`
10	`CHM_SP R = cholmod_transpose(L, /values/ 1, &c);`	10	`CHM_SP R = cholmod_transpose(L, /values/ 1, &c);`
11	`SEXP ans = PROTECT(chm_sparse_to_SEXP(R, 1/do_free/, 1/uploT/, 0/Rkind/,`	11	`SEXP ans = PROTECT(chm_sparse_to_SEXP(R, 1/do_free/, 1/uploT/, 0/Rkind/,`
12	`"N"/diag/, GET_SLOT(x, Matrix_DimNamesSym)));`	12	`"N"/diag/, GET_SLOT(x, Matrix_DimNamesSym)));`
-		13	`R_CheckStack();`
-		14
13	`cholmod_free_factor(&Ncp, &c);`	15	`cholmod_free_factor(&Ncp, &c);`
14	`cholmod_free_sparse(&L, &c);`	16	`cholmod_free_sparse(&L, &c);`
15	`if (asLogical(pivot)) {`	17	`if (asLogical(pivot)) {`
16	`SEXP piv = PROTECT(allocVector(INTSXP, N->n));`	18	`SEXP piv = PROTECT(allocVector(INTSXP, N->n));`
17	`int dest = INTEGER(piv), src = (int*)N->Perm, i;`	19	`int dest = INTEGER(piv), src = (int*)N->Perm, i;`
Line 48...		Line 50...
48	`if (LDL) fname[2] = 'D';`	50	`if (LDL) fname[2] = 'D';`
49	`Chol = get_factors(Ap, fname);`	51	`Chol = get_factors(Ap, fname);`
50	`if (Chol != R_NilValue) return Chol; /* return a cached factor */`	52	`if (Chol != R_NilValue) return Chol; /* return a cached factor */`
51		53
52	`A = AS_CHM_SP(Ap);`	54	`A = AS_CHM_SP(Ap);`
-		55	`R_CheckStack();`
53	`if (!A->stype)`	56	`if (!A->stype)`
54	`error("Non-symmetric matrix passed to dsCMatrix_chol");`	57	`error("Non-symmetric matrix passed to dsCMatrix_chol");`
55		58
56	`sup = c.supernodal;`	59	`sup = c.supernodal;`
57	`ll = c.final_ll;`	60	`ll = c.final_ll;`
Line 99...		Line 102...
99	`SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b)`	102	`SEXP dsCMatrix_Csparse_solve(SEXP a, SEXP b)`
100	`{`	103	`{`
101	`SEXP Chol = get_factor_pattern(a, "...Cholesky", 3);`	104	`SEXP Chol = get_factor_pattern(a, "...Cholesky", 3);`
102	`CHM_FR L;`	105	`CHM_FR L;`
103	`CHM_SP cx, cb = AS_CHM_SP(b);`	106	`CHM_SP cx, cb = AS_CHM_SP(b);`
-		107	`R_CheckStack();`
104		108
105	`if (Chol == R_NilValue) /* compute (and cache) "sPDCholesky" */`	109	`if (Chol == R_NilValue) /* compute (and cache) "sPDCholesky" */`
106	`Chol = dsCMatrix_Cholesky(a,`	110	`Chol = dsCMatrix_Cholesky(a,`
107	`ScalarLogical(1), /* permuted : "P" */`	111	`ScalarLogical(1), /* permuted : "P" */`
108	`ScalarLogical(1), /* LDL' : "D" */`	112	`ScalarLogical(1), /* LDL' : "D" */`
109	`ScalarLogical(0)); /* simplicial: "s" */`	113	`ScalarLogical(0)); /* simplicial: "s" */`
110	`L = AS_CHM_FR(Chol);`	114	`L = AS_CHM_FR(Chol);`
-		115	`R_CheckStack();`
111	`cx = cholmod_spsolve(CHOLMOD_A, L, cb, &c);`	116	`cx = cholmod_spsolve(CHOLMOD_A, L, cb, &c);`
112	`return chm_sparse_to_SEXP(cx, /do_free/ 1, /uploT/ 0,`	117	`return chm_sparse_to_SEXP(cx, /do_free/ 1, /uploT/ 0,`
113	`/Rkind/ 0, /diag/ "N",`	118	`/Rkind/ 0, /diag/ "N",`
114	`/dimnames = / R_NilValue);`	119	`/dimnames = / R_NilValue);`
115	`}`	120	`}`
Line 117...		Line 122...
117	`SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b)`	122	`SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b)`
118	`{`	123	`{`
119	`SEXP Chol = get_factor_pattern(a, "...Cholesky", 3);`	124	`SEXP Chol = get_factor_pattern(a, "...Cholesky", 3);`
120	`CHM_FR L;`	125	`CHM_FR L;`
121	`CHM_DN cx, cb = AS_CHM_DN(PROTECT(mMatrix_as_dgeMatrix(b)));`	126	`CHM_DN cx, cb = AS_CHM_DN(PROTECT(mMatrix_as_dgeMatrix(b)));`
-		127	`R_CheckStack();`
122		128
123	`if (Chol == R_NilValue) /* compute (and cache) "sPDCholesky" */`	129	`if (Chol == R_NilValue) /* compute (and cache) "sPDCholesky" */`
124	`Chol = dsCMatrix_Cholesky(a,`	130	`Chol = dsCMatrix_Cholesky(a,`
125	`ScalarLogical(1), /* permuted : "P" */`	131	`ScalarLogical(1), /* permuted : "P" */`
126	`ScalarLogical(1), /* LDL' : "D" */`	132	`ScalarLogical(1), /* LDL' : "D" */`
127	`ScalarLogical(0)); /* simplicial: "s" */`	133	`ScalarLogical(0)); /* simplicial: "s" */`
128	`L = AS_CHM_FR(Chol);`	134	`L = AS_CHM_FR(Chol);`
-		135	`R_CheckStack();`
129	`cx = cholmod_solve(CHOLMOD_A, L, cb, &c);`	136	`cx = cholmod_solve(CHOLMOD_A, L, cb, &c);`
130	`UNPROTECT(1);`	137	`UNPROTECT(1);`
131	`return chm_dense_to_SEXP(cx, 1, 0, /dimnames = / R_NilValue);`	138	`return chm_dense_to_SEXP(cx, 1, 0, /dimnames = / R_NilValue);`
132	`}`	139	`}`
133		140
Line 137...		Line 144...
137	`SEXP dsCMatrix_to_dgTMatrix(SEXP x)`	144	`SEXP dsCMatrix_to_dgTMatrix(SEXP x)`
138	`{`	145	`{`
139	`CHM_SP A = AS_CHM_SP(x);`	146	`CHM_SP A = AS_CHM_SP(x);`
140	`CHM_SP Afull = cholmod_copy(A, /stype/ 0, /mode/ 1, &c);`	147	`CHM_SP Afull = cholmod_copy(A, /stype/ 0, /mode/ 1, &c);`
141	`CHM_TR At = cholmod_sparse_to_triplet(Afull, &c);`	148	`CHM_TR At = cholmod_sparse_to_triplet(Afull, &c);`
-		149	`R_CheckStack();`
142		150
143	`if (!A->stype)`	151	`if (!A->stype)`
144	`error("Non-symmetric matrix passed to dsCMatrix_to_dgTMatrix");`	152	`error("Non-symmetric matrix passed to dsCMatrix_to_dgTMatrix");`
145	`cholmod_free_sparse(&Afull, &c);`	153	`cholmod_free_sparse(&Afull, &c);`
146	`return chm_triplet_to_SEXP(At, 1, /uploT/ 0, /Rkind/ 0, "",`	154	`return chm_triplet_to_SEXP(At, 1, /uploT/ 0, /Rkind/ 0, "",`