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