SCM

SCM Repository

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

Diff of /pkg/src/cscMatrix.c

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

revision 310, Wed Oct 27 21:40:08 2004 UTC revision 332, Fri Nov 12 21:04:36 2004 UTC
# Line 49  Line 49 
49      int j, *iVal, ncol = length(pslot) - 1, maxnz, nnz = 0, *pVal;      int j, *iVal, ncol = length(pslot) - 1, maxnz, nnz = 0, *pVal;
50      double *xVal;      double *xVal;
51    
52        SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));
53      maxnz = (ncol * (ncol + 1))/2;      maxnz = (ncol * (ncol + 1))/2;
54      iVal = Calloc(maxnz, int); xVal = Calloc(maxnz, double);      iVal = Calloc(maxnz, int); xVal = Calloc(maxnz, double);
55      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, ncol + 1));      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, ncol + 1));
# Line 100  Line 101 
101      return cscMatrix_set_Dim(ans, ncol);      return cscMatrix_set_Dim(ans, ncol);
102  }  }
103    
104    SEXP csc_tcrossprod(SEXP x)
105    {
106        SEXP pslot = GET_SLOT(x, Matrix_pSym),
107            ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))), tmp;
108        int *xp = INTEGER(pslot),
109            *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
110            *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
111        double *xx = REAL(GET_SLOT(x, Matrix_xSym));
112    
113        int j, ntrip, *iVal, nrow = dims[0], ncol = dims[1], *jVal, nnz, pos;
114        int *itmp, *ansp;
115        double *xVal, *xtmp;
116    
117        SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));
118        ntrip = nrow;               /* number of triplets */
119        for (j = 0; j < ncol; j++) {
120            int nzj = xp[j+1] - xp[j];
121            ntrip += (nzj * (nzj - 1))/2;
122        }
123        iVal = Calloc(ntrip, int); jVal = Calloc(ntrip, int);
124        xVal = Calloc(ntrip, double);
125        for (j = 0; j < nrow; j++) {
126            iVal[j] = jVal[j] = j;
127            xVal[j] = 0.;
128        }
129        pos = nrow;
130        for (j = 0; j < ncol; j++) {
131            int k, kk, k2 = xp[j+1];
132            for (k = xp[j]; k < k2; k++) {
133                int r1 = xi[k];
134                double x1 = xx[k];
135                xVal[r1] += x1 * x1;
136                for (kk = k + 1; kk < k2; kk++) {
137                    int r2 = xi[kk];
138                    double x2 = xx[kk];
139                    jVal[pos] = r1;
140                    iVal[pos] = r2;
141                    xVal[pos] = x1 * x2;
142                    pos++;
143                }
144            }
145        }
146        SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, nrow + 1));
147        ansp = INTEGER(GET_SLOT(ans, Matrix_pSym));
148        itmp = Calloc(ntrip, int); xtmp = Calloc(ntrip, double);
149        triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,
150                       ansp, itmp, xtmp);
151        nnz = ansp[nrow];
152        SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
153        SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
154        SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
155        Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);
156        Memcpy(REAL(GET_SLOT(ans, Matrix_xSym)), xtmp, nnz);
157        SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
158        dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
159        dims[0] = dims[1] = nrow;
160        Free(itmp); Free(xtmp); Free(iVal); Free(jVal); Free(xVal);
161        UNPROTECT(1);
162        return ans;
163    }
164    
165  SEXP csc_matrix_crossprod(SEXP x, SEXP y)  SEXP csc_matrix_crossprod(SEXP x, SEXP y)
166  {  {
167      SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;      SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;
# Line 188  Line 250 
250    
251      SET_SLOT(ans, Matrix_DimSym, duplicate(Dimslot));      SET_SLOT(ans, Matrix_DimSym, duplicate(Dimslot));
252      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nrow*ncol));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nrow*ncol));
253        SET_SLOT(ans, Matrix_rcondSym, allocVector(REALSXP, 0));
254        SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));
255      ax = REAL(GET_SLOT(ans, Matrix_xSym));      ax = REAL(GET_SLOT(ans, Matrix_xSym));
256      for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;      for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;
257      for (j = 0; j < ncol; j++) {      for (j = 0; j < ncol; j++) {
# Line 200  Line 264 
264      return ans;      return ans;
265  }  }
266    
 /* FIXME: This function does not appear to be used. */  
 static                          /* added to check if it was used */  
 SEXP csc_to_imagemat(SEXP x)  
 {  
     SEXP ans, pslot = GET_SLOT(x, Matrix_pSym);  
     int j, ncol = length(pslot) - 1,  
         nrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],  
         *xp = INTEGER(pslot),  
         *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),  
         *ax;  
                                 /* ans is the transpose of the indicator  
                                    of non-zero */  
     ax = INTEGER(ans = PROTECT(allocMatrix(INTSXP, ncol, nrow)));  
     for (j = 0; j < (ncol * nrow); j++) ax[j] = 0;  
     for (j = 0; j < ncol; j++) {  
         int ind;  
   
         for (ind = xp[j]; ind < xp[j+1]; ind++) {  
                                 /* reverse rows of transpose */  
             ax[j + ncol*(nrow - 1 - xi[ind])] = 1;  
         }  
     }  
     UNPROTECT(1);  
     return ans;  
 }  
   
267  SEXP matrix_to_csc(SEXP A)  SEXP matrix_to_csc(SEXP A)
268  {  {
269      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix")));
# Line 237  Line 275 
275      if (!(isMatrix(A) && isReal(A)))      if (!(isMatrix(A) && isReal(A)))
276          error("A must be a numeric matrix");          error("A must be a numeric matrix");
277      nrow = adims[0]; ncol = adims[1];      nrow = adims[0]; ncol = adims[1];
278        SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
279      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
280      vp = INTEGER(GET_SLOT(val, Matrix_pSym));      vp = INTEGER(GET_SLOT(val, Matrix_pSym));
281      maxnz = nrow * ncol;      maxnz = nrow * ncol;
# Line 318  Line 357 
357          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
358    
359      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
360        SET_SLOT(ans, Matrix_factorization, allocVector(VECSXP, 0));
361      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));
362      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
363      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
# Line 362  Line 402 
402      return val;      return val;
403  }  }
404    
405    SEXP csc_col_permute(SEXP x, SEXP perm)
406    {
407        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix"))), tmp;
408        int *iperm, *prm, *vi, *vp, *xi, *xp, j, k, ncol, pos;
409        double *vx, *xx;
410    
411        SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
412        tmp = GET_SLOT(x, Matrix_DimSym);
413        SET_SLOT(val, Matrix_DimSym, duplicate(tmp));
414        ncol = INTEGER(tmp)[1];
415        if (!(isInteger(perm) && length(perm) == ncol))
416            error("perm must be an integer vector of length %d",
417                  ncol);
418        prm = INTEGER(perm);
419        iperm = Calloc(ncol, int);
420        if (!ldl_valid_perm(ncol, prm, iperm))
421            error("perm is not a valid 0-based permutation");
422        for (j = 0; j < ncol; j++) iperm[prm[j]] = j;
423        tmp = GET_SLOT(x, Matrix_pSym);
424        xp = INTEGER(tmp);
425        SET_SLOT(val, Matrix_pSym, duplicate(tmp));
426        vp = INTEGER(GET_SLOT(val, Matrix_pSym));
427        tmp = GET_SLOT(x, Matrix_iSym);
428        xi = INTEGER(tmp);
429        SET_SLOT(val, Matrix_iSym, duplicate(tmp));
430        vi = INTEGER(GET_SLOT(val, Matrix_iSym));
431        tmp = GET_SLOT(x, Matrix_xSym);
432        xx = REAL(tmp);
433        SET_SLOT(val, Matrix_xSym, duplicate(tmp));
434        vx = REAL(GET_SLOT(val, Matrix_xSym));
435    
436        pos = vp[0] = 0;
437        for (j = 0; j < ncol; j++) {
438            int jj = iperm[j];
439            int j1 = xp[jj], j2 = xp[jj+1];
440            vp[j + 1] = vp[j] + (j2 - j1);
441            for (k = j1; k < j2; k++) {
442                vi[pos] = xi[k];
443                vx[pos] = xx[k];
444                pos++;
445            }
446        }
447        Free(iperm);
448        UNPROTECT(1);
449        return val;
450    }
451    
452    
453    

Legend:
Removed from v.310  
changed lines
  Added in v.332

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