SCM

SCM Repository

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

Annotation of /pkg/src/HBMM.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 825 #include "HBMM.h"
2 :     #include "iohb.h"
3 :     #include "mmio.h"
4 :    
5 : bates 852 #if 0
6 : bates 825 SEXP Matrix_readHarwellBoeing(SEXP filename)
7 :     {
8 :     char *fnm = CHAR(asChar(filename)), *Type = Calloc(4, char);
9 :     int M, N, nz, Nrhs;
10 :     SEXP ans = R_NilValue;
11 :    
12 :     readHB_info(fnm, &M, &N, &nz, &Type, &Nrhs);
13 :     Rprintf("Filename: %s, M=%d, N=%d, nz=%d, Type=\"%s\", Nrhs=%d\n",
14 :     fnm, M, N, nz, Type, Nrhs);
15 :    
16 :     if (toupper(Type[0]) == 'R') { /* Real (double precision) matrix */
17 :     char upT1 = toupper(Type[1]);
18 :     int *dims;
19 :    
20 :     if (upT1 == 'S') {
21 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));
22 :     SET_SLOT(ans, Matrix_uploSym, mkString("L"));
23 :     }
24 :     if (upT1 == 'R' || upT1 == 'U')
25 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
26 :     if (ans == R_NilValue) {Free(Type); return ans;}
27 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, N + 1));
28 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nz));
29 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nz));
30 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
31 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
32 :     dims[0] = M; dims[1] = N;
33 :     readHB_mat_double(fnm, INTEGER(GET_SLOT(ans, Matrix_pSym)),
34 :     INTEGER(GET_SLOT(ans, Matrix_iSym)),
35 :     REAL(GET_SLOT(ans, Matrix_xSym)));
36 :     }
37 :     Free(Type);
38 :     UNPROTECT(1);
39 :     return ans;
40 :     }
41 :    
42 :     SEXP Matrix_readMatrixMarket(SEXP filename)
43 :     {
44 : bates 852 FILE *conn = (FILE *) NULL; /* -Wall */
45 : bates 825 MM_typecode code;
46 :     int *dims, M, N, i, nz;
47 :     SEXP ans = R_NilValue;
48 :    
49 :     if (isString(filename)) {
50 :     conn = fopen(CHAR(asChar(filename)), "r");
51 :     if (conn == NULL) {
52 :     error("Unable to open file: %s", CHAR(asChar(filename)));
53 :     return R_NilValue;
54 :     }
55 :     } else {
56 :     error("non-string values not presently accepted");
57 :     }
58 :    
59 : bates 841 if ((nz = mm_read_banner(conn, &code))) {
60 : bates 838 fclose(conn);
61 : bates 825 error("mm_read_banner returned code %d", nz);
62 :     }
63 :     if (!mm_is_valid(code)) {
64 : bates 838 fclose(conn);
65 : bates 825 error("Invalid code: %s", mm_typecode_to_str(code));
66 :     }
67 :    
68 :     if (mm_is_sparse(code)) mm_read_mtx_crd_size(conn, &M, &N, &nz);
69 :     if (mm_is_dense(code)) mm_read_mtx_array_size(conn, &M, &N);
70 : maechler 951
71 : bates 825 if (mm_is_real(code)) {
72 :     if (mm_is_sparse(code)) {
73 :     int *ii, *jj;
74 :     double *vv;
75 :    
76 :     if (mm_is_general(code))
77 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));
78 :     if (mm_is_symmetric(code)) {
79 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsTMatrix")));
80 :     SET_SLOT(ans, Matrix_uploSym, mkString("L"));
81 :     }
82 :     if (ans == R_NilValue)
83 :     error("Unrecognized matrix type: %s",
84 :     mm_typecode_to_str(code));
85 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nz));
86 :     ii = INTEGER(GET_SLOT(ans, Matrix_iSym));
87 :     SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nz));
88 :     jj = INTEGER(GET_SLOT(ans, Matrix_jSym));
89 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nz));
90 :     vv = REAL(GET_SLOT(ans, Matrix_xSym));
91 :     for (i = 0; i < nz; i++) {
92 :     if (fscanf(conn, "%d %d %lg", &ii[i], &jj[i], &vv[i]) != 3)
93 :     error("Premature end of file");
94 :     ii[i] -= 1; /* change indices to zero-based */
95 :     jj[i] -= 1;
96 :     }
97 :     }
98 :     if (mm_is_dense(code)) {
99 :     int j;
100 :     double *vv;
101 :    
102 :     if (mm_is_general(code))
103 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
104 :     /* Need to adjust reading loop for symmetric matrices */
105 :     /* if (mm_is_symmetric(code)) { */
106 :     /* ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))); */
107 :     /* SET_SLOT(ans, Matrix_uploSym, mkString("L")); */
108 :     /* } */
109 :     if (ans == R_NilValue)
110 :     error("Unrecognized matrix type: %s",
111 :     mm_typecode_to_str(code));
112 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, M * N));
113 :     vv = REAL(GET_SLOT(ans, Matrix_xSym));
114 :     for (j = 0; j < N; j++) {
115 :     for (i = 0; i < M; i++) {
116 :     if (fscanf(conn, "%lg", &vv[i + j * M]) != 1)
117 :     error("Premature end of file");
118 :     }
119 :     }
120 :     }
121 :     } else {
122 :     error("Only real matrices handled at this point");
123 :     }
124 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
125 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
126 :     dims[0] = M; dims[1] = N;
127 :    
128 : bates 838 fclose(conn);
129 : bates 825 UNPROTECT(1);
130 :     return ans;
131 :     }
132 : maechler 951 #endif
133 : bates 825
134 : bates 835 SEXP Matrix_writeHarwellBoeing(SEXP obj, SEXP file, SEXP typep)
135 :     {
136 :     char *type = CHAR(asChar(typep)), Type[4] = "RUA";
137 : bates 852 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)),
138 :     *ii = (int *) NULL, *pp = (int *) NULL;
139 :     int M = dims[0], N = dims[1], nz = -1;
140 :     double *xx = (double *) NULL;
141 : bates 835
142 :     if (type[2] == 'C' || type[2] == 'T') {
143 :     SEXP islot = GET_SLOT(obj, Matrix_iSym);
144 :     nz = LENGTH(islot);
145 :     ii = INTEGER(islot);
146 : bates 841 if (type[2] == 'T') { /* create column pointers */
147 :     int *i1 = Calloc(nz, int);
148 :     double *x1 = Calloc(nz, double);
149 : maechler 951
150 : bates 841 pp = Calloc(N + 1, int);
151 :     triplet_to_col(M, N, nz, ii,
152 :     INTEGER(GET_SLOT(obj, Matrix_jSym)), xx,
153 :     pp, i1, x1);
154 :     nz = pp[N];
155 :     xx = x1;
156 :     ii = i1;
157 :     } else pp = INTEGER(GET_SLOT(obj, Matrix_pSym));
158 : bates 835 } else error("Only types 'C' and 'T' allowed");
159 :    
160 :     if (type[0] == 'D') {
161 :     xx = REAL(GET_SLOT(obj, Matrix_xSym));
162 :     } else error("Only real matrices allowed");
163 :    
164 : maechler 951 if (!isString(file))
165 : bates 835 error("non-string values for file not presently accepted");
166 :    
167 :     if (type[1] == 'S') {
168 : maechler 951 if (*uplo_P(obj) != 'L')
169 : bates 835 error("Symmetric matrices must be stored in lower triangle");
170 :     Type[1] = 'S';
171 :     }
172 :    
173 :     writeHB_mat_double(CHAR(asChar(file)), M, N, nz, pp, ii, xx, 0,
174 : maechler 951 (double *)NULL, (double *)NULL, (double *)NULL,
175 : bates 840 "", "", Type, (char*)NULL, (char*)NULL,
176 : maechler 951 (char*)NULL, (char*)NULL, "RUA");
177 : bates 835
178 :     if (type[2] == 'T') {Free(ii); Free(pp); Free(xx);}
179 :     return R_NilValue;
180 :     }
181 :    
182 :     SEXP Matrix_writeMatrixMarket(SEXP obj, SEXP file, SEXP typep)
183 :     {
184 :     char *type = CHAR(asChar(typep));
185 : bates 852 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)),
186 :     *ii = (int *) NULL, *jj = (int *) NULL;
187 :     int M = dims[0], N = dims[1], nz = -1;
188 : bates 835 MM_typecode matcode;
189 : bates 852 double *xx = (double *) NULL;
190 : bates 835
191 :     mm_set_matrix(&matcode);
192 :     if (type[2] == 'C' || type[2] == 'T') {
193 :     SEXP islot = GET_SLOT(obj, Matrix_iSym);
194 :     nz = LENGTH(islot);
195 :     ii = INTEGER(islot);
196 :     mm_set_coordinate(&matcode);
197 :     } else error("Only types 'C' and 'T' allowed");
198 :    
199 :     if (type[0] == 'D') {
200 :     xx = REAL(GET_SLOT(obj, Matrix_xSym));
201 :     mm_set_real(&matcode);
202 :     } else error("Only real matrices allowed");
203 :    
204 : maechler 951 if (!isString(file))
205 : bates 835 error("non-string values for file not currently allowed");
206 :    
207 :     if (type[1] == 'S') {
208 : maechler 951 if (*uplo_P(obj) != 'L')
209 : bates 835 error("Symmetric matrices must be stored in lower triangle");
210 :     mm_set_symmetric(&matcode);
211 :     }
212 :     if (type[1] == 'G') mm_set_general(&matcode);
213 :    
214 :     if (type[2] == 'C')
215 :     jj = expand_cmprPt(N, INTEGER(GET_SLOT(obj, Matrix_pSym)),
216 :     Calloc(nz, int));
217 :     if (type[2] == 'T')
218 :     jj = INTEGER(GET_SLOT(obj, Matrix_jSym));
219 :     if (!jj) error("storage mode must be T or C");
220 :    
221 :     mm_write_mtx_crd(CHAR(STRING_ELT(file, 0)), M, N, nz, ii, jj, xx,
222 :     matcode);
223 :    
224 :     if (type[2] == 'C') Free(jj);
225 :     return R_NilValue;
226 : maechler 951
227 : bates 835 }

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