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 718, Mon May 9 01:24:35 2005 UTC revision 723, Tue May 10 01:42:05 2005 UTC
# Line 58  Line 58 
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  {  {
 /*     int annz = ap[tra ? m : k], bnnz = bp[trb ? k : n]; */  
61      int cnnz = cp[n], extra = 0;      int cnnz = cp[n], extra = 0;
62      int *ci, i, j, prot = 0;    /* prot is the number of PROTECTs to UNPROTECT */      int *ci, i, j, prot = 0;    /* prot is the number of PROTECTs to UNPROTECT */
63    
# Line 69  Line 68 
68          cnnz = 0;          cnnz = 0;
69          ci = (int *) NULL;          ci = (int *) NULL;
70      }      }
71      if (tra) {  
72          if (trb) {              /* t(A) %*% t(B) */      if (tra) {                  /* replace ai and ap by els for transpose */
73          } else {                /* t(A) %*% B */          int nz = ap[m];
74          }          int *Ai = Calloc(nz, int),
75      } else {              *Ti = Calloc(nz, int),
76          if (trb) {              /* A %*% t(B) */              *Ap = Calloc(k + 1, int);
77          } else {                /* A %*% B */  
78            triplet_to_col(k, m, nz, expand_cmprPt(m, ap, Ti), ai,
79                           (double *) NULL, Ap, Ai, (double *) NULL);
80            Free(Ti);
81            ai = Ai; ap = Ap;
82        }
83    
84        if (trb) {                  /* replace bi and bp by els for transpose */
85            int nz = bp[k];
86            int *Bi = Calloc(nz, int),
87                *Ti = Calloc(nz, int),
88                *Bp = Calloc(n + 1, int);
89    
90            triplet_to_col(n, k, nz, expand_cmprPt(k, bp, Ti), bi,
91                           (double *) NULL, Bp, Bi, (double *) NULL);
92            Free(Ti);
93            bi = Bi; bp = Bp;
94        }
95    
96              for (j = 0; j < n; j++) { /* col index for B and C */              for (j = 0; j < n; j++) { /* col index for B and C */
97                  int ii, ii2 = bp[j + 1];                  int ii, ii2 = bp[j + 1];
98                  for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */                  for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
# Line 85  Line 102 
102                          if (check_csc_index(cp, ci, ai[i], j, -1) < 0) extra++;                          if (check_csc_index(cp, ci, ai[i], j, -1) < 0) extra++;
103                  }                  }
104              }              }
105    
106              if (extra) {              if (extra) {
107                  int ntot = cnnz + extra;                  int ntot = cnnz + extra;
108                  int *Ti = Calloc(ntot, int),                  int *Ti = Calloc(ntot, int),
# Line 108  Line 126 
126                  Memcpy(INTEGER(CIP), Ti, cp[n]);                  Memcpy(INTEGER(CIP), Ti, cp[n]);
127                  Free(Ti); Free(rwInd);                  Free(Ti); Free(rwInd);
128              }              }
129          }  
130      }      if (tra) {Free(ai); Free(ap);}
131        if (trb) {Free(bi); Free(bp);}
132      UNPROTECT(prot);      UNPROTECT(prot);
133      return CIP;      return CIP;
134  }  }
# Line 158  Line 177 
177      UNPROTECT(1);      UNPROTECT(1);
178      return ans;      return ans;
179  }  }
180    
181    /**
182     * Replace C by AA' + beta*C or A'A + beta*C
183     *
184     * @param up Indicator of upper/lower triangle in the symmetric sparse matrix
185     * @param tra Transpose, in the sense of dsyrk.  That is, tra TRUE indicates A'A
186     * @param n size of the product matrix
187     * @param k number of columns in A if tra is FALSE, otherwise the number of rows
188     * @param ai row indices for A
189     * @param ap column pointers for A
190     * @param beta TRUE if existing elements in C are to be preserved
191     * @param CIP SEXP whose INTEGER part is the row indices of C (not used if beta is FALSE)
192     * @param cp column pointers for C
193     *
194     * @return SEXP whose INTEGER part is the updated row indices of C
195     */
196    SEXP Matrix_lgCsyrk(int up, int tra, int n, int k, const int ai[], const int ap[],
197                        int beta, SEXP CIP, int cp[])
198    {
199        int extra = 0, i, ii, j, prot = 0;
200        int *ci, cnnz = cp[n];
201    
202        if (beta) {
203            ci = INTEGER(CIP);
204        } else {                    /* blank the C matrix */
205            for (j = 0; j <= n; j++) cp[j] = 0;
206            cnnz = 0;
207            ci = (int *) NULL;
208        }
209    
210        if (tra) {                  /* replace ai and ap by els for transpose */
211            int nz = ap[n];
212            int *Ai = Calloc(nz, int),
213                *Ti = Calloc(nz, int),
214                *Ap = Calloc(k + 1, int);
215    
216            triplet_to_col(k, n, nz, expand_cmprPt(k, ap, Ti), ai,
217                           (double *) NULL, Ap, Ai, (double *) NULL);
218            Free(Ti);
219            ai = Ai; ap = Ap;
220        }
221        for (j = 0; j < k; j++) {
222            int i2 = ap[j + 1];
223            for (i = ap[j]; i < i2; i++) {
224                int r1 = ai[i];
225                for (ii = i; ii < i2; ii++) {
226                    int r2 = ai[ii];
227                    if (check_csc_index(cp, ci, up?r1:r2, up?r2:r1, -1) < 0)
228                        extra++;
229                }
230            }
231        }
232    
233        if (extra) {
234            int ntot = cnnz + extra;
235            int *Ti = Memcpy(Calloc(ntot, int), ci, cnnz),
236                *Tj = expand_cmprPt(n, cp, Calloc(ntot, int)),
237                *Ci = Calloc(ntot, int),
238                pos = cnnz;
239    
240            for (j = 0; j < k; j++) {
241                int i2 = ap[j + 1];
242                for (i = ap[j]; i < i2; i++) {
243                    int r1 = ai[i];
244                    for (ii = i; ii < i2; ii++) {
245                        int r2 = ai[ii];
246                        int row = up ? r1 : r2, col = up ? r2 : r1;
247                        if (r2 < r1) error("[j,i,ii,r1,r2] = [%d,%d,%d,%d,%d]",
248                                           j,i,ii,r1,r2);
249                        if (check_csc_index(cp, ci, row, col, -1) < 0) {
250                            Ti[pos] = row;
251                            Tj[pos] = col;
252                            pos++;
253                        }
254                    }
255                }
256            }
257    
258            triplet_to_col(n, n, pos, Ti, Tj, (double *) NULL,
259                           cp, Ci, (double *) NULL);
260            PROTECT(CIP = allocVector(INTSXP, cp[n])); prot++;
261            Memcpy(INTEGER(CIP), Ci, cp[n]);
262            Free(Ti); Free(Tj); Free(Ci);
263        }
264    
265        if (tra) {Free(ai); Free(ap);}
266        UNPROTECT(prot);
267        return CIP;
268    }
269    
270    /**
271     * Create the cross-product or transpose cross-product of a logical
272     * sparse matrix in column-oriented compressed storage mode.
273     *
274     * @param x Pointer to a lgCMatrix
275     * @param trans logical indicator of transpose, in the sense of dsyrk.
276     * That is, trans TRUE is used for crossprod.
277     *
278     * @return An lsCMatrix of the form if(trans) X'X else XX'
279     */
280    SEXP lgCMatrix_crossprod(SEXP x, SEXP trans)
281    {
282        int tra = asLogical(trans);
283        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));
284        int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
285            *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
286        int k = xdims[tra ? 0 : 1], n = xdims[tra ? 1 : 0];
287    
288        adims[0] = adims[1] = n;
289        SET_SLOT(ans, Matrix_uploSym, mkString("U"));
290        SET_SLOT(ans, Matrix_iSym,
291                 Matrix_lgCsyrk(1, tra, n, k,
292                                INTEGER(GET_SLOT(x, Matrix_iSym)),
293                                INTEGER(GET_SLOT(x, Matrix_pSym)),
294                                0, R_NilValue,
295                                INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1))));
296        UNPROTECT(1);
297        return ans;
298    }

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

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