SCM Repository

[matrix] Diff of /pkg/src/dense.c
 [matrix] / pkg / src / dense.c

Diff of /pkg/src/dense.c

revision 357, Sat Nov 27 17:56:13 2004 UTC revision 582, Mon Feb 28 18:15:21 2005 UTC
# Line 25  Line 25
25      int i, jj;      int i, jj;
26
27      if (j >= k)      if (j >= k)
28          error("incorrect left cyclic shift, j (%d) >= k (%d)", j, k);          error(_("incorrect left cyclic shift, j (%d) >= k (%d)"), j, k);
29      if (j < 0)      if (j < 0)
30          error("incorrect left cyclic shift, j (%d) < 0", j, k);          error(_("incorrect left cyclic shift, j (%d) < 0"), j, k);
31      if (ldx < k)      if (ldx < k)
32          error("incorrect left cyclic shift, k (%d) > ldx (%d)", k, ldx);          error(_("incorrect left cyclic shift, k (%d) > ldx (%d)"), k, ldx);
33      lastcol = (double*) R_alloc(k+1, sizeof(double));      lastcol = (double*) R_alloc(k+1, sizeof(double));
34                                  /* keep a copy of column j */                                  /* keep a copy of column j */
35      for(i = 0; i <= j; i++) lastcol[i] = x[i + j*ldx];      for(i = 0; i <= j; i++) lastcol[i] = x[i + j*ldx];
# Line 74  Line 74
74      SET_STRING_ELT(nms, 2, mkChar("cosines"));      SET_STRING_ELT(nms, 2, mkChar("cosines"));
75      SET_STRING_ELT(nms, 3, mkChar("sines"));      SET_STRING_ELT(nms, 3, mkChar("sines"));
76      if (left_cyclic(x, ldx, jmin, rank - 1, REAL(cosines), REAL(sines)))      if (left_cyclic(x, ldx, jmin, rank - 1, REAL(cosines), REAL(sines)))
77          error("Unknown error in getGivens");          error(_("Unknown error in getGivens"));
78      UNPROTECT(1);      UNPROTECT(1);
79      return ans;      return ans;
80  }  }
# Line 86  Line 86
86      int  *Xdims;      int  *Xdims;
87
88      if (!(isReal(X) & isMatrix(X)))      if (!(isReal(X) & isMatrix(X)))
89          error("X must be a numeric (double precision) matrix");          error(_("X must be a numeric (double precision) matrix"));
90      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
91      SET_VECTOR_ELT(ans, 1, getGivens(REAL(Xcp), Xdims[0],      SET_VECTOR_ELT(ans, 1, getGivens(REAL(Xcp), Xdims[0],
92                                       asInteger(jmin), asInteger(rank)));                                       asInteger(jmin), asInteger(rank)));
# Line 102  Line 102
102      double *xpx, d_one = 1., d_zero = 0.;      double *xpx, d_one = 1., d_zero = 0.;
103
104      if (!(isReal(X) & isMatrix(X)))      if (!(isReal(X) & isMatrix(X)))
105          error("X must be a numeric (double precision) matrix");          error(_("X must be a numeric (double precision) matrix"));
106      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
107      n = Xdims[0];      n = Xdims[0];
108      p = Xdims[1];      p = Xdims[1];
109      if (!(isReal(y) & isMatrix(y)))      if (!(isReal(y) & isMatrix(y)))
110          error("y must be a numeric (double precision) matrix");          error(_("y must be a numeric (double precision) matrix"));
111      ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));      ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
112      if (ydims[0] != n)      if (ydims[0] != n)
113          error(          error(_(
114              "number of rows in y (%d) does not match number of rows in X (%d)",              "number of rows in y (%d) does not match number of rows in X (%d)"),
115              ydims[0], n);              ydims[0], n);
116      k = ydims[1];      k = ydims[1];
117      if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);      if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
# Line 122  Line 122
122      F77_CALL(dsyrk)("U", "T", &p, &n, &d_one, REAL(X), &n, &d_zero,      F77_CALL(dsyrk)("U", "T", &p, &n, &d_one, REAL(X), &n, &d_zero,
123                      xpx, &p);                      xpx, &p);
124      F77_CALL(dposv)("U", &p, &k, xpx, &p, REAL(ans), &p, &info);      F77_CALL(dposv)("U", &p, &k, xpx, &p, REAL(ans), &p, &info);
125      if (info) error("Lapack routine dposv returned error code %d", info);      if (info) error(_("Lapack routine dposv returned error code %d"), info);
126      UNPROTECT(1);      UNPROTECT(1);
127      return ans;      return ans;
128  }  }
# Line 135  Line 135
135      double *work, tmp, *xvals;      double *work, tmp, *xvals;
136
137      if (!(isReal(X) & isMatrix(X)))      if (!(isReal(X) & isMatrix(X)))
138          error("X must be a numeric (double precision) matrix");          error(_("X must be a numeric (double precision) matrix"));
139      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
140      n = Xdims[0];      n = Xdims[0];
141      p = Xdims[1];      p = Xdims[1];
142      if (!(isReal(y) & isMatrix(y)))      if (!(isReal(y) & isMatrix(y)))
143          error("y must be a numeric (double precision) matrix");          error(_("y must be a numeric (double precision) matrix"));
144      ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));      ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
145      if (ydims[0] != n)      if (ydims[0] != n)
146          error(          error(_(
147              "number of rows in y (%d) does not match number of rows in X (%d)",              "number of rows in y (%d) does not match number of rows in X (%d)"),
148              ydims[0], n);              ydims[0], n);
149      k = ydims[1];      k = ydims[1];
150      if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);      if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
# Line 155  Line 155
155      F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,      F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
156                      &tmp, &lwork, &info);                      &tmp, &lwork, &info);
157      if (info)      if (info)
158          error("First call to Lapack routine dgels returned error code %d",          error(_("First call to Lapack routine dgels returned error code %d"),
159                info);                info);
160      lwork = (int) tmp;      lwork = (int) tmp;
161      work = (double *) R_alloc(lwork, sizeof(double));      work = (double *) R_alloc(lwork, sizeof(double));
162      F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,      F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
163                      work, &lwork, &info);                      work, &lwork, &info);
164      if (info)      if (info)
165          error("Second call to Lapack routine dgels returned error code %d",          error(_("Second call to Lapack routine dgels returned error code %d"),
166                info);                info);
167      UNPROTECT(1);      UNPROTECT(1);
168      return ans;      return ans;
# Line 175  Line 175
175      double rcond = 0., tol = asReal(tl), *work;      double rcond = 0., tol = asReal(tl), *work;
176
177      if (!(isReal(Xin) & isMatrix(Xin)))      if (!(isReal(Xin) & isMatrix(Xin)))
178          error("X must be a real (numeric) matrix");          error(_("X must be a real (numeric) matrix"));
179      if (tol < 0.) error("tol, given as %g, must be non-negative", tol);      if (tol < 0.) error(_("tol, given as %g, must be non-negative"), tol);
180      if (tol > 1.) error("tol, given as %g, must be <= 1", tol);      if (tol > 1.) error(_("tol, given as %g, must be <= 1"), tol);
181      ans = PROTECT(allocVector(VECSXP,5));      ans = PROTECT(allocVector(VECSXP,5));
182      SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));      SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));
183      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));      Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
# Line 201  Line 201
201          lwork = -1;          lwork = -1;
202          F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork, &info);          F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork, &info);
203          if (info)          if (info)
204              error("First call to dgeqrf returned error code %d", info);              error(_("First call to dgeqrf returned error code %d"), info);
205          lwork = (int) tmp;          lwork = (int) tmp;
206          work = (double *) R_alloc((lwork < 3*trsz) ? 3*trsz : lwork,          work = (double *) R_alloc((lwork < 3*trsz) ? 3*trsz : lwork,
207                                    sizeof(double));                                    sizeof(double));
208          F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork, &info);          F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork, &info);
209          if (info)          if (info)
210              error("Second call to dgeqrf returned error code %d", info);              error(_("Second call to dgeqrf returned error code %d"), info);
211          iwork = (int *) R_alloc(trsz, sizeof(int));          iwork = (int *) R_alloc(trsz, sizeof(int));
212          F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,          F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
213                           work, iwork, &info);                           work, iwork, &info);
214          if (info)          if (info)
215              error("Lapack routine dtrcon returned error code %d", info);              error(_("Lapack routine dtrcon returned error code %d"), info);
216          while (rcond < tol) {   /* check diagonal elements */          while (rcond < tol) {   /* check diagonal elements */
217              double minabs = (xpt[0] < 0.) ? -xpt[0]: xpt[0];              double minabs = (xpt[0] < 0.) ? -xpt[0]: xpt[0];
218              int jmin = 0;              int jmin = 0;
# Line 232  Line 232
232              F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,              F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
233                               work, iwork, &info);                               work, iwork, &info);
234              if (info)              if (info)
235                  error("Lapack routine dtrcon returned error code %d", info);                  error(_("Lapack routine dtrcon returned error code %d"), info);
236          }          }
237      }      }
238      SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));      SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));

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