SCM

SCM Repository

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

Diff of /pkg/src/dsyMatrix.c

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

revision 1462, Tue Aug 29 16:34:49 2006 UTC revision 1463, Tue Aug 29 22:30:57 2006 UTC
# Line 75  Line 75 
75      return val;      return val;
76  }  }
77    
78  SEXP dsyMatrix_dgeMatrix_solve(SEXP a, SEXP b)  SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)
79  {  {
80      SEXP trf = dsyMatrix_trf(a),      SEXP trf = dsyMatrix_trf(a),
81          val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));          val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
82      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
83          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),          *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym)),
84          info;          info;
85    
86      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
87          error(_("Dimensions of system to be solved are inconsistent"));          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)));  
88      F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,      F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,
89                       REAL(GET_SLOT(trf, Matrix_xSym)), adims,                       REAL(GET_SLOT(trf, Matrix_xSym)), adims,
90                       INTEGER(GET_SLOT(trf, Matrix_permSym)),                       INTEGER(GET_SLOT(trf, Matrix_permSym)),
# Line 96  Line 94 
94      return val;      return val;
95  }  }
96    
 SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)  
 {  
     SEXP trf = dsyMatrix_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)))  
         error(_("Argument b must be a numeric matrix"));  
     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)  
         error(_("Dimensions of system to be solved are inconsistent"));  
     F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,  
                      REAL(GET_SLOT(trf, Matrix_xSym)), adims,  
                      INTEGER(GET_SLOT(trf, Matrix_permSym)),  
                      REAL(val), bdims, &info);  
     UNPROTECT(1);  
     return val;  
 }  
   
 #if 0                           /* no longer used */  
 SEXP dsyMatrix_as_dgeMatrix(SEXP from)  
 {  
     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));  
   
     SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));  
     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(from, Matrix_xSym)));  
     SET_SLOT(val, Matrix_DimSym,  
              duplicate(GET_SLOT(from, Matrix_DimSym)));  
     SET_SLOT(val, Matrix_DimNamesSym,  
              duplicate(GET_SLOT(from, Matrix_DimNamesSym)));  
     make_d_matrix_symmetric(REAL(GET_SLOT(val, Matrix_xSym)), from);  
     UNPROTECT(1);  
     return val;  
 }  
 #endif  
   
97  SEXP dsyMatrix_as_matrix(SEXP from)  SEXP dsyMatrix_as_matrix(SEXP from)
98  {  {
99      int n = INTEGER(GET_SLOT(from, Matrix_DimSym))[0];      int n = INTEGER(GET_SLOT(from, Matrix_DimSym))[0];
# Line 146  Line 107 
107      return val;      return val;
108  }  }
109    
110  SEXP dsyMatrix_dgeMatrix_mm(SEXP a, SEXP b)  SEXP dsyMatrix_matrix_mm(SEXP a, SEXP b, SEXP rtP)
 {  
     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),  
         *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),  
         *cdims,  
         m = adims[0], n = bdims[1], k = adims[1];  
     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));  
     double one = 1., zero = 0.;  
   
     if (bdims[0] != k)  
         error(_("Matrices are not conformable for multiplication"));  
     if (m < 1 || n < 1 || k < 1)  
         error(_("Matrices with zero extents cannot be multiplied"));  
     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));  
     cdims[0] = m; cdims[1] = n;  
     F77_CALL(dsymm)("L", uplo_P(a), adims, bdims+1, &one,  
                     REAL(GET_SLOT(a, Matrix_xSym)), adims,  
                     REAL(GET_SLOT(b, Matrix_xSym)), bdims,  
                     &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);  
     UNPROTECT(1);  
     return val;  
 }  
   
 SEXP dsyMatrix_dgeMatrix_mm_R(SEXP a, SEXP b)  
111  {  {
112        SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
113      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
114          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
115          *cdims,          m = bdims[0], n = bdims[1], rt = asLogical(rtP);
         m = adims[0], n = bdims[1], k = adims[1];  
     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));  
116      double one = 1., zero = 0.;      double one = 1., zero = 0.;
117    
118      if (bdims[0] != k)      if ((rt && n != adims[0]) || (!rt && m != adims[0]))
119          error(_("Matrices are not conformable for multiplication"));          error(_("Matrices are not conformable for multiplication"));
120      if (m < 1 || n < 1 || k < 1)      if (m < 1 || n < 1)
121          error(_("Matrices with zero extents cannot be multiplied"));          error(_("Matrices with zero extents cannot be multiplied"));
122      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      F77_CALL(dsymm)(rt ? "R" :"L", uplo_P(a), &m, &n, &one,
     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));  
     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));  
     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));  
     cdims[0] = m; cdims[1] = n;  
     F77_CALL(dsymm)("R", uplo_P(a), adims, bdims+1, &one,  
123                      REAL(GET_SLOT(a, Matrix_xSym)), adims,                      REAL(GET_SLOT(a, Matrix_xSym)), adims,
124                      REAL(GET_SLOT(b, Matrix_xSym)), bdims,                      REAL(GET_SLOT(b, Matrix_xSym)), &m,
125                      &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);                      &zero, REAL(GET_SLOT(val, Matrix_xSym)), &m);
126      UNPROTECT(1);      UNPROTECT(1);
127      return val;      return val;
128  }  }

Legend:
Removed from v.1462  
changed lines
  Added in v.1463

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