![[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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/mat.c,v 1.6 2002/04/10 02:59:10 saito Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
1.4 saito 51: #include "../parse/parse.h"
52:
53: extern int StrassenSize;
1.1 noro 54:
55: void addmat(vl,a,b,c)
56: VL vl;
57: MAT a,b,*c;
58: {
59: int row,col,i,j;
1.4 saito 60: MAT t;
61: pointer *ab,*bb,*tb;
62:
63: if ( !a )
64: *c = b;
65: else if ( !b )
66: *c = a;
67: else if ( (a->row != b->row) || (a->col != b->col) ) {
68: *c = 0; error("addmat : size mismatch add");
69: } else {
70: row = a->row; col = a->col;
71: MKMAT(t,row,col);
72: for ( i = 0; i < row; i++ )
73: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
74: j < col; j++ )
75: addr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
76: *c = t;
77: }
1.1 noro 78: }
79:
80: void submat(vl,a,b,c)
81: VL vl;
82: MAT a,b,*c;
83: {
84: int row,col,i,j;
1.4 saito 85: MAT t;
86: pointer *ab,*bb,*tb;
1.1 noro 87:
1.4 saito 88: if ( !a )
89: chsgnmat(b,c);
90: else if ( !b )
91: *c = a;
92: else if ( (a->row != b->row) || (a->col != b->col) ) {
93: *c = 0; error("submat : size mismatch sub");
94: } else {
95: row = a->row; col = a->col;
96: MKMAT(t,row,col);
97: for ( i = 0; i < row; i++ )
98: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
99: j < col; j++ )
100: subr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
101: *c = t;
102: }
1.1 noro 103: }
104:
105: void mulmat(vl,a,b,c)
106: VL vl;
107: Obj a,b,*c;
108: {
109: if ( !a || !b )
110: *c = 0;
111: else if ( OID(a) <= O_R )
112: mulrmat(vl,(Obj)a,(MAT)b,(MAT *)c);
113: else if ( OID(b) <= O_R )
114: mulrmat(vl,(Obj)b,(MAT)a,(MAT *)c);
115: else
116: switch ( OID(a) ) {
117: case O_VECT:
118: switch ( OID(b) ) {
119: case O_MAT:
120: mulvectmat(vl,(VECT)a,(MAT)b,(VECT *)c); break;
121: case O_VECT: default:
122: notdef(vl,a,b,c); break;
123: }
124: break;
125: case O_MAT:
126: switch ( OID(b) ) {
127: case O_VECT:
128: mulmatvect(vl,(MAT)a,(VECT)b,(VECT *)c); break;
129: case O_MAT:
130: mulmatmat(vl,(MAT)a,(MAT)b,(MAT *)c); break;
131: default:
132: notdef(vl,a,b,c); break;
133: }
134: break;
135: default:
136: notdef(vl,a,b,c); break;
137: }
138: }
139:
140: void divmat(vl,a,b,c)
141: VL vl;
142: Obj a,b,*c;
143: {
144: Obj t;
145:
146: if ( !b )
147: error("divmat : division by 0");
148: else if ( !a )
149: *c = 0;
150: else if ( OID(b) > O_R )
151: notdef(vl,a,b,c);
152: else {
153: divr(vl,(Obj)ONE,b,&t); mulrmat(vl,t,(MAT)a,(MAT *)c);
154: }
155: }
156:
157: void chsgnmat(a,b)
158: MAT a,*b;
159: {
160: MAT t;
161: int row,col,i,j;
162: pointer *ab,*tb;
163:
164: if ( !a )
165: *b = 0;
166: else {
167: row = a->row; col = a->col;
168: MKMAT(t,row,col);
169: for ( i = 0; i < row; i++ )
170: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i];
171: j < col; j++ )
172: chsgnr((Obj)ab[j],(Obj *)&tb[j]);
173: *b = t;
174: }
175: }
176:
177: void pwrmat(vl,a,r,c)
178: VL vl;
179: MAT a;
180: Obj r;
181: MAT *c;
182: {
183: if ( !a )
184: *c = 0;
185: else if ( !r || !NUM(r) || !RATN(r) ||
186: !INT(r) || (SGN((Q)r)<0) || (PL(NM((Q)r))>1) ) {
187: *c = 0; error("pwrmat : invalid exponent");
188: } else if ( a->row != a->col ) {
189: *c = 0; error("pwrmat : non square matrix");
190: } else
191: pwrmatmain(vl,a,QTOS((Q)r),c);
192: }
193:
194: void pwrmatmain(vl,a,e,c)
195: VL vl;
196: MAT a;
197: int e;
198: MAT *c;
199: {
200: MAT t,s;
201:
202: if ( e == 1 ) {
203: *c = a;
204: return;
205: }
206:
207: pwrmatmain(vl,a,e/2,&t);
208: mulmat(vl,(Obj)t,(Obj)t,(Obj *)&s);
209: if ( e % 2 )
210: mulmat(vl,(Obj)s,(Obj)a,(Obj *)c);
211: else
212: *c = s;
213: }
214:
215: void mulrmat(vl,a,b,c)
216: VL vl;
217: Obj a;
218: MAT b,*c;
219: {
220: int row,col,i,j;
221: MAT t;
222: pointer *bb,*tb;
223:
224: if ( !a || !b )
225: *c = 0;
226: else {
227: row = b->row; col = b->col;
228: MKMAT(t,row,col);
229: for ( i = 0; i < row; i++ )
230: for ( j = 0, bb = BDY(b)[i], tb = BDY(t)[i];
231: j < col; j++ )
232: mulr(vl,(Obj)a,(Obj)bb[j],(Obj *)&tb[j]);
233: *c = t;
234: }
235: }
236:
237: void mulmatmat(vl,a,b,c)
238: VL vl;
239: MAT a,b,*c;
240: {
1.4 saito 241: #if 0
1.1 noro 242: int arow,bcol,i,j,k,m;
243: MAT t;
244: pointer s,u,v;
245: pointer *ab,*tb;
246:
1.5 saito 247: /* Mismach col and row */
1.1 noro 248: if ( a->col != b->row ) {
249: *c = 0; error("mulmat : size mismatch");
250: } else {
251: arow = a->row; m = a->col; bcol = b->col;
1.4 saito 252: MKMAt(t,arow,bcol);
1.1 noro 253: for ( i = 0; i < arow; i++ )
254: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i]; j < bcol; j++ ) {
255: for ( k = 0, s = 0; k < m; k++ ) {
1.4 saito 256: mulr(vl,(Obj)ab[k],(Obj)BDY(b)[k][j],(Obj *)&u);
257: addr(vl,(Obj)s,(Obj)u,(Obj *)&v);
258: s = v;
1.1 noro 259: }
260: tb[j] = s;
261: }
262: *c = t;
263: }
264: }
1.4 saito 265:
266: void Strassen(arg, c)
267: NODE arg;
268: Obj *c;
269: {
270: MAT a,b;
271: VL vl;
272:
273: /* tomo */
274: a = (MAT)ARG0(arg);
275: b = (MAT)ARG1(arg);
276: vl = CO;
277: strassen(CO, a, b, c);
278: }
279:
280: void strassen(vl,a,b,c)
281: VL vl;
282: MAT a,b,*c;
283: {
284: #endif
285: int arow,bcol,i,j,k,m, h, arowh, bcolh;
286: MAT t, a11, a12, a21, a22;
287: MAT p, b11, b12, b21, b22;
288: MAT ans1, ans2, ans3, c11, c12, c21, c22;
289: MAT s1, s2, t1, t2, u1, v1, w1, aa, bb;
290: pointer s,u,v;
291: pointer *ab,*tb;
292: int a1row,a2row, a3row,a4row, a1col, a2col, a3col, a4col;
293: int b1row,b2row, b3row,b4row, b1col, b2col, b3col, b4col;
294: int pflag1, pflag2;
1.5 saito 295: /* mismach col and row */
1.4 saito 296: if ( a->col != b->row ) {
297: *c = 0; error("mulmat : size mismatch");
298: }
299: else {
300: pflag1 = 0; pflag2 = 0;
301: arow = a->row; m = a->col; bcol = b->col;
302: arowh = arow/2; bcolh = bcol/2;
303: MKMAT(t,arow,bcol);
304: /* StrassenSize == 0 or matrix size less then StrassenSize,
1.5 saito 305: then calc cannonical algorithm. */
306: if((StrassenSize == 0)||(a->row<=StrassenSize || a->col <= StrassenSize) || (b->row<=StrassenSize || b->col <= StrassenSize)) {
1.4 saito 307: for ( i = 0; i < arow; i++ )
308: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i]; j < bcol; j++ ) {
309: for ( k = 0, s = 0; k < m; k++ ) {
310: mulr(vl,(Obj)ab[k],(Obj)BDY(b)[k][j],(Obj *)&u);
311: addr(vl,(Obj)s,(Obj)u,(Obj *)&v);
312: s = v;
313: }
314: tb[j] = s;
315: }
316: *c = t;
317: return;
318: }
1.5 saito 319: /* padding odd col and row to even number for zero */
1.4 saito 320: i = arow/2;
321: j = arow - i;
322: if (i != j) {
323: arow++;
324: pflag1 = 1;
325: }
326: i = m/2;
327: j = m - i;
328: if (i != j) {
329: m++;
330: pflag2 = 1;
331: }
332: MKMAT(aa, arow, m);
333: for (i = 0; i < a->row; i++) {
334: for (j = 0; j < a->col; j++) {
335: aa->body[i][j] = a->body[i][j];
336: }
337: }
338: i = bcol/2;
339: j = bcol - i;
340: if (i != j) {
341: bcol++;
342: }
343: MKMAT(bb, m, bcol);
344: for (i = 0; i < b->row; i++) {
345: for ( j = 0; j < b->col; j++) {
346: bb->body[i][j] = b->body[i][j];
347: }
348: }
349:
1.5 saito 350: /* split matrix A and B */
1.4 saito 351: a1row = aa->row/2; a1col = aa->col/2;
352: MKMAT(a11,a1row,a1col);
353: MKMAT(a21,a1row,a1col);
354: MKMAT(a12,a1row,a1col);
355: MKMAT(a22,a1row,a1col);
356:
357: b1row = bb->row/2; b1col = bb->col/2;
358: MKMAT(b11,b1row,b1col);
359: MKMAT(b21,b1row,b1col);
360: MKMAT(b12,b1row,b1col);
361: MKMAT(b22,b1row,b1col);
362:
1.5 saito 363: /* make a11 matrix */
1.4 saito 364: for (i = 0; i < a1row; i++) {
365: for (j = 0; j < a1col; j++) {
366: a11->body[i][j] = aa->body[i][j];
367: }
368: }
369:
1.5 saito 370: /* make a21 matrix */
1.4 saito 371: for (i = a1row; i < aa->row; i++) {
372: for (j = 0; j < a1col; j++) {
373: a21->body[i-a1row][j] = aa->body[i][j];
374: }
375: }
376:
1.5 saito 377: /* create a12 matrix */
1.4 saito 378: for (i = 0; i < a1row; i++) {
379: for (j = a1col; j < aa->col; j++) {
380: a12->body[i][j-a1col] = aa->body[i][j];
381: }
382: }
383:
1.5 saito 384: /* create a22 matrix */
1.4 saito 385: for (i = a1row; i < aa->row; i++) {
386: for (j = a1col; j < aa->col; j++) {
387: a22->body[i-a1row][j-a1col] = aa->body[i][j];
388: }
389: }
390:
391:
1.5 saito 392: /* create b11 submatrix */
1.4 saito 393: for (i = 0; i < b1row; i++) {
394: for (j = 0; j < b1col; j++) {
395: b11->body[i][j] = bb->body[i][j];
396: }
397: }
398:
1.5 saito 399: /* create b21 submatrix */
1.4 saito 400: for (i = b1row; i < bb->row; i++) {
401: for (j = 0; j < b1col; j++) {
402: b21->body[i-b1row][j] = bb->body[i][j];
403: }
404: }
405:
1.5 saito 406: /* create b12 submatrix */
1.4 saito 407: for (i = 0; i < b1row; i++) {
408: for (j = b1col; j < bb->col; j++) {
409: b12->body[i][j-b1col] = bb->body[i][j];
410: }
411: }
412:
1.5 saito 413: /* create b22 submatrix */
1.4 saito 414: for (i = b1row; i < bb->row; i++) {
415: for (j = b1col; j < bb->col; j++) {
416: b22->body[i-b1row][j-b1col] = bb->body[i][j];
417: }
418: }
1.5 saito 419: /* expand matrix by Strassen-Winograd algorithm */
1.4 saito 420: /* s1=A21+A22 */
421: addmat(vl,a21,a22,&s1);
422:
423: /* s2=s1-A11 */
424: submat(vl,s1,a11,&s2);
425:
426: /* t1=B12-B11 */
427: submat(vl, b12, b11, &t1);
428:
429: /* t2=B22-t1 */
430: submat(vl, b22, t1, &t2);
431:
432: /* u=(A11-A21)*(B22-B12) */
433: submat(vl, a11, a21, &ans1);
434: submat(vl, b22, b12, &ans2);
435: mulmatmat(vl, ans1, ans2, &u1);
436:
437: /* v=s1*t1 */
438: mulmatmat(vl, s1, t1, &v1);
439:
440: /* w=A11*B11+s2*t2 */
441: mulmatmat(vl, a11, b11, &ans1);
442: mulmatmat(vl, s2, t2, &ans2);
443: addmat(vl, ans1, ans2, &w1);
444:
445: /* C11 = A11*B11+A12*B21 */
446: mulmatmat(vl, a12, b21, &ans2);
447: addmat(vl, ans1, ans2, &c11);
448:
449: /* C12 = w1+v1+(A12-s2)*B22 */
450: submat(vl, a12, s2, &ans1);
451: mulmatmat(vl, ans1, b22, &ans2);
452: addmat(vl, w1, v1, &ans1);
453: addmat(vl, ans1, ans2, &c12);
454:
455: /* C21 = w1+u1+A22*(B21-t2) */
456: submat(vl, b21, t2, &ans1);
457: mulmatmat(vl, a22, ans1, &ans2);
1.6 saito 458: addmat(vl, w1, u1, &ans1);
459: addmat(vl, ans1, ans2, &c21);
460:
461: /* C22 = w1 + u1 + v1 */
462: addmat(vl, ans1, v1, &c22);
463: }
464:
465: for(i =0; i<c11->row; i++) {
466: for ( j=0; j < c11->col; j++) {
467: t->body[i][j] = c11->body[i][j];
468: }
469: }
470: if (pflag1 == 0) {
471: k = c21->row;
472: } else {
473: k = c21->row - 1;
474: }
475: for(i =0; i<k; i++) {
476: for ( j=0; j < c21->col; j++) {
477: t->body[i+c11->row][j] = c21->body[i][j];
478: }
479: }
480: if (pflag2 == 0) {
481: h = c12->col;
482: } else {
483: h = c12->col -1;
484: }
485: for(i =0; i<c12->row; i++) {
1.4 saito 486: for ( j=0; j < k; j++) {
487: t->body[i][j+c11->col] = c12->body[i][j];
488: }
489: }
490: if (pflag1 == 0) {
491: k = c22->row;
492: } else {
493: k = c22->row -1;
494: }
495: if (pflag2 == 0) {
496: h = c22->col;
497: } else {
498: h = c22->col - 1;
499: }
500: for(i =0; i<k; i++) {
501: for ( j=0; j < h; j++) {
502: t->body[i+c11->row][j+c11->col] = c22->body[i][j];
503: }
504: }
505: *c = t;
506: }
507:
508:
1.1 noro 509:
510: void mulmatvect(vl,a,b,c)
511: VL vl;
512: MAT a;
513: VECT b;
514: VECT *c;
515: {
516: int arow,i,j,m;
517: VECT t;
518: pointer s,u,v;
519: pointer *ab;
520:
521: if ( !a || !b )
522: *c = 0;
523: else if ( a->col != b->len ) {
524: *c = 0; error("mulmatvect : size mismatch");
525: } else {
1.7 ! noro 526: for ( i = 0; i < b->len; i++ )
! 527: if ( BDY(b)[i] && OID((Obj)BDY(b)[i]) > O_R )
! 528: error("mulmatvect : invalid argument");
1.1 noro 529: arow = a->row; m = a->col;
530: MKVECT(t,arow);
531: for ( i = 0; i < arow; i++ ) {
532: for ( j = 0, s = 0, ab = BDY(a)[i]; j < m; j++ ) {
533: mulr(vl,(Obj)ab[j],(Obj)BDY(b)[j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
534: }
535: BDY(t)[i] = s;
536: }
537: *c = t;
538: }
539: }
540:
541: void mulvectmat(vl,a,b,c)
542: VL vl;
543: VECT a;
544: MAT b;
545: VECT *c;
546: {
547: int bcol,i,j,m;
548: VECT t;
549: pointer s,u,v;
550:
551: if ( !a || !b )
552: *c = 0;
553: else if ( a->len != b->row ) {
554: *c = 0; error("mulvectmat : size mismatch");
555: } else {
1.7 ! noro 556: for ( i = 0; i < a->len; i++ )
! 557: if ( BDY(a)[i] && OID((Obj)BDY(a)[i]) > O_R )
! 558: error("mulvectmat : invalid argument");
1.1 noro 559: bcol = b->col; m = a->len;
560: MKVECT(t,bcol);
561: for ( j = 0; j < bcol; j++ ) {
562: for ( i = 0, s = 0; i < m; i++ ) {
563: mulr(vl,(Obj)BDY(a)[i],(Obj)BDY(b)[i][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
564: }
565: BDY(t)[j] = s;
566: }
567: *c = t;
568: }
569: }
570:
571: int compmat(vl,a,b)
572: VL vl;
573: MAT a,b;
574: {
575: int i,j,t,row,col;
576:
577: if ( !a )
578: return b?-1:0;
579: else if ( !b )
580: return 1;
581: else if ( a->row != b->row )
582: return a->row>b->row ? 1 : -1;
583: else if (a->col != b->col )
584: return a->col > b->col ? 1 : -1;
585: else {
586: row = a->row; col = a->col;
587: for ( i = 0; i < row; i++ )
588: for ( j = 0; j < col; j++ )
589: if ( t = compr(vl,(Obj)BDY(a)[i][j],(Obj)BDY(b)[i][j]) )
590: return t;
591: return 0;
592: }
593: }
594:
595: pointer **almat_pointer(n,m)
596: int n,m;
597: {
598: pointer **mat;
599: int i;
600:
601: mat = (pointer **)MALLOC(n*sizeof(pointer *));
602: for ( i = 0; i < n; i++ )
603: mat[i] = (pointer *)CALLOC(m,sizeof(pointer));
604: return mat;
605: }