Line data Source code
1 : /*
2 : * File: MatrixFactory.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on July 12, 2015, 7:50 PM
6 : *
7 : * ForBES is free software: you can redistribute it and/or modify
8 : * it under the terms of the GNU Lesser General Public License as published by
9 : * the Free Software Foundation, either version 3 of the License, or
10 : * (at your option) any later version.
11 : *
12 : * ForBES is distributed in the hope that it will be useful,
13 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : * GNU Lesser General Public License for more details.
16 : *
17 : * You should have received a copy of the GNU Lesser General Public License
18 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
19 : */
20 :
21 : #include <set>
22 :
23 : #include "MatrixFactory.h"
24 : #include "Matrix.h"
25 :
26 : #include <vector> // std::vector
27 : #include <algorithm> // std::random_shuffle
28 : #include <ctime> // std::time
29 : #include <functional>
30 :
31 : #include <cmath>
32 : #include <sstream>
33 :
34 : typedef std::pair<size_t, size_t> nice_pair;
35 :
36 7 : Matrix MatrixFactory::MakeIdentity(size_t n, double alpha) {
37 7 : Matrix mat(n, n, Matrix::MATRIX_DIAGONAL);
38 1344 : for (size_t i = 0; i < n; i++) {
39 1337 : mat[i] = alpha;
40 : }
41 7 : return mat;
42 : }
43 :
44 6661 : Matrix MatrixFactory::MakeRandomSparse(size_t nrows, size_t ncols, size_t nnz, float offset, float scale) {
45 6661 : Matrix R = MakeSparse(nrows, ncols, nnz, Matrix::SPARSE_UNSYMMETRIC);
46 13320 : std::set<nice_pair> s;
47 6660 : nice_pair p;
48 : while (true) { // construct pairs
49 512383 : p.first = (std::rand() % nrows);
50 512383 : p.second = (std::rand() % ncols);
51 512383 : s.insert(p);
52 512383 : if (s.size() == nnz) {
53 6660 : break;
54 : }
55 : }
56 365832 : for (std::set<nice_pair>::iterator it = s.begin(); it != s.end(); ++it) {
57 : double rand;
58 359172 : rand = offset + (scale * std::rand()) / RAND_MAX;
59 359172 : R.set(it->first, it->second, rand);
60 : }
61 13320 : return R;
62 : }
63 :
64 6235 : Matrix MatrixFactory::MakeRandomMatrix(size_t nrows, size_t ncols, float offset, float scale, Matrix::MatrixType type) {
65 6235 : size_t len = 0;
66 6235 : switch (type) {
67 : case Matrix::MATRIX_DENSE:
68 5460 : len = nrows * ncols;
69 5460 : break;
70 : case Matrix::MATRIX_LOWERTR:
71 : case Matrix::MATRIX_SYMMETRIC:
72 162 : len = nrows * (nrows + 1) / 2;
73 162 : break;
74 : case Matrix::MATRIX_DIAGONAL:
75 612 : len = nrows;
76 612 : break;
77 : case Matrix::MATRIX_SPARSE:
78 1 : return MakeRandomSparse(nrows, ncols, std::ceil((nrows * ncols) / 3), offset, scale);
79 : break;
80 : default:
81 0 : throw std::logic_error("unsupported");
82 : }
83 6234 : Matrix mat(nrows, ncols, type);
84 1603276 : for (size_t j = 0; j < len; j++) {
85 1597042 : mat[j] = static_cast<double> (offset + (scale * std::rand()) / RAND_MAX);
86 : }
87 6234 : return mat;
88 : }
89 :
90 3947 : Matrix MatrixFactory::MakeRandomMatrix(size_t nrows, size_t ncols, float offset, float scale) {
91 3947 : return MatrixFactory::MakeRandomMatrix(nrows, ncols, offset, scale, Matrix::MATRIX_DENSE);
92 : }
93 :
94 6678 : Matrix MatrixFactory::MakeSparse(size_t nrows, size_t ncols, size_t max_nnz, Matrix::SparseMatrixType stype) {
95 6678 : if (max_nnz > nrows * ncols) {
96 1 : std::ostringstream oss;
97 1 : oss << "Matrix " << nrows << "x" << ncols << "(max_size=" << (nrows * ncols)
98 1 : << " cannot allocate " << max_nnz << "non-zeros";
99 1 : throw std::invalid_argument(oss.str().c_str());
100 : }
101 6677 : Matrix matrix(nrows, ncols, Matrix::MATRIX_SPARSE);
102 6677 : matrix.m_triplet = cholmod_allocate_triplet(nrows, ncols, max_nnz, stype, CHOLMOD_REAL, Matrix::cholmod_handle());
103 6677 : matrix.m_sparseStorageType = Matrix::CHOLMOD_TYPE_TRIPLET;
104 6677 : return matrix;
105 : }
106 :
107 2 : Matrix MatrixFactory::MakeSparseSymmetric(size_t n, size_t max_nnz) {
108 2 : return MakeSparse(n, n, max_nnz, Matrix::SPARSE_SYMMETRIC_L);
109 : }
110 :
111 3 : Matrix MatrixFactory::ReadSparse(FILE* fp) {
112 : cholmod_sparse *sp;
113 3 : sp = cholmod_read_sparse(fp, Matrix::cholmod_handle());
114 :
115 3 : Matrix mat(sp->nrow, sp->ncol, Matrix::MATRIX_SPARSE);
116 3 : mat.m_sparse = sp;
117 3 : mat.m_sparseStorageType = Matrix::CHOLMOD_TYPE_SPARSE;
118 3 : mat.m_triplet = cholmod_sparse_to_triplet(sp, Matrix::cholmod_handle());
119 3 : return mat;
120 : }
121 :
122 67 : Matrix MatrixFactory::ShallowMatrix(const Matrix& orig) {
123 67 : Matrix mat_shallow = Matrix(true);
124 67 : mat_shallow.m_transpose = orig.m_transpose;
125 67 : mat_shallow.m_nrows = orig.m_nrows;
126 67 : mat_shallow.m_ncols = orig.m_ncols;
127 67 : mat_shallow.m_dataLength = orig.m_dataLength;
128 67 : mat_shallow.m_delete_data = false;
129 67 : mat_shallow.m_data = orig.m_data;
130 67 : mat_shallow.m_type = orig.m_type;
131 67 : mat_shallow.m_sparseStorageType = orig.m_sparseStorageType;
132 67 : mat_shallow.m_dense = orig.m_dense;
133 67 : mat_shallow.m_triplet = orig.m_triplet;
134 67 : mat_shallow.m_sparse = orig.m_sparse;
135 67 : return mat_shallow;
136 : }
137 :
138 242 : Matrix MatrixFactory::ShallowVector(const Matrix& orig, size_t size, size_t offset) {
139 242 : if (orig.getType() != Matrix::MATRIX_DENSE) {
140 1 : throw std::invalid_argument("This method can only be applied to dense vectors");
141 : }
142 241 : if (orig.getNcols() != 1 && orig.getNrows() != 1) {
143 1 : throw std::invalid_argument("This method can only be applied to vectors");
144 : }
145 240 : Matrix v_shallow = Matrix(true);
146 240 : v_shallow.m_transpose = orig.m_transpose;
147 240 : v_shallow.m_nrows = orig.m_transpose ? 1 : size;
148 240 : v_shallow.m_ncols = orig.m_transpose ? size : 1;
149 240 : v_shallow.m_dataLength = v_shallow.m_nrows * v_shallow.m_ncols;
150 240 : v_shallow.m_delete_data = false;
151 240 : v_shallow.m_data = orig.m_data + offset;
152 240 : return v_shallow;
153 : }
154 :
155 10 : Matrix MatrixFactory::ShallowVector(const Matrix& orig, size_t offset) {
156 10 : return MatrixFactory::ShallowVector(orig, orig.length() - offset, offset);
157 : }
158 :
159 167 : Matrix MatrixFactory::ShallowVector(double* data, size_t size, size_t offset) {
160 167 : Matrix v_shallow = Matrix(true);
161 167 : v_shallow.m_transpose = false;
162 167 : v_shallow.m_nrows = size;
163 167 : v_shallow.m_ncols = 1;
164 167 : v_shallow.m_dataLength = v_shallow.m_nrows;
165 167 : v_shallow.m_delete_data = false;
166 167 : v_shallow.m_data = data + offset;
167 167 : return v_shallow;
168 : }
169 :
170 3 : Matrix MatrixFactory::ShallowVector() {
171 3 : Matrix v_shallow = Matrix(true);
172 3 : return v_shallow;
173 : }
174 :
175 96 : Matrix MatrixFactory::ShallowSubVector(Matrix& A, size_t j) {
176 96 : double * data = A.getData();
177 96 : data += A.getNrows() * j;
178 96 : return MatrixFactory::ShallowVector(data, A.getNrows(), 0);
179 : }
180 :
181 :
182 :
183 :
184 :
|
