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 477, Wed Feb 2 14:14:59 2005 UTC revision 478, Wed Feb 2 14:33:51 2005 UTC
# Line 1  Line 1 
1  #include "geMatrix.h"  #include "dgeMatrix.h"
2    
3  SEXP geMatrix_validate(SEXP obj)  SEXP dgeMatrix_validate(SEXP obj)
4  {  {
5      SEXP x = GET_SLOT(obj, Matrix_xSym),      SEXP x = GET_SLOT(obj, Matrix_xSym),
6          Dim = GET_SLOT(obj, Matrix_DimSym),          Dim = GET_SLOT(obj, Matrix_DimSym),
# Line 38  Line 38 
38                              dims, work);                              dims, work);
39  }  }
40    
41  SEXP geMatrix_norm(SEXP obj, SEXP type)  SEXP dgeMatrix_norm(SEXP obj, SEXP type)
42  {  {
43      return ScalarReal(get_norm(obj, CHAR(asChar(type))));      return ScalarReal(get_norm(obj, CHAR(asChar(type))));
44  }  }
# Line 52  Line 52 
52    
53      typnm[0] = rcond_type(typstr);      typnm[0] = rcond_type(typstr);
54      if (R_IsNA(rcond)) {      if (R_IsNA(rcond)) {
55          SEXP LU = geMatrix_LU(obj);          SEXP LU = dgeMatrix_LU(obj);
56          int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;          int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;
57          double anorm = get_norm(obj, typstr);          double anorm = get_norm(obj, typstr);
58    
# Line 69  Line 69 
69      return rcond;      return rcond;
70  }  }
71    
72  SEXP geMatrix_rcond(SEXP obj, SEXP type)  SEXP dgeMatrix_rcond(SEXP obj, SEXP type)
73  {  {
74    return ScalarReal(set_rcond(obj, CHAR(asChar(type))));    return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
75  }  }
76    
77  SEXP geMatrix_crossprod(SEXP x)  SEXP dgeMatrix_crossprod(SEXP x)
78  {  {
79      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("poMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
80      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
81          *vDims;          *vDims;
82      int i, n = Dims[1];      int i, n = Dims[1];
# Line 99  Line 99 
99      return val;      return val;
100  }  }
101    
102  SEXP geMatrix_geMatrix_crossprod(SEXP x, SEXP y)  SEXP dgeMatrix_dgeMatrix_crossprod(SEXP x, SEXP y)
103  {  {
104      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
105      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
106          *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),          *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),
107          *vDims;          *vDims;
# Line 127  Line 127 
127      return val;      return val;
128  }  }
129    
130  SEXP geMatrix_matrix_crossprod(SEXP x, SEXP y)  SEXP dgeMatrix_matrix_crossprod(SEXP x, SEXP y)
131  {  {
132      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
133      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
134          *yDims = INTEGER(getAttrib(y, R_DimSymbol)),          *yDims = INTEGER(getAttrib(y, R_DimSymbol)),
135          *vDims;          *vDims;
# Line 157  Line 157 
157      return val;      return val;
158  }  }
159    
160  SEXP geMatrix_getDiag(SEXP x)  SEXP dgeMatrix_getDiag(SEXP x)
161  {  {
162      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
163      int i, m = dims[0], nret = (m < dims[1]) ? m : dims[1];      int i, m = dims[0], nret = (m < dims[1]) ? m : dims[1];
# Line 171  Line 171 
171      return ret;      return ret;
172  }  }
173    
174  SEXP geMatrix_LU(SEXP x)  SEXP dgeMatrix_LU(SEXP x)
175  {  {
176      SEXP val = get_factors(x, "LU");      SEXP val = get_factors(x, "LU");
177      int *dims, npiv, info;      int *dims, npiv, info;
# Line 193  Line 193 
193      return set_factors(x, val, "LU");      return set_factors(x, val, "LU");
194  }  }
195    
196  SEXP geMatrix_determinant(SEXP x, SEXP logarithm)  SEXP dgeMatrix_determinant(SEXP x, SEXP logarithm)
197  {  {
198      int lg = asLogical(logarithm);      int lg = asLogical(logarithm);
199      SEXP lu = geMatrix_LU(x);      SEXP lu = dgeMatrix_LU(x);
200      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
201          *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),          *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),
202          i, n = dims[0], sign = 1;          i, n = dims[0], sign = 1;
# Line 224  Line 224 
224      return as_det_obj(modulus, lg, sign);      return as_det_obj(modulus, lg, sign);
225  }  }
226    
227  SEXP geMatrix_solve(SEXP a)  SEXP dgeMatrix_solve(SEXP a)
228  {  {
229      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
230          lu = geMatrix_LU(a);          lu = dgeMatrix_LU(a);
231      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
232          *pivot = INTEGER(GET_SLOT(lu, install("pivot")));          *pivot = INTEGER(GET_SLOT(lu, install("pivot")));
233      double *x, tmp;      double *x, tmp;
# Line 249  Line 249 
249      return val;      return val;
250  }  }
251    
252  SEXP geMatrix_geMatrix_mm(SEXP a, SEXP b)  SEXP dgeMatrix_dgeMatrix_mm(SEXP a, SEXP b)
253  {  {
254      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
255          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
256          *cdims,          *cdims,
257          m = adims[0], n = bdims[1], k = adims[1];          m = adims[0], n = bdims[1], k = adims[1];
258      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
259      char *trans = "N";      char *trans = "N";
260      double one = 1., zero = 0.;      double one = 1., zero = 0.;
261    
# Line 277  Line 277 
277      return val;      return val;
278  }  }
279    
280  SEXP geMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)  SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
281  {  {
282      int /* nu = asInteger(nnu),      int /* nu = asInteger(nnu),
283             nv = asInteger(nnv), */             nv = asInteger(nnv), */
# Line 342  Line 342 
342   *   *
343   * @return matrix exponential of x   * @return matrix exponential of x
344   */   */
345  SEXP geMatrix_exp(SEXP x)  SEXP dgeMatrix_exp(SEXP x)
346  {  {
347      SEXP val = PROTECT(duplicate(x));      SEXP val = PROTECT(duplicate(x));
348      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
# Line 373  Line 373 
373    
374      /* Preconditioning 2. Balancing with dgebal. */      /* Preconditioning 2. Balancing with dgebal. */
375      F77_CALL(dgebal)("P", &nc, v, &nc, &ilo, &ihi, perm, &j);      F77_CALL(dgebal)("P", &nc, v, &nc, &ilo, &ihi, perm, &j);
376      if (j) error("geMatrix_exp: LAPACK routine dgebal returned %d", j);      if (j) error("dgeMatrix_exp: LAPACK routine dgebal returned %d", j);
377      F77_CALL(dgebal)("S", &nc, v, &nc, &ilos, &ihis, scale, &j);      F77_CALL(dgebal)("S", &nc, v, &nc, &ilos, &ihis, scale, &j);
378      if (j) error("geMatrix_exp: LAPACK routine dgebal returned %d", j);      if (j) error("dgeMatrix_exp: LAPACK routine dgebal returned %d", j);
379    
380      /* Preconditioning 3. Scaling according to infinity norm */      /* Preconditioning 3. Scaling according to infinity norm */
381      inf_norm = F77_CALL(dlange)("I", &nc, &nc, v, &nc, work);      inf_norm = F77_CALL(dlange)("I", &nc, &nc, v, &nc, work);
# Line 413  Line 413 
413    
414      /* Pade' approximation is solve(dpp, npp) */      /* Pade' approximation is solve(dpp, npp) */
415      F77_CALL(dgetrf)(&nc, &nc, dpp, &nc, pivot, &j);      F77_CALL(dgetrf)(&nc, &nc, dpp, &nc, pivot, &j);
416      if (j) error("geMatrix_exp: dgetrf returned error code %d", j);      if (j) error("dgeMatrix_exp: dgetrf returned error code %d", j);
417      F77_CALL(dgetrs)("N", &nc, &nc, dpp, &nc, pivot, npp, &nc, &j);      F77_CALL(dgetrs)("N", &nc, &nc, dpp, &nc, pivot, npp, &nc, &j);
418      if (j) error("geMatrix_exp: dgetrs returned error code %d", j);      if (j) error("dgeMatrix_exp: dgetrs returned error code %d", j);
419      Memcpy(v, npp, ncsqr);      Memcpy(v, npp, ncsqr);
420    
421      /* Now undo all of the preconditioning */      /* Now undo all of the preconditioning */

Legend:
Removed from v.477  
changed lines
  Added in v.478

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