SCM

SCM Repository

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

Diff of /pkg/src/dgeMatrix.c

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

revision 692, Mon Apr 18 14:34:33 2005 UTC revision 714, Fri May 6 23:52:34 2005 UTC
# Line 82  Line 82 
82      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
83      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
84          *vDims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));          *vDims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
85      int i, k = tr ? Dims[1] : Dims[0], n = tr ? Dims[0] : Dims[1];      int k = tr ? Dims[1] : Dims[0], n = tr ? Dims[0] : Dims[1];
86      double one = 1.0, *xvals, zero = 0.0;      double one = 1.0, zero = 0.0;
87    
88      SET_SLOT(val, Matrix_uploSym, mkString("U"));      SET_SLOT(val, Matrix_uploSym, mkString("U"));
89      vDims[0] = vDims[1] = n;      vDims[0] = vDims[1] = n;
# Line 183  Line 183 
183                       dims,                       dims,
184                       INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, npiv)),                       INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, npiv)),
185                       &info);                       &info);
186      if (info)      if (info < 0)
187          error(_("Lapack routine %s returned error code %d"), "dgetrf", info);          error(_("Lapack routine %s returned error code %d"), "dgetrf", info);
188        else if (info > 0)
189            warning(_("Exact singularity detected during LU decomposition."));
190      UNPROTECT(1);      UNPROTECT(1);
191      return set_factors(x, val, "LU");      return set_factors(x, val, "LU");
192  }  }
# Line 525  Line 527 
527      UNPROTECT(1);      UNPROTECT(1);
528      return val;      return val;
529  }  }
530    
531    SEXP dgeMatrix_colsums(SEXP x, SEXP naRmP, SEXP cols, SEXP mean)
532    {
533        int keepNA = !asLogical(naRmP);
534        int doMean = asLogical(mean);
535        int useCols = asLogical(cols);
536        int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
537        int cnt = 0, i, j, n = dims[0], p = dims[1];
538        SEXP ans = PROTECT(allocVector(REALSXP, (useCols) ? p : n));
539        double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *rx, sum;
540    
541        if (useCols) {
542            cnt = n;
543            for (j = 0; j < p; j++) {
544                rx = xx + n*j;
545                if (keepNA)
546                    for (sum = 0., i = 0; i < n; i++) sum += *rx++;
547                else {
548                    for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
549                        if (!ISNAN(*rx)) {cnt++; sum += *rx;}
550                }
551                if (doMean) {
552                    if (cnt > 0) sum /= cnt; else sum = NA_REAL;
553                }
554                REAL(ans)[j] = sum;
555            }
556        } else {
557            double *rans = REAL(ans), *ra = rans, *rx = xx, *Cnt = NULL, *c;
558            cnt = p;
559            if (!keepNA && doMean) Cnt = Calloc(n, double);
560            for (ra = rans, i = 0; i < n; i++) *ra++ = 0.0;
561            for (j = 0; j < p; j++) {
562                ra = rans;
563                if (keepNA)
564                    for (i = 0; i < n; i++) *ra++ += *rx++;
565                else
566                    for (c = Cnt, i = 0; i < n; i++, ra++, rx++, c++)
567                        if (!ISNAN(*rx)) {
568                            *ra += *rx;
569                            if (doMean) (*c)++;
570                        }
571            }
572            if (doMean) {
573                if (keepNA)
574                    for (ra = rans, i = 0; i < n; i++)
575                        *ra++ /= p;
576                else {
577                    for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
578                        if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
579                    Free(Cnt);
580                }
581            }
582        }
583    
584        UNPROTECT(1);
585        return ans;
586    }

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

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