SCM

SCM Repository

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

Diff of /pkg/Matrix/src/factorizations.c

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

revision 3114, Tue May 19 16:33:13 2015 UTC revision 3161, Thu Feb 25 09:58:12 2016 UTC
# Line 2  Line 2 
2    
3  SEXP MatrixFactorization_validate(SEXP obj)  SEXP MatrixFactorization_validate(SEXP obj)
4  {  {
5      SEXP Dim = GET_SLOT(obj, Matrix_DimSym), val;      SEXP val;
6      if (isString(val = dim_validate(Dim, "MatrixFactorization")))      if (isString(val = dim_validate(GET_SLOT(obj, Matrix_DimSym),
7                                        "MatrixFactorization")))
8          return(val);          return(val);
9      return ScalarLogical(1);      return ScalarLogical(1);
10  }  }
# Line 60  Line 61 
61          dd = GET_SLOT(x, Matrix_DimSym);          dd = GET_SLOT(x, Matrix_DimSym);
62      int *iperm, *perm, *pivot = INTEGER(GET_SLOT(x, Matrix_permSym)),      int *iperm, *perm, *pivot = INTEGER(GET_SLOT(x, Matrix_permSym)),
63          *dim = INTEGER(dd), m = dim[0], n = dim[1], nn = m, i;          *dim = INTEGER(dd), m = dim[0], n = dim[1], nn = m, i;
64        size_t m_ = (size_t) m; // to prevent integer (multiplication) overflow
65      Rboolean is_sq = (n == m), L_is_tri = TRUE, U_is_tri = TRUE;      Rboolean is_sq = (n == m), L_is_tri = TRUE, U_is_tri = TRUE;
66    
67      // nn :=  min(n,m) ==  length(pivot[])      // nn :=  min(n,m) ==  length(pivot[])
# Line 82  Line 84 
84          SET_SLOT(L, Matrix_xSym, duplicate(lux));          SET_SLOT(L, Matrix_xSym, duplicate(lux));
85          SET_SLOT(L, Matrix_DimSym, duplicate(dd));          SET_SLOT(L, Matrix_DimSym, duplicate(dd));
86      } else { // !is_sq && L_is_tri -- m < n -- "wide" -- L is  m x m      } else { // !is_sq && L_is_tri -- m < n -- "wide" -- L is  m x m
87          double *Lx = REAL(ALLOC_SLOT(L, Matrix_xSym, REALSXP, m * m));          size_t m2 = m_ * m;
88            double *Lx = REAL(ALLOC_SLOT(L, Matrix_xSym, REALSXP, m2));
89          int *dL = INTEGER(ALLOC_SLOT(L, Matrix_DimSym, INTSXP, 2));          int *dL = INTEGER(ALLOC_SLOT(L, Matrix_DimSym, INTSXP, 2));
90          dL[0] = dL[1] = m;          dL[0] = dL[1] = m;
91          // fill lower-diagonal (non-{0,1}) part -- remainder by make_d_matrix*() below:          // fill lower-diagonal (non-{0,1}) part -- remainder by make_d_matrix*() below:
92          Memcpy(Lx, REAL(lux), m * m);          Memcpy(Lx, REAL(lux), m2);
93      }      }
94      if(is_sq || !U_is_tri) {      if(is_sq || !U_is_tri) {
95          SET_SLOT(U, Matrix_xSym, duplicate(lux));          SET_SLOT(U, Matrix_xSym, duplicate(lux));
96          SET_SLOT(U, Matrix_DimSym, duplicate(dd));          SET_SLOT(U, Matrix_DimSym, duplicate(dd));
97      } else { // !is_sq && U_is_tri -- m > n -- "long" -- U is  n x n      } else { // !is_sq && U_is_tri -- m > n -- "long" -- U is  n x n
98          double *Ux = REAL(ALLOC_SLOT(U, Matrix_xSym, REALSXP, n * n)),          double *Ux = REAL(ALLOC_SLOT(U, Matrix_xSym, REALSXP, ((size_t) n) * n)),
99                 *xx = REAL(lux);                 *xx = REAL(lux);
100          int *dU = INTEGER(ALLOC_SLOT(U, Matrix_DimSym, INTSXP, 2));          int *dU = INTEGER(ALLOC_SLOT(U, Matrix_DimSym, INTSXP, 2));
101          dU[0] = dU[1] = n;          dU[0] = dU[1] = n;
102          /* fill upper-diagonal (non-0) part -- remainder by make_d_matrix*() below:          /* fill upper-diagonal (non-0) part -- remainder by make_d_matrix*() below:
103           * this is more complicated than in the L case, as the x / lux part we need           * this is more complicated than in the L case, as the x / lux part we need
104           * is  *not*  continguous:  Memcpy(Ux, REAL(lux), n * n); -- is  WRONG */           * is  *not*  continguous:  Memcpy(Ux, REAL(lux), n * n); -- is  WRONG */
105          for (int j = 0; j < n; j++) {          for (size_t j = 0; j < n; j++) {
106              Memcpy(Ux+j*n, xx+j*m, j+1);              Memcpy(Ux+j*n, xx+j*m, j+1);
107              // for (i = 0; i <= j; i++)              // for (i = 0; i <= j; i++)
108              //   Ux[i + j*n] = xx[i + j*m];              //   Ux[i + j*n] = xx[i + j*m];
# Line 112  Line 115 
115      } else { // L is "unit-diagonal" trapezoidal -- m > n -- "long"      } else { // L is "unit-diagonal" trapezoidal -- m > n -- "long"
116          // fill the upper right part with 0  *and* the diagonal with 1          // fill the upper right part with 0  *and* the diagonal with 1
117          double *Lx = REAL(GET_SLOT(L, Matrix_xSym));          double *Lx = REAL(GET_SLOT(L, Matrix_xSym));
118          int ii;          size_t ii;
119          for (i = 0, ii = 0; i < n; i++, ii+=(m+1)) { // ii = i*(m+1)          for (i = 0, ii = 0; i < n; i++, ii+=(m+1)) { // ii = i*(m+1)
120              Lx[ii] = 1.;              Lx[ii] = 1.;
121              for (int j = i*m; j < ii; j++)              for (size_t j = i*m_; j < ii; j++)
122                  Lx[j] = 0.;                  Lx[j] = 0.;
123          }          }
124      }      }
# Line 128  Line 131 
131          // fill the lower left part with 0          // fill the lower left part with 0
132          double *Ux = REAL(GET_SLOT(U, Matrix_xSym));          double *Ux = REAL(GET_SLOT(U, Matrix_xSym));
133          for (i = 0; i < m; i++)          for (i = 0; i < m; i++)
134              for (int j = i*(m+1)+1; j < (i+1)*m; j++)              for (size_t j = i*(m_+1) +1; j < (i+1)*m_; j++)
135                  Ux[j] = 0.;                  Ux[j] = 0.;
136      }      }
137    

Legend:
Removed from v.3114  
changed lines
  Added in v.3161

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