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