SCM

SCM Repository

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

Annotation of /pkg/src/dsyMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1463 - (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 : maechler 890 SEXP dsyMatrix_validate(SEXP obj)
16 :     {
17 : maechler 1164 /* since "dsy" inherits from "symmetric", and "dMatrix", only need this:*/
18 :     return dense_nonpacked_validate(obj);
19 : maechler 890 }
20 :    
21 : bates 296 double get_norm_sy(SEXP obj, char *typstr)
22 : maechler 534 {
23 : bates 296 char typnm[] = {'\0', '\0'};
24 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
25 :     double *work = (double *) NULL;
26 :    
27 :     typnm[0] = norm_type(typstr);
28 :     if (*typnm == 'I' || *typnm == 'O') {
29 :     work = (double *) R_alloc(dims[0], sizeof(double));
30 :     }
31 : maechler 951 return F77_CALL(dlansy)(typnm, uplo_P(obj),
32 : bates 296 dims, REAL(GET_SLOT(obj, Matrix_xSym)),
33 :     dims, work);
34 :     }
35 :    
36 : bates 478 SEXP dsyMatrix_norm(SEXP obj, SEXP type)
37 : bates 296 {
38 :     return ScalarReal(get_norm_sy(obj, CHAR(asChar(type))));
39 :     }
40 :    
41 : bates 1167
42 :     SEXP dsyMatrix_rcond(SEXP obj, SEXP type)
43 : bates 296 {
44 : bates 1167 SEXP trf = dsyMatrix_trf(obj);
45 : bates 296 char typnm[] = {'\0', '\0'};
46 : bates 1167 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
47 :     double anorm = get_norm_sy(obj, "O");
48 : bates 296 double rcond;
49 :    
50 : bates 1167 typnm[0] = rcond_type(CHAR(asChar(type)));
51 :     F77_CALL(dsycon)(uplo_P(trf), dims,
52 :     REAL (GET_SLOT(trf, Matrix_xSym)), dims,
53 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
54 :     &anorm, &rcond,
55 :     (double *) R_alloc(2*dims[0], sizeof(double)),
56 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
57 :     return ScalarReal(rcond);
58 : bates 296 }
59 :    
60 : bates 478 SEXP dsyMatrix_solve(SEXP a)
61 : bates 296 {
62 : bates 643 SEXP trf = dsyMatrix_trf(a);
63 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix")));
64 :     int *dims = INTEGER(GET_SLOT(trf, Matrix_DimSym)), info;
65 :    
66 :     SET_SLOT(val, Matrix_uploSym, duplicate(GET_SLOT(trf, Matrix_uploSym)));
67 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(trf, Matrix_xSym)));
68 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(trf, Matrix_DimSym)));
69 : maechler 951 F77_CALL(dsytri)(uplo_P(val), dims,
70 :     REAL(GET_SLOT(val, Matrix_xSym)), dims,
71 : bates 643 INTEGER(GET_SLOT(trf, Matrix_permSym)),
72 :     (double *) R_alloc((long) dims[0], sizeof(double)),
73 :     &info);
74 :     UNPROTECT(1);
75 :     return val;
76 : bates 296 }
77 :    
78 : bates 1463 SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)
79 : bates 643 {
80 :     SEXP trf = dsyMatrix_trf(a),
81 : bates 1463 val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
82 : bates 643 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
83 : bates 1463 *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym)),
84 : bates 643 info;
85 :    
86 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
87 :     error(_("Dimensions of system to be solved are inconsistent"));
88 : maechler 951 F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,
89 : bates 643 REAL(GET_SLOT(trf, Matrix_xSym)), adims,
90 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
91 :     REAL(GET_SLOT(val, Matrix_xSym)),
92 :     bdims, &info);
93 :     UNPROTECT(1);
94 :     return val;
95 :     }
96 :    
97 : bates 478 SEXP dsyMatrix_as_matrix(SEXP from)
98 : bates 10 {
99 :     int n = INTEGER(GET_SLOT(from, Matrix_DimSym))[0];
100 :     SEXP val = PROTECT(allocMatrix(REALSXP, n, n));
101 : maechler 534
102 : maechler 1200 make_d_matrix_symmetric(Memcpy(REAL(val),
103 :     REAL(GET_SLOT(from, Matrix_xSym)), n * n),
104 :     from);
105 : maechler 628 setAttrib(val, R_DimNamesSymbol, GET_SLOT(from, Matrix_DimNamesSym));
106 : bates 10 UNPROTECT(1);
107 :     return val;
108 :     }
109 :    
110 : bates 1463 SEXP dsyMatrix_matrix_mm(SEXP a, SEXP b, SEXP rtP)
111 : bates 10 {
112 : bates 1463 SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
113 : bates 10 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
114 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
115 : bates 1463 m = bdims[0], n = bdims[1], rt = asLogical(rtP);
116 : bates 10 double one = 1., zero = 0.;
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 : bates 1463 if (m < 1 || n < 1)
121 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
122 : bates 1463 F77_CALL(dsymm)(rt ? "R" :"L", uplo_P(a), &m, &n, &one,
123 : bates 10 REAL(GET_SLOT(a, Matrix_xSym)), adims,
124 : bates 1463 REAL(GET_SLOT(b, Matrix_xSym)), &m,
125 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), &m);
126 : bates 10 UNPROTECT(1);
127 :     return val;
128 :     }
129 :    
130 : bates 631 SEXP dsyMatrix_trf(SEXP x)
131 :     {
132 :     SEXP val = get_factors(x, "BunchKaufman"),
133 :     dimP = GET_SLOT(x, Matrix_DimSym),
134 :     uploP = GET_SLOT(x, Matrix_uploSym);
135 :     int *dims = INTEGER(dimP), *perm, info;
136 :     int lwork = -1, n = dims[0];
137 :     char *uplo = CHAR(STRING_ELT(uploP, 0));
138 :     double tmp, *vx, *work;
139 :    
140 :     if (val != R_NilValue) return val;
141 :     dims = INTEGER(dimP);
142 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("BunchKaufman")));
143 :     SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
144 :     SET_SLOT(val, Matrix_diagSym, mkString("N"));
145 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
146 :     vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n));
147 :     AZERO(vx, n * n);
148 :     F77_CALL(dlacpy)(uplo, &n, &n, REAL(GET_SLOT(x, Matrix_xSym)), &n, vx, &n);
149 :     perm = INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, n));
150 :     F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, &tmp, &lwork, &info);
151 :     lwork = (int) tmp;
152 :     work = Calloc(lwork, double);
153 :     F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, work, &lwork, &info);
154 :     if (info) error(_("Lapack routine dsytrf returned error code %d"), info);
155 :     UNPROTECT(1);
156 :     Free(work);
157 :     return set_factors(x, val, "BunchKaufman");
158 :     }
159 :    
160 : bates 643 SEXP dsyMatrix_as_dspMatrix(SEXP from)
161 :     {
162 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dspMatrix"))),
163 :     uplo = GET_SLOT(from, Matrix_uploSym),
164 :     dimP = GET_SLOT(from, Matrix_DimSym);
165 :     int n = *INTEGER(dimP);
166 :    
167 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
168 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
169 : maechler 952 full_to_packed_double(
170 :     REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
171 :     REAL(GET_SLOT(from, Matrix_xSym)), n,
172 :     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW, NUN);
173 : bates 1236 SET_SLOT(val, Matrix_DimNamesSym,
174 :     duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
175 : bates 643 UNPROTECT(1);
176 :     return val;
177 :     }

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