SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/src/dtCMatrix.c
ViewVC logotype

Diff of /pkg/Matrix/src/dtCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

pkg/src/tscMatrix.c revision 70, Mon Apr 12 12:10:01 2004 UTC pkg/Matrix/src/dtCMatrix.c revision 2889, Thu Aug 8 21:06:22 2013 UTC
# Line 1  Line 1 
1                                  /* Sparse triangular matrices */                                  /* Sparse triangular numeric matrices */
2  #include "tscMatrix.h"  #include "dtCMatrix.h"
3    #include "cs_utils.h"
4    
5  SEXP tsc_validate(SEXP x)  #define RETURN(_CH_)   UNPROTECT(1); return (_CH_);
6    
7    /* This is used for *BOTH* triangular and symmetric Csparse: */
8    SEXP tCMatrix_validate(SEXP x)
9  {  {
10      return ScalarLogical(1);      SEXP val = xCMatrix_validate(x);/* checks x slot */
11        if(isString(val))
12            return(val);
13        else {
14            SEXP
15                islot = GET_SLOT(x, Matrix_iSym),
16                pslot = GET_SLOT(x, Matrix_pSym);
17            int uploT = (*uplo_P(x) == 'U'),
18                k, nnz = length(islot),
19                *xi = INTEGER(islot),
20                *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
21    
22            expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23    
24            /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
25            if(uploT) {
26                for (k = 0; k < nnz; k++)
27                    if(xi[k] > xj[k]) {
28                        RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
29                    }
30            }
31            else {
32                for (k = 0; k < nnz; k++)
33                    if(xi[k] < xj[k]) {
34                        RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
35                    }
36  }  }
37    
38  SEXP tsc_transpose(SEXP x)          RETURN(ScalarLogical(1));
39        }
40    }
41    
42    /* This is used for *BOTH* triangular and symmetric Rsparse: */
43    SEXP tRMatrix_validate(SEXP x)
44  {  {
45        SEXP val = xRMatrix_validate(x);/* checks x slot */
46        if(isString(val))
47            return(val);
48        else {
49      SEXP      SEXP
50          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix"))),              jslot = GET_SLOT(x, Matrix_jSym),
51          islot = GET_SLOT(x, Matrix_iSym);              pslot = GET_SLOT(x, Matrix_pSym);
52      int nnz = length(islot),          int uploT = (*uplo_P(x) == 'U'),
53          *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),              k, nnz = length(jslot),
54          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));              *xj = INTEGER(jslot),
55                *xi = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
56      adims[0] = xdims[1]; adims[1] = xdims[0];  
57      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')          expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xi);
58          SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));  
59      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
60      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));          if(uploT) {
61      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));              for (k = 0; k < nnz; k++)
62      csc_components_transpose(xdims[0], xdims[1], nnz,                  if(xi[k] > xj[k]) {
63                               INTEGER(GET_SLOT(x, Matrix_pSym)),                      RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
64                               INTEGER(islot),                  }
65                               REAL(GET_SLOT(x, Matrix_xSym)),          }
66                               INTEGER(GET_SLOT(ans, Matrix_pSym)),          else {
67                               INTEGER(GET_SLOT(ans, Matrix_iSym)),              for (k = 0; k < nnz; k++)
68                               REAL(GET_SLOT(ans, Matrix_xSym)));                  if(xi[k] < xj[k]) {
69      UNPROTECT(1);                      RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
70      return ans;                  }
71  }  }
72    
73  SEXP tsc_to_triplet(SEXP x)          RETURN(ScalarLogical(1));
74        }
75    }
76    
77    SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
78  {  {
79      SEXP ans;      int cl = asLogical(classed);
80      if (toupper(CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0]) != 'U')      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
81          ans = csc_to_triplet(x);      CSP A = AS_CSP(a);
82      else {                      /* unit triangular matrix */      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
83          SEXP islot = GET_SLOT(x, Matrix_iSym),          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
84              pslot = GET_SLOT(x, Matrix_pSym);                           getAttrib(b, R_DimSymbol));
85          int *ai, *aj, j,      int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
86              n = length(pslot) - 1,      double *bx;
87              nod = length(islot),      R_CheckStack();
88              nout = n + nod,  
89              *p = INTEGER(pslot);      if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
90          double *ax;          error(_("Dimensions of system to be solved are inconsistent"));
91        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
92          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix")));      /* FIXME: copy dimnames or Dimnames as well */
93          SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));      bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
94          SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));                  REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
95          ai = INTEGER(GET_SLOT(ans, Matrix_iSym));      for (j = 0; j < nrhs; j++)
96          Memcpy(ai, INTEGER(islot), nod);          lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
97          SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));      RETURN(ans);
98          aj = INTEGER(GET_SLOT(ans, Matrix_jSym));  }
99          SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));  
100          ax = REAL(GET_SLOT(ans, Matrix_xSym));  SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
101          Memcpy(ax, REAL(GET_SLOT(x, Matrix_xSym)), nod);  {
102          for (j = 0; j < n; j++) {      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
103              int jj, npj = nod + j,  p2 = p[j+1];      CSP A = AS_CSP(a), B = AS_CSP(b);
104              aj[npj] = ai[npj] = j;      R_CheckStack();
105              ax[npj] = 1.;      if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
106              for (jj = p[j]; jj < p2; jj++) aj[jj] = j;          error(_("Dimensions of system to be solved are inconsistent"));
107        // *before* Calloc()ing below [memory leak]!
108    
109        int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
110            xnz = 10 * B->p[B->n];  /* initial estimate of nnz in x */
111        int k, lo = uplo_P(a)[0] == 'L', pos = 0;
112        int    *ti = Calloc(xnz, int),     *xi = Calloc(2*A->n, int); /* for cs_reach */
113        double *tx = Calloc(xnz, double), *wrk = Calloc(  A->n, double);
114    
115        slot_dup(ans, b, Matrix_DimSym);
116        SET_DimNames(ans, b);
117        xp[0] = 0;
118        for (k = 0; k < B->n; k++) {
119            int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
120            int nz = A->n - top;
121    
122            xp[k + 1] = nz + xp[k];
123            if (xp[k + 1] > xnz) {
124                while (xp[k + 1] > xnz) xnz *= 2;
125                ti = Realloc(ti, xnz, int);
126                tx = Realloc(tx, xnz, double);
127            }
128            if (lo)                 /* increasing row order */
129                for(int p = top; p < A->n; p++, pos++) {
130                    ti[pos] = xi[p];
131                    tx[pos] = wrk[xi[p]];
132          }          }
133          UNPROTECT(1);          else                    /* decreasing order, reverse copy */
134                for(int p = A->n - 1; p >= top; p--, pos++) {
135                    ti[pos] = xi[p];
136                    tx[pos] = wrk[xi[p]];
137      }      }
     return ans;  
138  }  }
139        xnz = xp[B->n];
140        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  xnz)), ti, xnz);
141        Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);
142    
143        Free(ti);  Free(tx);
144        Free(wrk); Free(xi);
145    
146        RETURN(ans);
147    }
148    #undef RETURN
149    

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

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge