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 1167 - (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 : bates 1167 SEXP dtrMatrix_rcond(SEXP obj, SEXP type)
49 : bates 10 {
50 :     char typnm[] = {'\0', '\0'};
51 : bates 1167 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
52 :     double rcond;
53 : bates 10
54 : bates 1167 typnm[0] = rcond_type(CHAR(asChar(type)));
55 :     F77_CALL(dtrcon)(typnm, uplo_P(obj), diag_P(obj), dims,
56 :     REAL(GET_SLOT(obj, Matrix_xSym)), dims, &rcond,
57 :     (double *) R_alloc(3*dims[0], sizeof(double)),
58 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
59 :     return ScalarReal(rcond);
60 : bates 10 }
61 :    
62 : bates 478 SEXP dtrMatrix_solve(SEXP a)
63 : bates 10 {
64 :     SEXP val = PROTECT(duplicate(a));
65 :     int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));
66 : maechler 951 F77_CALL(dtrtri)(uplo_P(val), diag_P(val), Dim,
67 :     REAL(GET_SLOT(val, Matrix_xSym)), Dim, &info);
68 : bates 10 UNPROTECT(1);
69 :     return val;
70 :     }
71 :    
72 : bates 654 SEXP dtrMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
73 : maechler 655 {
74 : bates 654 int cl = asLogical(classed);
75 : bates 846 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
76 : bates 654 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
77 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
78 :     getAttrib(b, R_DimSymbol));
79 : maechler 655 int n = bdims[0], nrhs = bdims[1];
80 : bates 654 int sz = n * nrhs;
81 : bates 296 double one = 1.0;
82 :    
83 : bates 654 if (*adims != *bdims || bdims[1] < 1 || *adims < 1 || *adims != adims[1])
84 :     error(_("Dimensions of system to be solved are inconsistent"));
85 : bates 846 Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
86 : maechler 951 F77_CALL(dtrsm)("L", uplo_P(a),
87 :     "N", diag_P(a),
88 : bates 654 &n, &nrhs, &one, REAL(GET_SLOT(a, Matrix_xSym)), &n,
89 : bates 846 Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, sz)),
90 : bates 654 REAL(cl ? GET_SLOT(b, Matrix_xSym):b), sz), &n);
91 : bates 296 UNPROTECT(1);
92 : bates 846 return ans;
93 : bates 296 }
94 :    
95 : bates 654 SEXP dtrMatrix_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
96 :     {
97 :     int cl = asLogical(classed), rt = asLogical(right);
98 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
99 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
100 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
101 :     getAttrib(b, R_DimSymbol)),
102 :     *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
103 :     int m, n, sz;
104 :     double one = 1.;
105 :    
106 :     if (!cl && !(isReal(b) && isMatrix(b)))
107 :     error(_("Argument b must be a numeric matrix"));
108 : maechler 688 if (adims[0] != adims[1]) error(_("dtrMatrix in %*% must be square"));
109 : bates 654 m = rt ? bdims[0] : adims[0];
110 :     n = rt ? adims[1] : bdims[1];
111 :     if ((rt && (adims[0] != m)) || (!rt && (bdims[0] != m)))
112 :     error(_("Matrices are not conformable for multiplication"));
113 :     if (m < 1 || n < 1)
114 :     error(_("Matrices with zero extents cannot be multiplied"));
115 :     cdims[0] = m; cdims[1] = n; sz = m * n;
116 : maechler 951 F77_CALL(dtrmm)(rt ? "R" : "L", uplo_P(a),
117 :     "N", diag_P(a), &m, &n,
118 : bates 654 &one, REAL(GET_SLOT(a, Matrix_xSym)), rt ? &n : &m,
119 :     Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
120 :     REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), sz),
121 :     rt ? &m : &n);
122 :     UNPROTECT(1);
123 :     return val;
124 :     }
125 :    
126 : bates 478 SEXP dtrMatrix_as_dgeMatrix(SEXP from)
127 : bates 10 {
128 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
129 : maechler 534
130 : bates 296 SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(from, Matrix_xSym)));
131 : maechler 628 SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(from, Matrix_DimSym)));
132 :     SET_SLOT(val, Matrix_DimNamesSym, duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
133 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
134 : bates 296 make_array_triangular(REAL(GET_SLOT(val, Matrix_xSym)), from);
135 : bates 10 UNPROTECT(1);
136 :     return val;
137 :     }
138 :    
139 : bates 478 SEXP dtrMatrix_as_matrix(SEXP from)
140 : bates 10 {
141 :     int *Dim = INTEGER(GET_SLOT(from, Matrix_DimSym));
142 :     int m = Dim[0], n = Dim[1];
143 :     SEXP val = PROTECT(allocMatrix(REALSXP, m, n));
144 :     make_array_triangular(Memcpy(REAL(val),
145 : bates 296 REAL(GET_SLOT(from, Matrix_xSym)), m * n),
146 : bates 10 from);
147 : maechler 628 setAttrib(val, R_DimNamesSymbol, GET_SLOT(from, Matrix_DimNamesSym));
148 : bates 10 UNPROTECT(1);
149 :     return val;
150 :     }
151 :    
152 : bates 478 SEXP dtrMatrix_getDiag(SEXP x)
153 : bates 10 {
154 :     int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
155 :     SEXP ret = PROTECT(allocVector(REALSXP, n)),
156 : bates 296 xv = GET_SLOT(x, Matrix_xSym);
157 : bates 10
158 : maechler 951 if ('U' == diag_P(x)[0]) {
159 : bates 10 for (i = 0; i < n; i++) REAL(ret)[i] = 1.;
160 :     } else {
161 :     for (i = 0; i < n; i++) {
162 :     REAL(ret)[i] = REAL(xv)[i * (n + 1)];
163 :     }
164 :     }
165 :     UNPROTECT(1);
166 :     return ret;
167 :     }
168 : maechler 534
169 : bates 478 SEXP dtrMatrix_dgeMatrix_mm_R(SEXP a, SEXP b)
170 : bates 10 {
171 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
172 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
173 :     m = adims[0], n = bdims[1], k = adims[1];
174 :     SEXP val = PROTECT(duplicate(b));
175 :     double one = 1.;
176 : maechler 534
177 : bates 10 if (bdims[0] != k)
178 : bates 582 error(_("Matrices are not conformable for multiplication"));
179 : bates 10 if (m < 1 || n < 1 || k < 1)
180 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
181 : maechler 951 F77_CALL(dtrmm)("R", uplo_P(a), "N", diag_P(a), adims, bdims+1, &one,
182 : bates 10 REAL(GET_SLOT(a, Matrix_xSym)), adims,
183 :     REAL(GET_SLOT(val, Matrix_xSym)), bdims);
184 :     UNPROTECT(1);
185 :     return val;
186 :     }
187 : bates 597
188 :     SEXP dtrMatrix_as_dtpMatrix(SEXP from)
189 :     {
190 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtpMatrix"))),
191 :     uplo = GET_SLOT(from, Matrix_uploSym),
192 :     diag = GET_SLOT(from, Matrix_diagSym),
193 :     dimP = GET_SLOT(from, Matrix_DimSym);
194 :     int n = *INTEGER(dimP);
195 :    
196 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
197 :     SET_SLOT(val, Matrix_diagSym, duplicate(diag));
198 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
199 : maechler 952 full_to_packed_double(
200 :     REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
201 :     REAL(GET_SLOT(from, Matrix_xSym)), n,
202 :     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW,
203 :     *CHAR(STRING_ELT(diag, 0)) == 'U' ? UNT : NUN);
204 : bates 597 UNPROTECT(1);
205 :     return val;
206 :     }

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