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 651, Tue Mar 15 01:54:50 2005 UTC revision 652, Wed Mar 16 01:08:36 2005 UTC
# Line 186  Line 186 
186      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
187      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
188      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));  
189      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
190                       dims, INTEGER(GET_SLOT(val, install("pivot"))),                       dims,
191                         INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, npiv)),
192                       &info);                       &info);
193      if (info) error(_("Lapack routine dgetrf returned error code %d"), info);      if (info)
194            error(_("Lapack routine %s returned error code %d"), "dgetrf", info);
195      UNPROTECT(1);      UNPROTECT(1);
196      return set_factors(x, val, "LU");      return set_factors(x, val, "LU");
197  }  }
# Line 200  Line 201 
201      int lg = asLogical(logarithm);      int lg = asLogical(logarithm);
202      SEXP lu = dgeMatrix_LU(x);      SEXP lu = dgeMatrix_LU(x);
203      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
204          *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),          *jpvt = INTEGER(GET_SLOT(lu, Matrix_permSym)),
205          i, n = dims[0], sign = 1;          i, n = dims[0], sign = 1;
206      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;
207    
# Line 231  Line 232 
232      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
233          lu = dgeMatrix_LU(a);          lu = dgeMatrix_LU(a);
234      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
235          *pivot = INTEGER(GET_SLOT(lu, install("pivot")));          *pivot = INTEGER(GET_SLOT(lu, Matrix_permSym));
236      double *x, tmp;      double *x, tmp;
237      int info, lwork = -1;      int info, lwork = -1;
238    
239    
240      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));  
241      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
242      x = REAL(GET_SLOT(val, Matrix_xSym));      x = REAL(GET_SLOT(val, Matrix_xSym));
243      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 250 
250      return val;      return val;
251  }  }
252    
253  SEXP dgeMatrix_dgeMatrix_mm(SEXP a, SEXP b)  SEXP dgeMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
254  {  {
255      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int cl = asLogical(classed);
256          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
257          *cdims,          lu = dgeMatrix_LU(a);
258          m = adims[0], n = bdims[1], k = adims[1];      int *adims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
259            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
260                             getAttrib(b, R_DimSymbol));
261        int info, n = bdims[0], nrhs = bdims[1];
262        int sz = n * nrhs;
263    
264        if (*adims != *bdims || bdims[1] < 1 || *adims < 1 || *adims != adims[1])
265            error(_("Dimensions of system to be solved are inconsistent"));
266        Memcpy(INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)), bdims, 2);
267        F77_CALL(dgetrs)("N", &n, &nrhs, REAL(GET_SLOT(val, Matrix_xSym)), &n,
268                         INTEGER(GET_SLOT(lu, Matrix_permSym)),
269                         Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
270                                REAL(cl ? GET_SLOT(b, Matrix_xSym):b), sz),
271                         &n, &info);
272        UNPROTECT(1);
273        return val;
274    }
275    
276    SEXP dgeMatrix_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
277    {
278        int cl = asLogical(classed), rt = asLogical(right);
279      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
280      char *trans = "N";      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
281            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
282                             getAttrib(b, R_DimSymbol)),
283            *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
284      double one = 1., zero = 0.;      double one = 1., zero = 0.;
285    
286        if (asLogical(right)) {
287            int m = bdims[0], n = adims[1], k = bdims[1];
288            if (adims[0] != k)
289                error(_("Matrices are not conformable for multiplication"));
290            if (m < 1 || n < 1 || k < 1)
291                error(_("Matrices with zero extents cannot be multiplied"));
292            cdims[0] = m; cdims[1] = n;
293            F77_CALL(dgemm) ("N", "N", &m, &n, &k, &one,
294                             REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &m,
295                             REAL(GET_SLOT(a, Matrix_xSym)), &k, &zero,
296                             REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)),
297                             &m);
298        } else {
299            int m = adims[0], n = bdims[1], k = adims[1];
300    
301      if (bdims[0] != k)      if (bdims[0] != k)
302          error(_("Matrices are not conformable for multiplication"));          error(_("Matrices are not conformable for multiplication"));
303      if (m < 1 || n < 1 || k < 1)      if (m < 1 || n < 1 || k < 1)
304          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));  
305      cdims[0] = m; cdims[1] = n;      cdims[0] = m; cdims[1] = n;
306      F77_CALL(dgemm)(trans, trans, adims, bdims+1, bdims, &one,          F77_CALL(dgemm)
307                      REAL(GET_SLOT(a, Matrix_xSym)), adims,              ("N", "N", &m, &n, &k, &one, REAL(GET_SLOT(a, Matrix_xSym)),
308                      REAL(GET_SLOT(b, Matrix_xSym)), bdims,               &m, REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &k, &zero,
309                      &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);               REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)), &m);
310        }
311      UNPROTECT(1);      UNPROTECT(1);
312      return val;      return val;
313  }  }
314    
315    
316  SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)  SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
317  {  {
318      int /* nu = asInteger(nnu),      int /* nu = asInteger(nnu),

Legend:
Removed from v.651  
changed lines
  Added in v.652

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