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

revision 492, Thu Feb 3 14:24:03 2005 UTC revision 683, Thu Mar 31 17:10:16 2005 UTC
# Line 14  Line 14 
14    
15      nrow = dims[0];      nrow = dims[0];
16      if (length(islot) != length(xslot))      if (length(islot) != length(xslot))
17          return ScalarString(mkChar("lengths of slots i and x must match"));          return mkString(_("lengths of slots i and x must match"));
18      if (length(pslot) <= 0)      if (length(pslot) <= 0)
19          return ScalarString(mkChar("slot p must have length > 0"));          return mkString(_("slot p must have length > 0"));
20      if (xp[0] != 0)      if (xp[0] != 0)
21          return ScalarString(mkChar("first element of slot p must be zero"));          return mkString(_("first element of slot p must be zero"));
22      if (length(islot) != xp[ncol])      if (length(islot) != xp[ncol])
23          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"));  
24      for (j = 0; j < ncol; j++) {      for (j = 0; j < ncol; j++) {
25          if (xp[j] > xp[j+1])          if (xp[j] > xp[j+1])
26              return ScalarString(mkChar("slot p must be non-decreasing"));              return mkString(_("slot p must be non-decreasing"));
27      }      }
28      for (j = 0; j < length(islot); j++) {      for (j = 0; j < length(islot); j++) {
29          if (xi[j] < 0 || xi[j] >= nrow)          if (xi[j] < 0 || xi[j] >= nrow)
30              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"));  
31      }      }
32      if (csc_unsorted_columns(ncol, xp, xi)) {      if (csc_unsorted_columns(ncol, xp, xi)) {
33          csc_sort_columns(ncol, xp, xi, REAL(xslot));          csc_sort_columns(ncol, xp, xi, REAL(xslot));
# Line 152  Line 149 
149      triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,      triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,
150                     ansp, itmp, xtmp);                     ansp, itmp, xtmp);
151      nnz = ansp[nrow];      nnz = ansp[nrow];
152      SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));      SET_SLOT(ans, Matrix_uploSym, mkString("L"));
153      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
154      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
155      Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);      Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);
# Line 164  Line 161 
161      return ans;      return ans;
162  }  }
163    
164  SEXP csc_matrix_crossprod(SEXP x, SEXP y)  SEXP csc_matrix_crossprod(SEXP x, SEXP y, SEXP classed)
165  {  {
166      SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;      int cl = asLogical(classed);
167      int j,      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
168          *xp = INTEGER(pslot),      int *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
169          *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),          *ydims = INTEGER(cl ? GET_SLOT(y, Matrix_DimSym) :
170          xncol = length(pslot) - 1,                           getAttrib(y, R_DimSymbol)),
171          xnrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],          *vdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
172          *ydims;      int *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
173      double *xx = REAL(GET_SLOT(x, Matrix_xSym));          *xp = INTEGER(GET_SLOT(x, Matrix_pSym));
174        int j, k = xdims[0], m = xdims[1], n = ydims[1];
175      if (!(isMatrix(y) && isReal(y))) error("y must be a numeric matrix");      double *vx, *xx = REAL(GET_SLOT(x, Matrix_xSym)),
176      ydims = INTEGER(getAttrib(y, R_DimSymbol));          *yx = REAL(cl ? GET_SLOT(y, Matrix_xSym) : y);
177      if (xnrow != ydims[0]) error("x and y must have the same number of rows");  
178      ans = PROTECT(allocMatrix(REALSXP, xncol, ydims[1]));      if (!cl && !(isMatrix(y) && isReal(y)))
179      for (j = 0; j < ydims[1]; j++) {          error(_("y must be a numeric matrix"));
180          int i; double *ypt = REAL(y) + j * ydims[0];      if (ydims[0] != k)
181          for(i = 0; i < xncol; i++) {          error(_("x and y must have the same number of rows"));
182        if (m < 1 || n < 1 || k < 1)
183            error(_("Matrices with zero extents cannot be multiplied"));
184        vdims[0] = m; vdims[1] = n;
185        vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
186        for (j = 0; j < n; j++) {
187            int i; double *ypt = yx + j * k;
188            for(i = 0; i < m; i++) {
189              int ii; double accum = 0.;              int ii; double accum = 0.;
190              for (ii = xp[i]; ii < xp[i+1]; ii++) {              for (ii = xp[i]; ii < xp[i+1]; ii++) {
191                  accum += xx[ii] * ypt[xi[ii]];                  accum += xx[ii] * ypt[xi[ii]];
192              }              }
193              REAL(ans)[i + j * xncol] = accum;              vx[i + j * m] = accum;
194          }          }
195      }      }
196      UNPROTECT(1);      UNPROTECT(1);
197      return ans;      return val;
198  }  }
199    
200  SEXP csc_to_dgTMatrix(SEXP x)  SEXP compressed_to_dgTMatrix(SEXP x, SEXP colP)
201  {  {
202      SEXP      int col = asLogical(colP);
203          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),      SEXP indSym = col ? Matrix_iSym : Matrix_jSym;
204          dimslot = GET_SLOT(x, Matrix_DimSym),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
205          islot = GET_SLOT(x, Matrix_iSym),          indP = GET_SLOT(x, indSym),
206          pslot = GET_SLOT(x, Matrix_pSym);          pP = GET_SLOT(x, Matrix_pSym);
207      int *dims = INTEGER(dimslot), j, jj,      int npt = length(pP) - 1;
         *xp = INTEGER(pslot), *yj;  
208    
209      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));  
210      SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
211      SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, length(islot)));      SET_SLOT(ans, indSym, duplicate(indP));
212      yj = INTEGER(GET_SLOT(ans, Matrix_jSym));      expand_cmprPt(npt, INTEGER(pP),
213      jj = 0;                    INTEGER(ALLOC_SLOT(ans, col ? Matrix_jSym : Matrix_iSym,
214      for (j = 0; j < dims[1]; j++) {                                       INTSXP, length(indP))));
         while (jj < xp[j + 1]) {  
             yj[jj] = j;  
             jj++;  
         }  
     }  
215      UNPROTECT(1);      UNPROTECT(1);
216      return ans;      return ans;
217  }  }
# Line 275  Line 272 
272      double *vx;      double *vx;
273    
274      if (!(isMatrix(A) && isReal(A)))      if (!(isMatrix(A) && isReal(A)))
275          error("A must be a numeric matrix");          error(_("A must be a numeric matrix"));
276      nrow = adims[0]; ncol = adims[1];      nrow = adims[0]; ncol = adims[1];
277      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
278      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
# Line 357  Line 354 
354      int *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
355          nnz = length(islot);          nnz = length(islot);
356    
357      SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));      adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
     adims = INTEGER(GET_SLOT(ans, Matrix_DimSym));  
358      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
359      SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));      csc_compTr(xdims[0], xdims[1], nnz,
360      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));  
     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));  
     csc_components_transpose(xdims[0], xdims[1], nnz,  
                              INTEGER(GET_SLOT(x, Matrix_pSym)),  
                              INTEGER(islot),  
361                               REAL(GET_SLOT(x, Matrix_xSym)),                               REAL(GET_SLOT(x, Matrix_xSym)),
362                               INTEGER(GET_SLOT(ans, Matrix_pSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
363                               INTEGER(GET_SLOT(ans, Matrix_iSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
364                               REAL(GET_SLOT(ans, Matrix_xSym)));                 REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
365      UNPROTECT(1);      UNPROTECT(1);
366      return ans;      return ans;
367  }  }
# Line 386  Line 377 
377      SEXP val;      SEXP val;
378    
379      if (bdim[0] != r)      if (bdim[0] != r)
380          error("Matrices of sizes (%d,%d) and (%d,%d) cannot be multiplied",          error(_("Matrices of sizes (%d,%d) and (%d,%d) cannot be multiplied"),
381                m, r, bdim[0], n);                m, r, bdim[0], n);
382      val = PROTECT(allocMatrix(REALSXP, m, n));      val = PROTECT(allocMatrix(REALSXP, m, n));
383      for (j = 0; j < n; j++) {   /* across columns of b */      for (j = 0; j < n; j++) {   /* across columns of b */
# Line 416  Line 407 
407      SET_SLOT(val, Matrix_DimSym, duplicate(tmp));      SET_SLOT(val, Matrix_DimSym, duplicate(tmp));
408      ncol = INTEGER(tmp)[1];      ncol = INTEGER(tmp)[1];
409      if (!(isInteger(perm) && length(perm) == ncol))      if (!(isInteger(perm) && length(perm) == ncol))
410          error("perm must be an integer vector of length %d",          error(_("perm must be an integer vector of length %d"),
411                ncol);                ncol);
412      prm = INTEGER(perm);      prm = INTEGER(perm);
413      if (!R_ldl_valid_perm(ncol, prm))      if (!R_ldl_valid_perm(ncol, prm))
414          error("perm is not a valid 0-based permutation");          error(_("perm is not a valid 0-based permutation"));
415      iperm = Calloc(ncol, int);      iperm = Calloc(ncol, int);
416      for (j = 0; j < ncol; j++) iperm[prm[j]] = j;      for (j = 0; j < ncol; j++) iperm[prm[j]] = j;
417      tmp = GET_SLOT(x, Matrix_pSym);      tmp = GET_SLOT(x, Matrix_pSym);

Legend:
Removed from v.492  
changed lines
  Added in v.683

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