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

pkg/src/geMatrix.c revision 301, Sun Oct 17 05:06:04 2004 UTC pkg/src/dgeMatrix.c revision 652, Wed Mar 16 01:08:36 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),
7            fact = GET_SLOT(obj, Matrix_factorSym),
8            rc = GET_SLOT(obj, Matrix_rcondSym);
9      int m, n;      int m, n;
10    
11      if (length(Dim) != 2)      if (length(Dim) != 2)
12          return ScalarString(mkChar("Dim slot must have length 2"));          return mkString(_("Dim slot must have length 2"));
13      m = INTEGER(Dim)[0]; n = INTEGER(Dim)[1];      m = INTEGER(Dim)[0]; n = INTEGER(Dim)[1];
14      if (m < 0 || n < 0)      if (m < 0 || n < 0)
15          return          return mkString(_("Negative value(s) in Dim"));
             ScalarString(mkChar("Negative value(s) in Dim"));  
16      if (length(x) != m * n)      if (length(x) != m * n)
17          return          return mkString(_("length of x slot != prod(Dim)"));
18              ScalarString(mkChar("x slot must have length that is prod(Dim)"));      if (!isReal(x))
19            return mkString(_("x slot must be numeric \"double\""));
20        if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)
21            return mkString(_("factors slot must be named list"));
22        if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)
23            return mkString(_("rcond slot must be named numeric vector"));
24      return ScalarLogical(1);      return ScalarLogical(1);
25  }  }
26    
# Line 34  Line 40 
40                              dims, work);                              dims, work);
41  }  }
42    
43  SEXP geMatrix_norm(SEXP obj, SEXP type)  SEXP dgeMatrix_norm(SEXP obj, SEXP type)
44  {  {
45      return ScalarReal(get_norm(obj, CHAR(asChar(type))));      return ScalarReal(get_norm(obj, CHAR(asChar(type))));
46  }  }
# Line 43  Line 49 
49  double set_rcond(SEXP obj, char *typstr)  double set_rcond(SEXP obj, char *typstr)
50  {  {
51      char typnm[] = {'\0', '\0'};      char typnm[] = {'\0', '\0'};
52      SEXP rcv = GET_SLOT(obj, install("rcond"));      SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
53      double rcond = get_double_by_name(rcv, typnm);      double rcond = get_double_by_name(rcv, typnm);
54    
55      typnm[0] = rcond_type(typstr);      typnm[0] = rcond_type(typstr);
56      if (R_IsNA(rcond)) {      if (R_IsNA(rcond)) {
57          SEXP LU = geMatrix_LU(obj);          SEXP LU = dgeMatrix_LU(obj);
58          int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;          int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;
59          double anorm = get_norm(obj, typstr);          double anorm = get_norm(obj, typstr);
60    
61          if (dims[0] != dims[1] || dims[0] < 1)          if (dims[0] != dims[1] || dims[0] < 1)
62              error("rcond requires a square, non-empty matrix");              error(_("rcond requires a square, non-empty matrix"));
63          F77_CALL(dgecon)(typnm,          F77_CALL(dgecon)(typnm,
64                           dims, REAL(GET_SLOT(LU, Matrix_xSym)),                           dims, REAL(GET_SLOT(LU, Matrix_xSym)),
65                           dims, &anorm, &rcond,                           dims, &anorm, &rcond,
66                           (double *) R_alloc(4*dims[0], sizeof(double)),                           (double *) R_alloc(4*dims[0], sizeof(double)),
67                           (int *) R_alloc(dims[0], sizeof(int)), &info);                           (int *) R_alloc(dims[0], sizeof(int)), &info);
68          SET_SLOT(obj, install("rcond"),          SET_SLOT(obj, Matrix_rcondSym,
69                   set_double_by_name(rcv, rcond, typnm));                   set_double_by_name(rcv, rcond, typnm));
70      }      }
71      return rcond;      return rcond;
72  }  }
73    
74  SEXP geMatrix_rcond(SEXP obj, SEXP type)  SEXP dgeMatrix_rcond(SEXP obj, SEXP type)
75  {  {
76    return ScalarReal(set_rcond(obj, CHAR(asChar(type))));    return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
77  }  }
78    
79  SEXP geMatrix_crossprod(SEXP x)  SEXP dgeMatrix_crossprod(SEXP x)
80  {  {
81      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("poMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
82      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
83          *vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));          *vDims;
84      int n = Dims[1];      int i, n = Dims[1];
85      double one = 1.0, zero = 0.0;      int nsqr = n * n;
86        double one = 1.0, *xvals, zero = 0.0;
87      if ((*Dims) > 0 && n > 0) {  
88        SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
89        SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
90        SET_SLOT(val, Matrix_uploSym, mkString("U"));
91        SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
92        vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
93          vDims[0] = vDims[1] = n;          vDims[0] = vDims[1] = n;
94          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, n * n));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nsqr));
95          F77_CALL(dsyrk)(CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),      xvals = REAL(GET_SLOT(val, Matrix_xSym));
96                         "T", Dims + 1, Dims, &one,      for (i = 0; i < nsqr; i++) xvals[i] = 0.; /* keep valgrind happy */
97                         REAL(GET_SLOT(x, Matrix_xSym)),      F77_CALL(dsyrk)("U", "T", vDims, Dims,
98                         Dims, &zero,                      &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims,
99                         REAL(GET_SLOT(val, Matrix_xSym)), &n);                      &zero, xvals, vDims);
     }  
100      UNPROTECT(1);      UNPROTECT(1);
101      return val;      return val;
102  }  }
103    
104  SEXP geMatrix_geMatrix_crossprod(SEXP x, SEXP y)  SEXP dgeMatrix_dgeMatrix_crossprod(SEXP x, SEXP y)
105  {  {
106      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
107      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
108          *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),          *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),
109          *vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));          *vDims;
110      int m = xDims[1], n = yDims[1];      int m = xDims[1], n = yDims[1];
111      double one = 1.0, zero = 0.0;      double one = 1.0, zero = 0.0;
112    
113        SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
114        SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
115        SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
116        vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
117      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
118          if (*xDims != *yDims)          if (*xDims != *yDims)
119              error("Dimensions of x and y are not compatible for crossprod");              error(_("Dimensions of x and y are not compatible for crossprod"));
120          vDims[0] = m; vDims[1] = n;          vDims[0] = m; vDims[1] = n;
121          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
122          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
# Line 115  Line 129 
129      return val;      return val;
130  }  }
131    
132  SEXP geMatrix_matrix_crossprod(SEXP x, SEXP y)  SEXP dgeMatrix_matrix_crossprod(SEXP x, SEXP y)
133  {  {
134      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
135      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
136          *yDims = INTEGER(getAttrib(y, R_DimSymbol)),          *yDims = INTEGER(getAttrib(y, R_DimSymbol)),
137          *vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));          *vDims;
138      int m = xDims[1], n = yDims[1];      int m = xDims[1], n = yDims[1];
139      double one = 1.0, zero = 0.0;      double one = 1.0, zero = 0.0;
140    
141      if (!(isMatrix(y) && isReal(y)))      if (!(isMatrix(y) && isReal(y)))
142          error("Argument y must be a numeric (real) matrix");          error(_("Argument y must be a numeric (real) matrix"));
143        SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
144        SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
145        SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
146        vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
147      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
148          if (*xDims != *yDims)          if (*xDims != *yDims)
149              error("Dimensions of x and y are not compatible for crossprod");              error(_("Dimensions of x and y are not compatible for crossprod"));
150          vDims[0] = m; vDims[1] = n;          vDims[0] = m; vDims[1] = n;
151          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
152          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
# Line 141  Line 159 
159      return val;      return val;
160  }  }
161    
162  SEXP geMatrix_getDiag(SEXP x)  SEXP dgeMatrix_getDiag(SEXP x)
163  {  {
164      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
165      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 155  Line 173 
173      return ret;      return ret;
174  }  }
175    
176  SEXP geMatrix_LU(SEXP x)  SEXP dgeMatrix_LU(SEXP x)
177  {  {
178      SEXP val = get_factorization(x, "LU");      SEXP val = get_factors(x, "LU");
179      int *dims, npiv, info;      int *dims, npiv, info;
180    
181      if (val != R_NilValue) return val;      if (val != R_NilValue) return val;
182      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
183      if (dims[0] < 1 || dims[1] < 1)      if (dims[0] < 1 || dims[1] < 1)
184          error("Cannot factor a matrix with zero extents");          error(_("Cannot factor a matrix with zero extents"));
185      npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];      npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];
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_factorization(x, val, "LU");      return set_factors(x, val, "LU");
197  }  }
198    
199  SEXP geMatrix_determinant(SEXP x, SEXP logarithm)  SEXP dgeMatrix_determinant(SEXP x, SEXP logarithm)
200  {  {
201      int lg = asLogical(logarithm);      int lg = asLogical(logarithm);
202      SEXP lu = geMatrix_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    
208      if (n != dims[1])      if (n != dims[1])
209          error("Determinant requires a square matrix");          error(_("Determinant requires a square matrix"));
210      for (i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign;      for (i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign;
211      if (lg) {      if (lg) {
212          modulus = 0.0;          modulus = 0.0;
# Line 208  Line 227 
227      return as_det_obj(modulus, lg, sign);      return as_det_obj(modulus, lg, sign);
228  }  }
229    
230  SEXP geMatrix_solve(SEXP a)  SEXP dgeMatrix_solve(SEXP a)
231  {  {
232      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
233          lu = geMatrix_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"));
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 231  Line 250 
250      return val;      return val;
251  }  }
252    
253  SEXP geMatrix_geMatrix_mm(SEXP a, SEXP b)  SEXP dgeMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
254  {  {
255        int cl = asLogical(classed);
256        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
257            lu = dgeMatrix_LU(a);
258        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")));
280      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
281          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
282          *cdims,                           getAttrib(b, R_DimSymbol)),
283          m = adims[0], n = bdims[1], k = adims[1];          *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));  
     char *trans = "N";  
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_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  SEXP geMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)  
316    SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
317  {  {
318      int nu = asInteger(nnu),      int /* nu = asInteger(nnu),
319          nv = asInteger(nnv),             nv = asInteger(nnv), */
320          *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));          *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
321      double *xx = REAL(GET_SLOT(x, Matrix_xSym));      double *xx = REAL(GET_SLOT(x, Matrix_xSym));
322      SEXP val = PROTECT(allocVector(VECSXP, 3));      SEXP val = PROTECT(allocVector(VECSXP, 3));
# Line 270  Line 326 
326              lwork = -1, info;              lwork = -1, info;
327          int *iwork = Calloc(8 * mm, int);          int *iwork = Calloc(8 * mm, int);
328          double tmp, *work;          double tmp, *work;
329          int bdspac = 3*m*m + 4*m,  /*      int bdspac = 3*m*m + 4*m, */
330              wrkbl, maxwrk, minwrk, itmp,  /*          wrkbl, maxwrk, minwrk, itmp, */
331              ione = 1, iminus1 = -1;  /*          ione = 1, iminus1 = -1; */
332          int i1, i2, i3;  /*      int i1, i2, i3; */
333    
334          SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));          SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
335          SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));          SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
# Line 284  Line 340 
340                           REAL(VECTOR_ELT(val, 2)), &mm,                           REAL(VECTOR_ELT(val, 2)), &mm,
341                           &tmp, &lwork, iwork, &info);                           &tmp, &lwork, iwork, &info);
342          lwork = (int) tmp;          lwork = (int) tmp;
343          F77_CALL(foo)(&i1, &i2, &i3);  /*      F77_CALL(foo)(&i1, &i2, &i3); */
344          wrkbl = 3*m+(m+n)*i1;  /*      wrkbl = 3*m+(m+n)*i1; */
345          if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp;  /*      if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp; */
346          if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp;  /*      if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp; */
347          itmp = bdspac+3*m;  /*      itmp = bdspac+3*m; */
348          maxwrk = (wrkbl > itmp) ? wrkbl : itmp;  /*      maxwrk = (wrkbl > itmp) ? wrkbl : itmp; */
349          minwrk = 3*m + ((bdspac > n) ?  bdspac : n);  /*      minwrk = 3*m + ((bdspac > n) ?  bdspac : n); */
350          work = Calloc(lwork, double);          work = Calloc(lwork, double);
351          F77_CALL(dgesdd)("S", &m, &n, xx, &m,          F77_CALL(dgesdd)("S", &m, &n, xx, &m,
352                           REAL(VECTOR_ELT(val, 0)),                           REAL(VECTOR_ELT(val, 0)),
# Line 302  Line 358 
358      UNPROTECT(1);      UNPROTECT(1);
359      return val;      return val;
360  }  }
361    
362    static double padec [] =   /*  Constants for matrix exponential calculation. */
363    {
364      5.0000000000000000e-1,
365      1.1666666666666667e-1,
366      1.6666666666666667e-2,
367      1.6025641025641026e-3,
368      1.0683760683760684e-4,
369      4.8562548562548563e-6,
370      1.3875013875013875e-7,
371      1.9270852604185938e-9,
372    };
373    
374    /**
375     * Matrix exponential - based on the code for Octave's expm function.
376     *
377     * @param x real square matrix to exponentiate
378     *
379     * @return matrix exponential of x
380     */
381    SEXP dgeMatrix_exp(SEXP x)
382    {
383        SEXP val = PROTECT(duplicate(x));
384        int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
385        int i, ilo, ilos, ihi, ihis, j, nc = Dims[1], sqpow;
386        int ncp1 = Dims[1] + 1, ncsqr = nc * nc;
387        int *pivot = Calloc(nc, int);
388        int *iperm = Calloc(nc, int);
389        double *dpp = Calloc(ncsqr, double), /* denominator power Pade' */
390            *npp = Calloc(ncsqr, double), /* numerator power Pade' */
391            *perm = Calloc(nc, double),
392            *scale = Calloc(nc, double),
393            *v = REAL(GET_SLOT(val, Matrix_xSym)),
394            *work = Calloc(ncsqr, double), inf_norm, m1_j, /* (-1)^j */
395            one = 1., trshift, zero = 0.;
396    
397        if (nc < 1 || Dims[0] != nc)
398            error(_("Matrix exponential requires square, non-null matrix"));
399    
400        /* FIXME: Add special treatment for nc == 1 */
401    
402        /* Preconditioning 1.  Shift diagonal by average diagonal if positive. */
403        trshift = 0;                /* determine average diagonal element */
404        for (i = 0; i < nc; i++) trshift += v[i * ncp1];
405        trshift /= nc;
406        if (trshift > 0.) {         /* shift diagonal by -trshift */
407            for (i = 0; i < nc; i++) v[i * ncp1] -= trshift;
408        }
409    
410        /* Preconditioning 2. Balancing with dgebal. */
411        F77_CALL(dgebal)("P", &nc, v, &nc, &ilo, &ihi, perm, &j);
412        if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
413        F77_CALL(dgebal)("S", &nc, v, &nc, &ilos, &ihis, scale, &j);
414        if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
415    
416        /* Preconditioning 3. Scaling according to infinity norm */
417        inf_norm = F77_CALL(dlange)("I", &nc, &nc, v, &nc, work);
418        sqpow = (inf_norm > 0) ? (int) (1 + log(inf_norm)/log(2.)) : 0;
419        if (sqpow < 0) sqpow = 0;
420        if (sqpow > 0) {
421            double scale_factor = 1.0;
422            for (i = 0; i < sqpow; i++) scale_factor *= 2.;
423            for (i = 0; i < ncsqr; i++) v[i] /= scale_factor;
424        }
425    
426        /* Pade' approximation. Powers v^8, v^7, ..., v^1 */
427        AZERO(npp, ncsqr);
428        AZERO(dpp, ncsqr);
429        m1_j = -1;
430        for (j = 7; j >=0; j--) {
431            double mult = padec[j];
432            /* npp = m * npp + padec[j] *m */
433            F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, npp, &nc,
434                            &zero, work, &nc);
435            for (i = 0; i < ncsqr; i++) npp[i] = work[i] + mult * v[i];
436            /* dpp = m * dpp * (m1_j * padec[j]) * m */
437            mult *= m1_j;
438            F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, dpp, &nc,
439                            &zero, work, &nc);
440            for (i = 0; i < ncsqr; i++) dpp[i] = work[i] + mult * v[i];
441            m1_j *= -1;
442        }
443        /* Zero power */
444        for (i = 0; i < ncsqr; i++) dpp[i] *= -1.;
445        for (j = 0; j < nc; j++) {
446            npp[j * ncp1] += 1.;
447            dpp[j * ncp1] += 1.;
448        }
449    
450        /* Pade' approximation is solve(dpp, npp) */
451        F77_CALL(dgetrf)(&nc, &nc, dpp, &nc, pivot, &j);
452        if (j) error(_("dgeMatrix_exp: dgetrf returned error code %d"), j);
453        F77_CALL(dgetrs)("N", &nc, &nc, dpp, &nc, pivot, npp, &nc, &j);
454        if (j) error(_("dgeMatrix_exp: dgetrs returned error code %d"), j);
455        Memcpy(v, npp, ncsqr);
456    
457        /* Now undo all of the preconditioning */
458        /* Preconditioning 3: square the result for every power of 2 */
459        while (sqpow--) {
460            F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, v, &nc,
461                            &zero, work, &nc);
462            Memcpy(v, work, ncsqr);
463        }
464        /* Preconditioning 2: apply inverse scaling */
465        for (j = 0; j < nc; j++)
466            for (i = 0; i < nc; i++)
467                v[i + j * nc] *= scale[i]/scale[j];
468        /* Construct balancing permutation vector */
469        for (i = 0; i < nc; i++) iperm[i] = i; /* identity permutation */
470        /* Leading permutations applied in forward order */
471        for (i = 0; i < (ilo - 1); i++) {
472            int swapidx = (int) (perm[i]) - 1;
473            int tmp = iperm[i];
474            iperm[i] = iperm[swapidx];
475            iperm[swapidx] = tmp;
476        }
477        /* Trailing permutations applied in reverse order */
478        for (i = nc - 1; i >= ihi; i--) {
479            int swapidx = (int) (perm[i]) - 1;
480            int tmp = iperm[i];
481            iperm[i] = iperm[swapidx];
482            iperm[swapidx] = tmp;
483        }
484        /* Construct inverse balancing permutation vector */
485        Memcpy(pivot, iperm, nc);
486        for (i = 0; i < nc; i++) iperm[pivot[i]] = i;
487        /* Apply inverse permutation */
488        Memcpy(work, v, ncsqr);
489        for (j = 0; j < nc; j++)
490            for (i = 0; i < nc; i++)
491                v[i + j * nc] = work[iperm[i] + iperm[j] * nc];
492    
493        /* Preconditioning 1: Trace normalization */
494        if (trshift > 0.) {
495            double mult = exp(trshift);
496            for (i = 0; i < ncsqr; i++) v[i] *= mult;
497        }
498    
499        /* Clean up */
500        Free(dpp); Free(npp); Free(perm); Free(iperm); Free(pivot); Free(scale); Free(work);
501        UNPROTECT(1);
502        return val;
503    }
504    
505    SEXP dgeMatrix_Schur(SEXP x, SEXP vectors)
506    {
507        int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
508        int vecs = asLogical(vectors), info, izero = 0, lwork = -1, n = dims[0];
509        double *work, tmp;
510        char *nms[] = {"WR", "WI", "T", "Z", ""};
511        SEXP val = PROTECT(Matrix_make_named(VECSXP, nms));
512    
513        if (n != dims[1] || n < 1)
514            error(_("dgeMatrix_Schur: argument x must be a non-null square matrix"));
515        SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));
516        SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));
517        SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));
518        Memcpy(REAL(VECTOR_ELT(val, 2)), REAL(GET_SLOT(x, Matrix_xSym)), n * n);
519        SET_VECTOR_ELT(val, 3, allocMatrix(REALSXP, vecs ? n : 0, vecs ? n : 0));
520        F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
521                        (double *) NULL, (double *) NULL, (double *) NULL, dims,
522                        &tmp, &lwork, (int *) NULL, &info);
523        if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
524        lwork = (int) tmp;
525        work = Calloc(lwork, double);
526        F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
527                        &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
528                        REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
529                        (int *) NULL, &info);
530        if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
531        Free(work);
532        UNPROTECT(1);
533        return val;
534    }

Legend:
Removed from v.301  
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