~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/matvec/mat.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /*
  2  *
  3  * mat.c
  4  *
  5  * Dynamically allocated matrix of doubles.
  6  *
  7  */
  8 
  9 #include <tina/sys.h>
 10 #include <tina/sysfuncs.h>
 11 #include <tina/math.h>
 12 #include <tina/mathfuncs.h>
 13 
 14 Mat *mat_make(int m, int n)
 15 {
 16     Mat *a = talloc(Mat);
 17     a->m = m;
 18     a->n = n;
 19     a->el = tarray_alloc(0, 0, m, n, double);
 20     return(a);
 21 }
 22 
 23 void mat_free(Mat *a)
 24 {
 25     if(a == NULL)
 26         return;
 27     tarray_free(a->el, 0, 0, a->m, a->n, double);
 28     rfree(a);
 29 }
 30 
 31 void mat_copy(Mat *a, Mat *b)
 32 {
 33     if(a == NULL)
 34         return;
 35     tarray_copy_inplace(a->el, b->el, 0, 0, a->m, a->n, double);
 36 }
 37 
 38 void mat_zero(Mat *a)
 39 {
 40     int i, j;
 41     if(a == NULL)
 42         return;
 43     for(i = 0; i < a->m; i++)
 44         for(j = 0; j < a->n; j++)
 45             a->el[i][j] = 0;
 46 }
 47 
 48 void mat_unit(Mat *a)
 49 {
 50     int i, j;
 51     if(a == NULL)
 52         return;
 53     for(i = 0; i < a->m; i++)
 54         for(j = 0; j < a->n; j++)
 55             a->el[i][j] = (i == j)? 1: 0;
 56 }
 57 
 58 void mat_rand_unif(Mat *a, double p, double q)
 59 {
 60     int i, j;
 61     if(a == NULL)
 62         return;
 63     for(i = 0; i < a->m; i++)
 64         for(j = 0; j < a->n; j++)
 65             a->el[i][j] = rand_unif(p, q);
 66 }
 67 
 68 void mat_rand_normal(Mat *a, double m, double s)
 69 {
 70     int i, j;
 71     if(a == NULL)
 72         return;
 73     for(i = 0; i < a->m; i++)
 74         for(j = 0; j < a->n; j++)
 75             a->el[i][j] = rand_normal(m, s);
 76 }
 77 
 78 void mat_transp(Mat *at, Mat *a)
 79 {
 80     int i, j;
 81     if(at == NULL || a == NULL)
 82         return;
 83     for(i = 0; i < a->m; i++)
 84         for(j = 0; j < a->n; j++)
 85             at->el[j][i] = a->el[i][j];
 86 }
 87 
 88 void mat_stransp(Mat *a)
 89 {
 90     int i, j;
 91     if(a == NULL)
 92         return;
 93     for(i = 0; i < a->m; i++)
 94         for(j = 0; j < i; j++)
 95             SWAP(double, a->el[i][j], a->el[j][i]);
 96 }
 97 
 98 void mat_prod(Mat *ab, Mat *a, Mat *b)
 99 {
100     int i, j, k;
101     if(ab == NULL || a == NULL || b == NULL)
102         return;
103     for(i = 0; i < a->m; i++)
104         for(j = 0; j < b->n; j++)
105         {
106             double sum = 0;
107             for(k = 0; k < a->n; k++)
108                 sum += a->el[i][k]*b->el[k][j];
109             ab->el[i][j] = sum;
110         }
111 }
112 
113 void mat_sqr(Mat *aa, Mat *a)
114 {
115     int i, j, k;
116     if(aa == NULL || a == NULL)
117         return;
118     for(i = 0; i < a->m; i++)
119         for(j = 0; j < a->m; j++)
120         {
121             double sum = 0;
122             for(k = 0; k < a->n; k++)
123                 sum += a->el[i][k]*a->el[j][k];
124             aa->el[i][j] = sum;
125         }
126 }
127 
128 void mat_vprod(Vec *av, Mat *a, Vec *v)
129 {
130     int i, j;
131     if(av == NULL || a == NULL || v == NULL)
132         return;
133     for(i = 0; i < a->m; i++)
134     {
135         double sum = 0;
136         for(j = 0; j < a->n; j++)
137             sum += a->el[i][j]*v->el[j];
138         av->el[i] = sum;
139     }
140 }
141 
142 void mat_dprod(Vec *d, Mat *a)
143 {
144     int i, j;
145     if(a == NULL || d == NULL)
146         return;
147     for(i = 0; i < a->m; i++)
148         for(j = 0; j < a->n; j++)
149             a->el[i][j] *= d->el[i];
150 }
151 
152 double mat_sprod(Vec *v, Mat *a, Vec *w)
153 {
154     int i, j;
155     double sum = 0;
156     if(a == NULL || v == NULL)
157         return(0);
158     for(i = 0; i < a->m; i++)
159         for(j = 0; j < a->n; j++)
160             sum += v->el[i]*a->el[i][j]*w->el[j];
161     return(sum);
162 }
163 
164 void mat_sum(Mat *a, Mat *b)
165 {
166     int i, j;
167     if(a == NULL || b == NULL)
168         return;
169     for(i = 0; i < a->m; i++)
170         for(j = 0; j < a->n; j++)
171             a->el[i][j] += b->el[i][j];
172 }
173 
174 void mat_diff(Mat *a, Mat *b)
175 {
176     int i, j;
177     if(a == NULL || b == NULL)
178         return;
179     for(i = 0; i < a->m; i++)
180         for(j = 0; j < a->n; j++)
181             a->el[i][j] -= b->el[i][j];
182 }
183 
184 void mat_times(double k, Mat *a)
185 {
186     int i, j;
187     if(a == NULL)
188         return;
189     for(i = 0; i < a->m; i++)
190         for(j = 0; j < a->n; j++)
191             a->el[i][j] *= k;
192 }
193 
194 void mat_minus(Mat *a)
195 {
196     int i, j;
197     if(a == NULL)
198         return;
199     for(i = 0; i < a->m; i++)
200         for(j = 0; j < a->n; j++)
201             a->el[i][j] = -a->el[i][j];
202 }
203 
204 void mat_accum(Mat *a, double k, Mat *b)
205 {
206     int i, j;
207     if(a == NULL || b == NULL)
208         return;
209     for(i = 0; i < a->m; i++)
210         for(j = 0; j < a->n; j++)
211             a->el[i][j] += k*b->el[i][j];
212 }
213 
214 void mat_sum_tensor(Mat *a, Vec *v, Vec *w)
215 {
216     int i, j;
217     if(a == NULL || v == NULL || w == NULL)
218         return;
219     for(i = 0; i < a->m; i++)
220         for(j = 0; j < a->n; j++)
221             a->el[i][j] += v->el[i]*w->el[j];
222 }
223 
224 void mat_accum_tensor(Mat *a, double k, Vec *v, Vec *w)
225 {
226     int i, j;
227     if(a == NULL || v == NULL || w == NULL)
228         return;
229     for(i = 0; i < a->m; i++)
230         for(j = 0; j < a->n; j++)
231             a->el[i][j] += k*v->el[i]*w->el[j];
232 }
233 
234 void mat_print(FILE *fp, char *fmt, Mat *a)
235 {
236     int i, j;
237     if(a == NULL)
238     {
239         fprintf(fp, "0 0 :\n");
240         return;
241     }
242     fprintf(fp, "%d %d :\n", a->m, a->n);
243     for(i = 0; i < a->m; i++)
244     {
245         for(j = 0; j < a->n; j++)
246             fprintf(fp, fmt, a->el[i][j]);
247         fprintf(fp, "\n");
248     }
249 }
250 
251 Mat *mat_read(FILE *fp)
252 {
253     Mat *a;
254     int i, j, m, n;
255     fscanf(fp, "%d %d %*s", &m, &n);
256     a = mat_make(m, n);
257     for(i = 0; i < a->m; i++)
258         for(j = 0; j < a->n; j++)
259             fscanf(fp, "%lf", &a->el[i][j]);
260     return(a);
261 }
262 
263 void mat_row(Vec *row, Mat *a, int i)
264 {
265     int j;
266     if(row == NULL || a == NULL)
267         return;
268     for(j = 0; j < a->n; j++)
269         row->el[j] = a->el[i][j];
270 }
271 
272 void mat_row_get(Mat *a, int i, Vec *row)
273 {
274     int j;
275     if(row == NULL || a == NULL)
276         return;
277     for(j = 0; j < row->n; j++)
278         row->el[j] = a->el[i][j];
279 }
280 
281 void mat_row_set(Mat *a, int i, Vec *row)
282 {
283     int j;
284     if(row == NULL || a == NULL)
285         return;
286     for(j = 0; j < row->n; j++)
287         a->el[i][j] = row->el[j];
288 }
289 
290 void mat_col_get(Mat *a, int j, Vec *col)
291 {
292     int i;
293     if(col == NULL || a == NULL)
294         return;
295     for(i = 0; i < col->n; i++)
296         col->el[i] = a->el[i][j];
297 }
298 
299 void mat_col_set(Mat *a, int j, Vec *col)
300 {
301     int i;
302     if(col == NULL || a == NULL)
303         return;
304     for(i = 0; i < col->n; i++)
305         a->el[i][j] = col->el[i];
306 }
307 
308 void mat_block_get(Mat *a, int li, int lj, Mat *b)
309 {
310     int i, j;
311     for(i = 0; i < b->m; i++)
312         for(j = 0; j < b->n; j++)
313             b->el[i][j] = a->el[i+li][j+lj];
314 }
315 
316 void mat_block_set(Mat *a, int li, int lj, Mat *b)
317 {
318     int i, j;
319     for(i = 0; i < b->m; i++)
320         for(j = 0; j < b->n; j++)
321             a->el[i+li][j+lj] = b->el[i][j];
322 }
323 
324 void mat_index_get(Mat *a, Ivec *indi,  Ivec *indj, Mat *b)
325 {
326     int i, j;
327     for(i = 0; i < b->m; i++)
328         for(j = 0; j < b->n; j++)
329             b->el[i][j] = a->el[indi->el[i]][indj->el[j]];
330 }
331 
332 void mat_index_set(Mat *a, Ivec *indi,  Ivec *indj, Mat *b)
333 {
334     int i, j;
335     for(i = 0; i < b->m; i++)
336         for(j = 0; j < b->n; j++)
337             a->el[indi->el[i]][indj->el[j]] = b->el[i][j];
338 }
339 
340 void mat_mat3(Mat *a, Mat3 a3)
341 {
342     a->el[0][0] = mat3_xx(a3);
343     a->el[0][1] = mat3_xy(a3);
344     a->el[0][2] = mat3_xz(a3);
345     a->el[1][0] = mat3_yx(a3);
346     a->el[1][1] = mat3_yy(a3);
347     a->el[1][2] = mat3_yz(a3);
348     a->el[2][0] = mat3_zx(a3);
349     a->el[2][1] = mat3_zy(a3);
350     a->el[2][2] = mat3_zz(a3);
351 }
352 
353 Mat3 mat3_mat(Mat *a)
354 {
355     Mat3 a3;
356     mat3_xx(a3) = a->el[0][0];
357     mat3_xy(a3) = a->el[0][1];
358     mat3_xz(a3) = a->el[0][2];
359     mat3_yx(a3) = a->el[1][0];
360     mat3_yy(a3) = a->el[1][1];
361     mat3_yz(a3) = a->el[1][2];
362     mat3_zx(a3) = a->el[2][0];
363     mat3_zy(a3) = a->el[2][1];
364     mat3_zz(a3) = a->el[2][2];
365     return(a3);
366 }
367 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.