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 642, Tue Mar 15 00:49:32 2005 UTC revision 643, Tue Mar 15 00:50:48 2005 UTC
# Line 43  Line 43 
43      typnm[0] = rcond_type(typstr);      typnm[0] = rcond_type(typstr);
44      rcond = get_double_by_name(rcv, typnm);      rcond = get_double_by_name(rcv, typnm);
45    
 /* FIXME: Need a factorization here. */  
46      if (R_IsNA(rcond)) {      if (R_IsNA(rcond)) {
47            SEXP trf = dsyMatrix_trf(obj);
48          int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;          int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
49          double anorm = get_norm_sy(obj, "O");          double anorm = get_norm_sy(obj, "O");
50    
51          error(_("Code for set_rcond_sy not yet written"));          F77_CALL(dsycon)(CHAR(asChar(GET_SLOT(trf, Matrix_uploSym))),
52          F77_CALL(dsycon)(CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),                           dims, REAL(GET_SLOT(trf, Matrix_xSym)),
53                           dims, REAL(GET_SLOT(obj, Matrix_xSym)),                           dims, INTEGER(GET_SLOT(trf, Matrix_permSym)),
                          dims, INTEGER(GET_SLOT(obj, install("pivot"))),  
54                           &anorm, &rcond,                           &anorm, &rcond,
55                           (double *) R_alloc(2*dims[0], sizeof(double)),                           (double *) R_alloc(2*dims[0], sizeof(double)),
56                           (int *) R_alloc(dims[0], sizeof(int)), &info);                           (int *) R_alloc(dims[0], sizeof(int)), &info);
# Line 63  Line 62 
62    
63  SEXP dsyMatrix_rcond(SEXP obj, SEXP type)  SEXP dsyMatrix_rcond(SEXP obj, SEXP type)
64  {  {
65  /* FIXME: This is a stub */      return ScalarReal(set_rcond_sy(obj, CHAR(asChar(type))));
 /*     return ScalarReal(set_rcond_sy(obj, CHAR(asChar(type)))); */  
     return ScalarReal(NA_REAL);  
66  }  }
67    
68  static  static
# Line 89  Line 86 
86    
87  SEXP dsyMatrix_solve(SEXP a)  SEXP dsyMatrix_solve(SEXP a)
88  {  {
89  /* FIXME: Write the code */      SEXP trf = dsyMatrix_trf(a);
90      error(_("code for dsyMatrix_solve not yet written"));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix")));
91      return R_NilValue;      int *dims = INTEGER(GET_SLOT(trf, Matrix_DimSym)), info;
92    
93        SET_SLOT(val, Matrix_uploSym, duplicate(GET_SLOT(trf, Matrix_uploSym)));
94        SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(trf, Matrix_xSym)));
95        SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(trf, Matrix_DimSym)));
96        SET_SLOT(val, Matrix_rcondSym, duplicate(GET_SLOT(a, Matrix_rcondSym)));
97        F77_CALL(dsytri)(CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),
98                         dims, REAL(GET_SLOT(val, Matrix_xSym)), dims,
99                         INTEGER(GET_SLOT(trf, Matrix_permSym)),
100                         (double *) R_alloc((long) dims[0], sizeof(double)),
101                         &info);
102        UNPROTECT(1);
103        return val;
104    }
105    
106    SEXP dsyMatrix_dgeMatrix_solve(SEXP a, SEXP b)
107    {
108        SEXP trf = dsyMatrix_trf(a),
109            val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
110        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
111            *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
112            info;
113    
114        if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
115            error(_("Dimensions of system to be solved are inconsistent"));
116        SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(b, Matrix_DimSym)));
117        SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(b, Matrix_xSym)));
118        F77_CALL(dsytrs)(CHAR(asChar(GET_SLOT(trf, Matrix_uploSym))),
119                         adims, bdims + 1,
120                         REAL(GET_SLOT(trf, Matrix_xSym)), adims,
121                         INTEGER(GET_SLOT(trf, Matrix_permSym)),
122                         REAL(GET_SLOT(val, Matrix_xSym)),
123                         bdims, &info);
124        UNPROTECT(1);
125        return val;
126  }  }
127    
128  SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)  SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)
129  {  {
130  /* FIXME: Write the code */      SEXP trf = dsyMatrix_trf(a),
131      error(_("code for dsyMatrix_matrix_solve not yet written"));          val = PROTECT(duplicate(b));
132      return R_NilValue;      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
133            *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
134            info;
135    
136        if (!(isReal(b) && isMatrix(b)))
137            error(_("Argument b must be a numeric matrix"));
138        if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
139            error(_("Dimensions of system to be solved are inconsistent"));
140        F77_CALL(dsytrs)(CHAR(asChar(GET_SLOT(trf, Matrix_uploSym))),
141                         adims, bdims + 1,
142                         REAL(GET_SLOT(trf, Matrix_xSym)), adims,
143                         INTEGER(GET_SLOT(trf, Matrix_permSym)),
144                         REAL(val), bdims, &info);
145        UNPROTECT(1);
146        return val;
147  }  }
148    
149  SEXP dsyMatrix_as_dgeMatrix(SEXP from)  SEXP dsyMatrix_as_dgeMatrix(SEXP from)
# Line 217  Line 262 
262      return set_factors(x, val, "BunchKaufman");      return set_factors(x, val, "BunchKaufman");
263  }  }
264    
265    SEXP dsyMatrix_as_dspMatrix(SEXP from)
266    {
267        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dspMatrix"))),
268            uplo = GET_SLOT(from, Matrix_uploSym),
269            dimP = GET_SLOT(from, Matrix_DimSym);
270        int n = *INTEGER(dimP);
271    
272        SET_SLOT(val, Matrix_rcondSym,
273                 duplicate(GET_SLOT(from, Matrix_rcondSym)));
274        SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
275        SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
276        full_to_packed(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
277                       REAL(GET_SLOT(from, Matrix_xSym)), n,
278                       *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW, NUN);
279        UNPROTECT(1);
280        return val;
281    }

Legend:
Removed from v.642  
changed lines
  Added in v.643

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