SCM

SCM Repository

[matrix] Annotation of /pkg/src/dense.c
ViewVC logotype

Annotation of /pkg/src/dense.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (view) (download) (as text)

1 : bates 10 #include "dense.h"
2 :    
3 :    
4 :     /**
5 :     * Perform a left cyclic shift of columns j to k in the upper triangular
6 :     * matrix x, then restore it to upper triangular form with Givens rotations.
7 :     * The algorithm is based on the Fortran routine DCHEX from Linpack.
8 :     *
9 :     * The lower triangle of x is not modified.
10 :     *
11 :     * @param x Matrix stored in column-major order
12 :     * @param ldx leading dimension of x
13 :     * @param j column number (0-based) that will be shifted to position k
14 :     * @param k last column number (0-based) to be shifted
15 :     * @param cosines cosines of the Givens rotations
16 :     * @param sines sines of the Givens rotations
17 :     *
18 :     * @return 0 for success
19 :     */
20 :     static
21 :     int left_cyclic(double x[], int ldx, int j, int k,
22 :     double cosines[], double sines[])
23 :     {
24 :     double *lastcol;
25 :     int i, jj;
26 :    
27 :     if (j >= k)
28 :     error("incorrect left cyclic shift, j (%d) >= k (%d)", j, k);
29 :     if (j < 0)
30 :     error("incorrect left cyclic shift, j (%d) < 0", j, k);
31 :     if (ldx < k)
32 :     error("incorrect left cyclic shift, k (%d) > ldx (%d)", k, ldx);
33 :     lastcol = (double*) R_alloc(k+1, sizeof(double));
34 :     /* keep a copy of column j */
35 :     for(i = 0; i <= j; i++) lastcol[i] = x[i + j*ldx];
36 :     /* For safety, zero the rest */
37 :     for(i = j+1; i <= k; i++) lastcol[i] = 0.;
38 :     for(jj = j+1; jj <= k; jj++) { /* columns to be shifted */
39 :     int diagind = jj*(ldx+1), ind = (jj-j) - 1;
40 :     double tmp = x[diagind], cc, ss;
41 :     /* Calculate the Givens rotation. */
42 :     /* This modified the super-diagonal element */
43 :     F77_CALL(drotg)(x + diagind-1, &tmp, cosines + ind, sines + ind);
44 :     cc = cosines[ind]; ss = sines[ind];
45 :     /* Copy column jj+1 to column jj. */
46 :     for(i = 0; i < jj; i++) x[i + (jj-1)*ldx] = x[i+jj*ldx];
47 :     /* Apply rotation to columns up to k */
48 :     for(i = jj; i < k; i++) {
49 :     tmp = cc*x[(jj-1)+i*ldx] + ss*x[jj+i*ldx];
50 :     x[jj+i*ldx] = cc*x[jj+i*ldx] - ss*x[(jj-1)+i*ldx];
51 :     x[(jj-1)+i*ldx] = tmp;
52 :     }
53 :     /* Apply rotation to lastcol */
54 :     lastcol[jj] = -ss*lastcol[jj-1]; lastcol[jj-1] *= cc;
55 :     }
56 :     /* Copy lastcol to column k */
57 :     for(i = 0; i <= k; i++) x[i+k*ldx] = lastcol[i];
58 :     return 0;
59 :     }
60 :    
61 :     static
62 :     SEXP getGivens(double x[], int ldx, int jmin, int rank)
63 :     {
64 :     int shiftlen = (rank - jmin) - 1;
65 :     SEXP ans = PROTECT(allocVector(VECSXP, 4)), nms, cosines, sines;
66 :    
67 :     SET_VECTOR_ELT(ans, 0, ScalarInteger(jmin));
68 :     SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
69 :     SET_VECTOR_ELT(ans, 2, cosines = allocVector(REALSXP, shiftlen));
70 :     SET_VECTOR_ELT(ans, 3, sines = allocVector(REALSXP, shiftlen));
71 :     setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 4));
72 :     SET_STRING_ELT(nms, 0, mkChar("jmin"));
73 :     SET_STRING_ELT(nms, 1, mkChar("rank"));
74 :     SET_STRING_ELT(nms, 2, mkChar("cosines"));
75 :     SET_STRING_ELT(nms, 3, mkChar("sines"));
76 :     if (left_cyclic(x, ldx, jmin, rank - 1, REAL(cosines), REAL(sines)))
77 :     error("Unknown error in getGivens");
78 :     UNPROTECT(1);
79 :     return ans;
80 :     }
81 :    
82 :     SEXP checkGivens(SEXP X, SEXP jmin, SEXP rank)
83 :     {
84 :     SEXP ans = PROTECT(allocVector(VECSXP, 2)),
85 :     Xcp = PROTECT(duplicate(X));
86 :     int *Xdims;
87 :    
88 :     if (!(isReal(X) & isMatrix(X)))
89 :     error("X must be a numeric (double precision) matrix");
90 :     Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
91 :     SET_VECTOR_ELT(ans, 1, getGivens(REAL(Xcp), Xdims[0],
92 :     asInteger(jmin), asInteger(rank)));
93 :     SET_VECTOR_ELT(ans, 0, Xcp);
94 :     UNPROTECT(2);
95 :     return ans;
96 :     }
97 :    
98 :     SEXP lsq_dense_Chol(SEXP X, SEXP y)
99 :     {
100 :     SEXP ans;
101 :     int info, n, p, k, *Xdims, *ydims;
102 :     double *xpx, d_one = 1., d_zero = 0.;
103 :    
104 :     if (!(isReal(X) & isMatrix(X)))
105 :     error("X must be a numeric (double precision) matrix");
106 :     Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
107 :     n = Xdims[0];
108 :     p = Xdims[1];
109 :     if (!(isReal(y) & isMatrix(y)))
110 :     error("y must be a numeric (double precision) matrix");
111 :     ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
112 :     if (ydims[0] != n)
113 :     error(
114 :     "number of rows in y (%d) does not match number of rows in X (%d)",
115 :     ydims[0], n);
116 :     k = ydims[1];
117 :     if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
118 :     ans = PROTECT(allocMatrix(REALSXP, p, k));
119 :     F77_CALL(dgemm)("T", "N", &p, &k, &n, &d_one, REAL(X), &n, REAL(y), &n,
120 :     &d_zero, REAL(ans), &p);
121 :     xpx = (double *) R_alloc(p * p, sizeof(double));
122 :     F77_CALL(dsyrk)("U", "T", &p, &n, &d_one, REAL(X), &n, &d_zero,
123 :     xpx, &p);
124 :     F77_CALL(dposv)("U", &p, &k, xpx, &p, REAL(ans), &p, &info);
125 :     if (info) error("Lapack routine dposv returned error code %d", info);
126 :     UNPROTECT(1);
127 :     return ans;
128 :     }
129 :    
130 :    
131 :     SEXP lsq_dense_QR(SEXP X, SEXP y)
132 :     {
133 :     SEXP ans;
134 :     int info, n, p, k, *Xdims, *ydims, lwork;
135 :     double *work, tmp, *xvals;
136 :    
137 :     if (!(isReal(X) & isMatrix(X)))
138 :     error("X must be a numeric (double precision) matrix");
139 :     Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
140 :     n = Xdims[0];
141 :     p = Xdims[1];
142 :     if (!(isReal(y) & isMatrix(y)))
143 :     error("y must be a numeric (double precision) matrix");
144 :     ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
145 :     if (ydims[0] != n)
146 :     error(
147 :     "number of rows in y (%d) does not match number of rows in X (%d)",
148 :     ydims[0], n);
149 :     k = ydims[1];
150 :     if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
151 :     xvals = (double *) R_alloc(n * p, sizeof(double));
152 :     Memcpy(xvals, REAL(X), n * p);
153 :     ans = PROTECT(duplicate(y));
154 :     lwork = -1;
155 :     F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
156 :     &tmp, &lwork, &info);
157 :     if (info)
158 :     error("First call to Lapack routine dgels returned error code %d",
159 :     info);
160 :     lwork = (int) tmp;
161 :     work = (double *) R_alloc(lwork, sizeof(double));
162 :     F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
163 :     work, &lwork, &info);
164 :     if (info)
165 :     error("Second call to Lapack routine dgels returned error code %d",
166 :     info);
167 :     UNPROTECT(1);
168 :     return ans;
169 :     }
170 :    
171 :     SEXP lapack_qr(SEXP Xin, SEXP tl)
172 :     {
173 :     SEXP ans, Givens, Gcpy, nms, pivot, qraux, X;
174 :     int i, n, nGivens = 0, p, trsz, *Xdims, rank;
175 :     double rcond = 0., tol = asReal(tl), *work;
176 :    
177 :     if (!(isReal(Xin) & isMatrix(Xin)))
178 :     error("X must be a real (numeric) matrix");
179 :     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);
181 :     ans = PROTECT(allocVector(VECSXP,5));
182 :     SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));
183 :     Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
184 :     n = Xdims[0]; p = Xdims[1];
185 :     SET_VECTOR_ELT(ans, 2, qraux = allocVector(REALSXP, (n < p) ? n : p));
186 :     SET_VECTOR_ELT(ans, 3, pivot = allocVector(INTSXP, p));
187 :     for (i = 0; i < p; i++) INTEGER(pivot)[i] = i + 1;
188 :     trsz = (n < p) ? n : p; /* size of triangular part of decomposition */
189 :     rank = trsz;
190 :     Givens = PROTECT(allocVector(VECSXP, rank - 1));
191 :     setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 5));
192 :     SET_STRING_ELT(nms, 0, mkChar("qr"));
193 :     SET_STRING_ELT(nms, 1, mkChar("rank"));
194 :     SET_STRING_ELT(nms, 2, mkChar("qraux"));
195 :     SET_STRING_ELT(nms, 3, mkChar("pivot"));
196 :     SET_STRING_ELT(nms, 4, mkChar("Givens"));
197 :     if (n > 0 && p > 0) {
198 :     int info, *iwork, lwork;
199 :     double *xpt = REAL(X), tmp;
200 :    
201 :     lwork = -1;
202 :     F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork, &info);
203 :     if (info)
204 :     error("First call to dgeqrf returned error code %d", info);
205 :     lwork = (int) tmp;
206 :     work = (double *) R_alloc((lwork < 3*trsz) ? 3*trsz : lwork,
207 :     sizeof(double));
208 :     F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork, &info);
209 :     if (info)
210 :     error("Second call to dgeqrf returned error code %d", info);
211 :     iwork = (int *) R_alloc(trsz, sizeof(double));
212 :     F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
213 :     work, iwork, &info);
214 :     if (info)
215 :     error("Lapack routine dtrcon returned error code %d", info);
216 :     while (rcond < tol) { /* check diagonal elements */
217 :     double minabs = (xpt[0] < 0.) ? -xpt[0]: xpt[0];
218 :     int jmin = 0;
219 :     for (i = 1; i < rank; i++) {
220 :     double el = xpt[i*(n+1)];
221 :     el = (el < 0.) ? -el: el;
222 :     if (el < minabs) {
223 :     jmin = i;
224 :     minabs = el;
225 :     }
226 :     }
227 :     if (jmin < (rank - 1)) {
228 :     SET_VECTOR_ELT(Givens, nGivens, getGivens(xpt, n, jmin, rank));
229 :     nGivens++;
230 :     }
231 :     rank--;
232 :     F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
233 :     work, iwork, &info);
234 :     if (info)
235 :     error("Lapack routine dtrcon returned error code %d", info);
236 :     }
237 :     }
238 :     SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));
239 :     for (i = 0; i < nGivens; i++)
240 :     SET_VECTOR_ELT(Gcpy, i, VECTOR_ELT(Givens, i));
241 :     SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
242 :     setAttrib(ans, install("useLAPACK"), ScalarLogical(1));
243 :     setAttrib(ans, install("rcond"), ScalarReal(rcond));
244 :     UNPROTECT(2);
245 :     return ans;
246 :     }

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