SCM

SCM Repository

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

Diff of /pkg/src/dspMatrix.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 84  Line 84 
84      return val;      return val;
85  }  }
86    
87  SEXP dspMatrix_dgeMatrix_solve(SEXP a, SEXP b)  SEXP dspMatrix_matrix_solve(SEXP a, SEXP b, SEXP classedP)
88  {  {
89        int classed = asLogical(classedP);
90      SEXP trf = dspMatrix_trf(a),      SEXP trf = dspMatrix_trf(a),
91          val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));          val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
92      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
93          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),          *bdims = (classed ? INTEGER(GET_SLOT(b, Matrix_DimSym)) :
94          info;                    INTEGER(getAttrib(b, R_DimSymbol)));
95        int n = bdims[0], nrhs = bdims[1], info;
96        int sz = n * nrhs;
97        double *bx = (classed ? REAL(GET_SLOT(b, Matrix_xSym)) : REAL(b));
98    
99      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)      if (!classed && !(isReal(b) && isMatrix(b)))
         error(_("Dimensions of system to be solved are inconsistent"));  
     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(b, Matrix_DimSym)));  
     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(b, Matrix_xSym)));  
     F77_CALL(dsptrs)(CHAR(asChar(GET_SLOT(trf, Matrix_uploSym))),  
                      adims, bdims + 1,  
                      REAL(GET_SLOT(trf, Matrix_xSym)),  
                      INTEGER(GET_SLOT(trf, Matrix_permSym)),  
                      REAL(GET_SLOT(val, Matrix_xSym)),  
                      bdims, &info);  
     UNPROTECT(1);  
     return val;  
 }  
   
 SEXP dspMatrix_matrix_solve(SEXP a, SEXP b)  
 {  
     SEXP trf = dspMatrix_trf(a),  
         val = PROTECT(duplicate(b));  
     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),  
         *bdims = INTEGER(getAttrib(b, R_DimSymbol)),  
         info;  
   
     if (!(isReal(b) && isMatrix(b)))  
100          error(_("Argument b must be a numeric matrix"));          error(_("Argument b must be a numeric matrix"));
101      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
102          error(_("Dimensions of system to be solved are inconsistent"));          error(_("Dimensions of system to be solved are inconsistent"));
103        Memcpy(INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)), bdims, 2);
104      F77_CALL(dsptrs)(CHAR(asChar(GET_SLOT(trf, Matrix_uploSym))),      F77_CALL(dsptrs)(CHAR(asChar(GET_SLOT(trf, Matrix_uploSym))),
105                       adims, bdims + 1,                       &n, &nrhs, REAL(GET_SLOT(trf, Matrix_xSym)),
                      REAL(GET_SLOT(trf, Matrix_xSym)),  
106                       INTEGER(GET_SLOT(trf, Matrix_permSym)),                       INTEGER(GET_SLOT(trf, Matrix_permSym)),
107                       REAL(val), bdims, &info);                       Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
108                                bx, sz), &n, &info);
109      UNPROTECT(1);      UNPROTECT(1);
110      return val;      return val;
111  }  }
# Line 147  Line 130 
130      return val;      return val;
131  }  }
132    
133  SEXP dspMatrix_dgeMatrix_mm(SEXP a, SEXP b)  SEXP dspMatrix_matrix_mm(SEXP a, SEXP b, SEXP classedP)
134  {  {
135        int classed = asLogical(classedP);
136        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
137            bdimP = (classed ? GET_SLOT(b, Matrix_DimSym) :
138                     getAttrib(b, R_DimSymbol));
139      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
140          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym));          *bdims = INTEGER(bdimP);
141      int i, ione = 1, n = adims[0], nrhs = bdims[1];      int i, ione = 1, n = bdims[0], nrhs = bdims[1], info;
142      SEXP val = PROTECT(duplicate(b));      int sz = n * nrhs;
143      char *uplo = CHAR(STRING_ELT(GET_SLOT(a, Matrix_uploSym), 0));      char *uplo = CHAR(STRING_ELT(GET_SLOT(a, Matrix_uploSym), 0));
144      double *ax = REAL(GET_SLOT(a, Matrix_xSym)),      double *ax = REAL(GET_SLOT(a, Matrix_xSym)), one = 1., zero = 0.,
145          *vx = REAL(GET_SLOT(val, Matrix_xSym)), one = 1., zero = 0.;          *bx = (classed ? REAL(GET_SLOT(b, Matrix_xSym)) : REAL(b)),
146            *vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz));
147    
148      if (bdims[0] != n)      if (bdims[0] != n)
149          error(_("Matrices are not conformable for multiplication"));          error(_("Matrices are not conformable for multiplication"));
150      if (nrhs < 1 || n < 1)      if (nrhs < 1 || n < 1)
151          error(_("Matrices with zero extents cannot be multiplied"));          error(_("Matrices with zero extents cannot be multiplied"));
     for (i = 0; i < nrhs; i++)  
         F77_CALL(dspmv)(uplo, &n, &one, ax, vx + i * n, &ione,  
                         &zero, vx + i * n, &ione);  
     UNPROTECT(1);  
     return val;  
 }  
   
 SEXP dspMatrix_matrix_mm(SEXP a, SEXP b)  
 {  
     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),  
         *bdims = INTEGER(getAttrib(b, R_DimSymbol));  
     int i, ione = 1, n = adims[0], nrhs = bdims[1];  
     SEXP val = PROTECT(duplicate(b));  
     char *uplo = CHAR(STRING_ELT(GET_SLOT(a, Matrix_uploSym), 0));  
     double *ax = REAL(GET_SLOT(a, Matrix_xSym)),  
         *vx = REAL(val), one = 1., zero = 0.;  
152    
153      if (bdims[0] != n)      SET_SLOT(val, Matrix_DimSym, duplicate(bdimP));
         error(_("Matrices are not conformable for multiplication"));  
     if (nrhs < 1 || n < 1)  
         error(_("Matrices with zero extents cannot be multiplied"));  
154      for (i = 0; i < nrhs; i++)      for (i = 0; i < nrhs; i++)
155          F77_CALL(dspmv)(uplo, &n, &one, ax, vx + i * n, &ione,          F77_CALL(dspmv)(uplo, &n, &one, ax, bx + i * n, &ione,
156                          &zero, vx + i * n, &ione);                          &zero, vx + i * n, &ione);
157      UNPROTECT(1);      UNPROTECT(1);
158      return val;      return val;

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