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 534, Tue Feb 8 08:59:31 2005 UTC revision 582, Mon Feb 28 18:15:21 2005 UTC
# Line 9  Line 9 
9      int m, n;      int m, n;
10    
11      if (length(Dim) != 2)      if (length(Dim) != 2)
12          return mkString("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 mkString("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 mkString("length of x slot != prod(Dim)");          return mkString(_("length of x slot != prod(Dim)"));
18      if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)      if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)
19          return mkString("factors slot must be named list");          return mkString(_("factors slot must be named list"));
20      if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)      if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)
21          return mkString("rcond slot must be named numeric vector");          return mkString(_("rcond slot must be named numeric vector"));
22      return ScalarLogical(1);      return ScalarLogical(1);
23  }  }
24    
# Line 57  Line 57 
57          double anorm = get_norm(obj, typstr);          double anorm = get_norm(obj, typstr);
58    
59          if (dims[0] != dims[1] || dims[0] < 1)          if (dims[0] != dims[1] || dims[0] < 1)
60              error("rcond requires a square, non-empty matrix");              error(_("rcond requires a square, non-empty matrix"));
61          F77_CALL(dgecon)(typnm,          F77_CALL(dgecon)(typnm,
62                           dims, REAL(GET_SLOT(LU, Matrix_xSym)),                           dims, REAL(GET_SLOT(LU, Matrix_xSym)),
63                           dims, &anorm, &rcond,                           dims, &anorm, &rcond,
# Line 114  Line 114 
114      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
115      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
116          if (*xDims != *yDims)          if (*xDims != *yDims)
117              error("Dimensions of x and y are not compatible for crossprod");              error(_("Dimensions of x and y are not compatible for crossprod"));
118          vDims[0] = m; vDims[1] = n;          vDims[0] = m; vDims[1] = n;
119          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
120          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
# Line 137  Line 137 
137      double one = 1.0, zero = 0.0;      double one = 1.0, zero = 0.0;
138    
139      if (!(isMatrix(y) && isReal(y)))      if (!(isMatrix(y) && isReal(y)))
140          error("Argument y must be a numeric (real) matrix");          error(_("Argument y must be a numeric (real) matrix"));
141      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
142      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
143      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
144      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));      vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
145      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {      if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
146          if (*xDims != *yDims)          if (*xDims != *yDims)
147              error("Dimensions of x and y are not compatible for crossprod");              error(_("Dimensions of x and y are not compatible for crossprod"));
148          vDims[0] = m; vDims[1] = n;          vDims[0] = m; vDims[1] = n;
149          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));          SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
150          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,          F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
# Line 179  Line 179 
179      if (val != R_NilValue) return val;      if (val != R_NilValue) return val;
180      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
181      if (dims[0] < 1 || dims[1] < 1)      if (dims[0] < 1 || dims[1] < 1)
182          error("Cannot factor a matrix with zero extents");          error(_("Cannot factor a matrix with zero extents"));
183      npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];      npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];
184      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));      val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
185      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
# Line 188  Line 188 
188      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),      F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
189                       dims, INTEGER(GET_SLOT(val, install("pivot"))),                       dims, INTEGER(GET_SLOT(val, install("pivot"))),
190                       &info);                       &info);
191      if (info) error("Lapack routine dgetrf returned error code %d", info);      if (info) error(_("Lapack routine dgetrf returned error code %d"), info);
192      UNPROTECT(1);      UNPROTECT(1);
193      return set_factors(x, val, "LU");      return set_factors(x, val, "LU");
194  }  }
# Line 203  Line 203 
203      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;      double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;
204    
205      if (n != dims[1])      if (n != dims[1])
206          error("Determinant requires a square matrix");          error(_("Determinant requires a square matrix"));
207      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;
208      if (lg) {      if (lg) {
209          modulus = 0.0;          modulus = 0.0;
# Line 234  Line 234 
234      int info, lwork = -1;      int info, lwork = -1;
235    
236    
237      if (dims[0] != dims[1]) error("Solve requires a square matrix");      if (dims[0] != dims[1]) error(_("Solve requires a square matrix"));
238      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
239      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
240      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
# Line 260  Line 260 
260      double one = 1., zero = 0.;      double one = 1., zero = 0.;
261    
262      if (bdims[0] != k)      if (bdims[0] != k)
263          error("Matrices are not conformable for multiplication");          error(_("Matrices are not conformable for multiplication"));
264      if (m < 1 || n < 1 || k < 1)      if (m < 1 || n < 1 || k < 1)
265          error("Matrices with zero extents cannot be multiplied");          error(_("Matrices with zero extents cannot be multiplied"));
266      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));      SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
267      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
268      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
# Line 359  Line 359 
359          one = 1., trshift, zero = 0.;          one = 1., trshift, zero = 0.;
360    
361      if (nc < 1 || Dims[0] != nc)      if (nc < 1 || Dims[0] != nc)
362          error("Matrix exponential requires square, non-null matrix");          error(_("Matrix exponential requires square, non-null matrix"));
363    
364      /* FIXME: Add special treatment for nc == 1 */      /* FIXME: Add special treatment for nc == 1 */
365    
# 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("dgeMatrix_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("dgeMatrix_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("dgeMatrix_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("dgeMatrix_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 */
# Line 475  Line 475 
475      SEXP val = PROTECT(Matrix_make_named(VECSXP, nms));      SEXP val = PROTECT(Matrix_make_named(VECSXP, nms));
476    
477      if (n != dims[1] || n < 1)      if (n != dims[1] || n < 1)
478          error("dgeMatrix_Schur: argument x must be a non-null square matrix");          error(_("dgeMatrix_Schur: argument x must be a non-null square matrix"));
479      SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));      SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));
480      SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));      SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));
481      SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));      SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));
# Line 484  Line 484 
484      F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,      F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
485                      (double *) NULL, (double *) NULL, (double *) NULL, dims,                      (double *) NULL, (double *) NULL, (double *) NULL, dims,
486                      &tmp, &lwork, (int *) NULL, &info);                      &tmp, &lwork, (int *) NULL, &info);
487      if (info) error("dgeMatrix_Schur: first call to dgees failed");      if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
488      lwork = (int) tmp;      lwork = (int) tmp;
489      work = Calloc(lwork, double);      work = Calloc(lwork, double);
490      F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,      F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
491                      &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),                      &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
492                      REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,                      REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
493                      (int *) NULL, &info);                      (int *) NULL, &info);
494      if (info) error("dgeMatrix_Schur: dgees returned code %d", info);      if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
495      Free(work);      Free(work);
496      UNPROTECT(1);      UNPROTECT(1);
497      return val;      return val;

Legend:
Removed from v.534  
changed lines
  Added in v.582

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