SCM

SCM Repository

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

Annotation of /pkg/src/dspMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 635 #include "dspMatrix.h"
2 :    
3 : maechler 945 /* Note: also used for lspMatrix */
4 : bates 635 SEXP dspMatrix_validate(SEXP obj)
5 :     {
6 : maechler 890 SEXP val = symmetricMatrix_validate(obj);
7 :     if(isString(val))
8 :     return(val);
9 :     else { /* identical to the test in dtpMatrix_validate() : */
10 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
11 :     if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))
12 :     return(mkString(_("Incorrect length of 'x' slot")));
13 :     return ScalarLogical(1);
14 :     }
15 : bates 635 }
16 :    
17 :     double get_norm_sp(SEXP obj, char *typstr)
18 :     {
19 :     char typnm[] = {'\0', '\0'};
20 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
21 :     double *work = (double *) NULL;
22 :    
23 :     typnm[0] = norm_type(typstr);
24 :     if (*typnm == 'I' || *typnm == 'O') {
25 : maechler 951 work = (double *) R_alloc(dims[0], sizeof(double));
26 : bates 635 }
27 : maechler 951 return F77_CALL(dlansp)(typnm, uplo_P(obj), dims,
28 :     REAL(GET_SLOT(obj, Matrix_xSym)), work);
29 : bates 635 }
30 :    
31 :     SEXP dspMatrix_norm(SEXP obj, SEXP type)
32 :     {
33 :     return ScalarReal(get_norm_sp(obj, CHAR(asChar(type))));
34 :     }
35 :    
36 : bates 1167 SEXP dspMatrix_rcond(SEXP obj, SEXP type)
37 : bates 635 {
38 : bates 1167 SEXP trf = dspMatrix_trf(obj);
39 : bates 635 char typnm[] = {'\0', '\0'};
40 : bates 1167 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
41 :     double anorm = get_norm_sp(obj, "O"), rcond;
42 : bates 635
43 : bates 1167 typnm[0] = rcond_type(CHAR(asChar(type)));
44 :     F77_CALL(dspcon)(uplo_P(trf), dims,
45 :     REAL (GET_SLOT(trf, Matrix_xSym)),
46 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
47 :     &anorm, &rcond,
48 :     (double *) R_alloc(2*dims[0], sizeof(double)),
49 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
50 :     return ScalarReal(rcond);
51 : bates 635 }
52 :    
53 :     SEXP dspMatrix_solve(SEXP a)
54 :     {
55 :     SEXP trf = dspMatrix_trf(a);
56 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dspMatrix")));
57 :     int *dims = INTEGER(GET_SLOT(trf, Matrix_DimSym)), info;
58 :    
59 :     SET_SLOT(val, Matrix_uploSym, duplicate(GET_SLOT(trf, Matrix_uploSym)));
60 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(trf, Matrix_xSym)));
61 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(trf, Matrix_DimSym)));
62 : maechler 951 F77_CALL(dsptri)(uplo_P(val), dims, REAL(GET_SLOT(val, Matrix_xSym)),
63 : bates 635 INTEGER(GET_SLOT(trf, Matrix_permSym)),
64 :     (double *) R_alloc((long) dims[0], sizeof(double)),
65 :     &info);
66 :     UNPROTECT(1);
67 :     return val;
68 :     }
69 :    
70 : maechler 1432 SEXP dspMatrix_matrix_solve(SEXP a, SEXP b)
71 : bates 635 {
72 :     SEXP trf = dspMatrix_trf(a),
73 : maechler 1432 val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
74 : bates 635 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
75 : bates 1410 *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
76 : bates 652 int n = bdims[0], nrhs = bdims[1], info;
77 : bates 635
78 : bates 1463 if (adims[0] != n || nrhs < 1 || n < 1)
79 : bates 635 error(_("Dimensions of system to be solved are inconsistent"));
80 : maechler 951 F77_CALL(dsptrs)(uplo_P(trf),
81 : bates 652 &n, &nrhs, REAL(GET_SLOT(trf, Matrix_xSym)),
82 : bates 635 INTEGER(GET_SLOT(trf, Matrix_permSym)),
83 : bates 1410 REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);
84 : bates 635 UNPROTECT(1);
85 :     return val;
86 :     }
87 :    
88 :     SEXP dspMatrix_as_dsyMatrix(SEXP from)
89 :     {
90 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))),
91 :     uplo = GET_SLOT(from, Matrix_uploSym),
92 :     dimP = GET_SLOT(from, Matrix_DimSym),
93 :     dmnP = GET_SLOT(from, Matrix_DimNamesSym);
94 :     int n = *INTEGER(dimP);
95 :    
96 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
97 :     SET_SLOT(val, Matrix_DimNamesSym, duplicate(dmnP));
98 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
99 : maechler 952 packed_to_full_double(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),
100 :     REAL(GET_SLOT(from, Matrix_xSym)), n,
101 :     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);
102 : bates 635 UNPROTECT(1);
103 :     return val;
104 :     }
105 :    
106 : bates 1463 SEXP dspMatrix_matrix_mm(SEXP a, SEXP b)
107 : bates 635 {
108 : bates 1463 SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
109 :     int *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
110 : maechler 655 int i, ione = 1, n = bdims[0], nrhs = bdims[1];
111 : maechler 951 char *uplo = uplo_P(a);
112 : bates 652 double *ax = REAL(GET_SLOT(a, Matrix_xSym)), one = 1., zero = 0.,
113 : bates 1463 *vx = REAL(GET_SLOT(val, Matrix_xSym));
114 :     double *bx = Memcpy(Calloc(n * nrhs, double), vx, n * nrhs);
115 : bates 635
116 :     if (bdims[0] != n)
117 :     error(_("Matrices are not conformable for multiplication"));
118 :     if (nrhs < 1 || n < 1)
119 :     error(_("Matrices with zero extents cannot be multiplied"));
120 :     for (i = 0; i < nrhs; i++)
121 : bates 652 F77_CALL(dspmv)(uplo, &n, &one, ax, bx + i * n, &ione,
122 : bates 635 &zero, vx + i * n, &ione);
123 : bates 1463 Free(bx);
124 : bates 635 UNPROTECT(1);
125 :     return val;
126 :     }
127 :    
128 :     SEXP dspMatrix_trf(SEXP x)
129 :     {
130 :     SEXP val = get_factors(x, "pBunchKaufman"),
131 :     dimP = GET_SLOT(x, Matrix_DimSym),
132 :     uploP = GET_SLOT(x, Matrix_uploSym);
133 :     int *dims = INTEGER(dimP), *perm, info;
134 :     int n = dims[0];
135 :     char *uplo = CHAR(STRING_ELT(uploP, 0));
136 :    
137 :     if (val != R_NilValue) return val;
138 :     dims = INTEGER(dimP);
139 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("pBunchKaufman")));
140 :     SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
141 :     SET_SLOT(val, Matrix_diagSym, mkString("N"));
142 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
143 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
144 :     perm = INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, n));
145 :     F77_CALL(dsptrf)(uplo, dims, REAL(GET_SLOT(val, Matrix_xSym)), perm, &info);
146 :     if (info) error(_("Lapack routine %s returned error code %d"), "dsptrf", info);
147 :     UNPROTECT(1);
148 :     return set_factors(x, val, "pBunchKaufman");
149 :     }
150 :    

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