![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to iarray.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin|
Return to iarray.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.1 ! saito 1: /*
! 2: * $OpenXM$
! 3: */
! 4: #include "ca.h"
! 5: #include "base.h"
! 6: #include "parse.h"
! 7: #include "inline.h"
! 8:
! 9: extern int StrassenSize;
! 10:
! 11: struct ftab imat_tab[] = {
! 12: {"newimat", Pnewimat,2},
! 13: {"m2im", Pm2Im,1},
! 14: {"im2m", PIm2m,1},
! 15: {0,0,0},
! 16: };
! 17:
! 18: const IENT zent = { -1, -1, -1, 0}; /* last ient const */
! 19:
! 20: void MEnt(cr, row, col, trg, ent)
! 21: int cr, row, col;
! 22: Obj trg;
! 23: IENT *ent;
! 24: {
! 25: /* make Index Entry IENT */
! 26: /* set:cr, row, trag */
! 27: /* return: value ent */
! 28: /* MEnt(cr, row, col, trg, &ent); */
! 29:
! 30: ent->cr = cr;
! 31: ent->row = row;
! 32: ent->col = col;
! 33: ent->body = (pointer)trg;
! 34: }
! 35:
! 36: void GetNextIent(Im, ent, c)
! 37: IMATC *Im;
! 38: IENT *ent;
! 39: int *c;
! 40: {
! 41: /* get next IENT */
! 42: /* set: Im, c */
! 43: /* return: ent, c */
! 44: /* GetNextIent(&Im, &ent, &c); */
! 45:
! 46: if( (++*c) >= IMATCH ){
! 47: if ( ! (IMATC)(*Im)->next ) {
! 48: /* no more IENT */
! 49: *c = -1;
! 50: return;
! 51: }
! 52: *Im = (IMATC)(*Im)->next;
! 53: *c = 0;
! 54: }
! 55: *ent = (*Im)->ient[*c];
! 56: }
! 57:
! 58: void GetForeIent(Im, ent, c)
! 59: IMATC *Im;
! 60: IENT *ent;
! 61: int *c;
! 62: {
! 63: /* GetForeIent(&Im, &ent, &c); */
! 64:
! 65: /* find last IENT */
! 66: if( (--*c) < 0 ){
! 67: if ( !(IMATC)(*Im)->fore ){
! 68: /* no more IENT */
! 69: *c = -1;
! 70: return;
! 71: }
! 72: /* next IENT block */
! 73: *Im = (IMATC)(*Im)->fore;
! 74: *c = IMATCH - 1;
! 75: }
! 76: *ent = (*Im)->ient[*c];
! 77: }
! 78:
! 79: void AppendIent(m, row, col, body)
! 80: IMAT m;
! 81: int row, col;
! 82: Obj body;
! 83: {
! 84: /* append row, col, body to matrix m */
! 85: IMATC Im, Imt;
! 86: int d;
! 87:
! 88: if ( ! body ) return;
! 89: if ( ! m->clen ) {
! 90: NEWIENT(Im);
! 91: m->root = (pointer)Im;
! 92: m->toor = (pointer)Im;
! 93: m->clen = 1;
! 94: MEnt(row * m->row + col, row, col, body, &Im->ient[0]);
! 95: Im->ient[1] = zent;
! 96: } else {
! 97: Im = (pointer)m->toor;
! 98: for( d = 0; Im->ient[d].cr != -1; d++);
! 99: if ( d == IMATCH - 1 ) {
! 100: /* append point is chank end */
! 101: NEWIENT(Imt);
! 102: Im->next = (pointer)Imt;
! 103: m->toor = (pointer)Imt;
! 104: ++(m->clen);
! 105: Imt->fore = (pointer)Im;
! 106: Imt->next = 0;
! 107: MEnt(row * m->row + col, row, col, body, &Im->ient[d]);
! 108: Imt->ient[0] = zent;
! 109: } else {
! 110: MEnt(row * m->row + col, row, col, body, &Im->ient[d]);
! 111: Im->ient[d + 1] = zent;
! 112: }
! 113: }
! 114: }
! 115:
! 116: void PutIent(m, row, col, trg)
! 117: IMAT m;
! 118: int row, col;
! 119: Obj trg;
! 120: {
! 121: /* insert or replace IENT */
! 122: IMATC Im, Imn;
! 123: IENT ent, tent;
! 124: int cr, c, d;
! 125:
! 126: if ( m->row <= row || m->col <= col || row < 0 || col < 0 )
! 127: error("putim : Out of rage");
! 128: cr = row * m->row + col;
! 129: if ( ! m->clen ) {
! 130: if( trg == 0 ) return;
! 131: AppendIent(m, row, col, trg);
! 132: return;
! 133: }
! 134: MEnt(cr, row, col, trg, &tent);
! 135: Im = (pointer)m->root;
! 136: c = -1;
! 137: GetNextIent(&Im, &ent, &c);
! 138: while (ent.cr < cr) {
! 139: GetNextIent(&Im, &ent, &c);
! 140: if( c == -1 ) {
! 141: /* zero insert case */
! 142: if ( ! trg ) return;
! 143: /* last append case */
! 144: AppendIent(m, row, col, trg);
! 145: return;
! 146: }
! 147: }
! 148: if ( ent.cr == cr ) {
! 149: /* same row and col (replace) case */
! 150: if ( ! trg ) {
! 151: /* IENT remove case */
! 152: Imn = Im;
! 153: while ( c != -1 ) {
! 154: GetNextIent( &Im, &ent, &c );
! 155: if ( ! c ) {
! 156: Imn->ient[IMATCH-1] = ent;
! 157: Imn = Im;
! 158: } else {
! 159: Im->ient[c - 1] = ent;
! 160: }
! 161: }
! 162: } else {
! 163: /* replase case */
! 164: Im->ient[c] = tent;
! 165: }
! 166: } else {
! 167: /* IENT insert case */
! 168: while( c != -1 ) {
! 169: Im->ient[c] = tent;
! 170: tent = ent;
! 171: GetNextIent(&Im, &ent, &c);
! 172: }
! 173:
! 174: }
! 175: }
! 176:
! 177: void Pnewimat(arg, rp)
! 178: NODE arg;
! 179: IMAT *rp;
! 180: {
! 181: /* make new index type matrix for parser */
! 182: int row,col;
! 183: IMAT m;
! 184:
! 185: asir_assert(ARG0(arg),O_N,"newimat");
! 186: asir_assert(ARG1(arg),O_N,"newimat");
! 187: row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));
! 188: if ( row <= 0 || col <= 0 ) error("newimat : invalid size");
! 189: NEWIMAT(m);
! 190: m->col = col;
! 191: m->row = row;
! 192: *rp = m;
! 193: }
! 194:
! 195: void GetIbody(m, row, col, trg)
! 196: IMAT m;
! 197: int col, row;
! 198: Obj *trg;
! 199: {
! 200: /* get entry body from m for parser */
! 201: IMATC Im;
! 202: IENT ent;
! 203: int cr;
! 204: int c;
! 205:
! 206: if ( m->row <= row || m->col <= col || row < 0 || col < 0 )
! 207: error("putim : Out of rage");
! 208: if ( ! m->clen ) {
! 209: /* zero matrix case */
! 210: *trg = (Obj)0;
! 211: } else {
! 212: cr = row * m->row + col;
! 213: c = -1;
! 214: Im = (pointer)m->root;
! 215: GetNextIent( &Im, &ent, &c);
! 216: while( ent.cr < cr ) {
! 217: if ( ent.cr == -1 ) {
! 218: /* not fined ent to last case */
! 219: *trg = (Obj)0;
! 220: return;
! 221: }
! 222: GetNextIent( &Im, &ent, &c);
! 223: }
! 224: if ( ent.cr == cr ) *trg = (Obj)ent.body;
! 225: else *trg = (Obj)0; /* not fined end to mid case */
! 226: }
! 227: }
! 228:
! 229: void PChsgnI(arg, rp)
! 230: NODE arg;
! 231: IMAT *rp;
! 232: {
! 233: /* index type matrix all entry sgn cheng for parser */
! 234: VL vl;
! 235: IMAT m, n;
! 236:
! 237: asir_assert(ARG0(arg),O_IMAT,"chsgnind");
! 238: vl = CO;
! 239: m = (IMAT)ARG0(arg);
! 240: ChsgnI(m, &n);
! 241: *rp = n;
! 242: }
! 243:
! 244: void ChsgnI(a, c)
! 245: IMAT a;
! 246: IMAT *c;
! 247: {
! 248: /* index type matrix all entry sgn chg */
! 249: IMAT b;
! 250: IMATC ac;
! 251: IENT aent;
! 252: Obj r;
! 253: int ai;
! 254:
! 255: if( ! a->root ){
! 256: /* zero matrix case */
! 257: *c = a;
! 258: } else {
! 259: NEWIMAT(b);
! 260: b->col = a->col;
! 261: b->row = a->row;
! 262: ai = -1;
! 263: ac = (pointer)a->root;
! 264: GetNextIent(&ac, &aent, &ai);
! 265: while(aent.cr != -1){
! 266: arf_chsgn((Obj)aent.body, (Obj *)&r);
! 267: AppendIent(b, aent.row, aent.col, r);
! 268: GetNextIent(&ac, &aent, &ai);
! 269: }
! 270: *c = b;
! 271: }
! 272: }
! 273:
! 274: void Pm2Im(arg, rp)
! 275: NODE arg;
! 276: IMAT *rp;
! 277: {
! 278: /* matrix type convert from MAT to IMAT */
! 279: int i,j, row, col;
! 280: MAT m;
! 281: IMAT Im;
! 282: pointer *a;
! 283:
! 284: m = (MAT)ARG0(arg);
! 285: row = m->row;
! 286: col = m->col;
! 287: NEWIMAT(Im);
! 288: Im->col = m->col;
! 289: Im->row = m->row;
! 290: for (i = 0; i < row; i++)
! 291: for (j = 0, a = BDY(m)[i]; j < col; j++)
! 292: if (a[j]) AppendIent(Im, i, j, a[j]);
! 293: *rp = Im;
! 294: }
! 295:
! 296: void PIm2m(arg, rp)
! 297: NODE arg;
! 298: MAT *rp;
! 299: {
! 300: /* matrix type convert from IMAT to MAT */
! 301: IMAT Im;
! 302: MAT m;
! 303: int c, i, j;
! 304: IMATC Imc;
! 305: IENT ent;
! 306:
! 307: Im = (IMAT)ARG0(arg);
! 308: MKMAT(m,Im->row,Im->col);
! 309: if ( Im->root ) {
! 310: /* non zero matrix case */
! 311: Imc = (pointer)Im->root;
! 312: c = -1;
! 313: GetNextIent(&Imc, &ent, &c);
! 314: while (ent.cr != -1) {
! 315: BDY(m)[ent.row][ent.col] = ent.body;
! 316: GetNextIent(&Imc, &ent, &c);
! 317: }
! 318: }
! 319: *rp = m;
! 320: }
! 321:
! 322: void AddMatI(vl, a, b, c)
! 323: VL vl;
! 324: IMAT a, b, *c;
! 325: {
! 326: /* add index type matrix */
! 327: int ai, bi;
! 328: IMAT m;
! 329: IMATC ac, bc;
! 330: IENT aent, bent;
! 331: Obj pc;
! 332:
! 333: if ( a->col != b->col ) {
! 334: error("addmati : colum size mismatch");
! 335: *c = 0;
! 336: } else if ( a->row != b->row ) {
! 337: error("addmati : row size mismatch");
! 338: *c = 0;
! 339: } else if ( ! a->root ) {
! 340: *c = b; /* a : zero matrox */
! 341: } else if ( ! b->root ) {
! 342: *c = a; /* b : zero matrix */
! 343: } else {
! 344: NEWIMAT(m);
! 345: m->col = a->col;
! 346: m->row = a->row;
! 347: ai = -1;
! 348: ac = (pointer)a->root;
! 349: bi = -1;
! 350: bc = (pointer)b->root;
! 351: GetNextIent(&ac, &aent, &ai);
! 352: GetNextIent(&bc, &bent, &bi);
! 353: while( aent.cr != -1 ) {
! 354: if( aent.cr == bent.cr ) {
! 355: arf_add(vl, (Obj)aent.body, (Obj)bent.body, (Obj *)&pc);
! 356: AppendIent( m, aent.row, aent.col, pc);
! 357: GetNextIent( &ac, &aent, &ai );
! 358: GetNextIent( &bc, &bent, &bi );
! 359: } else if ( aent.cr < bent.cr ) {
! 360: AppendIent( m, aent.row, aent.col, (Obj)aent.body);
! 361: GetNextIent( &ac, &aent, &ai);
! 362: } else if ( bent.cr == -1 ) {
! 363: AppendIent( m, aent.row, aent.col, (Obj)aent.body);
! 364: GetNextIent( &ac, &aent, &ai);
! 365: } else {
! 366: AppendIent( m, bent.row, bent.col, (Obj)bent.body);
! 367: GetNextIent( &bc, &bent, &bi);
! 368: if ( bent.cr == -1 ){
! 369: while(aent.cr != -1){
! 370: AppendIent(m, aent.row, aent.col, (Obj)aent.body);
! 371: GetNextIent( &ac, &aent, &ai );
! 372: }
! 373: break;
! 374: }
! 375: }
! 376: }
! 377: while(bent.cr != -1) {
! 378: AppendIent( m, bent.row, bent.col, (Obj)bent.body);
! 379: GetNextIent( &bc, &bent, &bi );
! 380: }
! 381: *c = m;
! 382: }
! 383: return;
! 384: }
! 385:
! 386: void SubMatI(vl, a, b, c)
! 387: VL vl;
! 388: IMAT a, b, *c;
! 389: {
! 390: /* subtract index type matrix */
! 391: int ai, bi;
! 392: IMAT m;
! 393: IMATC ac, bc;
! 394: IENT aent, bent;
! 395: Obj obj;
! 396:
! 397: if ( a->col != b->col ) {
! 398: error("submati : colum size mismatch");
! 399: *c = 0;
! 400: } else if ( a->row != b->row ) {
! 401: error("submati : row size mismatch");
! 402: *c = 0;
! 403: } else if ( ! b->root ) {
! 404: *c = a;
! 405: } else if ( ! a->root ) {
! 406: ChsgnI(b, &m);
! 407: *c = m;
! 408: } else {
! 409: vl = CO;
! 410: NEWIMAT(m);
! 411: m->col = a->col;
! 412: m->row = a->row;
! 413: ai = -1;
! 414: ac = (pointer)a->root;
! 415: bi = -1;
! 416: bc = (pointer)b->root;
! 417:
! 418: GetNextIent(&ac, &aent, &ai);
! 419: GetNextIent(&bc, &bent, &bi);
! 420: while(aent.cr != -1) {
! 421: if( aent.cr == bent.cr ) {
! 422: arf_sub(vl, (Obj)aent.body, (Obj)bent.body, (Obj *)&obj);
! 423: AppendIent(m, aent.row, aent.col, obj);
! 424: GetNextIent(&ac, &aent, &ai);
! 425: GetNextIent(&bc, &bent, &bi);
! 426: } else if ( aent.cr < bent.cr ) {
! 427: AppendIent( m, aent.row, aent.col, (Obj)aent.body);
! 428: GetNextIent(&ac, &aent, &ai);
! 429: } else if ( bent.cr == -1 ) {
! 430: AppendIent( m, aent.row, aent.col, (Obj)aent.body);
! 431: GetNextIent(&ac, &aent, &ai);
! 432: } else {
! 433: arf_chsgn((Obj)bent.body, (Obj *)&obj);
! 434: AppendIent( m, bent.row, bent.col, (Obj)obj);
! 435: GetNextIent(&bc, &bent, &bi);
! 436: if (bent.cr == -1){
! 437: while(aent.cr != -1){
! 438: AppendIent(m, aent.row, aent.col, (Obj)aent.body);
! 439: GetNextIent(&ac, &aent, &ai);
! 440: }
! 441: break;
! 442: }
! 443: }
! 444: }
! 445: while(bent.cr != -1) {
! 446: arf_chsgn((Obj)bent.body, (Obj *)&obj);
! 447: AppendIent( m, bent.row, bent.col, (Obj)obj);
! 448: GetNextIent(&bc, &bent, &bi);
! 449: }
! 450: *c = m;
! 451: }
! 452: }
! 453:
! 454: void MulrMatI(vl, a, b, rp)
! 455: VL vl;
! 456: Obj a, b, *rp;
! 457: {
! 458: /* multiply a expression and a index type matrix */
! 459: IMAT m;
! 460: IENT ent;
! 461: IMATC Im;
! 462: Obj p;
! 463: int ia;
! 464:
! 465: NEWIMAT(m);
! 466: m->col = ((IMAT)b)->col;
! 467: m->row = ((IMAT)b)->row;
! 468: Im = (IMATC)(((IMAT)b)->root);
! 469: ia = -1;
! 470: GetNextIent(&Im, &ent, &ia);
! 471: while(ent.cr != -1) {
! 472: arf_mul(vl, (Obj)a, (Obj)ent.body, (Obj *)&p);
! 473: AppendIent(m, ent.row, ent.col, (Obj)p);
! 474: GetNextIent(&Im, &ent, &ia);
! 475: }
! 476: *rp = (Obj)m;
! 477: }
! 478:
! 479: void MulMatG(vl, a, b, c)
! 480: VL vl;
! 481: Obj a, b, *c;
! 482: {
! 483: /* multiply index type matrix general procedure */
! 484: IMAT m;
! 485:
! 486: if ( !a ) {
! 487: NEWIMAT(m);
! 488: m->col = ((IMAT)b)->col;
! 489: m->row = ((IMAT)b)->row;
! 490: *c = (Obj)m;
! 491: } else if ( !b ) {
! 492: NEWIMAT(m);
! 493: m->col = ((IMAT)a)->col;
! 494: m->row = ((IMAT)a)->row;
! 495: *c = (Obj)m;
! 496: } else if ( OID(a) <= O_R ) MulrMatI(vl, (Obj)a, (Obj)b, (Obj *)c);
! 497: else if ( OID(b) <= O_R ) MulrMatI(vl, (Obj)b, (Obj)a, (Obj *)c);
! 498: else MulMatS(vl, (IMAT)a, (IMAT)b, (IMAT *)c);
! 499: }
! 500:
! 501: void MulMatS(vl, m,n,rp)
! 502: VL vl;
! 503: IMAT m, n, *rp;
! 504: {
! 505: /* multiply index type matrix and index type matrix */
! 506: IMAT r, a11, a12, a21, a22, b11, b12, b21, b22;
! 507: IMAT ans1, ans2, c11, c12, c21, c22, s1, s2, t1, t2, u1, v1, w1;
! 508: pointer u, v, *ab;
! 509: int hmcol, hmrow, hncol, hnrow;
! 510: IENT ent, entn, entm;
! 511: IMATC in,im;
! 512: int *c, ai, bi, cr, row, col, ca11,ca12,ca21,ca22;
! 513: if ( m->col != n->row ) {
! 514: *rp =0; error("mulmati : size mismatch");
! 515: }
! 516: vl = CO;
! 517: NEWIMAT(r);
! 518: r->row = m->row;
! 519: r->col = n->col;
! 520: if( m->clen == 0 || n->clen == 0){
! 521: /* zero matrix case */
! 522: *rp = (pointer)r;
! 523: return;
! 524: } else if(StrassenSize == 0) {
! 525: /* not use Strassen algorithm */
! 526: MulMatI(vl, m, n, &r);
! 527: *rp = (pointer)r;
! 528: return;
! 529: #if 0
! 530: } else if(m->row == 1){
! 531: ab = (pointer)MALLOC((n->col +1)*sizeof(Obj));
! 532: ai = -1;
! 533: im = (pointer)m->root;
! 534: GetNextIent(&im, &entm, &ai);
! 535: bi = -1;
! 536: in = (pointer)n->root;
! 537: GetNextIent(&in, &entn, &bi);
! 538: while( entn.cr != -1 ){
! 539: if( entn.row == entm.col ){
! 540: arf_mul(vl,(Obj)entm.body,(Obj)entm.body,(Obj *)&u);
! 541: arf_add(vl,(Obj)u,(Obj)ab[entn.col],(Obj *)&v);
! 542: ab[entn.col] = (pointer)v;
! 543: GetNextIent(&in, &entn, &bi);
! 544: } else if( entn.row < entm.col ) {
! 545: GetNextIent(&in, &entn, &bi);
! 546: } else {
! 547: GetNextIent(&im, &entm, &ai);
! 548: }
! 549: }
! 550: for(ai=0; ai< m->col; ai++){
! 551: if(ab[ai] != 0) AppendIent(r, 1, ai, (Obj)&ab[ai]);
! 552: }
! 553: *rp=r;
! 554: return;
! 555: } else if(n->col == 1){
! 556: /* not yet */
! 557: *rp=0;
! 558: return;
! 559: } else if(m->col == 1){
! 560: ai = -1;
! 561: im = (pointer)m->root;
! 562: GetNextIent(&im, &entm, &ai);
! 563: while( entm.cr != -1 ){
! 564: bi = -1;
! 565: in = (pointer)n->root;
! 566: GetNextIent(&in, &entn, &bi);
! 567: while( entn.cr != -1 ){
! 568: arf_mul(vl,(Obj)entm.body,(Obj)entn.body,(Obj *)&u);
! 569: AppendIent(r, entm.row, entn.col, (Obj)&u);
! 570: GetNextIent(&in, &entn, &bi);
! 571: }
! 572: GetNextIent(&im, &entm, &ai);
! 573: }
! 574: *rp = r;
! 575: return;
! 576: #endif
! 577: } else if((m->row<=StrassenSize)||(m->col<=StrassenSize)|| (n->col<=StrassenSize)) {
! 578: /* not use Strassen algorithm */
! 579: MulMatI(vl, m, n, &r);
! 580: *rp = (pointer)r;
! 581: return;
! 582: }
! 583: /* Strassen Algorithm */
! 584: /* calcrate matrix half size */
! 585: hmrow = (m->row + (m->row & 1)) >> 1;
! 586: hmcol = (m->col + (m->col & 1)) >> 1;
! 587: hnrow = (n->row + (n->row & 1)) >> 1;
! 588: hncol = (n->col + (n->col & 1)) >> 1;
! 589: NEWIMAT(a11); a11->col = hmcol; a11->row = hmrow;
! 590: NEWIMAT(a12); a12->col = hmcol; a12->row = hmrow;
! 591: NEWIMAT(a21); a21->col = hmcol; a21->row = hmrow;
! 592: NEWIMAT(a22); a22->col = hmcol; a22->row = hmrow;
! 593: NEWIMAT(b11); b11->col = hncol; b11->row = hnrow;
! 594: NEWIMAT(b12); b12->col = hncol; b12->row = hnrow;
! 595: NEWIMAT(b21); b21->col = hncol; b21->row = hnrow;
! 596: NEWIMAT(b22); b22->col = hncol; b22->row = hnrow;
! 597:
! 598: /* copy matrix n to 4 small matrix a11,a12,a21,a22 */
! 599: im = (pointer)m->root;
! 600: ai = -1;
! 601: GetNextIent(&im, &ent, &ai);
! 602: while( ent.cr != -1 ){
! 603: if(ent.col < hmcol){
! 604: if(ent.row < hmrow){
! 605: AppendIent(a11,ent.row,ent.col,(Obj)ent.body);
! 606: } else {
! 607: AppendIent(a21,ent.row-hmrow,ent.col,(Obj)ent.body);
! 608: }
! 609: } else {
! 610: if(ent.row < hmrow){
! 611: AppendIent(a12,ent.row,ent.col-hmcol,(Obj)ent.body);
! 612: } else {
! 613: AppendIent(a22,ent.row-hmrow,ent.col-hmcol,(Obj)ent.body);
! 614: }
! 615: }
! 616: GetNextIent(&im, &ent, &ai);
! 617: }
! 618: /* copy matrix m to 4 small matrix b11,b12,b21,b22 */
! 619: im = (pointer)n->root;
! 620: ai = -1;
! 621: GetNextIent(&im, &ent, &ai);
! 622: while( ent.cr != -1 ){
! 623: if(ent.col < hmcol){
! 624: if(ent.row < hnrow){
! 625: AppendIent(b11,ent.row,ent.col,(Obj)ent.body);
! 626: } else {
! 627: AppendIent(b21,ent.row-hnrow,ent.col,(Obj)ent.body);
! 628: }
! 629: } else {
! 630: if(ent.row < hnrow){
! 631: AppendIent(b12,ent.row,ent.col-hncol,(Obj)ent.body);
! 632: } else {
! 633: AppendIent(b22,ent.row-hnrow,ent.col-hncol,(Obj)ent.body);
! 634: }
! 635: }
! 636: GetNextIent(&im, &ent, &ai);
! 637: }
! 638: /* expand matrix by Strassen-Winograd algorithm */
! 639: /* s1 = A21 + A22 */
! 640: AddMatI(vl,a21,a22, &s1);
! 641:
! 642: /* s2=s1-A11 */
! 643: SubMatI(vl,s1,a11,&s2);
! 644:
! 645: /* t1=B12-B11 */
! 646: SubMatI(vl,b12,b11,&t1);
! 647:
! 648: /* t2=B22-t1 */
! 649: SubMatI(vl,b22,t1,&t2);
! 650:
! 651: /* u=(A11-A21)*(B22-B12) */
! 652: SubMatI(vl,a11,a21,&ans1);
! 653: SubMatI(vl,b22,b12,&ans2);
! 654: MulMatS(vl, ans1,ans2,&u1);
! 655:
! 656: /* v=s1*t1 */
! 657: MulMatS(vl, s1,t1,&v1);
! 658:
! 659: /* w=A11*B11+s2*t2 */
! 660: MulMatS(vl, a11,b11,&ans1);
! 661: MulMatS(vl, s2,t2,&ans2);
! 662: AddMatI(vl,ans1,ans2,&w1);
! 663:
! 664: /* C11 = A11*B11+A12*B21 */
! 665: MulMatS(vl, a12,b21,&ans2);
! 666: AddMatI(vl,ans1,ans2,&c11);
! 667:
! 668: /* C12 = w1+v1+(A12-s2)*B22 */
! 669: SubMatI(vl,a12,s2,&ans1);
! 670: MulMatS(vl, ans1, b22, &ans2);
! 671: AddMatI(vl, w1, v1, &ans1);
! 672: AddMatI(vl, ans1, ans2, &c12);
! 673:
! 674: /* C21 = w1+u1+A22*(B21-t2) */
! 675: SubMatI(vl, b21, t2, &ans1);
! 676: MulMatS(vl, a22, ans1, &ans2);
! 677: AddMatI(vl, w1, u1, &ans1);
! 678: AddMatI(vl, ans1, ans2, &c21);
! 679:
! 680: /* C22 = w1 + u1 + v1 */
! 681: AddMatI(vl, ans1, v1, &c22);
! 682:
! 683: /* create return i matrix */
! 684: r->col = m->col;
! 685: r->row = n->row;
! 686:
! 687: /* copy c11 to r */
! 688: ca11 = -1;
! 689: if( c11->root ){
! 690: im = (pointer)c11->root;
! 691: ai = -1;
! 692: GetNextIent(&im, &ent, &ai);
! 693: while( ai != -1 ){
! 694: AppendIent(r,ent.row,ent.col,(Obj)ent.body);
! 695: GetNextIent(&im, &ent, &ai);
! 696: }
! 697: }
! 698: /* copy c21 to r */
! 699: if( c21->root ){
! 700: im = (pointer)c21->root;
! 701: ai = -1;
! 702: GetNextIent(&im, &ent, &ai);
! 703: while( ent.cr != -1 ){
! 704: PutIent(r, ent.row + hmrow, ent.col, (Obj)ent.body);
! 705: GetNextIent(&im, &ent, &ai);
! 706: }
! 707: }
! 708: /* copy c12 to r */
! 709: if( c12->root ){
! 710: im = (pointer)c12->root;
! 711: ai = -1;
! 712: GetNextIent(&im, &ent, &ai);
! 713: while( ent.cr != -1 ){
! 714: PutIent(r, ent.row, ent.col + hncol, (Obj)ent.body);
! 715: GetNextIent(&im, &ent, &ai);
! 716: }
! 717: }
! 718: /* copy c22 to r */
! 719: if( c22->root ){
! 720: im = (pointer)c22->root;
! 721: ai = -1;
! 722: GetNextIent(&im, &ent, &ai);
! 723: while( ent.cr != -1 ){
! 724: PutIent(r, ent.row + hmrow, ent.col + hncol, (Obj)ent.body);
! 725: GetNextIent(&im, &ent, &ai);
! 726: }
! 727: }
! 728: *rp = (pointer)r;
! 729: return;
! 730: }
! 731:
! 732: void MulMatI(vl,m,n,r)
! 733: VL vl;
! 734: IMAT m, n, *r;
! 735: {
! 736: /* inner prodocut algorithm for index type multiply procedure */
! 737: IMAT l;
! 738: IMATC im, imt;
! 739: IENT ent, tent;
! 740: int bc, bend, ac, aend;
! 741: int col, row, tcol, trow;
! 742: int i, j, k, ai, bi;
! 743: int bj, ik,kj;
! 744: pointer s, u, v;
! 745:
! 746: typedef struct oIND {
! 747: int row, col;
! 748: pointer body;
! 749: } *IND;
! 750: IND a,b;
! 751: /* make trans(n) */
! 752: a = (IND)MALLOC(IMATCH * m->clen *sizeof(struct oIND)+1);
! 753: b = (IND)MALLOC(IMATCH * n->clen *sizeof(struct oIND)+1);
! 754:
! 755: bend = -1;
! 756: row = n->row;
! 757: col = n->col;
! 758: bc = -1;
! 759: im = (pointer)n->root;
! 760: GetNextIent(&im, &ent, &bc);
! 761: b[++bend].body = ent.body;
! 762: b[bend].row = ent.col;
! 763: b[bend].col = ent.row;
! 764: GetNextIent(&im, &ent, &bc);
! 765: while ( ent.cr != -1 ) {
! 766: trow = ent.col;
! 767: tcol = ent.row;
! 768: for( i = bend; i >= 0; i--) {
! 769: if( b[i].row > trow ) {
! 770: b[i + 1] = b[i];
! 771: } else if ( b[i].col > tcol ) {
! 772: b[i + 1] = b[i];
! 773: } else {
! 774: b[i+1].col = tcol;
! 775: b[i+1].row = trow;
! 776: b[i+1].body = ent.body;
! 777: break;
! 778: }
! 779: }
! 780: if ( i == -1 ) {
! 781: b[0].col = tcol;
! 782: b[0].row = trow;
! 783: b[0].body = ent.body;
! 784: }
! 785: bend++;
! 786: GetNextIent(&im, &ent, &bc);
! 787: }
! 788:
! 789: im = (pointer)m->root;
! 790: ac = -1;
! 791: GetNextIent(&im, &ent, &ac);
! 792: aend = -1;
! 793: while( ent.cr != -1 ){
! 794: a[++aend].col = ent.col;
! 795: a[aend].row = ent.row;
! 796: a[aend].body = ent.body;
! 797: GetNextIent(&im, &ent, &ac);
! 798: }
! 799: /* make return matrix */
! 800: NEWIMAT(l);
! 801: l->row = m->row;
! 802: l->col = n->col;
! 803: ac = -1;
! 804: for( i = 0; i<= aend;) {
! 805: ai = a[i].row;
! 806: bj = b[0].row;
! 807: s = 0;
! 808: ik=i;
! 809: for(j=0; j<= bend;){
! 810: if( a[ik].col == b[j].col) {
! 811: arf_mul(CO,(Obj)a[ik].body, (Obj)b[j].body, (Obj *)&u);
! 812: arf_add(CO,(Obj)s,(Obj)u,(Obj *)&v);
! 813: s = v;
! 814: j++;
! 815: if ( bj != b[j].row ) {
! 816: AppendIent(l, ai, bj, (Obj)s);
! 817: ik = i;
! 818: bj = b[j].row;
! 819: s = 0;
! 820: } else {
! 821: ik++;
! 822: if ( ai != a[ik].row ) {
! 823: AppendIent(l, ai, bj, (Obj)s);
! 824: s = 0;
! 825: while ( bj == b[j].row ) {
! 826: j++;
! 827: if ( j > bend ) break;
! 828: }
! 829: ik = i;
! 830: bj = b[j].row;
! 831: }
! 832: }
! 833: } else if ( a[ik].col > b[j].col ) {
! 834: j++;
! 835: if ( bj != b[j].row ) {
! 836: AppendIent(l, ai, bj, (Obj)s);
! 837: s = 0;
! 838: ik = i;
! 839: bj = b[j].row;
! 840: }
! 841: } else /* if ( a[ik].col < b[j].col ) */ {
! 842: ik++;
! 843: if ( ai != a[ik].row ) {
! 844: AppendIent(l, ai, bj, (Obj)s);
! 845: s = 0;
! 846: while ( bj == b[j].row ) {
! 847: j++;
! 848: if ( j > bend ) break;
! 849: }
! 850: bj = b[j].row;
! 851: ik = i;
! 852: }
! 853: }
! 854: }
! 855: i = ik;
! 856: while ( ai == a[i].row ) {
! 857: i++;
! 858: if( i > aend) break;
! 859: }
! 860: }
! 861: *r = (IMAT)l;
! 862: FREE(a);
! 863: FREE(b);
! 864: return;
! 865: }