1 /* matutil.c
2
3 Written by Don Robert Maszle
4 18 September 1992
5
6 Copyright (c) 1992-2017 Free Software Foundation, Inc.
7
8 This file is part of GNU MCSim.
9
10 GNU MCSim is free software; you can redistribute it and/or
11 modify it under the terms of the GNU General Public License
12 as published by the Free Software Foundation; either version 3
13 of the License, or (at your option) any later version.
14
15 GNU MCSim is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
22
23 A bunch of matrix and vector utililties.
24 */
25
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30
31 #include "sim.h"
32
33 /* -----------------------------------------------------------------------------
34 LogTransformArray
35
36 log-transforms the elements of a Src array to a Destination array,
37 which is returned.
38
39 To transform in place, give the same pointer for both args.
40 */
41
LogTransformArray(long nElems,double * rgdSrc,double * rgdDes)42 double *LogTransformArray (long nElems, double *rgdSrc, double *rgdDes)
43 {
44 register long l;
45
46 for (l = 0; l < nElems; l++)
47 rgdDes[l] = log(rgdSrc[l]);
48
49 return rgdDes;
50
51 } /* LogTransformArray */
52
53
54 /* -----------------------------------------------------------------------------
55 InitdVector
56
57 initializes a vector of double.
58
59 The pointer to the vector is returned if no error occurs, otherwise NULL
60 is returned. It is an error to call it with cVectors = 0.
61 */
62
InitdVector(long cVectors)63 double *InitdVector (long cVectors)
64 {
65
66 if (cVectors == 0) {
67 printf ("Error: zero length array allocation in InitdVector - Exiting\n");
68 exit (0);
69 }
70 else
71 return (double *) malloc (cVectors * sizeof(double));
72
73 } /* InitdVector */
74
75
76 /* -----------------------------------------------------------------------------
77 InitiVector
78
79 initializes a vector of integers.
80
81 The pointer to the vector is returned if no error occurs, otherwise NULL
82 is returned. It is an error to call it with cVectors = 0.
83 */
84
InitiVector(long cVectors)85 int *InitiVector (long cVectors)
86 {
87
88 if (cVectors == 0) {
89 printf ("Error: zero length array allocation in InitiVector - Exiting\n");
90 exit (0);
91 }
92 else
93 return (int *) malloc (cVectors * sizeof(int));
94
95 } /* InitiVector */
96
97
98 /* ----------------------------------------------------------------------------
99 InitlVector
100
101 initializes a vector of longs.
102
103 The pointer to the vector is returned if no error occurs, otherwise NULL
104 is returned. It is an error to call it with cVectors = 0.
105 */
106
InitlVector(long cVectors)107 long *InitlVector (long cVectors)
108 {
109
110 if (cVectors == 0) {
111 printf ("Error: zero length array allocation in InitlVector - Exiting\n");
112 exit (0);
113 }
114 else
115 return (long *) malloc (cVectors * sizeof(long));
116
117 } /* InitlVector */
118
119
120 /* ----------------------------------------------------------------------------
121 InitpdVector
122
123 initializes a vector of pointers to doubles.
124
125 The pointer to the vector is returned if no error occurs, otherwise NULL
126 is returned. It is an error to call it with cVectors = 0.
127 */
128
InitpdVector(long cVectors)129 double **InitpdVector (long cVectors)
130 {
131
132 if (cVectors == 0) {
133 printf ("Error: zero length array allocation in InitpdVector - Exiting\n");
134 exit (0);
135 }
136 else
137 return (double **) malloc (cVectors * sizeof(double *));
138
139 } /* InitpdVector */
140
141
142 /* ----------------------------------------------------------------------------
143 InitdMatrix
144
145 initializes a 2 dimensional matrix of doubles.
146
147 The pointer to the vector of vectors is returned if no error occurs,
148 otherwise the NULL value is returned. It is an error to call it with
149 cVectors = 0 or cElemsEach = 0.
150 */
151
InitdMatrix(long cVectors,long cElemsEach)152 double **InitdMatrix (long cVectors, long cElemsEach)
153 {
154 register long i;
155 double **rgp;
156
157 if ((cVectors == 0) || (cElemsEach == 0)) {
158 printf ("Error: zero length array allocation in InitdMatrix - Exiting\n");
159 exit (0);
160 }
161
162 rgp = (double **) malloc (cVectors * sizeof(double *));
163
164 if (rgp) {
165 for (i = 0; i < cVectors; i++) {
166 rgp[i] = (double *) malloc (cElemsEach * sizeof(double));
167 if (!rgp[i]) {
168 rgp = 0;
169 break;
170 } /* if */
171 } /* for */
172 } /* if */
173
174 return (rgp);
175
176 } /* InitdMatrix */
177
178
179 /* ----------------------------------------------------------------------------
180 InitlMatrix
181
182 initializes a 2 dimensional matrix of longs.
183
184 The pointer to the vector of vectors is returned if no error occurs,
185 otherwise the NULL value is returned. It is an error to call it with
186 cVectors = 0 or cElemsEach = 0.
187 */
188
InitlMatrix(long cVectors,long cElemsEach)189 long **InitlMatrix (long cVectors, long cElemsEach)
190 {
191 register long i;
192 long **rgp;
193
194 if ((cVectors == 0) || (cElemsEach == 0)) {
195 printf ("Error: zero length array allocation in InitlMatrix - Exiting\n");
196 exit (0);
197 }
198
199 rgp = (long **) malloc (cVectors * sizeof(long *));
200
201 if (rgp) {
202 for (i = 0; i < cVectors; i++) {
203 rgp[i] = (long *) malloc (cElemsEach * sizeof(long));
204 if (!rgp[i]) {
205 rgp = 0;
206 break;
207 } /* if */
208 } /* for */
209 } /* if */
210
211 return (rgp);
212
213 } /* InitlMatrix */
214
215
216 /* ----------------------------------------------------------------------------
217 ColumnMeans
218 */
ColumnMeans(long cRows,long cCols,double ** x,double * x_bar)219 void ColumnMeans (long cRows, long cCols, double **x, double *x_bar)
220 {
221 register long i, l;
222
223 for (l = 0; l < cCols; l++) x_bar[l] = 0.0;
224
225 for (i = 0; i < cRows; i++)
226 for (l = 0; l < cCols; l++) x_bar[l] += x[i][l];
227
228 for (l = 0; l < cCols; l++) x_bar[l] /= cRows;
229
230 } /* ColumnMean */
231
232
233 /* -----------------------------------------------------------------------------
234 Cholesky
235
236 Does the Cholesky decomposition of a variance matrix:
237 Compute the matrix A such that AAt = VAR.
238
239 Returns 0 if successful, -1 otherwise.
240 The variance matrix is destroyed in the process !
241 */
242
Cholesky(PDOUBLE * prgdVariance,PDOUBLE * prgdComponent,long lNparams)243 int Cholesky (PDOUBLE *prgdVariance, PDOUBLE *prgdComponent, long lNparams)
244 {
245 register int i, j, k;
246 double dSum;
247
248 for (i = 0; i < lNparams; i++)
249 for (j = 0; j < lNparams; j++) prgdComponent[i][j] = 0.0;
250
251 for (i = 0; i < lNparams; i++)
252 for (j = i; j < lNparams ; j++) {
253 dSum = prgdVariance[i][j];
254 for (k = i - 1; k >= 0 ; k--)
255 dSum = dSum - prgdVariance[i][k] * prgdVariance[j][k];
256
257 if (i == j) {
258 if (dSum <= 0.0) {
259 printf ("Warning: input matrix for Cholesky is not "
260 "positive definite\n");
261 return 0;
262 }
263 else
264 prgdComponent[i][i] = sqrt(dSum);
265 }
266 else prgdVariance[j][i] = dSum / prgdComponent[i][i];
267 } /* fin de for j */
268
269 for (i = 0; i < lNparams ; i++)
270 for (j = i+1; j < lNparams ; j++)
271 prgdComponent[j][i] = prgdVariance[j][i];
272
273 /* success */
274 return 1;
275
276 } /* Cholesky */
277
278
279 /* End */
280