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 476, Wed Feb 2 11:51:24 2005 UTC pkg/src/dgeMatrix.c revision 692, Mon Apr 18 14:34:33 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 9  Line 9 
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 ScalarString(mkChar("Negative value(s) in Dim"));          return mkString(_("Negative value(s) in Dim"));
16      if (length(x) != m * n)      if (length(x) != m * n)
17          return ScalarString(mkChar("length of x slot != prod(Dim)"));          return mkString(_("length of x slot != prod(Dim)"));
18        if (!isReal(x))
19            return mkString(_("x slot must be numeric \"double\""));
20      if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)      if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)
21          return ScalarString(mkChar("factors slot must be named list"));          return mkString(_("factors slot must be named list"));
22      if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)      if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)
23          return ScalarString(mkChar("rcond slot must be named numeric vector"));          return mkString(_("rcond slot must be named numeric vector"));
24      return ScalarLogical(1);      return ScalarLogical(1);
25  }  }
26    
# Line 38  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 52  Line 54 
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,
# Line 69  Line 71 
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, SEXP trans)
80  {  {
81      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("poMatrix")));      int tr = asLogical(trans);
82        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
83      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
84          *vDims;          *vDims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
85      int i, n = Dims[1];      int i, k = tr ? Dims[1] : Dims[0], n = tr ? Dims[0] : Dims[1];
     int nsqr = n * n;  
86      double one = 1.0, *xvals, zero = 0.0;      double one = 1.0, *xvals, zero = 0.0;
87    
88      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_uploSym, mkString("U"));
     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));  
     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("U")));  
     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));  
     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));  
89      vDims[0] = vDims[1] = n;      vDims[0] = vDims[1] = n;
90      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nsqr));      F77_CALL(dsyrk)("U", tr ? "N" : "T", &n, &k,
91      xvals = REAL(GET_SLOT(val, Matrix_xSym));                      &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims, &zero,
92      for (i = 0; i < nsqr; i++) xvals[i] = 0.; /* keep valgrind happy */                      REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n)), &n);
     F77_CALL(dsyrk)("U", "T", vDims, Dims,  
                     &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims,  
                     &zero, xvals, vDims);  
93      UNPROTECT(1);      UNPROTECT(1);
94      return val;      return val;
95  }  }
96    
97  SEXP geMatrix_geMatrix_crossprod(SEXP x, SEXP y)  SEXP dgeMatrix_dgeMatrix_crossprod(SEXP x, SEXP y)
98  {  {
99      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
100      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
101          *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),          *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),
102          *vDims;          *vDims;
# Line 114  Line 109 
109      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
110      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
111          if (*xDims != *yDims)          if (*xDims != *yDims)
112              error("Dimensions of x and y are not compatible for crossprod");              error(_("Dimensions of x and y are not compatible for crossprod"));
113          vDims[0] = m; vDims[1] = n;          vDims[0] = m; vDims[1] = n;
114          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
115          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
# Line 127  Line 122 
122      return val;      return val;
123  }  }
124    
125  SEXP geMatrix_matrix_crossprod(SEXP x, SEXP y)  SEXP dgeMatrix_matrix_crossprod(SEXP x, SEXP y)
126  {  {
127      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
128      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
129          *yDims = INTEGER(getAttrib(y, R_DimSymbol)),          *yDims = INTEGER(getAttrib(y, R_DimSymbol)),
130          *vDims;          *vDims;
# Line 137  Line 132 
132      double one = 1.0, zero = 0.0;      double one = 1.0, zero = 0.0;
133    
134      if (!(isMatrix(y) && isReal(y)))      if (!(isMatrix(y) && isReal(y)))
135          error("Argument y must be a numeric (real) matrix");          error(_("Argument y must be a numeric (real) matrix"));
136      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
137      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
138      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
139      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
140      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
141          if (*xDims != *yDims)          if (*xDims != *yDims)
142              error("Dimensions of x and y are not compatible for crossprod");              error(_("Dimensions of x and y are not compatible for crossprod"));
143          vDims[0] = m; vDims[1] = n;          vDims[0] = m; vDims[1] = n;
144          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
145          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
# Line 157  Line 152 
152      return val;      return val;
153  }  }
154    
155  SEXP geMatrix_getDiag(SEXP x)  SEXP dgeMatrix_getDiag(SEXP x)
156  {  {
157      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
158      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 166 
166      return ret;      return ret;
167  }  }
168    
169  SEXP geMatrix_LU(SEXP x)  SEXP dgeMatrix_LU(SEXP x)
170  {  {
171      SEXP val = get_factors(x, "LU");      SEXP val = get_factors(x, "LU");
172      int *dims, npiv, info;      int *dims, npiv, info;
# Line 179  Line 174 
174      if (val != R_NilValue) return val;      if (val != R_NilValue) return val;
175      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
176      if (dims[0] < 1 || dims[1] < 1)      if (dims[0] < 1 || dims[1] < 1)
177          error("Cannot factor a matrix with zero extents");          error(_("Cannot factor a matrix with zero extents"));
178      npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];      npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];
179      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
180      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
181      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));  
182      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
183                       dims, INTEGER(GET_SLOT(val, install("pivot"))),                       dims,
184                         INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, npiv)),
185                       &info);                       &info);
186      if (info) error("Lapack routine dgetrf returned error code %d", info);      if (info)
187            error(_("Lapack routine %s returned error code %d"), "dgetrf", info);
188      UNPROTECT(1);      UNPROTECT(1);
189      return set_factors(x, val, "LU");      return set_factors(x, val, "LU");
190  }  }
191    
192  SEXP geMatrix_determinant(SEXP x, SEXP logarithm)  SEXP dgeMatrix_determinant(SEXP x, SEXP logarithm)
193  {  {
194      int lg = asLogical(logarithm);      int lg = asLogical(logarithm);
195      SEXP lu = geMatrix_LU(x);      SEXP lu = dgeMatrix_LU(x);
196      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
197          *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),          *jpvt = INTEGER(GET_SLOT(lu, Matrix_permSym)),
198          i, n = dims[0], sign = 1;          i, n = dims[0], sign = 1;
199      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;
200    
201      if (n != dims[1])      if (n != dims[1])
202          error("Determinant requires a square matrix");          error(_("Determinant requires a square matrix"));
203      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;
204      if (lg) {      if (lg) {
205          modulus = 0.0;          modulus = 0.0;
# Line 224  Line 220 
220      return as_det_obj(modulus, lg, sign);      return as_det_obj(modulus, lg, sign);
221  }  }
222    
223  SEXP geMatrix_solve(SEXP a)  SEXP dgeMatrix_solve(SEXP a)
224  {  {
225      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
226          lu = geMatrix_LU(a);          lu = dgeMatrix_LU(a);
227      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
228          *pivot = INTEGER(GET_SLOT(lu, install("pivot")));          *pivot = INTEGER(GET_SLOT(lu, Matrix_permSym));
229      double *x, tmp;      double *x, tmp;
230      int info, lwork = -1;      int info, lwork = -1;
231    
232    
233      if (dims[0] != dims[1]) error("Solve requires a square matrix");      if (dims[0] != dims[1]) error(_("Solve requires a square matrix"));
     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));  
     SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));  
234      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
235      x = REAL(GET_SLOT(val, Matrix_xSym));      x = REAL(GET_SLOT(val, Matrix_xSym));
236      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));
# Line 249  Line 243 
243      return val;      return val;
244  }  }
245    
246  SEXP geMatrix_geMatrix_mm(SEXP a, SEXP b)  SEXP dgeMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
247  {  {
248        int cl = asLogical(classed);
249        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
250            lu = dgeMatrix_LU(a);
251        int *adims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
252            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
253                             getAttrib(b, R_DimSymbol));
254        int info, n = bdims[0], nrhs = bdims[1];
255        int sz = n * nrhs;
256    
257        if (*adims != *bdims || bdims[1] < 1 || *adims < 1 || *adims != adims[1])
258            error(_("Dimensions of system to be solved are inconsistent"));
259        Memcpy(INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)), bdims, 2);
260        F77_CALL(dgetrs)("N", &n, &nrhs, REAL(GET_SLOT(lu, Matrix_xSym)), &n,
261                         INTEGER(GET_SLOT(lu, Matrix_permSym)),
262                         Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
263                                REAL(cl ? GET_SLOT(b, Matrix_xSym):b), sz),
264                         &n, &info);
265        UNPROTECT(1);
266        return val;
267    }
268    
269    SEXP dgeMatrix_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
270    {
271        int cl = asLogical(classed);
272        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
273      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
274          *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
275          *cdims,                           getAttrib(b, R_DimSymbol)),
276          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";  
277      double one = 1., zero = 0.;      double one = 1., zero = 0.;
278    
279        if (asLogical(right)) {
280            int m = bdims[0], n = adims[1], k = bdims[1];
281            if (adims[0] != k)
282                error(_("Matrices are not conformable for multiplication"));
283            if (m < 1 || n < 1 || k < 1)
284                error(_("Matrices with zero extents cannot be multiplied"));
285            cdims[0] = m; cdims[1] = n;
286            F77_CALL(dgemm) ("N", "N", &m, &n, &k, &one,
287                             REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &m,
288                             REAL(GET_SLOT(a, Matrix_xSym)), &k, &zero,
289                             REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)),
290                             &m);
291        } else {
292            int m = adims[0], n = bdims[1], k = adims[1];
293    
294      if (bdims[0] != k)      if (bdims[0] != k)
295          error("Matrices are not conformable for multiplication");              error(_("Matrices are not conformable for multiplication"));
296      if (m < 1 || n < 1 || k < 1)      if (m < 1 || n < 1 || k < 1)
297          error("Matrices with zero extents cannot be multiplied");              error(_("Matrices with zero extents cannot be multiplied"));
     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));  
     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));  
298      cdims[0] = m; cdims[1] = n;      cdims[0] = m; cdims[1] = n;
299      F77_CALL(dgemm)(trans, trans, adims, bdims+1, bdims, &one,          F77_CALL(dgemm)
300                      REAL(GET_SLOT(a, Matrix_xSym)), adims,              ("N", "N", &m, &n, &k, &one, REAL(GET_SLOT(a, Matrix_xSym)),
301                      REAL(GET_SLOT(b, Matrix_xSym)), bdims,               &m, REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &k, &zero,
302                      &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);               REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)), &m);
303        }
304      UNPROTECT(1);      UNPROTECT(1);
305      return val;      return val;
306  }  }
307    
308  SEXP geMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)  
309    SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
310  {  {
311      int /* nu = asInteger(nnu),      int /* nu = asInteger(nnu),
312             nv = asInteger(nnv), */             nv = asInteger(nnv), */
# Line 342  Line 371 
371   *   *
372   * @return matrix exponential of x   * @return matrix exponential of x
373   */   */
374  SEXP geMatrix_exp(SEXP x)  SEXP dgeMatrix_exp(SEXP x)
375  {  {
376      SEXP val = PROTECT(duplicate(x));      SEXP val = PROTECT(duplicate(x));
377      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
# Line 359  Line 388 
388          one = 1., trshift, zero = 0.;          one = 1., trshift, zero = 0.;
389    
390      if (nc < 1 || Dims[0] != nc)      if (nc < 1 || Dims[0] != nc)
391          error("Matrix exponential requires square, non-null matrix");          error(_("Matrix exponential requires square, non-null matrix"));
392    
393      /* FIXME: Add special treatment for nc == 1 */      /* FIXME: Add special treatment for nc == 1 */
394    
# Line 373  Line 402 
402    
403      /* Preconditioning 2. Balancing with dgebal. */      /* Preconditioning 2. Balancing with dgebal. */
404      F77_CALL(dgebal)("P", &nc, v, &nc, &ilo, &ihi, perm, &j);      F77_CALL(dgebal)("P", &nc, v, &nc, &ilo, &ihi, perm, &j);
405      if (j) error("geMatrix_exp: LAPACK routine dgebal returned %d", j);      if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
406      F77_CALL(dgebal)("S", &nc, v, &nc, &ilos, &ihis, scale, &j);      F77_CALL(dgebal)("S", &nc, v, &nc, &ilos, &ihis, scale, &j);
407      if (j) error("geMatrix_exp: LAPACK routine dgebal returned %d", j);      if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
408    
409      /* Preconditioning 3. Scaling according to infinity norm */      /* Preconditioning 3. Scaling according to infinity norm */
410      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 442 
442    
443      /* Pade' approximation is solve(dpp, npp) */      /* Pade' approximation is solve(dpp, npp) */
444      F77_CALL(dgetrf)(&nc, &nc, dpp, &nc, pivot, &j);      F77_CALL(dgetrf)(&nc, &nc, dpp, &nc, pivot, &j);
445      if (j) error("geMatrix_exp: dgetrf returned error code %d", j);      if (j) error(_("dgeMatrix_exp: dgetrf returned error code %d"), j);
446      F77_CALL(dgetrs)("N", &nc, &nc, dpp, &nc, pivot, npp, &nc, &j);      F77_CALL(dgetrs)("N", &nc, &nc, dpp, &nc, pivot, npp, &nc, &j);
447      if (j) error("geMatrix_exp: dgetrs returned error code %d", j);      if (j) error(_("dgeMatrix_exp: dgetrs returned error code %d"), j);
448      Memcpy(v, npp, ncsqr);      Memcpy(v, npp, ncsqr);
449    
450      /* Now undo all of the preconditioning */      /* Now undo all of the preconditioning */
# Line 465  Line 494 
494      UNPROTECT(1);      UNPROTECT(1);
495      return val;      return val;
496  }  }
497    
498    SEXP dgeMatrix_Schur(SEXP x, SEXP vectors)
499    {
500        int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
501        int vecs = asLogical(vectors), info, izero = 0, lwork = -1, n = dims[0];
502        double *work, tmp;
503        char *nms[] = {"WR", "WI", "T", "Z", ""};
504        SEXP val = PROTECT(Matrix_make_named(VECSXP, nms));
505    
506        if (n != dims[1] || n < 1)
507            error(_("dgeMatrix_Schur: argument x must be a non-null square matrix"));
508        SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));
509        SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));
510        SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));
511        Memcpy(REAL(VECTOR_ELT(val, 2)), REAL(GET_SLOT(x, Matrix_xSym)), n * n);
512        SET_VECTOR_ELT(val, 3, allocMatrix(REALSXP, vecs ? n : 0, vecs ? n : 0));
513        F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
514                        (double *) NULL, (double *) NULL, (double *) NULL, dims,
515                        &tmp, &lwork, (int *) NULL, &info);
516        if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
517        lwork = (int) tmp;
518        work = Calloc(lwork, double);
519        F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
520                        &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
521                        REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
522                        (int *) NULL, &info);
523        if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
524        Free(work);
525        UNPROTECT(1);
526        return val;
527    }

Legend:
Removed from v.476  
changed lines
  Added in v.692

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