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 operator routines.
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 | Oper.C - Last Edited 15. 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 intrest you as an alternative function or
64 | because they control this function somehow
65 =====================================================================*/
66
67
68 /*
69 * $Id: oper.c,v 1.1.1.1 2005/04/14 13:29:14 vierinen Exp $
70 *
71 * $Log: oper.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:51 jpr
76 *
77 * Added Id, started Log.
78 *
79 *
80 */
81
82 #include "elmer/matc.h"
83
84 #ifdef TYPE
85 #undef TYPE
86 #endif
87 #ifdef NROW
88 #undef NROW
89 #endif
90 #ifdef NCOL
91 #undef NCOL
92 #endif
93 #ifdef MATR
94 #undef MATR
95 #endif
96 #ifdef M
97 #undef M
98 #endif
99 #ifdef MATSIZE
100 #undef MATSIZE
101 #endif
102
103 #define TYPE(mat) (mat)->type
104 #define NROW(mat) (mat)->nrow
105 #define NCOL(mat) (mat)->ncol
106 #define MATR(mat) (mat)->data
107 #define MATSIZE(mat) (NROW(mat) * NCOL(mat) * sizeof(double))
108
109 #define MA(i,j) a[ncola * (i) + (j)]
110 #define MB(i,j) b[ncolb * (i) + (j)]
111 #define MC(i,j) c[ncolc * (i) + (j)]
112
mat_new(type,nrow,ncol)113 MATRIX *mat_new(type, nrow, ncol) int type, nrow, ncol;
114 {
115 MATRIX *res;
116
117 res = (MATRIX *)ALLOCMEM(MATRIXSIZE);
118 TYPE(res) = type;
119 NROW(res) = nrow;
120 NCOL(res) = ncol;
121 MATR(res) = (double *)ALLOCMEM(MATSIZE(res));
122
123 return res;
124 }
125
mat_copy(mat)126 MATRIX *mat_copy(mat) MATRIX *mat;
127 {
128 MATRIX *res;
129
130 if (mat == (MATRIX *)NULL) return NULL;
131
132 res = mat_new(TYPE(mat), NROW(mat), NCOL(mat));
133 memcpy((char *)MATR(res), (char *)MATR(mat), MATSIZE(mat));
134
135 return res;
136 }
137
mat_free(mat)138 void mat_free(mat) MATRIX *mat;
139 {
140 if (mat == (MATRIX *)NULL) return;
141
142 FREEMEM((char *)MATR(mat));
143 FREEMEM((char *)mat);
144 }
145
opr_vector(A,B)146 MATRIX *opr_vector(A, B)
147 MATRIX *A, *B;
148 {
149 MATRIX *C;
150
151 int i, n, inc;
152
153 double *a = MATR(A), *b = MATR(B), *c;
154
155 inc = ((int)*b > (int)*a) ? 1:(-1);
156 n = abs((int)*b - (int)*a) + 1;
157
158 C = mat_new(TYPE_DOUBLE, 1, n);
159 c = MATR(C);
160
161 for(i = 0; i < n; i++)
162 *c++ = (int)*a + i*inc;
163
164 return C;
165 }
166
opr_resize(A,B)167 MATRIX *opr_resize(A, B)
168 MATRIX *A, *B;
169 {
170 MATRIX *C;
171
172 double *a = MATR(A), *b = MATR(B), *c;
173 int i, j, n, m;
174
175 if (NCOL(B) >= 2)
176 {
177 i = *b++; j = *b++;
178 }
179 else
180 {
181 i = 1; j = *b;
182 }
183
184 if (i < 1 || j < 1)
185 error("resize: invalid size for and array");
186
187 C = mat_new(TYPE(A), i, j);
188 c = MATR(C);
189
190 n = i * j; m = NROW(A) * NCOL(A);
191
192 for(i = j = 0; i < n; i++)
193 {
194 *c++ = a[j++];
195 if (j == m) j = 0;
196 }
197
198 return C;
199 }
200
opr_apply(A)201 MATRIX *opr_apply(A)
202 MATRIX *A;
203 {
204 VARIABLE *store, *ptr;
205 MATRIX *C = NULL;
206
207 store = var_temp_new(TYPE_STRING, NROW(A), NCOL(A));
208 REFCNT(store) = 0;
209 mat_free(store->this);
210 store->this = A;
211 REFCNT(store)++;
212
213 ptr = (VARIABLE *)com_apply(store);
214
215 var_delete_temp(store);
216
217 if ( ptr ) C = mat_copy(ptr->this);
218
219 return C;
220 }
221
opr_add(A,B)222 MATRIX *opr_add(A, B)
223 MATRIX *A, *B;
224 {
225 MATRIX *C;
226
227 int i;
228
229 int nrowa = NROW(A), ncola = NCOL(A);
230 int nrowb = NROW(B), ncolb = NCOL(B);
231
232 double *a = MATR(A), *b = MATR(B), *c;
233
234 double value;
235
236 if (nrowa == nrowb && ncola == ncolb)
237 {
238 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
239 nrowa *= ncola;
240 for(i = 0; i < nrowa; i++) *c++ = *a++ + *b++;
241 }
242
243 else if (nrowa == 1 && ncola == 1)
244 {
245 C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
246 value = *a; nrowb *= ncolb;
247 for(i = 0; i < nrowb; i++) *c++ = value + *b++;
248 }
249
250 else if (nrowb == 1 && ncolb == 1)
251 {
252 C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);
253 value = *b; nrowa *= ncola;
254 for(i = 0; i < nrowa; i++) *c++ = value + *a++;
255 }
256
257 else
258 error("Add: Incompatible for addition.\n");
259
260 return C;
261 }
262
opr_minus(A)263 MATRIX *opr_minus(A)
264 MATRIX *A;
265 {
266 MATRIX *C;
267 int i;
268
269 int nrowa = NROW(A), ncola = NCOL(A);
270
271 double *a = MATR(A), *c;
272
273 C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);
274
275 nrowa *= ncola;
276 for(i = 0; i < nrowa; i++) *c++ = -*a++;
277
278 return C;
279 }
280
opr_subs(A,B)281 MATRIX *opr_subs(A, B)
282 MATRIX *A, *B;
283 {
284 MATRIX *C;
285
286 double value;
287
288 int i;
289
290 int nrowa = NROW(A), ncola = NCOL(A);
291 int nrowb = NROW(B), ncolb = NCOL(B);
292
293 double *a = MATR(A), *b = MATR(B), *c;
294
295 if (nrowa == nrowb && ncola == ncolb)
296 {
297 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
298 nrowa *= ncola;
299 for(i = 0; i < nrowa; i++) *c++ = *a++ - *b++;
300 }
301 else if (nrowa == 1 && ncola == 1)
302 {
303 C = mat_new(TYPE(B),nrowb, ncolb); c = MATR(C);
304 value = *a; nrowb *= ncolb;
305 for(i = 0; i < nrowb; i++) *c++ = value - *b++;
306 }
307 else if (nrowb == 1 && ncolb == 1)
308 {
309 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
310 value = *b; nrowa *= ncola;
311 for(i = 0; i < nrowa; i++) *c++ = *a++ - value;
312 }
313 else
314 error("Substr: Incompatible for addition.\n");
315
316 return C;
317 }
318
opr_mul(A,B)319 MATRIX *opr_mul(A, B)
320 MATRIX *A, *B;
321 {
322 MATRIX *C;
323
324 double value,s;
325 int i, j, k;
326
327 int nrowa = NROW(A), ncola = NCOL(A);
328 int nrowb = NROW(B), ncolb = NCOL(B);
329
330 double *a = MATR(A), *b = MATR(B), *c;
331
332 if (nrowa == 1 && ncola == 1)
333 {
334 C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
335 value = *a; nrowb *= ncolb;
336 for(i = 0; i < nrowb; i++) *c++ = value * *b++;
337 }
338 else if (nrowb == 1 && ncolb == 1)
339 {
340 C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);
341 value = *b; nrowa *= ncola;
342 for(i = 0; i < nrowa; i++) *c++ = value * *a++;
343 }
344 else if (ncola == nrowb)
345 {
346 C = mat_new(TYPE(A), nrowa,ncolb);
347 c = MATR(C);
348 for(i = 0; i < nrowa; i++)
349 {
350 for(j = 0; j < ncolb; j++)
351 {
352 s = 0;
353 for(k = 0; k < ncola; k++) s += a[k] * MB(k,j);
354 *c++ = s;
355 }
356 a += ncola;
357 }
358 }
359 else if ( ncola == ncolb && nrowa == nrowb )
360 {
361 C = mat_new(TYPE(A), nrowa,ncolb);
362 c = MATR(C);
363 k = 0;
364 for( i = 0; i < nrowa; i++ )
365 for( j = 0; j < ncolb; j++,k++ ) c[k] = a[k] * b[k];
366 }
367 else
368 {
369 error("Mul: Incompatible for multiplication.\n");
370 }
371
372 return C;
373 }
374
opr_pmul(A,B)375 MATRIX *opr_pmul(A, B)
376 MATRIX *A, *B;
377 {
378 MATRIX *C;
379
380 double value;
381 int i;
382
383 int nrowa = NROW(A), ncola = NCOL(A);
384 int nrowb = NROW(B), ncolb = NCOL(B);
385
386 double *a = MATR(A), *b = MATR(B), *c;
387
388 if (nrowa == nrowb && ncola == ncolb)
389 {
390 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
391 nrowa *= ncola;
392 for(i = 0; i < nrowa; i++) *c++ = *a++ * *b++;
393 }
394 else if (nrowa == 1 && ncola == 1)
395 {
396 C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
397 value = *a; nrowb *= ncolb;
398 for(i = 0; i < nrowb; i++) *c++ = value * *b++;
399 }
400 else if (nrowb == 1 && ncolb == 1)
401 {
402 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
403 value = *b; nrowa *= ncola;
404 for(i = 0; i < nrowa; i++) *c++ = *a++ * value;
405 }
406 else
407 error("PMul: Incompatible for pointwise multiplication.\n");
408
409 return C;
410 }
411
opr_div(A,B)412 MATRIX *opr_div(A, B)
413 MATRIX *A, *B;
414 {
415 MATRIX *C;
416
417 double value;
418 int i;
419
420 int nrowa = NROW(A), ncola = NCOL(A);
421 int nrowb = NROW(B), ncolb = NCOL(B);
422
423 double *a = MATR(A), *b = MATR(B), *c;
424
425 if (nrowa == nrowb && ncola == ncolb)
426 {
427 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
428 nrowa *= ncola;
429 for(i = 0; i < nrowa; i++) *c++ = *a++ / *b++;
430 }
431 else if (nrowa == 1 && ncola == 1)
432 {
433 C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
434 value = *a; nrowb *= ncolb;
435 for(i = 0; i < nrowb; i++) *c++ = value / *b++;
436 }
437 else if (nrowb == 1 && ncolb == 1)
438 {
439 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
440 value = *b; nrowa *= ncola;
441 for(i = 0; i < nrowa; i++) *c++ = *a++ / value;
442 }
443 else
444 error("Div: Incompatible for division.\n");
445
446 return C;
447 }
448
opr_pow(A,B)449 MATRIX *opr_pow(A,B)
450 MATRIX *A, *B;
451 {
452 MATRIX *C;
453
454 int i, j, k, l, power;
455
456 int nrowa = NROW(A), ncola = NCOL(A);
457 int nrowb = NROW(B), ncolb = NCOL(B);
458 int ncolc;
459
460 double *a = MATR(A), *b = MATR(B), *c;
461
462 double *v, value;
463
464 if (nrowb != 1 || ncolb != 1) error("Pow: Matrix ^ Matrix ?.\n");
465
466 if (nrowa == 1 || ncola != nrowa)
467 {
468 C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);
469 value = *b; nrowa *= ncola;
470 for(i = 0; i < nrowa; i++) *c++ = pow(*a++, value);
471
472 return C;
473 }
474
475 power = (int)*b;
476
477 if (power == 0) {
478 C = mat_new(TYPE(A), nrowa, ncola);
479 ncolc = ncola; c = MATR(C);
480 for(i = 0; i < nrowa; i++) MC(i,i) = 1;
481 }
482 else if (abs(power) == 1)
483 {
484 C = mat_copy(A);
485 }
486 else
487 {
488 v = (double *)ALLOCMEM(nrowa * sizeof(double));
489
490 C = mat_new(TYPE(A), nrowa, nrowa);
491 c = MATR(C); b = MATR(A);
492
493 for(l = 1; l < abs(power); l++) {
494 for(i = 0; i < nrowa; i++)
495 {
496 for(j = 0; j < nrowa; j++)
497 {
498 v[j] = 0.0;
499 for(k = 0; k < nrowa; k++) v[j] += b[k] * MA(k,j);
500 }
501 for(j = 0; j < nrowa; j++) *c++ = v[j];
502 b += nrowa;
503 }
504 a = MATR(A); b = c = MATR(C);
505 }
506 FREEMEM((char *)v);
507 }
508
509 if (power < 0)
510 {
511 VARIABLE *inv, *tmp;
512 tmp = (VARIABLE *)ALLOCMEM(VARIABLESIZE);
513 tmp -> this = C;
514 inv = mtr_inv(tmp);
515 mat_free(C);
516 FREEMEM((char *)tmp);
517 C = inv->this;
518 REFCNT(inv)++;
519 var_delete_temp(inv);
520 }
521
522 return C;
523 }
524
opr_trans(A)525 MATRIX *opr_trans(A)
526 MATRIX *A;
527 {
528 MATRIX *C;
529
530 int i, j;
531
532 int ncola = NCOL(A), nrowa = NROW(A);
533 int ncolc;
534
535 double *a = MATR(A), *c;
536
537 C = mat_new(TYPE(A),ncola,nrowa);
538 c = MATR(C); ncolc = nrowa;
539
540 for(i = 0; i < nrowa; i++)
541 for(j = 0; j < ncola; j++) MC(j,i) = *a++;
542
543 return C;
544 }
545
opr_reduction(A,B)546 MATRIX *opr_reduction(A, B)
547 MATRIX *A, *B;
548 {
549 MATRIX *C;
550
551 int i;
552
553 int nrowa = NROW(A), ncola = NCOL(A);
554 int nrowb = NROW(B), ncolb = NCOL(B);
555
556 double *a = MATR(A), *b = MATR(B), *c;
557
558 if ((nrowa == nrowb) && (ncola == ncolb))
559 {
560 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
561 nrowa *= ncola;
562 for(i = 0; i < nrowa; i++, a++) *c++ = (*b++) ? (*a) : (0);
563 }
564 else
565 error("Incompatible for reduction.\n");
566
567 return C;
568 }
569
570
opr_lt(A,B)571 MATRIX *opr_lt(A, B)
572 MATRIX *A, *B;
573 {
574 MATRIX *C;
575
576 int i;
577
578 int nrowa = NROW(A), ncola = NCOL(A);
579 int nrowb = NROW(B), ncolb = NCOL(B);
580
581 double *a = MATR(A), *b = MATR(B), *c;
582
583 if ((nrowa == 1) && (ncola == 1))
584 {
585 C = mat_new(TYPE(B),nrowb, ncolb); c = MATR(C);
586 nrowb *= ncolb;
587 for(i = 0; i < nrowb; i++,c++) if (*a < b[i]) *c = 1;
588 }
589
590 else if ((nrowb == 1) && (ncolb == 1))
591 {
592 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
593 nrowa *= ncola;
594 for(i = 0; i < nrowa; i++,c++) if (a[i] < *b) *c = 1;
595 }
596
597 else if ((nrowa == nrowb) && (ncola == ncolb))
598 {
599 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
600 nrowa *= ncola;
601 for(i = 0; i < nrowa; i++,c++) if (a[i] < b[i]) *c = 1;
602 }
603 else
604 error("lt: Incompatible for comparison.\n");
605
606 return C;
607 }
608
609
opr_le(A,B)610 MATRIX *opr_le(A, B)
611 MATRIX *A, *B;
612 {
613 MATRIX *C;
614
615 int i;
616
617 int nrowa = NROW(A), ncola = NCOL(A);
618 int nrowb = NROW(B), ncolb = NCOL(B);
619
620 double *a = MATR(A), *b = MATR(B), *c;
621
622 if ((nrowa == 1) && (ncola == 1))
623 {
624 C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
625 nrowb *= ncolb;
626 for(i = 0; i < nrowb; i++,c++) if (*a <= b[i]) *c = 1;
627 }
628
629 else if ((nrowb == 1) && (ncolb == 1))
630 {
631 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
632 nrowa *= ncola;
633 for(i = 0; i < nrowa; i++,c++) if (a[i] <= *b) *c = 1;
634 }
635
636 else if ((nrowa == nrowb) && (ncola == ncolb))
637 {
638 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
639 nrowa *= ncola;
640 for(i = 0; i < nrowa; i++,c++) if (a[i] <= b[i]) *c = 1;
641 }
642 else
643 error("le: Incompatible for comparison.\n");
644
645 return C;
646 }
647
648
opr_gt(A,B)649 MATRIX *opr_gt(A, B)
650 MATRIX *A, *B;
651 {
652 MATRIX *C;
653 int i;
654
655 int nrowa = NROW(A), ncola = NCOL(A);
656 int nrowb = NROW(B), ncolb = NCOL(B);
657
658 double *a = MATR(A), *b = MATR(B), *c;
659
660 if ((nrowa == 1) && (ncola == 1))
661 {
662 C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
663 nrowb *= ncolb;
664 for(i = 0; i < nrowb; i++,c++) if (*a > b[i]) *c = 1;
665 }
666
667 else if ((nrowb == 1) && (ncolb == 1))
668 {
669 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
670 nrowa *= ncola;
671 for(i = 0; i < nrowa; i++,c++) if (a[i] > *b) *c = 1;
672 }
673
674 else if ((nrowa == nrowb) && (ncola == ncolb))
675 {
676 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
677 nrowa *= ncola;
678 for(i = 0; i < nrowa; i++,c++) if (a[i] > b[i]) *c = 1;
679 }
680
681 else
682 error("gt: Incompatible for comparison.\n");
683
684 return C;
685 }
686
687
opr_ge(A,B)688 MATRIX *opr_ge(A, B)
689 MATRIX *A, *B;
690 {
691 MATRIX *C;
692 int i;
693
694 int nrowa = NROW(A), ncola = NCOL(A);
695 int nrowb = NROW(B), ncolb = NCOL(B);
696
697 double *a = MATR(A), *b = MATR(B), *c;
698
699 if ((nrowa == 1) && (ncola == 1))
700 {
701 C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
702 nrowb *= ncolb;
703 for(i = 0; i < nrowa; i++,c++) if (*a >= b[i]) *c = 1;
704 }
705
706 else if ((nrowb == 1) && (ncolb == 1))
707 {
708 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
709 nrowa *= ncola;
710 for(i = 0; i < nrowa; i++,c++) if (a[i] >= *b) *c = 1;
711 }
712
713 else if ((nrowa == nrowb) && (ncola == ncolb))
714 {
715 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
716 nrowa *= ncola;
717 for(i = 0; i < nrowa; i++,c++) if (a[i] >= b[i]) *c = 1;
718 }
719
720 else
721 error("ge: Incompatible for comparison.\n");
722
723 return C;
724 }
725
726
opr_eq(A,B)727 MATRIX *opr_eq(A, B)
728 MATRIX *A, *B;
729 {
730 MATRIX *C;
731 int i;
732
733 int nrowa = NROW(A), ncola = NCOL(A);
734 int nrowb = NROW(B), ncolb = NCOL(B);
735
736 double *a = MATR(A), *b = MATR(B), *c;
737
738 if ((nrowa == 1) && (ncola == 1))
739 {
740 C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
741 nrowb *= ncolb;
742 for(i = 0; i < nrowb; i++,c++) if (*a == b[i]) *c = 1;
743 }
744
745 else if ((nrowb == 1) && (ncolb == 1))
746 {
747 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
748 nrowa *= ncola;
749 for(i = 0; i < nrowa; i++,c++) if (a[i] == *b) *c = 1;
750 }
751
752 else if ((nrowa == nrowb) && (ncola == ncolb))
753 {
754 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
755 nrowa *= ncola;
756 for(i = 0; i < nrowa; i++,c++) if (a[i] == b[i]) *c = 1;
757 }
758
759 else
760 error("eq: Incompatible for comparison.\n");
761
762 return C;
763 }
764
765
opr_neq(A,B)766 MATRIX *opr_neq(A, B)
767 MATRIX *A, *B;
768 {
769 MATRIX *C;
770 int i;
771
772 int nrowa = NROW(A), ncola = NCOL(A);
773 int nrowb = NROW(B), ncolb = NCOL(B);
774
775 double *a = MATR(A), *b = MATR(B), *c;
776
777 if ((nrowa == 1) && (ncola == 1))
778 {
779 C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
780 nrowb *= ncolb;
781 for(i = 0; i < nrowb; i++,c++) if (*a != b[i]) *c = 1;
782 }
783
784 else if ((nrowb == 1) && (ncolb == 1))
785 {
786 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
787 nrowa *= ncola;
788 for(i = 0; i < nrowa; i++,c++) if (a[i] != *b) *c = 1;
789 }
790
791 else if ((nrowa == nrowb) && (ncola == ncolb))
792 {
793 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
794 nrowa *= ncola;
795 for(i = 0; i < nrowa; i++,c++) if (a[i] != b[i]) *c = 1;
796 }
797
798 else
799 error("neq: Incompatible for comparison.\n");
800
801 return C;
802 }
803
opr_or(A,B)804 MATRIX *opr_or(A, B)
805 MATRIX *A, *B;
806 {
807 MATRIX *C;
808 int i;
809
810 int nrowa = NROW(A), ncola = NCOL(A);
811 int nrowb = NROW(B), ncolb = NCOL(B);
812
813 double *a = MATR(A), *b = MATR(B), *c;
814
815 if ((nrowa == 1) && (ncola == 1))
816 {
817 C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
818 nrowb *= ncolb;
819 for(i = 0; i < nrowb; i++) *c++ = (*a) || (b[i]);
820 }
821
822 else if ((nrowb == 1) && (ncolb == 1))
823 {
824 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
825 nrowa *= ncola;
826 for(i = 0; i < nrowa; i++) *c++ = (a[i]) || (*b);
827 }
828
829 else if ((nrowa == nrowb) && (ncola == ncolb))
830 {
831 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
832 nrowa *= ncola;
833 for(i = 0; i < nrowa; i++) *c++ = (a[i]) || (b[i]);
834 }
835
836 else
837 error("or: Incompatible for comparison.\n");
838
839 return C;
840 }
841
opr_and(A,B)842 MATRIX *opr_and(A, B)
843 MATRIX *A, *B;
844 {
845 MATRIX *C;
846 int i;
847
848 int nrowa = NROW(A), ncola = NCOL(A);
849 int nrowb = NROW(B), ncolb = NCOL(B);
850
851 double *a = MATR(A), *b = MATR(B), *c;
852
853 if ((nrowa == 1) && (ncola == 1))
854 {
855 C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
856 nrowb *= ncolb;
857 for(i = 0; i < nrowb; i++) *c++ = (*a) && (b[i]);
858 }
859
860 else if ((nrowb == 1) && (ncolb == 1))
861 {
862 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
863 nrowa *= ncola;
864 for(i = 0; i < nrowa; i++) *c++ = (a[i]) && (*b);
865 }
866
867 else if ((nrowa == nrowb) && (ncola == ncolb))
868 {
869 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
870 nrowa *= ncola;
871 for(i = 0; i < nrowa; i++) *c++ = (a[i]) && (b[i]);
872 }
873
874 else
875 error("and: Incompatible for comparison.\n");
876
877 return C;
878 }
879
opr_not(A)880 MATRIX *opr_not(A)
881 MATRIX *A;
882 {
883 MATRIX *C;
884 int i;
885
886 int nrowa = NROW(A), ncola = NCOL(A);
887
888 double *a = MATR(A), *c;
889
890 C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
891 nrowa *= ncola;
892 for(i = 0; i < nrowa; i++,c++) if (*a == 0) *c = 1;
893
894 return C;
895 }
896