SCM

SCM Repository

[matrix] Diff of /branches/Matrix-mer2/src/dgCMatrix.c
ViewVC logotype

Diff of /branches/Matrix-mer2/src/dgCMatrix.c

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

pkg/src/cscMatrix.c revision 371, Sun Dec 5 13:52:25 2004 UTC pkg/src/dgCMatrix.c revision 956, Fri Sep 30 17:28:00 2005 UTC
# Line 1  Line 1 
1  #include "cscMatrix.h"  #include "dgCMatrix.h"
2    
3  SEXP csc_validate(SEXP x)  #ifdef USE_CHOLMOD
4    #include "chm_common.h"
5    #endif
6    
7    /* FIXME -- we "forget" about dimnames almost everywhere : */
8    
9    SEXP dgCMatrix_validate(SEXP x)
10  {  {
11      SEXP pslot = GET_SLOT(x, Matrix_pSym),      SEXP pslot = GET_SLOT(x, Matrix_pSym),
12          islot = GET_SLOT(x, Matrix_iSym),          islot = GET_SLOT(x, Matrix_iSym),
# Line 14  Line 20 
20    
21      nrow = dims[0];      nrow = dims[0];
22      if (length(islot) != length(xslot))      if (length(islot) != length(xslot))
23          return ScalarString(mkChar("lengths of slots i and x must match"));          return mkString(_("lengths of slots i and x must match"));
24      if (length(pslot) <= 0)      if (length(pslot) <= 0)
25          return ScalarString(mkChar("slot p must have length > 0"));          return mkString(_("slot p must have length > 0"));
26      if (xp[0] != 0)      if (xp[0] != 0)
27          return ScalarString(mkChar("first element of slot p must be zero"));          return mkString(_("first element of slot p must be zero"));
28      if (length(islot) != xp[ncol])      if (length(islot) != xp[ncol])
29          return ScalarString(          return mkString(_("last element of slot p must match length of slots i and x"));
             mkChar(  
                 "last element of slot p must match length of slots i and x"));  
30      for (j = 0; j < ncol; j++) {      for (j = 0; j < ncol; j++) {
31          if (xp[j] > xp[j+1])          if (xp[j] > xp[j+1])
32              return ScalarString(mkChar("slot p must be non-decreasing"));              return mkString(_("slot p must be non-decreasing"));
33      }      }
34      for (j = 0; j < length(islot); j++) {      for (j = 0; j < length(islot); j++) {
35          if (xi[j] < 0 || xi[j] >= nrow)          if (xi[j] < 0 || xi[j] >= nrow)
36              return ScalarString(              return mkString(_("all row indices must be between 0 and nrow-1"));
                 mkChar("all row indices must be between 0 and nrow-1"));  
37      }      }
38      if (csc_unsorted_columns(ncol, xp, xi)) {      if (csc_unsorted_columns(ncol, xp, xi))
39          csc_sort_columns(ncol, xp, xi, REAL(xslot));          csc_sort_columns(ncol, xp, xi, REAL(xslot));
40      }  
41      return ScalarLogical(1);      return ScalarLogical(1);
42  }  }
43    
44  SEXP csc_crossprod(SEXP x)  SEXP csc_crossprod(SEXP x)
45  {  {
46      SEXP pslot = GET_SLOT(x, Matrix_pSym),      SEXP pslot = GET_SLOT(x, Matrix_pSym),
47          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))), tmp;          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))), tmp;
48      int *xp = INTEGER(pslot),      int *xp = INTEGER(pslot),
49          *xi = INTEGER(GET_SLOT(x, Matrix_iSym));          *xi = INTEGER(GET_SLOT(x, Matrix_iSym));
50      double *xx = REAL(GET_SLOT(x, Matrix_xSym));      double *xx = REAL(GET_SLOT(x, Matrix_xSym));
# Line 49  Line 52 
52      int j, *iVal, ncol = length(pslot) - 1, maxnz, nnz = 0, *pVal;      int j, *iVal, ncol = length(pslot) - 1, maxnz, nnz = 0, *pVal;
53      double *xVal;      double *xVal;
54    
55      SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));      SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));
56        SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
57        SET_SLOT(ans, Matrix_uploSym, mkString("L"));
58      maxnz = (ncol * (ncol + 1))/2;      maxnz = (ncol * (ncol + 1))/2;
59      iVal = Calloc(maxnz, int); xVal = Calloc(maxnz, double);      iVal = Calloc(maxnz, int); xVal = Calloc(maxnz, double);
60      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, ncol + 1));      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, ncol + 1));
# Line 98  Line 103 
103      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
104      Memcpy(REAL(GET_SLOT(ans, Matrix_xSym)), xVal, nnz);      Memcpy(REAL(GET_SLOT(ans, Matrix_xSym)), xVal, nnz);
105      Free(iVal); Free(xVal); UNPROTECT(1);      Free(iVal); Free(xVal); UNPROTECT(1);
106      return cscMatrix_set_Dim(ans, ncol);      return dgCMatrix_set_Dim(ans, ncol);
107  }  }
108    
109  SEXP csc_tcrossprod(SEXP x)  SEXP csc_tcrossprod(SEXP x)
110  {  {
111    #ifdef USE_CHOLMOD
112        cholmod_sparse *cha = cholmod_aat(as_cholmod_sparse(x),
113            (int *) NULL, 0, 1, &c);
114    
115        cha->stype = -1;            /* set the symmetry */
116        cholmod_sort(cha, &c);      /* drop redundant entries */
117        return chm_sparse_to_SEXP(cha, -1);
118    #else
119    
120      SEXP pslot = GET_SLOT(x, Matrix_pSym),      SEXP pslot = GET_SLOT(x, Matrix_pSym),
121          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix")));          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));
122      int *xp = INTEGER(pslot),      int *xp = INTEGER(pslot),
123          *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),          *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
124          *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));          *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
# Line 114  Line 128 
128      int *itmp, *ansp;      int *itmp, *ansp;
129      double *xVal, *xtmp;      double *xVal, *xtmp;
130    
131      SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));      SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));
132        SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
133      ntrip = nrow;               /* number of triplets */      ntrip = nrow;               /* number of triplets */
134      for (j = 0; j < ncol; j++) {      for (j = 0; j < ncol; j++) {
135          int nzj = xp[j+1] - xp[j];          int nzj = xp[j+1] - xp[j];
# Line 149  Line 164 
164      triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,      triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,
165                     ansp, itmp, xtmp);                     ansp, itmp, xtmp);
166      nnz = ansp[nrow];      nnz = ansp[nrow];
167      SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));      SET_SLOT(ans, Matrix_uploSym, mkString("L"));
168      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
169      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
170      Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);      Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);
171      Memcpy(REAL(GET_SLOT(ans, Matrix_xSym)), xtmp, nnz);      Memcpy(REAL(GET_SLOT(ans, Matrix_xSym)), xtmp, nnz);
     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));  
172      dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));      dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
173      dims[0] = dims[1] = nrow;      dims[0] = dims[1] = nrow;
174      Free(itmp); Free(xtmp); Free(iVal); Free(jVal); Free(xVal);      Free(itmp); Free(xtmp); Free(iVal); Free(jVal); Free(xVal);
175      UNPROTECT(1);      UNPROTECT(1);
176      return ans;      return ans;
177    #endif /* USE_CHOLMOD */
178  }  }
179    
180  SEXP csc_matrix_crossprod(SEXP x, SEXP y)  SEXP csc_matrix_crossprod(SEXP x, SEXP y, SEXP classed)
181  {  {
182      SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;      int cl = asLogical(classed);
183      int j,      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
184          *xp = INTEGER(pslot),      int *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
185          *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),          *ydims = INTEGER(cl ? GET_SLOT(y, Matrix_DimSym) :
186          xncol = length(pslot) - 1,                           getAttrib(y, R_DimSymbol)),
187          xnrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],          *vdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
188          *ydims;      int *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
189      double *xx = REAL(GET_SLOT(x, Matrix_xSym));          *xp = INTEGER(GET_SLOT(x, Matrix_pSym));
190        int j, k = xdims[0], m = xdims[1], n = ydims[1];
191      if (!(isMatrix(y) && isReal(y))) error("y must be a numeric matrix");      double *vx, *xx = REAL(GET_SLOT(x, Matrix_xSym)),
192      ydims = INTEGER(getAttrib(y, R_DimSymbol));          *yx = REAL(cl ? GET_SLOT(y, Matrix_xSym) : y);
193      if (xnrow != ydims[0]) error("x and y must have the same number of rows");  
194      ans = PROTECT(allocMatrix(REALSXP, xncol, ydims[1]));      if (!cl && !(isMatrix(y) && isReal(y)))
195      for (j = 0; j < ydims[1]; j++) {          error(_("y must be a numeric matrix"));
196          int i; double *ypt = REAL(y) + j * ydims[0];      if (ydims[0] != k)
197          for(i = 0; i < xncol; i++) {          error(_("x and y must have the same number of rows"));
198        if (m < 1 || n < 1 || k < 1)
199            error(_("Matrices with zero extents cannot be multiplied"));
200        vdims[0] = m; vdims[1] = n;
201        vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
202        for (j = 0; j < n; j++) {
203            int i; double *ypt = yx + j * k;
204            for(i = 0; i < m; i++) {
205              int ii; double accum = 0.;              int ii; double accum = 0.;
206              for (ii = xp[i]; ii < xp[i+1]; ii++) {              for (ii = xp[i]; ii < xp[i+1]; ii++) {
207                  accum += xx[ii] * ypt[xi[ii]];                  accum += xx[ii] * ypt[xi[ii]];
208              }              }
209              REAL(ans)[i + j * xncol] = accum;              vx[i + j * m] = accum;
210          }          }
211      }      }
212      UNPROTECT(1);      UNPROTECT(1);
213      return ans;      return val;
214  }  }
215    
216  SEXP csc_to_triplet(SEXP x)  SEXP compressed_to_dgTMatrix(SEXP x, SEXP colP)
217  {  {
218      SEXP      int col = asLogical(colP); /* 1 if "C"olumn compressed;  0 if "R"ow */
219          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix"))),      SEXP indSym = col ? Matrix_iSym : Matrix_jSym;
220          dimslot = GET_SLOT(x, Matrix_DimSym),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
221          islot = GET_SLOT(x, Matrix_iSym),          indP = GET_SLOT(x, indSym),
222          pslot = GET_SLOT(x, Matrix_pSym);          pP = GET_SLOT(x, Matrix_pSym);
223      int *dims = INTEGER(dimslot), j, jj,      int npt = length(pP) - 1;
         *xp = INTEGER(pslot), *yj;  
224    
225      SET_SLOT(ans, Matrix_iSym, duplicate(islot));      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
     SET_SLOT(ans, Matrix_DimSym, duplicate(dimslot));  
226      SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
227      SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, length(islot)));      SET_SLOT(ans, indSym, duplicate(indP));
228      yj = INTEGER(GET_SLOT(ans, Matrix_jSym));      expand_cmprPt(npt, INTEGER(pP),
229      jj = 0;                    INTEGER(ALLOC_SLOT(ans, col ? Matrix_jSym : Matrix_iSym,
230      for (j = 0; j < dims[1]; j++) {                                       INTSXP, length(indP))));
231          while (jj < xp[j + 1]) {      UNPROTECT(1);
232              yj[jj] = j;      return ans;
             jj++;  
         }  
233      }      }
234    
235    SEXP compressed_non_0_ij(SEXP x, SEXP colP)
236    {
237        int col = asLogical(colP); /* 1 if "C"olumn compressed;  0 if "R"ow */
238        SEXP ans, indSym = col ? Matrix_iSym : Matrix_jSym;
239        SEXP indP = GET_SLOT(x, indSym),
240            pP = GET_SLOT(x, Matrix_pSym);
241        int n_el = length(indP), i, *ij;
242    
243        ij = INTEGER(ans = PROTECT(allocMatrix(INTSXP, n_el, 2)));
244        /* expand the compressed margin to 'i' or 'j' : */
245        expand_cmprPt(length(pP) - 1, INTEGER(pP), &ij[col ? n_el : 0]);
246        /* and copy the other one: */
247        if (col)
248            for(i = 0; i < n_el; i++)
249                ij[i] = INTEGER(indP)[i];
250        else /* row compressed */
251            for(i = 0; i < n_el; i++)
252                ij[i + n_el] = INTEGER(indP)[i];
253    
254      UNPROTECT(1);      UNPROTECT(1);
255      return ans;      return ans;
256  }  }
# Line 238  Line 276 
276      return ans;      return ans;
277  }  }
278    
279  SEXP csc_to_geMatrix(SEXP x)  SEXP csc_to_dgeMatrix(SEXP x)
280  {  {
281      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
282          Dimslot = GET_SLOT(x, Matrix_DimSym);          Dimslot = GET_SLOT(x, Matrix_DimSym);
283      int *dims = INTEGER(Dimslot),      int *dims = INTEGER(Dimslot),
284          *xp = INTEGER(GET_SLOT(x, Matrix_pSym)),          *xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
# Line 251  Line 289 
289      SET_SLOT(ans, Matrix_DimSym, duplicate(Dimslot));      SET_SLOT(ans, Matrix_DimSym, duplicate(Dimslot));
290      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nrow*ncol));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nrow*ncol));
291      SET_SLOT(ans, Matrix_rcondSym, allocVector(REALSXP, 0));      SET_SLOT(ans, Matrix_rcondSym, allocVector(REALSXP, 0));
292      SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));      SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));
293      ax = REAL(GET_SLOT(ans, Matrix_xSym));      ax = REAL(GET_SLOT(ans, Matrix_xSym));
294      for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;      for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;
295      for (j = 0; j < ncol; j++) {      for (j = 0; j < ncol; j++) {
# Line 264  Line 302 
302      return ans;      return ans;
303  }  }
304    
305  SEXP matrix_to_csc(SEXP A)  SEXP double_to_csc(double *a, int *dim_a)
306  {  {
307      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
308      int *adims = INTEGER(getAttrib(A, R_DimSymbol)), j,      int j, maxnz, nrow, ncol, nnz, *vp, *vi;
         maxnz, nrow, ncol, nnz, *vp, *vi;  
   
309      double *vx;      double *vx;
310    
311      if (!(isMatrix(A) && isReal(A)))      nrow = dim_a[0]; ncol = dim_a[1];
312          error("A must be a numeric matrix");      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
313      nrow = adims[0]; ncol = adims[1];      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));  
314      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
315      vp = INTEGER(GET_SLOT(val, Matrix_pSym));      vp = INTEGER(GET_SLOT(val, Matrix_pSym));
316      maxnz = nrow * ncol;      maxnz = nrow * ncol;
# Line 285  Line 320 
320          int i;          int i;
321          vp[j] = nnz;          vp[j] = nnz;
322          for (i = 0; i < nrow; i++) {          for (i = 0; i < nrow; i++) {
323              double val = REAL(A)[i + j * nrow];              double val = a[i + j * nrow];
324              if (val != 0.) {              if (val != 0.) {
325                  vi[nnz] = i;                  vi[nnz] = i;
326                  vx[nnz] = val;                  vx[nnz] = val;
# Line 300  Line 335 
335      Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), vx, nnz);      Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), vx, nnz);
336      Free(vi); Free(vx);      Free(vi); Free(vx);
337      UNPROTECT(1);      UNPROTECT(1);
338      return cscMatrix_set_Dim(val, nrow);      return dgCMatrix_set_Dim(val, nrow);
339  }  }
340    
341    SEXP matrix_to_csc(SEXP A)
342    {
343        if (!(isMatrix(A) && isReal(A)))
344            error(_("A must be a numeric matrix"));
345        return double_to_csc(REAL(A),
346                             INTEGER(getAttrib(A, R_DimSymbol)));
347    }
348    
349  SEXP triplet_to_csc(SEXP triplet)  SEXP dgeMatrix_to_csc(SEXP x)
350  {  {
351      SEXP Tisl = GET_SLOT(triplet, Matrix_iSym);      return double_to_csc(   REAL(GET_SLOT(x, Matrix_xSym)),
352                             INTEGER(GET_SLOT(x, Matrix_DimSym)));
353    }
354    
355    
356    
357    SEXP dgTMatrix_to_csc(SEXP dgTMatrix)
358    {
359        SEXP Tisl = GET_SLOT(dgTMatrix, Matrix_iSym);
360      int *Ti = INTEGER(Tisl),      int *Ti = INTEGER(Tisl),
361          *Tj = INTEGER(GET_SLOT(triplet, Matrix_jSym)),          *Tj = INTEGER(GET_SLOT(dgTMatrix, Matrix_jSym)),
362          i, nrow, ncol,          i, nrow, ncol,
363          nz = length(Tisl);          nz = length(Tisl);
364    
# Line 318  Line 368 
368          if (Tj[i] > ncol) ncol = Tj[i];          if (Tj[i] > ncol) ncol = Tj[i];
369      }      }
370      return triple_as_SEXP(nrow + 1, ncol + 1, nz, Ti, Tj,      return triple_as_SEXP(nrow + 1, ncol + 1, nz, Ti, Tj,
371                            REAL(GET_SLOT(triplet, Matrix_xSym)),                            REAL(GET_SLOT(dgTMatrix, Matrix_xSym)),
372                            "cscMatrix");                            "dgCMatrix");
373  }  }
374    
375  SEXP csc_getDiag(SEXP x)  SEXP csc_getDiag(SEXP x)
# Line 349  Line 399 
399    
400  SEXP csc_transpose(SEXP x)  SEXP csc_transpose(SEXP x)
401  {  {
402      SEXP  #ifdef USE_CHOLMOD
403          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix"))),      cholmod_sparse *chx = as_cholmod_sparse(x);
404          islot = GET_SLOT(x, Matrix_iSym);      SEXP ans =
405      int nnz = length(islot),          chm_sparse_to_SEXP(cholmod_transpose(chx, 1, &c), 1);
406          *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),      Free(chx);
407          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));      return ans;
408    #else
409        SEXP xi = GET_SLOT(x, Matrix_iSym);
410        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
411        int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
412            *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
413            nz = length(xi);
414        int *xj = Calloc(nz, int);
415        SEXP adn = ALLOC_SLOT(ans, Matrix_DimNamesSym, VECSXP, 2),
416            xdn = GET_SLOT(x, Matrix_DimNamesSym);
417    
418      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
419      SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));      SET_VECTOR_ELT(adn, 0, VECTOR_ELT(xdn, 1));
420      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));      SET_VECTOR_ELT(adn, 1, VECTOR_ELT(xdn, 0));
421      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      triplet_to_col(adims[0], adims[1], nz,
422      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));                     expand_cmprPt(xdims[1], INTEGER(GET_SLOT(x, Matrix_pSym)), xj),
423      csc_components_transpose(xdims[0], xdims[1], nnz,                     INTEGER(xi),
                              INTEGER(GET_SLOT(x, Matrix_pSym)),  
                              INTEGER(islot),  
424                               REAL(GET_SLOT(x, Matrix_xSym)),                               REAL(GET_SLOT(x, Matrix_xSym)),
425                               INTEGER(GET_SLOT(ans, Matrix_pSym)),                     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, adims[1] + 1)),
426                               INTEGER(GET_SLOT(ans, Matrix_iSym)),                     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
427                               REAL(GET_SLOT(ans, Matrix_xSym)));                     REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)));
428        Free(xj);
429      UNPROTECT(1);      UNPROTECT(1);
430      return ans;      return ans;
431    #endif /* USE_CHOLMOD */
432  }  }
433    
434  SEXP csc_matrix_mm(SEXP a, SEXP b)  SEXP csc_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
435  {  {
436      int *adim = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int cl = asLogical(classed), rt = asLogical(right);
437        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
438        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
439          *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),          *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),
440          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
441          *bdim = INTEGER(getAttrib(b, R_DimSymbol));          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
442      int j, k, m = adim[0], n = bdim[1], r = adim[1];                           getAttrib(b, R_DimSymbol)),
443      double *ax = REAL(GET_SLOT(a, Matrix_xSym));          *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)),
444      SEXP val;          chk, ione = 1, j, jj, k, m, n;
445        double *ax = REAL(GET_SLOT(a, Matrix_xSym)),
446      if (bdim[0] != r)          *bx = REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), *cx;
447          error("Matrices of sizes (%d,%d) and (%d,%d) cannot be multiplied",  
448                m, r, bdim[0], n);      if (rt) {
449      val = PROTECT(allocMatrix(REALSXP, m, n));          m = bdims[0]; n = adims[1]; k = bdims[1]; chk = adims[0];
450      for (j = 0; j < n; j++) {   /* across columns of b */      } else {
451          double *ccol = REAL(val) + j * m,          m = adims[0]; n = bdims[1]; k = adims[1]; chk = bdims[0];
452              *bcol = REAL(b) + j * r;      }
453        if (chk != k)
454          for (k = 0; k < m; k++) ccol[k] = 0.; /* zero the accumulators */          error(_("Matrices are not conformable for multiplication"));
455          for (k = 0; k < r; k++) { /* across columns of a */      if (m < 1 || n < 1 || k < 1)
456              int kk, k2 = ap[k + 1];          error(_("Matrices with zero extents cannot be multiplied"));
457              for (kk = ap[k]; kk < k2; kk++) {      cx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
458                  ccol[ai[kk]] += ax[kk] * bcol[k];      AZERO(cx, m * n); /* zero the accumulators */
459        for (j = 0; j < n; j++) { /* across columns of c */
460            if (rt) {
461                int kk, k2 = ap[j + 1];
462                for (kk = ap[j]; kk < k2; kk++) {
463                    F77_CALL(daxpy)(&m, &ax[kk], &bx[ai[kk]*m],
464                                    &ione, &cx[j*m], &ione);
465                }
466            } else {
467                double *ccol = cx + j * m,
468                    *bcol = bx + j * k;
469    
470                for (jj = 0; jj < k; jj++) { /* across columns of a */
471                    int kk, k2 = ap[jj + 1];
472                    for (kk = ap[jj]; kk < k2; kk++) {
473                        ccol[ai[kk]] += ax[kk] * bcol[jj];
474                    }
475              }              }
476          }          }
477      }      }
478        cdims[0] = m; cdims[1] = n;
479      UNPROTECT(1);      UNPROTECT(1);
480      return val;      return val;
481  }  }
482    
483  SEXP csc_col_permute(SEXP x, SEXP perm)  SEXP csc_col_permute(SEXP x, SEXP perm)
484  {  {
485      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix"))), tmp;      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix"))), tmp;
486      int *iperm, *prm, *vi, *vp, *xi, *xp, j, k, ncol, pos;      int *iperm, *prm, *vi, *vp, *xi, *xp, j, k, ncol, pos;
487      double *vx, *xx;      double *vx, *xx;
488    
489      SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
490      tmp = GET_SLOT(x, Matrix_DimSym);      tmp = GET_SLOT(x, Matrix_DimSym);
491      SET_SLOT(val, Matrix_DimSym, duplicate(tmp));      SET_SLOT(val, Matrix_DimSym, duplicate(tmp));
492      ncol = INTEGER(tmp)[1];      ncol = INTEGER(tmp)[1];
493      if (!(isInteger(perm) && length(perm) == ncol))      if (!(isInteger(perm) && length(perm) == ncol))
494          error("perm must be an integer vector of length %d",          error(_("perm must be an integer vector of length %d"),
495                ncol);                ncol);
496      prm = INTEGER(perm);      prm = INTEGER(perm);
497      if (!R_ldl_valid_perm(ncol, prm))      if (!R_ldl_valid_perm(ncol, prm))
498          error("perm is not a valid 0-based permutation");          error(_("perm is not a valid 0-based permutation"));
499      iperm = Calloc(ncol, int);      iperm = Calloc(ncol, int);
500      for (j = 0; j < ncol; j++) iperm[prm[j]] = j;      for (j = 0; j < ncol; j++) iperm[prm[j]] = j;
501      tmp = GET_SLOT(x, Matrix_pSym);      tmp = GET_SLOT(x, Matrix_pSym);
# Line 448  Line 526 
526      UNPROTECT(1);      UNPROTECT(1);
527      return val;      return val;
528  }  }
   
   
   

Legend:
Removed from v.371  
changed lines
  Added in v.956

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