![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mat.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to mat.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/mat.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3:
4: void addmat(vl,a,b,c)
5: VL vl;
6: MAT a,b,*c;
7: {
8: int row,col,i,j;
9: MAT t;
10: pointer *ab,*bb,*tb;
11:
12: if ( !a )
13: *c = b;
14: else if ( !b )
15: *c = a;
16: else if ( (a->row != b->row) || (a->col != b->col) ) {
17: *c = 0; error("addmat : size mismatch");
18: } else {
19: row = a->row; col = a->col;
20: MKMAT(t,row,col);
21: for ( i = 0; i < row; i++ )
22: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
23: j < col; j++ )
24: addr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
25: *c = t;
26: }
27: }
28:
29: void submat(vl,a,b,c)
30: VL vl;
31: MAT a,b,*c;
32: {
33: int row,col,i,j;
34: MAT t;
35: pointer *ab,*bb,*tb;
36:
37: if ( !a )
38: chsgnmat(b,c);
39: else if ( !b )
40: *c = a;
41: else if ( (a->row != b->row) || (a->col != b->col) ) {
42: *c = 0; error("submat : size mismatch");
43: } else {
44: row = a->row; col = a->col;
45: MKMAT(t,row,col);
46: for ( i = 0; i < row; i++ )
47: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
48: j < col; j++ )
49: subr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
50: *c = t;
51: }
52: }
53:
54: void mulmat(vl,a,b,c)
55: VL vl;
56: Obj a,b,*c;
57: {
58: if ( !a || !b )
59: *c = 0;
60: else if ( OID(a) <= O_R )
61: mulrmat(vl,(Obj)a,(MAT)b,(MAT *)c);
62: else if ( OID(b) <= O_R )
63: mulrmat(vl,(Obj)b,(MAT)a,(MAT *)c);
64: else
65: switch ( OID(a) ) {
66: case O_VECT:
67: switch ( OID(b) ) {
68: case O_MAT:
69: mulvectmat(vl,(VECT)a,(MAT)b,(VECT *)c); break;
70: case O_VECT: default:
71: notdef(vl,a,b,c); break;
72: }
73: break;
74: case O_MAT:
75: switch ( OID(b) ) {
76: case O_VECT:
77: mulmatvect(vl,(MAT)a,(VECT)b,(VECT *)c); break;
78: case O_MAT:
79: mulmatmat(vl,(MAT)a,(MAT)b,(MAT *)c); break;
80: default:
81: notdef(vl,a,b,c); break;
82: }
83: break;
84: default:
85: notdef(vl,a,b,c); break;
86: }
87: }
88:
89: void divmat(vl,a,b,c)
90: VL vl;
91: Obj a,b,*c;
92: {
93: Obj t;
94:
95: if ( !b )
96: error("divmat : division by 0");
97: else if ( !a )
98: *c = 0;
99: else if ( OID(b) > O_R )
100: notdef(vl,a,b,c);
101: else {
102: divr(vl,(Obj)ONE,b,&t); mulrmat(vl,t,(MAT)a,(MAT *)c);
103: }
104: }
105:
106: void chsgnmat(a,b)
107: MAT a,*b;
108: {
109: MAT t;
110: int row,col,i,j;
111: pointer *ab,*tb;
112:
113: if ( !a )
114: *b = 0;
115: else {
116: row = a->row; col = a->col;
117: MKMAT(t,row,col);
118: for ( i = 0; i < row; i++ )
119: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i];
120: j < col; j++ )
121: chsgnr((Obj)ab[j],(Obj *)&tb[j]);
122: *b = t;
123: }
124: }
125:
126: void pwrmat(vl,a,r,c)
127: VL vl;
128: MAT a;
129: Obj r;
130: MAT *c;
131: {
132: if ( !a )
133: *c = 0;
134: else if ( !r || !NUM(r) || !RATN(r) ||
135: !INT(r) || (SGN((Q)r)<0) || (PL(NM((Q)r))>1) ) {
136: *c = 0; error("pwrmat : invalid exponent");
137: } else if ( a->row != a->col ) {
138: *c = 0; error("pwrmat : non square matrix");
139: } else
140: pwrmatmain(vl,a,QTOS((Q)r),c);
141: }
142:
143: void pwrmatmain(vl,a,e,c)
144: VL vl;
145: MAT a;
146: int e;
147: MAT *c;
148: {
149: MAT t,s;
150:
151: if ( e == 1 ) {
152: *c = a;
153: return;
154: }
155:
156: pwrmatmain(vl,a,e/2,&t);
157: mulmat(vl,(Obj)t,(Obj)t,(Obj *)&s);
158: if ( e % 2 )
159: mulmat(vl,(Obj)s,(Obj)a,(Obj *)c);
160: else
161: *c = s;
162: }
163:
164: void mulrmat(vl,a,b,c)
165: VL vl;
166: Obj a;
167: MAT b,*c;
168: {
169: int row,col,i,j;
170: MAT t;
171: pointer *bb,*tb;
172:
173: if ( !a || !b )
174: *c = 0;
175: else {
176: row = b->row; col = b->col;
177: MKMAT(t,row,col);
178: for ( i = 0; i < row; i++ )
179: for ( j = 0, bb = BDY(b)[i], tb = BDY(t)[i];
180: j < col; j++ )
181: mulr(vl,(Obj)a,(Obj)bb[j],(Obj *)&tb[j]);
182: *c = t;
183: }
184: }
185:
186: void mulmatmat(vl,a,b,c)
187: VL vl;
188: MAT a,b,*c;
189: {
190: int arow,bcol,i,j,k,m;
191: MAT t;
192: pointer s,u,v;
193: pointer *ab,*tb;
194:
195: if ( a->col != b->row ) {
196: *c = 0; error("mulmat : size mismatch");
197: } else {
198: arow = a->row; m = a->col; bcol = b->col;
199: MKMAT(t,arow,bcol);
200: for ( i = 0; i < arow; i++ )
201: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i]; j < bcol; j++ ) {
202: for ( k = 0, s = 0; k < m; k++ ) {
203: mulr(vl,(Obj)ab[k],(Obj)BDY(b)[k][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
204: }
205: tb[j] = s;
206: }
207: *c = t;
208: }
209: }
210:
211: void mulmatvect(vl,a,b,c)
212: VL vl;
213: MAT a;
214: VECT b;
215: VECT *c;
216: {
217: int arow,i,j,m;
218: VECT t;
219: pointer s,u,v;
220: pointer *ab;
221:
222: if ( !a || !b )
223: *c = 0;
224: else if ( a->col != b->len ) {
225: *c = 0; error("mulmatvect : size mismatch");
226: } else {
227: arow = a->row; m = a->col;
228: MKVECT(t,arow);
229: for ( i = 0; i < arow; i++ ) {
230: for ( j = 0, s = 0, ab = BDY(a)[i]; j < m; j++ ) {
231: mulr(vl,(Obj)ab[j],(Obj)BDY(b)[j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
232: }
233: BDY(t)[i] = s;
234: }
235: *c = t;
236: }
237: }
238:
239: void mulvectmat(vl,a,b,c)
240: VL vl;
241: VECT a;
242: MAT b;
243: VECT *c;
244: {
245: int bcol,i,j,m;
246: VECT t;
247: pointer s,u,v;
248:
249: if ( !a || !b )
250: *c = 0;
251: else if ( a->len != b->row ) {
252: *c = 0; error("mulvectmat : size mismatch");
253: } else {
254: bcol = b->col; m = a->len;
255: MKVECT(t,bcol);
256: for ( j = 0; j < bcol; j++ ) {
257: for ( i = 0, s = 0; i < m; i++ ) {
258: mulr(vl,(Obj)BDY(a)[i],(Obj)BDY(b)[i][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
259: }
260: BDY(t)[j] = s;
261: }
262: *c = t;
263: }
264: }
265:
266: int compmat(vl,a,b)
267: VL vl;
268: MAT a,b;
269: {
270: int i,j,t,row,col;
271:
272: if ( !a )
273: return b?-1:0;
274: else if ( !b )
275: return 1;
276: else if ( a->row != b->row )
277: return a->row>b->row ? 1 : -1;
278: else if (a->col != b->col )
279: return a->col > b->col ? 1 : -1;
280: else {
281: row = a->row; col = a->col;
282: for ( i = 0; i < row; i++ )
283: for ( j = 0; j < col; j++ )
284: if ( t = compr(vl,(Obj)BDY(a)[i][j],(Obj)BDY(b)[i][j]) )
285: return t;
286: return 0;
287: }
288: }
289:
290: pointer **almat_pointer(n,m)
291: int n,m;
292: {
293: pointer **mat;
294: int i;
295:
296: mat = (pointer **)MALLOC(n*sizeof(pointer *));
297: for ( i = 0; i < n; i++ )
298: mat[i] = (pointer *)CALLOC(m,sizeof(pointer));
299: return mat;
300: }