![[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.4 ! saito 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/mat.c,v 1.3 2000/08/22 05:04:05 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: {
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.4 ! saito 247: /* 行列のcol,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;
! 295: /* 行列のcol,rowの数があわない場合 */
! 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,
! 305: then calc cannonical algorizm. */
! 306: if((StrassenSize == 0)||(a->row<=StrassenSize || a->col <= StrassenSize)) {
! 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: }
! 319: /* 行列が奇数次の場合は偶数次になるように0でpadding */
! 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:
! 350: /* 行列A,Bを分割 */
! 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:
! 363: /* a11の行列を作る */
! 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:
! 370: /* a21の行列を作る */
! 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:
! 377: /* a12の行列を作る */
! 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:
! 384: /* a22の行列を作る */
! 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:
! 392: /* b11の行列を作る */
! 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:
! 399: /* b21の行列を作る */
! 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:
! 406: /* b12の行列を作る */
! 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:
! 413: /* b22の行列を作る */
! 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: }
! 419: /* Strassen-Winogradの方法で展開 */
! 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: /* tomo
! 437: strassen(vl, ans1, ans2, &u1);
! 438: */
! 439:
! 440: /* v=s1*t1 */
! 441: mulmatmat(vl, s1, t1, &v1);
! 442: /* tomo
! 443: strassen(vl, s1, t1, &v1);
! 444: */
! 445:
! 446: /* w=A11*B11+s2*t2 */
! 447: mulmatmat(vl, a11, b11, &ans1);
! 448: mulmatmat(vl, s2, t2, &ans2);
! 449: /* tomo
! 450: strassen(vl, a11, b11, &ans1);
! 451: strassen(vl, s2, t2, &ans2);
! 452: */
! 453: addmat(vl, ans1, ans2, &w1);
! 454:
! 455: /* C11 = A11*B11+A12*B21 */
! 456: mulmatmat(vl, a12, b21, &ans2);
! 457: /* tomo
! 458: strassen(vl, a12, b21, &ans2);
! 459: */
! 460: addmat(vl, ans1, ans2, &c11);
! 461:
! 462: /* C12 = w1+v1+(A12-s2)*B22 */
! 463: submat(vl, a12, s2, &ans1);
! 464: mulmatmat(vl, ans1, b22, &ans2);
! 465: /* tomo
! 466: strassen(vl, ans1, b22, &ans2);
! 467: */
! 468: addmat(vl, w1, v1, &ans1);
! 469: addmat(vl, ans1, ans2, &c12);
! 470:
! 471: /* C21 = w1+u1+A22*(B21-t2) */
! 472: submat(vl, b21, t2, &ans1);
! 473: mulmatmat(vl, a22, ans1, &ans2);
! 474: /* tomo
! 475: strassen(vl, a22, ans1, &ans2);
! 476: */
! 477: addmat(vl, w1, u1, &ans1);
! 478: addmat(vl, ans1, ans2, &c21);
! 479:
! 480: /* C22 = w1 + u1 + v1 */
! 481: addmat(vl, ans1, v1, &c22);
! 482:
! 483: }
! 484:
! 485: /* 解の領域tに計算結果を戻す */
! 486: for(i =0; i<c11->row; i++) {
! 487: for ( j=0; j < c11->col; j++) {
! 488: t->body[i][j] = c11->body[i][j];
! 489: }
! 490: }
! 491: if (pflag1 == 0) {
! 492: k = c21->row;
! 493: } else {
! 494: k = c21->row - 1;
! 495: }
! 496: for(i =0; i<k; i++) {
! 497: for ( j=0; j < c21->col; j++) {
! 498: t->body[i+c11->row][j] = c21->body[i][j];
! 499: }
! 500: }
! 501: if (pflag2 == 0) {
! 502: h = c12->col;
! 503: } else {
! 504: h = c12->col -1;
! 505: }
! 506: for(i =0; i<c12->row; i++) {
! 507: for ( j=0; j < k; j++) {
! 508: t->body[i][j+c11->col] = c12->body[i][j];
! 509: }
! 510: }
! 511: if (pflag1 == 0) {
! 512: k = c22->row;
! 513: } else {
! 514: k = c22->row -1;
! 515: }
! 516: if (pflag2 == 0) {
! 517: h = c22->col;
! 518: } else {
! 519: h = c22->col - 1;
! 520: }
! 521: for(i =0; i<k; i++) {
! 522: for ( j=0; j < h; j++) {
! 523: t->body[i+c11->row][j+c11->col] = c22->body[i][j];
! 524: }
! 525: }
! 526: *c = t;
! 527: }
! 528:
! 529:
1.1 noro 530:
531: void mulmatvect(vl,a,b,c)
532: VL vl;
533: MAT a;
534: VECT b;
535: VECT *c;
536: {
537: int arow,i,j,m;
538: VECT t;
539: pointer s,u,v;
540: pointer *ab;
541:
542: if ( !a || !b )
543: *c = 0;
544: else if ( a->col != b->len ) {
545: *c = 0; error("mulmatvect : size mismatch");
546: } else {
547: arow = a->row; m = a->col;
548: MKVECT(t,arow);
549: for ( i = 0; i < arow; i++ ) {
550: for ( j = 0, s = 0, ab = BDY(a)[i]; j < m; j++ ) {
551: mulr(vl,(Obj)ab[j],(Obj)BDY(b)[j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
552: }
553: BDY(t)[i] = s;
554: }
555: *c = t;
556: }
557: }
558:
559: void mulvectmat(vl,a,b,c)
560: VL vl;
561: VECT a;
562: MAT b;
563: VECT *c;
564: {
565: int bcol,i,j,m;
566: VECT t;
567: pointer s,u,v;
568:
569: if ( !a || !b )
570: *c = 0;
571: else if ( a->len != b->row ) {
572: *c = 0; error("mulvectmat : size mismatch");
573: } else {
574: bcol = b->col; m = a->len;
575: MKVECT(t,bcol);
576: for ( j = 0; j < bcol; j++ ) {
577: for ( i = 0, s = 0; i < m; i++ ) {
578: mulr(vl,(Obj)BDY(a)[i],(Obj)BDY(b)[i][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
579: }
580: BDY(t)[j] = s;
581: }
582: *c = t;
583: }
584: }
585:
586: int compmat(vl,a,b)
587: VL vl;
588: MAT a,b;
589: {
590: int i,j,t,row,col;
591:
592: if ( !a )
593: return b?-1:0;
594: else if ( !b )
595: return 1;
596: else if ( a->row != b->row )
597: return a->row>b->row ? 1 : -1;
598: else if (a->col != b->col )
599: return a->col > b->col ? 1 : -1;
600: else {
601: row = a->row; col = a->col;
602: for ( i = 0; i < row; i++ )
603: for ( j = 0; j < col; j++ )
604: if ( t = compr(vl,(Obj)BDY(a)[i][j],(Obj)BDY(b)[i][j]) )
605: return t;
606: return 0;
607: }
608: }
609:
610: pointer **almat_pointer(n,m)
611: int n,m;
612: {
613: pointer **mat;
614: int i;
615:
616: mat = (pointer **)MALLOC(n*sizeof(pointer *));
617: for ( i = 0; i < n; i++ )
618: mat[i] = (pointer *)CALLOC(m,sizeof(pointer));
619: return mat;
620: }