1 /*****************************************************************************
2 *
3 * Elmer, A Finite Element Software for Multiphysical Problems
4 *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library 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 GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library (in file ../LGPL-2.1); if not, write
19 * to the Free Software Foundation, Inc., 51 Franklin Street,
20 * Fifth Floor, Boston, MA 02110-1301 USA
21 *
22 *****************************************************************************/
23
24 /*******************************************************************************
25 *
26 * MATC matrix utilities.
27 *
28 *******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: CSC - IT Center for Science Ltd.
33 * Keilaranta 14, P.O. BOX 405
34 * 02101 Espoo, Finland
35 * Tel. +358 0 457 2723
36 * Telefax: +358 0 457 2302
37 * EMail: Juha.Ruokolainen@csc.fi
38 *
39 * Date: 30 May 1996
40 *
41 * Modified by:
42 *
43 * Date of modification:
44 *
45 ******************************************************************************/
46 /***********************************************************************
47 |
48 | MATRIX.C - Last Edited 8. 8. 1988
49 |
50 ***********************************************************************/
51
52 /*======================================================================
53 |Syntax of the manual pages:
54 |
55 |FUNCTION NAME(...) params ...
56 |
57 $ usage of the function and type of the parameters
58 ? explane the effects of the function
59 = return value and the type of value if not of type int
60 @ globals effected directly by this routine
61 ! current known bugs or limitations
62 & functions called by this function
63 ~ these functions may interest you as an alternative function or
64 | because they control this function somehow
65 ^=====================================================================*/
66
67
68 /*
69 * $Id: matrix.c,v 1.1.1.1 2005/04/14 13:29:14 vierinen Exp $
70 *
71 * $Log: matrix.c,v $
72 * Revision 1.1.1.1 2005/04/14 13:29:14 vierinen
73 * initial matc automake package
74 *
75 * Revision 1.2 1998/08/01 12:34:50 jpr
76 *
77 * Added Id, started Log.
78 *
79 *
80 */
81
82 #include "elmer/matc.h"
83
84 #define MA(i,j) a[(i) * ncola + (j)]
85 #define MB(i,j) b[(i) * ncolb + (j)]
86 #define MC(i,j) c[(i) * ncolc + (j)]
87
func_abs(arg)88 double func_abs(arg)
89 double arg;
90 {
91 return abs(arg);
92 }
93
func_mod(x,y)94 double func_mod(x,y)
95 double x,y;
96 {
97 int ix, iy;
98
99 ix = x + 0.5;
100 iy = y + 0.5;
101 return (double)(ix % iy);
102 }
103
mtr_sum(A)104 VARIABLE *mtr_sum(A) VARIABLE *A;
105 {
106 VARIABLE *C;
107
108 int i, j;
109
110 int nrowa = NROW(A), ncola = NCOL(A);
111
112 double *a = MATR(A), *c;
113
114
115 if (nrowa == 1 || ncola == 1)
116 {
117 C = var_temp_new(TYPE_DOUBLE, 1, 1); c = MATR(C);
118 nrowa = (nrowa == 1) ? ncola : nrowa;
119 for(i = 0; i < nrowa; i++) *c += *a++;
120 }
121 else
122 {
123 C = var_temp_new(TYPE_DOUBLE, 1, ncola); c = MATR(C);
124 for(i = 0; i < ncola; i++)
125 for(j = 0; j < nrowa; j++) c[i] += MA(j, i);
126 }
127
128 return C;
129 }
130
mtr_trace(A)131 VARIABLE *mtr_trace(A) VARIABLE *A;
132 {
133 VARIABLE *C;
134
135 double temp = 0.0;
136 int i;
137
138 int nrowa = NROW(A), ncola = NCOL(A);
139 double *a = MATR(A);
140
141 if (nrowa != ncola) error("trace: not square.\n");
142
143 for(i = 0; i < nrowa; i++) temp += MA(i,i);
144
145 C = var_temp_new(TYPE(A), 1, 1); *MATR(C) = temp;
146
147 return C;
148 }
149
mtr_zeros(A)150 VARIABLE *mtr_zeros(A) VARIABLE *A;
151 {
152 VARIABLE *C;
153
154 int ind1 = 1, ind2 = 1;
155
156 if (NEXT(A) != NULL)
157 {
158 ind1 = (int)*MATR(A); ind2 = (int)*MATR(NEXT(A));
159 }
160 else
161 {
162 ind2 = (int)*MATR(A);
163 }
164
165 if (ind1 < 1 || ind2 < 1)
166 error("Zeros: invalid size for and array");
167
168 C = var_temp_new(TYPE_DOUBLE, ind1, ind2);
169
170 return C;
171 }
172
mtr_ones(A)173 VARIABLE *mtr_ones(A) VARIABLE *A;
174 {
175 VARIABLE *C;
176 double *c;
177 int i, n;
178
179 C = mtr_zeros(A); c = MATR(C);
180 n = NROW(C) * NCOL(C);
181
182 for(i = 0; i < n; i++) *c++ = 1.0;
183
184 return C;
185 }
186
mtr_rand(A)187 VARIABLE *mtr_rand(A) VARIABLE *A;
188 {
189 VARIABLE *C;
190
191 static int seed = 0;
192 #pragma omp threadprivate (seed)
193 int i, n;
194
195 double *c;
196
197 C = mtr_zeros(A); c = MATR(C);
198 n = NROW(C) * NCOL(C);
199
200 if (seed == 0) seed = time(NULL);
201 for(i = 0; i < n; i++) *c++ = urand(&seed);
202
203 return C;
204 }
205
mtr_resize(A)206 VARIABLE *mtr_resize(A) VARIABLE *A;
207 {
208 VARIABLE *C;
209
210 int i, j, n, m, ind1 = 1, ind2;
211
212 double *a = MATR(A), *c;
213
214 if (NEXT(NEXT(A)) != NULL)
215 {
216 ind1 = *MATR(NEXT(A)); ind2 = *MATR(NEXT(NEXT(A)));
217 }
218 else
219 {
220 ind2 = (int)*MATR(NEXT(A));;
221 }
222
223 if (ind1 < 1 || ind2 < 1)
224 error("resize: invalid size for and array");
225
226 C = var_temp_new(TYPE(A), ind1, ind2); c = MATR(C);
227 a = MATR(A); n = ind1 * ind2; m = NROW(A) * NCOL(A);
228
229 for(i = j = 0; i < n; i++)
230 {
231 *c++ = a[j++]; if (j == m) j = 0;
232 }
233
234 return C;
235 }
236
mtr_vector(A)237 VARIABLE *mtr_vector(A) VARIABLE *A;
238 {
239 VARIABLE *C;
240
241 double start, stop, incr, x, *c;
242
243 int i, eval;
244
245 start = *MATR(A); stop = *MATR(NEXT(A));
246
247 if (NEXT(NEXT(A)) != (VARIABLE *)NULL)
248 incr = *MATR(NEXT(NEXT(A)));
249 else
250 incr = (start < stop) ? (1) : (-1);
251
252 if (incr == 0)
253 incr = (start < stop) ? (1) : (-1);
254
255 eval = (int)(abs(stop-start) / abs(incr)) + 1;
256 if (eval < 1) return NULL;
257
258 C = var_temp_new(TYPE_DOUBLE, 1, eval); c = MATR(C);
259
260 x = start;
261 for(i = 0; i < eval; i++)
262 {
263 *c++ = x; x += incr;
264 }
265
266 return C;
267 }
268
mtr_eye(A)269 VARIABLE *mtr_eye(A) VARIABLE *A;
270 {
271 VARIABLE *C;
272 double *c;
273 int i, ind, ncolc;
274
275 if (*MATR(A) < 1)
276 {
277 error("eye: Invalid size for an array.\n");
278 }
279
280 ind = (int)*MATR(A);
281
282 C = var_temp_new(TYPE_DOUBLE, ind, ind);
283 c = MATR(C); ncolc = ind;
284
285 for (i = 0; i < ind; i++) MC(i,i) = 1.0;
286
287 return C;
288 }
289
mtr_size(A)290 VARIABLE *mtr_size(A) VARIABLE *A;
291 {
292 VARIABLE *C;
293 double *c;
294
295 C = var_temp_new(TYPE_DOUBLE, 1, 2); c = MATR(C);
296 *c++ = NROW(A); *c = NCOL(A);
297
298 return C;
299 }
300
mtr_min(A)301 VARIABLE *mtr_min(A) VARIABLE *A;
302 {
303 VARIABLE *C;
304
305 double *a = MATR(A), *c;
306
307 int nrowa = NROW(A), ncola = NCOL(A);
308 int i, j;
309
310 if (nrowa == 1 || ncola == 1)
311 {
312 C = var_temp_new(TYPE_DOUBLE, 1, 1); c = MATR(C);
313 *c = *a++; nrowa = max(ncola, nrowa);
314 for(i = 1; i < nrowa; i++, a++) *c = min(*c, *a);
315 }
316 else
317 {
318 C = var_temp_new(TYPE_DOUBLE, 1, ncola); c = MATR(C);
319 for(i = 0; i < ncola; i++, c++)
320 {
321 *c = MA(0, i);
322 for(j = 1; j < nrowa; j++) *c = min(*c, MA(j, i));
323 }
324 }
325
326 return C;
327 }
328
mtr_max(A)329 VARIABLE *mtr_max(A) VARIABLE *A;
330 {
331 VARIABLE *C;
332
333 double *a = MATR(A), *c;
334
335 int nrowa = NROW(A), ncola = NCOL(A);
336 int i, j;
337
338 if (nrowa == 1 || ncola == 1)
339 {
340 C = var_temp_new(TYPE_DOUBLE, 1, 1); c = MATR(C);
341 *c = *a++; nrowa = max(ncola, nrowa);
342 for(i = 1; i < nrowa; i++, a++) *c = max(*c, *a);
343 }
344 else
345 {
346 C = var_temp_new(TYPE_DOUBLE, 1, ncola); c = MATR(C);
347 for(i = 0; i < ncola; i++, c++)
348 {
349 *c = MA(0, i);
350 for(j = 1; j < nrowa; j++) *c = max(*c, MA(j, i));
351 }
352 }
353
354 return C;
355 }
356
mtr_diag(A)357 VARIABLE *mtr_diag(A) VARIABLE *A;
358 {
359 VARIABLE *C;
360
361 double *a = MATR(A), *c;
362
363 int nrowa = NROW(A), ncola = NCOL(A);
364 int ncolc;
365
366 int i;
367
368 if (nrowa == 1 || ncola == 1)
369 {
370 nrowa = max(nrowa, ncola); ncolc = nrowa;
371 C = var_temp_new(TYPE_DOUBLE, nrowa, nrowa); c = MATR(C);
372 for(i = 0; i < nrowa; i++) MC(i, i) = *a++;
373 }
374 else
375 {
376 C = var_temp_new(TYPE_DOUBLE, 1, nrowa); c = MATR(C);
377 for(i = 0; i < min(nrowa,ncola); i++) *c++ = MA(i, i);
378 }
379
380 return C;
381 }
382
mtr_pow(A)383 VARIABLE *mtr_pow(A) VARIABLE *A;
384 {
385 VARIABLE *B = NEXT(A), *C;
386 double *a = MATR(A), b = M(B,0,0), *c;
387
388 int nrowa = NROW(A), ncola = NCOL(A);
389
390 int i;
391
392 C = var_temp_new(TYPE_DOUBLE, nrowa,ncola );
393 c = MATR( C );
394 for(i = 0; i < nrowa*ncola; i++) *c++ = pow(*a++,b);
395
396 return C;
397 }
398
mtr_where(A)399 VARIABLE *mtr_where(A) VARIABLE *A;
400 {
401 VARIABLE *C;
402
403 double *a = MATR(A), *c;
404
405 int nrowa = NROW(A), ncola = NCOL(A);
406
407 int i,n=0;
408
409 for( i=0; i < nrowa*ncola; i++) if ( a[i] ) n++;
410
411 C = var_temp_new( TYPE_DOUBLE,1,n );
412 c = MATR(C);
413
414 for( i=0; i < nrowa*ncola; i++ ) if ( a[i] ) { *c++ = i; }
415
416 return C;
417 }
418
mtr_com_init()419 void mtr_com_init()
420 {
421 static char *minHelp =
422 {
423 "r = min(matrix)\n"
424 "Return value is a vector containing smallest element in columns of given matrix.\n"
425 "r=min(min(matrix) gives smallest element of the matrix.\n\n"
426 };
427
428 static char *maxHelp =
429 {
430 "r = max(matrix)\n"
431 "Return value is a vector containing largest element in columns of given matrix.\n"
432 "r=max(max(matrix)) gives largest element of the matrix.\n\n"
433 };
434
435 static char *sumHelp =
436 {
437 "r = sum(matrix)\n"
438 "Return vector is column sums of given matrix. r=sum(sum(matrix)) gives\n"
439 "the total sum of elements of the matrix.\n\n"
440 };
441
442 static char *traceHelp =
443 {
444 "r = trace(matrix)\n"
445 "Return value is sum of matrix diagonal elements.\n\n"
446 };
447
448 static char *detHelp =
449 {
450 "r = det(matrix)\n"
451 "Return value is determinant of given square matrix.\n\n"
452 };
453
454 static char *invHelp =
455 {
456 "r = inv(matrix)\n"
457 "Invert given square matrix. Computed also by r=matrix^(-1).\n\n"
458 };
459
460 static char *eigHelp =
461 {
462 "r = eig(matrix)\n"
463 "Return eigenvalues of given square matrix. r(n,0) is real part of the\n"
464 "n:th eigenvalue, r(n,1) is the imaginary part respectively\n\n"
465 };
466
467 static char *jacobHelp =
468 {
469 "r = jacob(a,b,eps)\n"
470 "Solve symmetric positive definite eigenvalue problem by Jacob iteration.\n"
471 "Return values are the eigenvalues. Also a variable eigv is created containing\n"
472 "eigenvectors.\n\n"
473 };
474
475 static char *ludHelp =
476 {
477 "r = lud(matrix)\n"
478 "Return value is lud decomposition of given matrix.\n\n"
479 };
480
481 static char *hesseHelp =
482 {
483 "r = hesse(matrix)\n"
484 "Return the upper hessenberg form of given matrix.\n\n"
485 };
486
487 static char *eyeHelp =
488 {
489 "r = eye(n)\n"
490 "Return n by n identity matrix.\n\n"
491 };
492
493 static char *zerosHelp =
494 {
495 "r = zeros(n,m)\n"
496 "Return n by m matrix with elements initialized to zero.\n"
497 };
498
499 static char *onesHelp =
500 {
501 "r = ones(n,m)\n"
502 "Return n by m matrix with elements initialized to one.\n"
503 };
504
505 static char *randHelp =
506 {
507 "r = rand(n,m)\n"
508 "Return n by m matrix with elements initialized to with random number from\n"
509 "zero to one.\n\n"
510 };
511
512 static char *diagHelp =
513 {
514 "r=diag(matrix) or r=diag(vector)\n"
515 "Given matrix return diagonal entries as a vector. Given vector return matrix\n"
516 "with diagonal elements from vector. r=diag(diag(a)) gives matrix with diagonal\n"
517 "elements from matrix a otherwise elements are zero.\n\n"
518 };
519
520 static char *vectorHelp =
521 {
522 "r=vector(start,end,inc)\n"
523 "Return vector of values going from start to end by inc.\n\n"
524 };
525
526 static char *sizeHelp =
527 {
528 "r = size(matrix)\n"
529 "Return size of given matrix.\n"
530 };
531
532 static char *resizeHelp =
533 {
534 "r = resize(matrix,n,m)\n"
535 "Make a matrix to look as a n by m matrix.\n\n"
536 };
537
538 static char *whereHelp =
539 {
540 "r = where(l)\n"
541 "Return linear indexes of where l is true.\n\n"
542 };
543
544 com_init( "sin" , TRUE, TRUE, (VARIABLE *(*)())sin , 1, 1, "r=sin(x)" );
545 com_init( "cos" , TRUE, TRUE, (VARIABLE *(*)())cos , 1, 1, "r=cos(x)" );
546 com_init( "tan" , TRUE, TRUE, (VARIABLE *(*)())tan , 1, 1, "r=tan(x)" );
547 com_init( "asin" , TRUE, TRUE, (VARIABLE *(*)())asin , 1, 1, "r=asin(x)" );
548 com_init( "acos" , TRUE, TRUE, (VARIABLE *(*)())acos , 1, 1, "r=acos(x)" );
549 com_init( "atan" , TRUE, TRUE, (VARIABLE *(*)())atan , 1, 1, "r=atan(x)" );
550 com_init( "atan2" , TRUE, TRUE, (VARIABLE *(*)())atan2 , 2, 2, "r=atan2(y,x)" );
551 com_init( "sinh" , TRUE, TRUE, (VARIABLE *(*)())sinh , 1, 1, "r=sinh(x)" );
552 com_init( "cosh" , TRUE, TRUE, (VARIABLE *(*)())cosh , 1, 1, "r=cosh(x)" );
553 com_init( "tanh" , TRUE, TRUE, (VARIABLE *(*)())tanh , 1, 1, "r=tanh(x)" );
554 com_init( "exp" , TRUE, TRUE, (VARIABLE *(*)())exp , 1, 1, "r=exp(x)" );
555 com_init( "ln" , TRUE, TRUE, (VARIABLE *(*)())log , 1, 1, "r=ln(x)\nNatural logarithm." );
556 com_init( "log" , TRUE, TRUE, (VARIABLE *(*)())log10 , 1, 1, "r=log(x)\nBase 10 logarithm." );
557 com_init( "sqrt" , TRUE, TRUE, (VARIABLE *(*)())sqrt , 1, 1, "r=sqrt(x)" );
558 com_init( "ceil" , TRUE, TRUE, (VARIABLE *(*)())ceil , 1, 1, "r=ceil(x)\nSmallest integer not less than x." );
559 com_init( "floor" , TRUE, TRUE, (VARIABLE *(*)())floor , 1, 1, "r=floor(x)\nLargest integer not more than x." );
560 com_init( "abs" , TRUE, TRUE, (VARIABLE *(*)())func_abs , 1, 1,"r=abs(x)");
561 com_init( "mod" , TRUE, TRUE, (VARIABLE *(*)())func_mod , 2, 2,"r=mod(x,y)");
562 com_init( "pow" , FALSE, TRUE, mtr_pow, 2, 2, "r=pow(x,y)" );
563 com_init( "min" , FALSE, TRUE, mtr_min, 1, 1, minHelp );
564 com_init( "max" , FALSE, TRUE, mtr_max, 1, 1, maxHelp );
565 com_init( "sum" , FALSE, TRUE, mtr_sum, 1, 1, sumHelp );
566 com_init( "trace" , FALSE, TRUE, mtr_trace, 1, 1, traceHelp );
567 com_init( "det" , FALSE, TRUE, mtr_det, 1, 1, detHelp );
568 com_init( "inv" , FALSE, TRUE, mtr_inv, 1, 1, invHelp );
569 com_init( "eig" , FALSE, TRUE, mtr_eig, 1, 1, eigHelp );
570 com_init( "jacob" , FALSE, TRUE, mtr_jacob, 3, 3, jacobHelp );
571 com_init( "lud" , FALSE, TRUE, mtr_LUD, 1, 1, ludHelp );
572 com_init( "hesse" , FALSE, TRUE, mtr_hesse, 1, 1, hesseHelp );
573 com_init( "eye" , FALSE, TRUE, mtr_eye, 1, 1, eyeHelp );
574 com_init( "zeros" , FALSE, TRUE, mtr_zeros, 1, 2, zerosHelp );
575 com_init( "ones" , FALSE, TRUE, mtr_ones, 1, 2, onesHelp );
576 com_init( "rand" , FALSE, FALSE, mtr_rand, 1, 2, randHelp );
577 com_init( "diag" , FALSE, TRUE, mtr_diag, 1, 1, diagHelp );
578 com_init( "vector" , FALSE, TRUE, mtr_vector, 2, 3, vectorHelp );
579 com_init( "size" , FALSE, TRUE, mtr_size, 1, 1, sizeHelp );
580 com_init( "resize" , FALSE, TRUE, mtr_resize, 2, 3, resizeHelp );
581 com_init( "where" , FALSE, FALSE, mtr_where, 1, 1, whereHelp );
582 }
583