SCM

SCM Repository

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

Annotation of /pkg/Matrix/src/dsyMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 478 #include "dsyMatrix.h"
2 : bates 10
3 : maechler 890 SEXP symmetricMatrix_validate(SEXP obj)
4 : bates 10 {
5 : maechler 890 SEXP val = GET_SLOT(obj, Matrix_DimSym);
6 :     if (LENGTH(val) < 2)
7 :     return mkString(_("'Dim' slot has length less than two"));
8 :     if (INTEGER(val)[0] != INTEGER(val)[1])
9 :     return mkString(_("Matrix is not square"));
10 : bates 592 if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
11 :     "LU", "uplo"))) return val;
12 : bates 10 return ScalarLogical(1);
13 :     }
14 :    
15 : bates 1867 double get_norm_sy(SEXP obj, const char *typstr)
16 : maechler 534 {
17 : bates 296 char typnm[] = {'\0', '\0'};
18 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
19 :     double *work = (double *) NULL;
20 :    
21 : maechler 2117 typnm[0] = La_norm_type(typstr);
22 : bates 296 if (*typnm == 'I' || *typnm == 'O') {
23 :     work = (double *) R_alloc(dims[0], sizeof(double));
24 :     }
25 : maechler 951 return F77_CALL(dlansy)(typnm, uplo_P(obj),
26 : bates 296 dims, REAL(GET_SLOT(obj, Matrix_xSym)),
27 :     dims, work);
28 :     }
29 :    
30 : bates 478 SEXP dsyMatrix_norm(SEXP obj, SEXP type)
31 : bates 296 {
32 :     return ScalarReal(get_norm_sy(obj, CHAR(asChar(type))));
33 :     }
34 :    
35 : bates 1167
36 :     SEXP dsyMatrix_rcond(SEXP obj, SEXP type)
37 : bates 296 {
38 : bates 1167 SEXP trf = dsyMatrix_trf(obj);
39 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
40 :     double anorm = get_norm_sy(obj, "O");
41 : bates 296 double rcond;
42 :    
43 : bates 1167 F77_CALL(dsycon)(uplo_P(trf), dims,
44 :     REAL (GET_SLOT(trf, Matrix_xSym)), dims,
45 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
46 :     &anorm, &rcond,
47 :     (double *) R_alloc(2*dims[0], sizeof(double)),
48 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
49 :     return ScalarReal(rcond);
50 : bates 296 }
51 :    
52 : bates 478 SEXP dsyMatrix_solve(SEXP a)
53 : bates 296 {
54 : bates 643 SEXP trf = dsyMatrix_trf(a);
55 : mmaechler 3273 SEXP val = PROTECT(NEW_OBJECT_OF_CLASS("dsyMatrix"));
56 : bates 643 int *dims = INTEGER(GET_SLOT(trf, Matrix_DimSym)), info;
57 :    
58 : maechler 1968 slot_dup(val, trf, Matrix_uploSym);
59 :     slot_dup(val, trf, Matrix_xSym);
60 :     slot_dup(val, trf, Matrix_DimSym);
61 : maechler 951 F77_CALL(dsytri)(uplo_P(val), dims,
62 :     REAL(GET_SLOT(val, Matrix_xSym)), dims,
63 : bates 643 INTEGER(GET_SLOT(trf, Matrix_permSym)),
64 :     (double *) R_alloc((long) dims[0], sizeof(double)),
65 :     &info);
66 :     UNPROTECT(1);
67 :     return val;
68 : bates 296 }
69 :    
70 : bates 1463 SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)
71 : bates 643 {
72 :     SEXP trf = dsyMatrix_trf(a),
73 : bates 1463 val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
74 : bates 643 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
75 : bates 1463 *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym)),
76 : bates 643 info;
77 :    
78 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
79 :     error(_("Dimensions of system to be solved are inconsistent"));
80 : maechler 951 F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,
81 : bates 643 REAL(GET_SLOT(trf, Matrix_xSym)), adims,
82 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
83 :     REAL(GET_SLOT(val, Matrix_xSym)),
84 :     bdims, &info);
85 :     UNPROTECT(1);
86 :     return val;
87 :     }
88 :    
89 : maechler 2005 SEXP dsyMatrix_as_matrix(SEXP from, SEXP keep_dimnames)
90 : bates 10 {
91 :     int n = INTEGER(GET_SLOT(from, Matrix_DimSym))[0];
92 :     SEXP val = PROTECT(allocMatrix(REALSXP, n, n));
93 : maechler 534
94 : maechler 1200 make_d_matrix_symmetric(Memcpy(REAL(val),
95 :     REAL(GET_SLOT(from, Matrix_xSym)), n * n),
96 :     from);
97 : maechler 2005 if(asLogical(keep_dimnames))
98 : mmaechler 3055 setAttrib(val, R_DimNamesSymbol, R_symmetric_Dimnames(from));
99 : bates 10 UNPROTECT(1);
100 :     return val;
101 :     }
102 :    
103 : bates 1463 SEXP dsyMatrix_matrix_mm(SEXP a, SEXP b, SEXP rtP)
104 : bates 10 {
105 : mmaechler 2787 SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));// incl. its dimnames
106 : maechler 2018 int rt = asLogical(rtP); /* if(rt), compute b %*% a, else a %*% b */
107 : bates 10 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
108 : maechler 2036 *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym)),
109 : maechler 2018 m = bdims[0], n = bdims[1];
110 : mmaechler 3105 double one = 1., zero = 0., mn = ((double) m) * ((double) n);
111 :     if (mn > INT_MAX)
112 :     error(_("Matrix dimension %d x %d (= %g) is too large"), m, n, mn);
113 :     // else: m * n will not overflow below
114 :     double *bcp, *vx = REAL(GET_SLOT(val, Matrix_xSym));
115 :     C_or_Alloca_TO(bcp, m * n, double);
116 :     Memcpy(bcp, vx, m * n);
117 : maechler 534
118 : bates 1463 if ((rt && n != adims[0]) || (!rt && m != adims[0]))
119 : bates 582 error(_("Matrices are not conformable for multiplication"));
120 : mmaechler 2990 if (m >=1 && n >= 1)
121 : maechler 2115 F77_CALL(dsymm)(rt ? "R" :"L", uplo_P(a), &m, &n, &one,
122 :     REAL(GET_SLOT(a, Matrix_xSym)), adims, bcp,
123 :     &m, &zero, vx, &m);
124 : mmaechler 2990 // add dimnames:
125 : mmaechler 3213 int nd = rt ?
126 :     1 : // v <- b %*% a : rownames(v) == rownames(b) are already there
127 :     0; // v <- a %*% b : colnames(v) == colnames(b) are already there
128 :     SEXP nms = PROTECT(duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), nd)));
129 :     SET_VECTOR_ELT(GET_SLOT(val, Matrix_DimNamesSym), nd, nms);
130 : mmaechler 3105 if(mn >= SMALL_4_Alloca) Free(bcp);
131 : mmaechler 3213 UNPROTECT(2);
132 : bates 10 return val;
133 :     }
134 :    
135 : bates 631 SEXP dsyMatrix_trf(SEXP x)
136 :     {
137 :     SEXP val = get_factors(x, "BunchKaufman"),
138 :     dimP = GET_SLOT(x, Matrix_DimSym),
139 :     uploP = GET_SLOT(x, Matrix_uploSym);
140 :     int *dims = INTEGER(dimP), *perm, info;
141 :     int lwork = -1, n = dims[0];
142 : bates 1867 const char *uplo = CHAR(STRING_ELT(uploP, 0));
143 : bates 631 double tmp, *vx, *work;
144 :    
145 :     if (val != R_NilValue) return val;
146 :     dims = INTEGER(dimP);
147 : mmaechler 3273 val = PROTECT(NEW_OBJECT_OF_CLASS("BunchKaufman"));
148 : bates 631 SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
149 :     SET_SLOT(val, Matrix_diagSym, mkString("N"));
150 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
151 :     vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n));
152 :     AZERO(vx, n * n);
153 :     F77_CALL(dlacpy)(uplo, &n, &n, REAL(GET_SLOT(x, Matrix_xSym)), &n, vx, &n);
154 :     perm = INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, n));
155 :     F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, &tmp, &lwork, &info);
156 :     lwork = (int) tmp;
157 : mmaechler 3105 C_or_Alloca_TO(work, lwork, double);
158 :    
159 : bates 631 F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, work, &lwork, &info);
160 : mmaechler 3105
161 :     if(lwork >= SMALL_4_Alloca) Free(work);
162 : bates 631 if (info) error(_("Lapack routine dsytrf returned error code %d"), info);
163 :     UNPROTECT(1);
164 :     return set_factors(x, val, "BunchKaufman");
165 :     }
166 :    
167 : mmaechler 2805 // this is very close to lsyMatrix_as_lsp*() in ./ldense.c -- keep synced !
168 : bates 643 SEXP dsyMatrix_as_dspMatrix(SEXP from)
169 :     {
170 : mmaechler 3273 SEXP val = PROTECT(NEW_OBJECT_OF_CLASS("dspMatrix")),
171 : bates 643 uplo = GET_SLOT(from, Matrix_uploSym),
172 :     dimP = GET_SLOT(from, Matrix_DimSym);
173 :     int n = *INTEGER(dimP);
174 :    
175 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
176 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
177 : maechler 952 full_to_packed_double(
178 :     REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
179 : mmaechler 2805 REAL( GET_SLOT(from, Matrix_xSym)), n,
180 : maechler 952 *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW, NUN);
181 : bates 1236 SET_SLOT(val, Matrix_DimNamesSym,
182 :     duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
183 : mmaechler 2271 SET_SLOT(val, Matrix_factorSym,
184 :     duplicate(GET_SLOT(from, Matrix_factorSym)));
185 : bates 643 UNPROTECT(1);
186 :     return val;
187 :     }

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge