SCM

SCM Repository

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

Annotation of /pkg/src/dtrMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : maechler 534 /* double (precision) TRiangular Matrices */
2 :    
3 : bates 478 #include "dtrMatrix.h"
4 : bates 10
5 : maechler 890 SEXP triangularMatrix_validate(SEXP obj)
6 : bates 10 {
7 : maechler 890 SEXP val = GET_SLOT(obj, Matrix_DimSym);
8 : maechler 534
9 : maechler 890 if (LENGTH(val) < 2)
10 :     return mkString(_("'Dim' slot has length less than two"));
11 :     if (INTEGER(val)[0] != INTEGER(val)[1])
12 : bates 846 return mkString(_("Matrix is not square"));
13 : bates 592 if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
14 :     "LU", "uplo"))) return val;
15 :     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),
16 :     "NU", "diag"))) return val;
17 : bates 10 return ScalarLogical(1);
18 :     }
19 :    
20 : maechler 890 SEXP dtrMatrix_validate(SEXP obj)
21 :     {
22 : maechler 1164 /* since "dtr" inherits from "triangular", and "dMatrix", only need this:*/
23 :     return dense_nonpacked_validate(obj);
24 : maechler 890 }
25 :    
26 :    
27 : bates 10 static
28 :     double get_norm(SEXP obj, char *typstr)
29 :     {
30 :     char typnm[] = {'\0', '\0'};
31 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
32 :     double *work = (double *) NULL;
33 :    
34 :     typnm[0] = norm_type(typstr);
35 :     if (*typnm == 'I') {
36 :     work = (double *) R_alloc(dims[0], sizeof(double));
37 :     }
38 : maechler 951 return F77_CALL(dlantr)(typnm, uplo_P(obj), diag_P(obj), dims, dims+1,
39 :     REAL(GET_SLOT(obj, Matrix_xSym)), dims, work);
40 : bates 10 }
41 :    
42 :    
43 : bates 478 SEXP dtrMatrix_norm(SEXP obj, SEXP type)
44 : bates 10 {
45 :     return ScalarReal(get_norm(obj, CHAR(asChar(type))));
46 :     }
47 :    
48 :     static
49 :     double set_rcond(SEXP obj, char *typstr)
50 :     {
51 :     char typnm[] = {'\0', '\0'};
52 : bates 342 SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
53 : bates 10 double rcond = get_double_by_name(rcv, typnm);
54 :    
55 :     typnm[0] = rcond_type(typstr);
56 :     if (R_IsNA(rcond)) {
57 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
58 : maechler 951 F77_CALL(dtrcon)(typnm, uplo_P(obj), diag_P(obj), dims,
59 :     REAL(GET_SLOT(obj, Matrix_xSym)), dims, &rcond,
60 : bates 10 (double *) R_alloc(3*dims[0], sizeof(double)),
61 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
62 : bates 342 SET_SLOT(obj, Matrix_rcondSym,
63 : bates 10 set_double_by_name(rcv, rcond, typnm));
64 :     }
65 :     return rcond;
66 :     }
67 :    
68 : bates 478 SEXP dtrMatrix_rcond(SEXP obj, SEXP type)
69 : bates 10 {
70 :     return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
71 :     }
72 :    
73 : bates 478 SEXP dtrMatrix_solve(SEXP a)
74 : bates 10 {
75 :     SEXP val = PROTECT(duplicate(a));
76 :     int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));
77 : maechler 951 F77_CALL(dtrtri)(uplo_P(val), diag_P(val), Dim,
78 :     REAL(GET_SLOT(val, Matrix_xSym)), Dim, &info);
79 : bates 10 UNPROTECT(1);
80 :     return val;
81 :     }
82 :    
83 : bates 654 SEXP dtrMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
84 : maechler 655 {
85 : bates 654 int cl = asLogical(classed);
86 : bates 846 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
87 : bates 654 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
88 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
89 :     getAttrib(b, R_DimSymbol));
90 : maechler 655 int n = bdims[0], nrhs = bdims[1];
91 : bates 654 int sz = n * nrhs;
92 : bates 296 double one = 1.0;
93 :    
94 : bates 654 if (*adims != *bdims || bdims[1] < 1 || *adims < 1 || *adims != adims[1])
95 :     error(_("Dimensions of system to be solved are inconsistent"));
96 : bates 846 Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
97 : maechler 951 F77_CALL(dtrsm)("L", uplo_P(a),
98 :     "N", diag_P(a),
99 : bates 654 &n, &nrhs, &one, REAL(GET_SLOT(a, Matrix_xSym)), &n,
100 : bates 846 Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, sz)),
101 : bates 654 REAL(cl ? GET_SLOT(b, Matrix_xSym):b), sz), &n);
102 : bates 296 UNPROTECT(1);
103 : bates 846 return ans;
104 : bates 296 }
105 :    
106 : bates 654 SEXP dtrMatrix_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
107 :     {
108 :     int cl = asLogical(classed), rt = asLogical(right);
109 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
110 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
111 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
112 :     getAttrib(b, R_DimSymbol)),
113 :     *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
114 :     int m, n, sz;
115 :     double one = 1.;
116 :    
117 :     if (!cl && !(isReal(b) && isMatrix(b)))
118 :     error(_("Argument b must be a numeric matrix"));
119 : maechler 688 if (adims[0] != adims[1]) error(_("dtrMatrix in %*% must be square"));
120 : bates 654 m = rt ? bdims[0] : adims[0];
121 :     n = rt ? adims[1] : bdims[1];
122 :     if ((rt && (adims[0] != m)) || (!rt && (bdims[0] != m)))
123 :     error(_("Matrices are not conformable for multiplication"));
124 :     if (m < 1 || n < 1)
125 :     error(_("Matrices with zero extents cannot be multiplied"));
126 :     cdims[0] = m; cdims[1] = n; sz = m * n;
127 : maechler 951 F77_CALL(dtrmm)(rt ? "R" : "L", uplo_P(a),
128 :     "N", diag_P(a), &m, &n,
129 : bates 654 &one, REAL(GET_SLOT(a, Matrix_xSym)), rt ? &n : &m,
130 :     Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
131 :     REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), sz),
132 :     rt ? &m : &n);
133 :     UNPROTECT(1);
134 :     return val;
135 :     }
136 :    
137 : bates 478 SEXP dtrMatrix_as_dgeMatrix(SEXP from)
138 : bates 10 {
139 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
140 : maechler 534
141 : maechler 628 SET_SLOT(val, Matrix_rcondSym, duplicate(GET_SLOT(from, Matrix_rcondSym)));
142 : bates 296 SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(from, Matrix_xSym)));
143 : maechler 628 SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(from, Matrix_DimSym)));
144 :     SET_SLOT(val, Matrix_DimNamesSym, duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
145 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
146 : bates 296 make_array_triangular(REAL(GET_SLOT(val, Matrix_xSym)), from);
147 : bates 10 UNPROTECT(1);
148 :     return val;
149 :     }
150 :    
151 : bates 478 SEXP dtrMatrix_as_matrix(SEXP from)
152 : bates 10 {
153 :     int *Dim = INTEGER(GET_SLOT(from, Matrix_DimSym));
154 :     int m = Dim[0], n = Dim[1];
155 :     SEXP val = PROTECT(allocMatrix(REALSXP, m, n));
156 :     make_array_triangular(Memcpy(REAL(val),
157 : bates 296 REAL(GET_SLOT(from, Matrix_xSym)), m * n),
158 : bates 10 from);
159 : maechler 628 setAttrib(val, R_DimNamesSymbol, GET_SLOT(from, Matrix_DimNamesSym));
160 : bates 10 UNPROTECT(1);
161 :     return val;
162 :     }
163 :    
164 : bates 478 SEXP dtrMatrix_getDiag(SEXP x)
165 : bates 10 {
166 :     int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
167 :     SEXP ret = PROTECT(allocVector(REALSXP, n)),
168 : bates 296 xv = GET_SLOT(x, Matrix_xSym);
169 : bates 10
170 : maechler 951 if ('U' == diag_P(x)[0]) {
171 : bates 10 for (i = 0; i < n; i++) REAL(ret)[i] = 1.;
172 :     } else {
173 :     for (i = 0; i < n; i++) {
174 :     REAL(ret)[i] = REAL(xv)[i * (n + 1)];
175 :     }
176 :     }
177 :     UNPROTECT(1);
178 :     return ret;
179 :     }
180 : maechler 534
181 : bates 478 SEXP dtrMatrix_dgeMatrix_mm_R(SEXP a, SEXP b)
182 : bates 10 {
183 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
184 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
185 :     m = adims[0], n = bdims[1], k = adims[1];
186 :     SEXP val = PROTECT(duplicate(b));
187 :     double one = 1.;
188 : maechler 534
189 : bates 10 if (bdims[0] != k)
190 : bates 582 error(_("Matrices are not conformable for multiplication"));
191 : bates 10 if (m < 1 || n < 1 || k < 1)
192 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
193 : maechler 951 F77_CALL(dtrmm)("R", uplo_P(a), "N", diag_P(a), adims, bdims+1, &one,
194 : bates 10 REAL(GET_SLOT(a, Matrix_xSym)), adims,
195 :     REAL(GET_SLOT(val, Matrix_xSym)), bdims);
196 :     UNPROTECT(1);
197 :     return val;
198 :     }
199 : bates 597
200 :     SEXP dtrMatrix_as_dtpMatrix(SEXP from)
201 :     {
202 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtpMatrix"))),
203 :     uplo = GET_SLOT(from, Matrix_uploSym),
204 :     diag = GET_SLOT(from, Matrix_diagSym),
205 :     dimP = GET_SLOT(from, Matrix_DimSym);
206 :     int n = *INTEGER(dimP);
207 :    
208 :     SET_SLOT(val, Matrix_rcondSym,
209 :     duplicate(GET_SLOT(from, Matrix_rcondSym)));
210 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
211 :     SET_SLOT(val, Matrix_diagSym, duplicate(diag));
212 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
213 : maechler 952 full_to_packed_double(
214 :     REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
215 :     REAL(GET_SLOT(from, Matrix_xSym)), n,
216 :     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW,
217 :     *CHAR(STRING_ELT(diag, 0)) == 'U' ? UNT : NUN);
218 : bates 597 UNPROTECT(1);
219 :     return val;
220 :     }

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