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 604, Fri Mar 4 17:27:51 2005 UTC revision 692, Mon Apr 18 14:34:33 2005 UTC
# Line 76  Line 76 
76    return ScalarReal(set_rcond(obj, CHAR(asChar(type))));    return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
77  }  }
78    
79  SEXP dgeMatrix_crossprod(SEXP x)  SEXP dgeMatrix_crossprod(SEXP x, SEXP trans)
80  {  {
81        int tr = asLogical(trans);
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;          *vDims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
85      int i, n = Dims[1];      int i, k = tr ? Dims[1] : Dims[0], n = tr ? Dims[0] : Dims[1];
     int nsqr = n * n;  
86      double one = 1.0, *xvals, zero = 0.0;      double one = 1.0, *xvals, zero = 0.0;
87    
     SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));  
     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));  
88      SET_SLOT(val, Matrix_uploSym, mkString("U"));      SET_SLOT(val, Matrix_uploSym, mkString("U"));
     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));  
     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));  
89      vDims[0] = vDims[1] = n;      vDims[0] = vDims[1] = n;
90      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nsqr));      F77_CALL(dsyrk)("U", tr ? "N" : "T", &n, &k,
91      xvals = REAL(GET_SLOT(val, Matrix_xSym));                      &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims, &zero,
92      for (i = 0; i < nsqr; i++) xvals[i] = 0.; /* keep valgrind happy */                      REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n)), &n);
     F77_CALL(dsyrk)("U", "T", vDims, Dims,  
                     &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims,  
                     &zero, xvals, vDims);  
93      UNPROTECT(1);      UNPROTECT(1);
94      return val;      return val;
95  }  }
# Line 186  Line 179 
179      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
180      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
181      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
     SET_SLOT(val, install("pivot"), allocVector(INTSXP, npiv));  
182      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
183                       dims, INTEGER(GET_SLOT(val, install("pivot"))),                       dims,
184                         INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, npiv)),
185                       &info);                       &info);
186      if (info) error(_("Lapack routine dgetrf returned error code %d"), info);      if (info)
187            error(_("Lapack routine %s returned error code %d"), "dgetrf", info);
188      UNPROTECT(1);      UNPROTECT(1);
189      return set_factors(x, val, "LU");      return set_factors(x, val, "LU");
190  }  }
# Line 200  Line 194 
194      int lg = asLogical(logarithm);      int lg = asLogical(logarithm);
195      SEXP lu = dgeMatrix_LU(x);      SEXP lu = dgeMatrix_LU(x);
196      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
197          *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),          *jpvt = INTEGER(GET_SLOT(lu, Matrix_permSym)),
198          i, n = dims[0], sign = 1;          i, n = dims[0], sign = 1;
199      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;
200    
# Line 231  Line 225 
225      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
226          lu = dgeMatrix_LU(a);          lu = dgeMatrix_LU(a);
227      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
228          *pivot = INTEGER(GET_SLOT(lu, install("pivot")));          *pivot = INTEGER(GET_SLOT(lu, Matrix_permSym));
229      double *x, tmp;      double *x, tmp;
230      int info, lwork = -1;      int info, lwork = -1;
231    
232    
233      if (dims[0] != dims[1]) error(_("Solve requires a square matrix"));      if (dims[0] != dims[1]) error(_("Solve requires a square matrix"));
     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));  
     SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));  
234      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
235      x = REAL(GET_SLOT(val, Matrix_xSym));      x = REAL(GET_SLOT(val, Matrix_xSym));
236      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));
# Line 251  Line 243 
243      return val;      return val;
244  }  }
245    
246  SEXP dgeMatrix_dgeMatrix_mm(SEXP a, SEXP b)  SEXP dgeMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
247  {  {
248      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int cl = asLogical(classed);
249          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
250          *cdims,          lu = dgeMatrix_LU(a);
251          m = adims[0], n = bdims[1], k = adims[1];      int *adims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
252            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
253                             getAttrib(b, R_DimSymbol));
254        int info, n = bdims[0], nrhs = bdims[1];
255        int sz = n * nrhs;
256    
257        if (*adims != *bdims || bdims[1] < 1 || *adims < 1 || *adims != adims[1])
258            error(_("Dimensions of system to be solved are inconsistent"));
259        Memcpy(INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)), bdims, 2);
260        F77_CALL(dgetrs)("N", &n, &nrhs, REAL(GET_SLOT(lu, Matrix_xSym)), &n,
261                         INTEGER(GET_SLOT(lu, Matrix_permSym)),
262                         Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
263                                REAL(cl ? GET_SLOT(b, Matrix_xSym):b), sz),
264                         &n, &info);
265        UNPROTECT(1);
266        return val;
267    }
268    
269    SEXP dgeMatrix_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
270    {
271        int cl = asLogical(classed);
272      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
273      char *trans = "N";      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
274            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
275                             getAttrib(b, R_DimSymbol)),
276            *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
277      double one = 1., zero = 0.;      double one = 1., zero = 0.;
278    
279        if (asLogical(right)) {
280            int m = bdims[0], n = adims[1], k = bdims[1];
281            if (adims[0] != k)
282                error(_("Matrices are not conformable for multiplication"));
283            if (m < 1 || n < 1 || k < 1)
284                error(_("Matrices with zero extents cannot be multiplied"));
285            cdims[0] = m; cdims[1] = n;
286            F77_CALL(dgemm) ("N", "N", &m, &n, &k, &one,
287                             REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &m,
288                             REAL(GET_SLOT(a, Matrix_xSym)), &k, &zero,
289                             REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)),
290                             &m);
291        } else {
292            int m = adims[0], n = bdims[1], k = adims[1];
293    
294      if (bdims[0] != k)      if (bdims[0] != k)
295          error(_("Matrices are not conformable for multiplication"));          error(_("Matrices are not conformable for multiplication"));
296      if (m < 1 || n < 1 || k < 1)      if (m < 1 || n < 1 || k < 1)
297          error(_("Matrices with zero extents cannot be multiplied"));          error(_("Matrices with zero extents cannot be multiplied"));
     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));  
     SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));  
     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));  
     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));  
     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));  
298      cdims[0] = m; cdims[1] = n;      cdims[0] = m; cdims[1] = n;
299      F77_CALL(dgemm)(trans, trans, adims, bdims+1, bdims, &one,          F77_CALL(dgemm)
300                      REAL(GET_SLOT(a, Matrix_xSym)), adims,              ("N", "N", &m, &n, &k, &one, REAL(GET_SLOT(a, Matrix_xSym)),
301                      REAL(GET_SLOT(b, Matrix_xSym)), bdims,               &m, REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &k, &zero,
302                      &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);               REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)), &m);
303        }
304      UNPROTECT(1);      UNPROTECT(1);
305      return val;      return val;
306  }  }
307    
308    
309  SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)  SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
310  {  {
311      int /* nu = asInteger(nnu),      int /* nu = asInteger(nnu),

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

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