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 692, Mon Apr 18 14:34:33 2005 UTC revision 718, Mon May 9 01:24:35 2005 UTC
# Line 46  Line 46 
46   * @param bi vector of row indices of TRUE elements in B   * @param bi vector of row indices of TRUE elements in B
47   * @param bp column pointers for B   * @param bp column pointers for B
48   * @param beta if non-zero existing TRUE elements in C are retained   * @param beta if non-zero existing TRUE elements in C are retained
49   * @param ciP pointer to the column indices of TRUE elements in C.   * @param ciP SEXP whose INTEGER part is the column indices of TRUE
50   *        This may be reallocated if the number of TRUE elements in C changes   * elements in C (not used if beta == 0).
51   * @param cp column pointers for C   * @param cp column pointers for C
52     *
53     * @return SEXP whose INTEGER part is the column indices of TRUE
54     * elements in the product.  Note that the contents of cp may be modified.
55   */   */
56  void Matrix_lgClgCmm(int tra, int trb, int m, int n, int k,  SEXP Matrix_lgClgCmm(int tra, int trb, int m, int n, int k,
57                       const int ai[], const int ap[],                       const int ai[], const int ap[],
58                       const int bi[], const int bp[],                       const int bi[], const int bp[],
59                       int beta, SEXP *CIP, int cp[])                       int beta, SEXP CIP, int cp[])
60  {  {
61      int annz = ap[tra ? m : k], bnnz = bp[trb ? k : n], cnnz = cp[n], extra = 0;  /*     int annz = ap[tra ? m : k], bnnz = bp[trb ? k : n]; */
62      int *ci = INTEGER(*CIP);      int cnnz = cp[n], extra = 0;
63        int *ci, i, j, prot = 0;    /* prot is the number of PROTECTs to UNPROTECT */
64    
65      if (beta) {      if (beta) {
66            ci = INTEGER(CIP);
67        } else {                    /* blank the C matrix */
68            for (j = 0; j <= n; j++) cp[j] = 0;
69            cnnz = 0;
70            ci = (int *) NULL;
71        }
72          if (tra) {          if (tra) {
73              if (trb) {          /* t(A) %*% t(B) */              if (trb) {          /* t(A) %*% t(B) */
74              } else {            /* t(A) %*% B */              } else {            /* t(A) %*% B */
   
75              }              }
76          } else {          } else {
77              if (trb) {          /* A %*% t(B) */              if (trb) {          /* A %*% t(B) */
78              } else {            /* A %*% B */              } else {            /* A %*% B */
79                for (j = 0; j < n; j++) { /* col index for B and C */
80                    int ii, ii2 = bp[j + 1];
81                    for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
82                        int jj = bi[ii]; /* row index of B; col index of A */
83                        int i, i2 = ap[jj + 1]; /* index into ai */
84                        for (i = ap[jj]; i < i2; i++)
85                            if (check_csc_index(cp, ci, ai[i], j, -1) < 0) extra++;
86                    }
87                }
88                if (extra) {
89                    int ntot = cnnz + extra;
90                    int *Ti = Calloc(ntot, int),
91                        *rwInd = Calloc(m, int), /* indicator of TRUE in column j */
92                        pos = 0;
93    
94                    cp[0] = 0;
95                    for (j = 0; j < n; j++) {
96                        int ii, ii2 = bp[j + 1];
97                        cp[j + 1] = cp[j];
98                        AZERO(rwInd, m); /* initialize column j of C to FALSE */
99                        for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
100                            int jj = bi[ii]; /* row index of B; col index of A */
101                            int i, i2 = ap[jj + 1]; /* index into ai */
102                            for (i = ap[jj]; i < i2; i++) rwInd[ai[i]] = 1;
103                        }
104                        for (i = 0; i < m; i++)
105                            if (rwInd[i]) {cp[j + 1]++; Ti[pos++] = i;}
106                    }
107                    PROTECT(CIP = allocVector(INTSXP, cp[n])); prot++;
108                    Memcpy(INTEGER(CIP), Ti, cp[n]);
109                    Free(Ti); Free(rwInd);
110              }              }
111          }          }
112      }      }
113        UNPROTECT(prot);
114        return CIP;
115  }  }
116    
117  SEXP lgCMatrix_lgCMatrix_mm(SEXP a, SEXP b, SEXP transa, SEXP transb)  SEXP lgCMatrix_lgCMatrix_mm(SEXP a, SEXP b)
118  {  {
119      return R_NilValue;      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lgCMatrix")));
120        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
121            *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
122            *cdims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
123        int k = adims[1], m = adims[0], n = bdims[1];
124        int *cp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
125    
126        if (bdims[0] != k)
127            error(_("Matrices are not conformable for multiplication"));
128        cdims[0] = m; cdims[1] = n;
129        SET_SLOT(ans, Matrix_iSym,
130                 Matrix_lgClgCmm(0, 0, m, n, k,
131                                 INTEGER(GET_SLOT(a, Matrix_iSym)),
132                                 INTEGER(GET_SLOT(a, Matrix_pSym)),
133                                 INTEGER(GET_SLOT(b, Matrix_iSym)),
134                                 INTEGER(GET_SLOT(b, Matrix_pSym)),
135                                 0, (SEXP) NULL, cp));
136        UNPROTECT(1);
137        return ans;
138    }
139    
140    SEXP lgCMatrix_trans(SEXP x)
141    {
142        SEXP xi = GET_SLOT(x, Matrix_iSym);
143        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lgCMatrix")));
144        int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
145            *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
146            nz = length(xi);
147        int *xj = Calloc(nz, int);
148    
149        adims[1] = xdims[0];
150        adims[0] = xdims[1];
151        triplet_to_col(adims[0], adims[1], nz,
152                       expand_cmprPt(xdims[1], INTEGER(GET_SLOT(x, Matrix_pSym)), xj),
153                       INTEGER(GET_SLOT(x, Matrix_iSym)), (double *) NULL,
154                       INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP,  adims[1] + 1)),
155                       INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)),
156                       (double *) NULL);
157        Free(xj);
158        UNPROTECT(1);
159        return ans;
160  }  }

Legend:
Removed from v.692  
changed lines
  Added in v.718

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