SCM

SCM Repository

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

Diff of /pkg/src/lgCMatrix.c

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

revision 726, Thu May 12 14:59:04 2005 UTC revision 732, Tue May 17 17:18:14 2005 UTC
# Line 72  Line 72 
72      if (tra) {                  /* replace ai and ap by els for transpose */      if (tra) {                  /* replace ai and ap by els for transpose */
73          int nz = ap[m];          int nz = ap[m];
74          int *Ai = Calloc(nz, int),          int *Ai = Calloc(nz, int),
75              *Ti = Calloc(nz, int),              *aj = expand_cmprPt(m, ap, Calloc(nz, int)),
76              *Ap = Calloc(k + 1, int);              *Ap = Calloc(k + 1, int);
77    
78          triplet_to_col(k, m, nz, expand_cmprPt(m, ap, Ti), ai,          triplet_to_col(m, k, nz, aj, ai, (double *) NULL,
79                         (double *) NULL, Ap, Ai, (double *) NULL);                         Ap, Ai, (double *) NULL);
80          Free(Ti);          Free(aj);
81          ai = Ai; ap = Ap;          ai = Ai; ap = Ap;
82      }      }
83    
84      if (trb) {                  /* replace bi and bp by els for transpose */      if (trb) {                  /* replace bi and bp by els for transpose */
85          int nz = bp[k];          int nz = bp[k];
86          int *Bi = Calloc(nz, int),          int *Bi = Calloc(nz, int),
87              *Ti = Calloc(nz, int),              *bj = expand_cmprPt(k, bp, Calloc(nz, int)),
88              *Bp = Calloc(n + 1, int);              *Bp = Calloc(n + 1, int);
89    
90          triplet_to_col(n, k, nz, expand_cmprPt(k, bp, Ti), bi,          triplet_to_col(k, n, nz, bj, bi, (double *) NULL,
91                         (double *) NULL, Bp, Bi, (double *) NULL);                         Bp, Bi, (double *) NULL);
92          Free(Ti);          Free(bj);
93          bi = Bi; bp = Bp;          bi = Bi; bp = Bp;
94      }      }
95    
# Line 105  Line 105 
105    
106      if (extra) {      if (extra) {
107          int ntot = cnnz + extra;          int ntot = cnnz + extra;
108          int *Ti = Calloc(ntot, int),          int *Cp = Calloc(n + 1, int),
109                *Ti = Calloc(ntot, int),
110              *rwInd = Calloc(m, int), /* indicator of TRUE in column j */              *rwInd = Calloc(m, int), /* indicator of TRUE in column j */
111              pos = 0;              pos = 0;
112    
113          cp[0] = 0;          Cp[0] = 0;
114          for (j = 0; j < n; j++) {          for (j = 0; j < n; j++) {
115              int ii, ii2 = bp[j + 1];              int ii, ii2 = bp[j + 1];
116              cp[j + 1] = cp[j];  
117              AZERO(rwInd, m); /* initialize column j of C to FALSE */              AZERO(rwInd, m);    /* initialize column j of C */
118                for (i = cp[j]; i < cp[j+1]; i++) rwInd[ci[i]] = 1;
119    
120                Cp[j + 1] = Cp[j];
121              for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */              for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
122                  int jj = bi[ii]; /* row index of B; col index of A */                  int jj = bi[ii]; /* row index of B; col index of A */
123                  int i, i2 = ap[jj + 1]; /* index into ai */                  int i, i2 = ap[jj + 1]; /* index into ai */
124                  for (i = ap[jj]; i < i2; i++) rwInd[ai[i]] = 1;                  for (i = ap[jj]; i < i2; i++) rwInd[ai[i]] = 1;
125              }              }
126              for (i = 0; i < m; i++)              for (i = 0; i < m; i++)
127                  if (rwInd[i]) {cp[j + 1]++; Ti[pos++] = i;}                  if (rwInd[i]) {Cp[j + 1]++; Ti[pos++] = i;}
128          }          }
129          PROTECT(CIP = allocVector(INTSXP, cp[n])); prot++;          PROTECT(CIP = allocVector(INTSXP, Cp[n])); prot++;
130          Memcpy(INTEGER(CIP), Ti, cp[n]);          Memcpy(INTEGER(CIP), Ti, Cp[n]);
131          Free(Ti); Free(rwInd);          Memcpy(cp, Cp, n + 1);
132            Free(Cp); Free(Ti); Free(rwInd);
133      }      }
134    
135      if (tra) {Free(ai); Free(ap);}      if (tra) {Free(ai); Free(ap);}
# Line 213  Line 218 
218      if (tra) {                  /* replace ai and ap by els for transpose */      if (tra) {                  /* replace ai and ap by els for transpose */
219          int nz = ap[n];          int nz = ap[n];
220          int *Ai = Calloc(nz, int),          int *Ai = Calloc(nz, int),
221              *Ti = Calloc(nz, int),              *aj = expand_cmprPt(n, ap, Calloc(nz, int)),
222              *Ap = Calloc(k + 1, int);              *Ap = Calloc(k + 1, int);
223    
224          triplet_to_col(k, n, nz, expand_cmprPt(k, ap, Ti), ai,          triplet_to_col(n, k, nz, aj, ai, (double *) NULL,
225                         (double *) NULL, Ap, Ai, (double *) NULL);                         Ap, Ai, (double *) NULL);
226          Free(Ti);          Free(aj);
227          ai = Ai; ap = Ap;          ai = Ai; ap = Ap;
228      }      }
229    
230      for (j = 0; j < k; j++) {      for (j = 0; j < k; j++) {
231          int i2 = ap[j + 1];          int i2 = ap[j + 1];
232          for (i = ap[j]; i < i2; i++) {          for (i = ap[j]; i < i2; i++) {
233              int r1 = ai[i];              int r1 = ai[i];
234                if (r1 < 0 || r1 >= n)
235                    error(_("row %d not in row range [0,%d]"), r1, n - 1);
236              for (ii = i; ii < i2; ii++) {              for (ii = i; ii < i2; ii++) {
237                  int r2 = ai[ii];                  int r2 = ai[ii];
238                    if (r2 < 0 || r2 >= n)
239                        error(_("row %d not in row range [0,%d]"), r2, n - 1);
240                  if (check_csc_index(cp, ci, up?r1:r2, up?r2:r1, -1) < 0)                  if (check_csc_index(cp, ci, up?r1:r2, up?r2:r1, -1) < 0)
241                      extra++;                      extra++;
242              }              }
# Line 276  Line 286 
286   *   *
287   * @param x Pointer to a lgCMatrix   * @param x Pointer to a lgCMatrix
288   * @param trans logical indicator of transpose, in the sense of dsyrk.   * @param trans logical indicator of transpose, in the sense of dsyrk.
289   * That is, trans TRUE is used for crossprod.   * That is, trans == TRUE is used for crossprod.
290     * @param C
291   *   *
292   * @return An lsCMatrix of the form if(trans) X'X else XX'   * @return An lsCMatrix of the form if(trans) X'X else XX'
293   */   */
294  SEXP lgCMatrix_crossprod(SEXP x, SEXP trans)  SEXP lgCMatrix_crossprod(SEXP x, SEXP trans, SEXP C)
295  {  {
296      int tra = asLogical(trans);      int tra = asLogical(trans);
297      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));      int *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
     int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),  
         *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));  
298      int k = xdims[tra ? 0 : 1], n = xdims[tra ? 1 : 0];      int k = xdims[tra ? 0 : 1], n = xdims[tra ? 1 : 0];
299    
300        if (C == R_NilValue) {
301            SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));
302    
303            adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
304      adims[0] = adims[1] = n;      adims[0] = adims[1] = n;
305      SET_SLOT(ans, Matrix_uploSym, mkString("U"));      SET_SLOT(ans, Matrix_uploSym, mkString("U"));
306      SET_SLOT(ans, Matrix_iSym,      SET_SLOT(ans, Matrix_iSym,
# Line 299  Line 312 
312      UNPROTECT(1);      UNPROTECT(1);
313      return ans;      return ans;
314  }  }
315        adims = INTEGER(GET_SLOT(C, Matrix_DimSym));
316        if (adims[0] != n || adims[1] != n)
317            error(_("Dimensions of x and y are not compatible for crossprod"));
318        SET_SLOT(C, Matrix_iSym,
319                 Matrix_lgCsyrk(CHAR(asChar(GET_SLOT(C, Matrix_uploSym)))[0] == 'U',
320                                tra, n, k,
321                                INTEGER(GET_SLOT(x, Matrix_iSym)),
322                                INTEGER(GET_SLOT(x, Matrix_pSym)),
323                                1, GET_SLOT(C, Matrix_iSym),
324                                INTEGER(GET_SLOT(C, Matrix_pSym))));
325        return C;
326    }
327    
328    /**
329     * Special-purpose function that returns a permutation of the columns
330     * of a lgTMatrix for which nrow(x) > ncol(x).  The ordering puts
331     * columns with fewer entries on the left.  Once a column has been
332     * moved to the left the rows in where that column is TRUE are removed
333     * from the counts.
334     *
335     * @param x Pointer to an lgTMatrix object
336     *
337     * @return 0-based permutation vector for the columns of x
338     */
339    SEXP lgCMatrix_picky_column(SEXP x)
340    {
341        int *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
342        int *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
343            *xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
344            m = xdims[0], n = xdims[1];
345        SEXP ans = PROTECT(allocVector(INTSXP, n));
346        int *actr = Calloc(m, int),
347            *actc = Calloc(n, int),
348            cj, i, j, mincount, minloc, pos;
349    
350        for (i = 0; i < m; i++) actr[i] = 1;
351        mincount = m + 1;
352        for (j = 0; j < n; j++) {
353            cj = xp[j + 1] - xp[j];
354            actc[j] = 1;
355            if (cj < mincount) {
356                mincount = cj;
357                minloc = j;
358            }
359        }
360    
361        pos = 0;
362        while (pos < n) {
363            INTEGER(ans)[pos++] = minloc;
364            actc[minloc] = 0;
365            for (i = xp[minloc]; i < xp[minloc + 1]; i++) actr[xi[i]] = 0;
366            mincount = m + 1;
367            for (j = 0; j < n; j++) {
368                if (actc[j]) {
369                    cj = 0;
370                    for (i = xp[j]; i < xp[j + 1]; i++) {
371                        if (actr[xi[i]]) cj++;
372                        if (cj < mincount) {
373                            mincount = cj;
374                            minloc = j;
375                        }
376                    }
377                }
378            }
379        }
380    
381        Free(actr); Free(actc);
382        UNPROTECT(1);
383        return ans;
384    }

Legend:
Removed from v.726  
changed lines
  Added in v.732

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