1 // Copyright (c) 2015-2016, Massachusetts Institute of Technology
2 // Copyright (c) 2016-2017 Sandia Corporation
3 // Copyright (c) 2017 NTESS, LLC.
4
5 // This file is part of the Compressed Continuous Computation (C3) Library
6 // Author: Alex A. Gorodetsky
7 // Contact: alex@alexgorodetsky.com
8
9 // All rights reserved.
10
11 // Redistribution and use in source and binary forms, with or without modification,
12 // are permitted provided that the following conditions are met:
13
14 // 1. Redistributions of source code must retain the above copyright notice,
15 // this list of conditions and the following disclaimer.
16
17 // 2. Redistributions in binary form must reproduce the above copyright notice,
18 // this list of conditions and the following disclaimer in the documentation
19 // and/or other materials provided with the distribution.
20
21 // 3. Neither the name of the copyright holder nor the names of its contributors
22 // may be used to endorse or promote products derived from this software
23 // without specific prior written permission.
24
25 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
28 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
29 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
31 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
32 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36 //Code
37
38
39
40
41
42
43
44
45 #include <stdlib.h>
46 #include <stdio.h>
47 #include <string.h>
48 #include <math.h>
49
50 #include "array.h"
51 #include "matrix_util.h"
52
53
54 /***********************************************************//**
55 Function v2m
56
57 Purpose: Convert a double array to a matrix
58
59 Parameters:
60 - N (IN) - Number of elements in array
61 - arr (IN) - array to convert to a matrix
62 - asrow (IN)
63 * 1 to create a row matrix (1 x N)
64 * 0 to create a column matrix (N x 1)
65
66 Returns: either row or column Matrix Structure *mat*
67 ***************************************************************/
v2m(size_t nVals,double * arr,int asrow)68 struct mat * v2m(size_t nVals,double * arr,int asrow)
69 {
70 struct mat * outmat;
71 if (NULL == ( outmat = malloc(sizeof(struct mat)))){
72 fprintf(stderr, "failed to allocate memory.\n");
73 exit(1);
74 }
75
76 outmat->vals = calloc_double(nVals);
77 memmove(outmat->vals, arr, nVals * sizeof(double));
78 if (asrow == 0){
79 outmat->nrows = nVals;
80 outmat->ncols = 1;
81 }
82 else if (asrow == 1){
83 outmat->nrows = 1;
84 outmat->ncols = nVals;
85 }
86 return outmat;
87 }
88
89 /***********************************************************//**
90 Function diagv2m
91
92 Purpose: Convert a double array to a diagonal matrix
93
94 Parameters:
95 - N (IN) - Number of elements in array
96 - arr (IN) - array to convert to a matrix
97
98 Returns: diagonal matrix *mat*
99 ***************************************************************/
diagv2m(size_t nVals,double * arr)100 struct mat * diagv2m(size_t nVals,double * arr)
101 {
102 struct mat * outmat;
103 if (NULL == ( outmat = malloc(sizeof(struct mat)))){
104 fprintf(stderr, "failed to allocate memory.\n");
105 exit(1);
106 }
107 outmat->vals = calloc_double(nVals * nVals);
108 outmat->nrows = nVals;
109 outmat->ncols = nVals;
110 size_t ii;
111 for (ii = 0; ii < nVals; ii++){
112 outmat->vals[ii*nVals+ii] = arr[ii];
113 }
114 return outmat;
115 }
116
117
118 /***********************************************************//**
119 Function horizcat
120
121 Purpose: Concatenate the columns of two matricies
122
123 Parameters:
124 - a (IN) - left matrix
125 - b (IN) - right matrix
126
127 Returns: horizontally stacked matrix *mat*
128 ***************************************************************/
horizcat(const struct mat * a,const struct mat * b)129 struct mat * horizcat(const struct mat * a, const struct mat * b)
130 {
131 if (a->nrows != b->nrows) {
132 printf("Should not be concatenating these matrices!!\n");
133 }
134
135 struct mat * m;
136 if (NULL == ( m = malloc(sizeof(struct mat)))){
137 fprintf(stderr, "failed to allocate memory.\n");
138 exit(1);
139 }
140
141 m->nrows = a->nrows;
142 m->ncols = a->ncols + b->ncols;
143
144 m->vals = calloc_double(m->nrows * m->ncols);
145
146 size_t siz = sizeof(double);
147 size_t ii;
148 for (ii = 0; ii < a->nrows; ii++){
149 memmove(m->vals + ii * m->ncols, a->vals + ii * a->ncols, a->ncols * siz);
150 memmove(m->vals + ii * m->ncols + a->ncols, b->vals + ii * b->ncols, b->ncols * siz);
151 }
152 return m;
153 }
154
155 /***********************************************************//**
156 Function block_diag
157
158 Purpose: Create a blog diagonal matrix
159
160 Parameters:
161 - mats (IN) - array of mats
162 - nMats (IN) - number of matrices
163
164 Returns: block diagonal matrix *m*
165 ***************************************************************/
block_diag(const struct mat ** mats,size_t nMats)166 struct mat * block_diag(const struct mat ** mats, size_t nMats){
167
168 struct mat * m;
169 if (NULL == ( m = malloc(sizeof(struct mat)))){
170 fprintf(stderr, "failed to allocate memory.\n");
171 exit(1);
172 }
173 m->nrows = 0;
174 m->ncols = 0;
175 size_t ii;
176 for (ii = 0; ii < nMats; ii++){
177 m->nrows += mats[ii]->nrows;
178 m->ncols += mats[ii]->ncols;
179 }
180 m->vals = calloc_double(m->nrows * m->ncols);
181
182 size_t jj;
183 size_t siz = sizeof(double);
184 size_t row = 0;
185 size_t cols = 0;
186 for (ii = 0; ii < nMats; ii++){
187 for (jj = 0; jj < mats[ii]->nrows; jj++){
188 memmove(m->vals + row * m->ncols + cols,
189 mats[ii]->vals + jj * mats[ii]->ncols,
190 siz * mats[ii]->ncols);
191 row += 1;
192 }
193 cols = cols + mats[ii]->ncols;
194 }
195 return m;
196 }
197
198 /***********************************************************//**
199 Function freemat
200
201 Purpose: free a matrix structure
202
203 Parameters:
204 - matin (IN) - matrix to free
205
206 Returns: Nothing
207 ***************************************************************/
freemat(struct mat * matin)208 void freemat(struct mat * matin){
209 free(matin->vals);
210 free(matin);
211 }
212
213 /***********************************************************//**
214 Function printmat
215
216 Purpose: print a matrix to screen
217
218 Parameters:
219 - matin (IN) - matrix to print
220
221 Returns: Nothing
222 ***************************************************************/
printmat(struct mat * matin)223 void printmat(struct mat * matin){
224
225 size_t ii, jj;
226 for (ii = 0; ii < matin->nrows; ii++){
227 for (jj = 0; jj < matin->ncols; jj++){
228 printf("%f ", matin->vals[ii*matin->ncols + jj]);
229 }
230 printf("\n");
231 }
232 }
233
mat_ones(struct mat * m,size_t nrows,size_t ncols)234 void mat_ones(struct mat * m, size_t nrows, size_t ncols)
235 {
236 size_t ii, jj;
237 m->nrows = nrows;
238 m->ncols = ncols;
239 m->vals = calloc_double(nrows*ncols);
240 for (ii = 0; ii < nrows; ii++){
241 for (jj = 0; jj < ncols; jj++){
242 m->vals[ii*ncols + jj] = 1.0;
243 }
244 }
245 }
246
247
248