# SCM Repository

[matrix] Diff of /pkg/Matrix/src/dgCMatrix.c
 [matrix] / pkg / Matrix / src / dgCMatrix.c

# Diff of /pkg/Matrix/src/dgCMatrix.c

revision 2888, Wed Aug 7 13:56:43 2013 UTC revision 2889, Thu Aug 8 21:06:22 2013 UTC
# Line 397  Line 397
397      return get_factors(Ap, "LU");      return get_factors(Ap, "LU");
398  }  }
399
400  SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b)  SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse)
401  {  {
402        Rboolean sparse = asLogical(give_sparse);
403        if(sparse) {
404            // FIXME: implement this
405            error(_("dgCMatrix_matrix_solve(.., sparse=TRUE) not yet implemented"));
406
407            /* Idea: in the  for(j = 0; j < nrhs ..) loop below, build the *sparse* result matrix
408             * ----- *column* wise -- which is perfect for dgCMatrix
409             * --> build (i,p,x) slots "increasingly" [well, allocate in batches ..]
410             *
411             * --> maybe first a protoype in R
412             */
413
414        }
415      SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)),      SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)),
416          lu, qslot;          lu, qslot;
417      CSP L, U;      CSP L, U;
# Line 426  Line 439
439          cs_pvec(p, ax + j * n, x, n);  /* x = b(p) */          cs_pvec(p, ax + j * n, x, n);  /* x = b(p) */
440          cs_lsolve(L, x);               /* x = L\x */          cs_lsolve(L, x);               /* x = L\x */
441          cs_usolve(U, x);               /* x = U\x */          cs_usolve(U, x);               /* x = U\x */
442          if (q)                         /* b(q) = x */          if (q)                         /* r(q) = x , hence r = Q' U{^-1} L{^-1} P b = A^{-1} b */
443              cs_ipvec(q, x, ax + j * n, n);              cs_ipvec(q, x, ax + j * n, n);
444          else          else
445              Memcpy(ax + j * n, x, n);              Memcpy(ax + j * n, x, n);

Legend:
 Removed from v.2888 changed lines Added in v.2889