![[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.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 ! noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.3 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/mat.c,v 1.2 2000/08/21 08:31:28 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51:
52: void addmat(vl,a,b,c)
53: VL vl;
54: MAT a,b,*c;
55: {
56: int row,col,i,j;
57: MAT t;
58: pointer *ab,*bb,*tb;
59:
60: if ( !a )
61: *c = b;
62: else if ( !b )
63: *c = a;
64: else if ( (a->row != b->row) || (a->col != b->col) ) {
65: *c = 0; error("addmat : size mismatch");
66: } else {
67: row = a->row; col = a->col;
68: MKMAT(t,row,col);
69: for ( i = 0; i < row; i++ )
70: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
71: j < col; j++ )
72: addr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
73: *c = t;
74: }
75: }
76:
77: void submat(vl,a,b,c)
78: VL vl;
79: MAT a,b,*c;
80: {
81: int row,col,i,j;
82: MAT t;
83: pointer *ab,*bb,*tb;
84:
85: if ( !a )
86: chsgnmat(b,c);
87: else if ( !b )
88: *c = a;
89: else if ( (a->row != b->row) || (a->col != b->col) ) {
90: *c = 0; error("submat : size mismatch");
91: } else {
92: row = a->row; col = a->col;
93: MKMAT(t,row,col);
94: for ( i = 0; i < row; i++ )
95: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
96: j < col; j++ )
97: subr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
98: *c = t;
99: }
100: }
101:
102: void mulmat(vl,a,b,c)
103: VL vl;
104: Obj a,b,*c;
105: {
106: if ( !a || !b )
107: *c = 0;
108: else if ( OID(a) <= O_R )
109: mulrmat(vl,(Obj)a,(MAT)b,(MAT *)c);
110: else if ( OID(b) <= O_R )
111: mulrmat(vl,(Obj)b,(MAT)a,(MAT *)c);
112: else
113: switch ( OID(a) ) {
114: case O_VECT:
115: switch ( OID(b) ) {
116: case O_MAT:
117: mulvectmat(vl,(VECT)a,(MAT)b,(VECT *)c); break;
118: case O_VECT: default:
119: notdef(vl,a,b,c); break;
120: }
121: break;
122: case O_MAT:
123: switch ( OID(b) ) {
124: case O_VECT:
125: mulmatvect(vl,(MAT)a,(VECT)b,(VECT *)c); break;
126: case O_MAT:
127: mulmatmat(vl,(MAT)a,(MAT)b,(MAT *)c); break;
128: default:
129: notdef(vl,a,b,c); break;
130: }
131: break;
132: default:
133: notdef(vl,a,b,c); break;
134: }
135: }
136:
137: void divmat(vl,a,b,c)
138: VL vl;
139: Obj a,b,*c;
140: {
141: Obj t;
142:
143: if ( !b )
144: error("divmat : division by 0");
145: else if ( !a )
146: *c = 0;
147: else if ( OID(b) > O_R )
148: notdef(vl,a,b,c);
149: else {
150: divr(vl,(Obj)ONE,b,&t); mulrmat(vl,t,(MAT)a,(MAT *)c);
151: }
152: }
153:
154: void chsgnmat(a,b)
155: MAT a,*b;
156: {
157: MAT t;
158: int row,col,i,j;
159: pointer *ab,*tb;
160:
161: if ( !a )
162: *b = 0;
163: else {
164: row = a->row; col = a->col;
165: MKMAT(t,row,col);
166: for ( i = 0; i < row; i++ )
167: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i];
168: j < col; j++ )
169: chsgnr((Obj)ab[j],(Obj *)&tb[j]);
170: *b = t;
171: }
172: }
173:
174: void pwrmat(vl,a,r,c)
175: VL vl;
176: MAT a;
177: Obj r;
178: MAT *c;
179: {
180: if ( !a )
181: *c = 0;
182: else if ( !r || !NUM(r) || !RATN(r) ||
183: !INT(r) || (SGN((Q)r)<0) || (PL(NM((Q)r))>1) ) {
184: *c = 0; error("pwrmat : invalid exponent");
185: } else if ( a->row != a->col ) {
186: *c = 0; error("pwrmat : non square matrix");
187: } else
188: pwrmatmain(vl,a,QTOS((Q)r),c);
189: }
190:
191: void pwrmatmain(vl,a,e,c)
192: VL vl;
193: MAT a;
194: int e;
195: MAT *c;
196: {
197: MAT t,s;
198:
199: if ( e == 1 ) {
200: *c = a;
201: return;
202: }
203:
204: pwrmatmain(vl,a,e/2,&t);
205: mulmat(vl,(Obj)t,(Obj)t,(Obj *)&s);
206: if ( e % 2 )
207: mulmat(vl,(Obj)s,(Obj)a,(Obj *)c);
208: else
209: *c = s;
210: }
211:
212: void mulrmat(vl,a,b,c)
213: VL vl;
214: Obj a;
215: MAT b,*c;
216: {
217: int row,col,i,j;
218: MAT t;
219: pointer *bb,*tb;
220:
221: if ( !a || !b )
222: *c = 0;
223: else {
224: row = b->row; col = b->col;
225: MKMAT(t,row,col);
226: for ( i = 0; i < row; i++ )
227: for ( j = 0, bb = BDY(b)[i], tb = BDY(t)[i];
228: j < col; j++ )
229: mulr(vl,(Obj)a,(Obj)bb[j],(Obj *)&tb[j]);
230: *c = t;
231: }
232: }
233:
234: void mulmatmat(vl,a,b,c)
235: VL vl;
236: MAT a,b,*c;
237: {
238: int arow,bcol,i,j,k,m;
239: MAT t;
240: pointer s,u,v;
241: pointer *ab,*tb;
242:
243: if ( a->col != b->row ) {
244: *c = 0; error("mulmat : size mismatch");
245: } else {
246: arow = a->row; m = a->col; bcol = b->col;
247: MKMAT(t,arow,bcol);
248: for ( i = 0; i < arow; i++ )
249: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i]; j < bcol; j++ ) {
250: for ( k = 0, s = 0; k < m; k++ ) {
251: mulr(vl,(Obj)ab[k],(Obj)BDY(b)[k][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
252: }
253: tb[j] = s;
254: }
255: *c = t;
256: }
257: }
258:
259: void mulmatvect(vl,a,b,c)
260: VL vl;
261: MAT a;
262: VECT b;
263: VECT *c;
264: {
265: int arow,i,j,m;
266: VECT t;
267: pointer s,u,v;
268: pointer *ab;
269:
270: if ( !a || !b )
271: *c = 0;
272: else if ( a->col != b->len ) {
273: *c = 0; error("mulmatvect : size mismatch");
274: } else {
275: arow = a->row; m = a->col;
276: MKVECT(t,arow);
277: for ( i = 0; i < arow; i++ ) {
278: for ( j = 0, s = 0, ab = BDY(a)[i]; j < m; j++ ) {
279: mulr(vl,(Obj)ab[j],(Obj)BDY(b)[j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
280: }
281: BDY(t)[i] = s;
282: }
283: *c = t;
284: }
285: }
286:
287: void mulvectmat(vl,a,b,c)
288: VL vl;
289: VECT a;
290: MAT b;
291: VECT *c;
292: {
293: int bcol,i,j,m;
294: VECT t;
295: pointer s,u,v;
296:
297: if ( !a || !b )
298: *c = 0;
299: else if ( a->len != b->row ) {
300: *c = 0; error("mulvectmat : size mismatch");
301: } else {
302: bcol = b->col; m = a->len;
303: MKVECT(t,bcol);
304: for ( j = 0; j < bcol; j++ ) {
305: for ( i = 0, s = 0; i < m; i++ ) {
306: mulr(vl,(Obj)BDY(a)[i],(Obj)BDY(b)[i][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
307: }
308: BDY(t)[j] = s;
309: }
310: *c = t;
311: }
312: }
313:
314: int compmat(vl,a,b)
315: VL vl;
316: MAT a,b;
317: {
318: int i,j,t,row,col;
319:
320: if ( !a )
321: return b?-1:0;
322: else if ( !b )
323: return 1;
324: else if ( a->row != b->row )
325: return a->row>b->row ? 1 : -1;
326: else if (a->col != b->col )
327: return a->col > b->col ? 1 : -1;
328: else {
329: row = a->row; col = a->col;
330: for ( i = 0; i < row; i++ )
331: for ( j = 0; j < col; j++ )
332: if ( t = compr(vl,(Obj)BDY(a)[i][j],(Obj)BDY(b)[i][j]) )
333: return t;
334: return 0;
335: }
336: }
337:
338: pointer **almat_pointer(n,m)
339: int n,m;
340: {
341: pointer **mat;
342: int i;
343:
344: mat = (pointer **)MALLOC(n*sizeof(pointer *));
345: for ( i = 0; i < n; i++ )
346: mat[i] = (pointer *)CALLOC(m,sizeof(pointer));
347: return mat;
348: }