![[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: /*
1.2 ! saito 2: * $OpenXM: OpenXM_contrib2/asir2000/builtin/iarray.c,v 1.1 2004/12/18 16:50:10 saito Exp $
1.1 saito 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;
1.2 ! saito 94: MEnt(row * m->col + col, row, col, body, &Im->ient[0]);
1.1 saito 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;
1.2 ! saito 107: MEnt(row * m->col + col, row, col, body, &Im->ient[d]);
1.1 saito 108: Imt->ient[0] = zent;
109: } else {
1.2 ! saito 110: MEnt(row * m->col + col, row, col, body, &Im->ient[d]);
1.1 saito 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");
1.2 ! saito 128: cr = row * m->col + col;
! 129: printf("cr = %d\n",cr);
1.1 saito 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 {
1.2 ! saito 213: cr = row * m->col + col;
1.1 saito 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: }