/*
 * $OpenXM: OpenXM_contrib2/asir2018/builtin/iarray.c,v 1.2 2018/09/28 08:20:27 noro Exp $
*/
#include "ca.h"
#include "base.h"
#include "parse.h"
#include "inline.h"

extern int StrassenSize;

struct ftab imat_tab[] = {
  {"newimat", Pnewimat,2},
  {"m2im",  Pm2Im,1},
  {"im2m",  PIm2m,1},
  {0,0,0},
};

const IENT zent = { -1, -1, -1, 0}; /* last ient const */

void MEnt(cr, row, col, trg, ent)
int cr, row, col;
Obj trg;
IENT *ent;
{
  /* make Index Entry IENT */
  /* set:cr, row, trag */
  /* return: value ent */
  /* MEnt(cr, row, col, trg, &ent); */

  ent->cr = cr;
  ent->row = row;
  ent->col = col;
  ent->body = (pointer)trg;
}

void GetNextIent(Im, ent, c)
IMATC *Im;
IENT *ent;
int *c;
{
  /* get next IENT */
  /* set: Im, c */
  /* return: ent, c */
  /* GetNextIent(&Im, &ent, &c); */

  if( (++*c) >= IMATCH ){
    if ( ! (IMATC)(*Im)->next ) {
    /* no more IENT */
      *c = -1;
      return;
    }
    *Im = (IMATC)(*Im)->next;
    *c = 0;
  }
  *ent = (*Im)->ient[*c];
}

void GetForeIent(Im, ent, c)
IMATC *Im;
IENT *ent;
int *c;
{
  /* GetForeIent(&Im, &ent, &c); */

  /* find last IENT */
  if( (--*c) < 0 ){
    if ( !(IMATC)(*Im)->fore ){
    /* no more IENT */
      *c = -1;
      return;
    }
    /* next IENT block */
    *Im = (IMATC)(*Im)->fore;
    *c = IMATCH - 1;
  }
  *ent = (*Im)->ient[*c];
}

void AppendIent(m, row, col, body)
IMAT m;
int row, col;
Obj body;
{
  /* append row, col, body to matrix m */
  IMATC Im, Imt;
  int d;

  if ( ! body ) return;
  if ( ! m->clen ) {
    NEWIENT(Im);
    m->root = (pointer)Im;
    m->toor = (pointer)Im;
    m->clen = 1;
    MEnt(row * m->col + col, row, col, body, &Im->ient[0]);
    Im->ient[1] = zent;
  } else {
    Im = (pointer)m->toor;
    for( d = 0; Im->ient[d].cr != -1; d++);
    if ( d == IMATCH - 1 ) {
      /* append point is chank end */
      NEWIENT(Imt);
      Im->next = (pointer)Imt;
      m->toor = (pointer)Imt;
      ++(m->clen);
      Imt->fore = (pointer)Im;
      Imt->next = 0;
      MEnt(row *  m->col + col, row, col, body, &Im->ient[d]);
      Imt->ient[0] = zent;
    } else {
      MEnt(row * m->col + col, row, col, body, &Im->ient[d]);
      Im->ient[d + 1] = zent;
    }
  }
}

void PutIent(m, row, col, trg)
IMAT m;
int row, col;
Obj trg;
{
  /* insert or replace IENT */
  IMATC Im, Imn;
  IENT ent, tent;
  int cr, c, d;

  if ( m->row <= row || m->col <= col || row < 0 || col < 0 )
    error("putim : Out of rage");
  cr = row * m->col + col;
printf("cr = %d\n",cr);
  if ( ! m->clen ) {
    if( trg == 0 ) return;
    AppendIent(m, row, col, trg);
    return;
  }
  MEnt(cr, row, col, trg, &tent);
  Im = (pointer)m->root;
  c = -1;
  GetNextIent(&Im, &ent, &c);
  while (ent.cr < cr) {
    GetNextIent(&Im, &ent, &c);
    if( c == -1 ) {
      /* zero insert case */
      if ( ! trg ) return;
      /* last append case */
      AppendIent(m, row, col, trg);
      return;
    }
  }
  if ( ent.cr == cr ) {
    /* same row and col (replace) case */
    if ( ! trg ) {
      /* IENT remove case */
      Imn = Im;
      while ( c != -1 ) {
        GetNextIent( &Im, &ent, &c );
        if ( ! c ) {
          Imn->ient[IMATCH-1] = ent;
          Imn = Im;
        } else {
          Im->ient[c - 1] = ent;
        }
      }
    } else {
      /* replase case */
      Im->ient[c] = tent;
    }
  } else {
    /* IENT insert case */
    while( c != -1 ) {
      Im->ient[c] = tent;
      tent = ent;
      GetNextIent(&Im, &ent, &c);
    }

  }
}

void Pnewimat(arg, rp)
NODE arg;
IMAT *rp;
{
  /* make new index type matrix for parser */
  int row,col;
  IMAT m;

  asir_assert(ARG0(arg),O_N,"newimat");
  asir_assert(ARG1(arg),O_N,"newimat");
  row = ZTOS((Q)ARG0(arg)); col = ZTOS((Q)ARG1(arg));
  if ( row <= 0 || col <= 0 ) error("newimat : invalid size");
  NEWIMAT(m);
  m->col = col;
  m->row = row;
  *rp = m;
}

void GetIbody(m, row, col, trg)
IMAT m;
int col, row;
Obj *trg;
{
  /* get entry body from m for parser */
  IMATC Im;
  IENT ent;
  int cr;
  int c;

  if ( m->row <= row || m->col <= col || row < 0 || col < 0 )
    error("putim : Out of rage");
  if ( ! m->clen ) {
  /* zero matrix case */
    *trg = (Obj)0;
  } else {
    cr = row * m->col + col;
    c = -1;
    Im = (pointer)m->root;
    GetNextIent( &Im, &ent, &c);
    while( ent.cr < cr ) {
      if ( ent.cr == -1 ) {
      /* not fined ent to last case */
        *trg = (Obj)0;
        return;
      }
      GetNextIent( &Im, &ent, &c);
    }
    if ( ent.cr == cr ) *trg = (Obj)ent.body;
    else *trg = (Obj)0; /* not fined end to mid case */
  }
}

void PChsgnI(arg, rp)
NODE arg;
IMAT *rp;
{
  /* index type matrix all entry sgn cheng for parser */
  VL vl;
  IMAT m, n;

  asir_assert(ARG0(arg),O_IMAT,"chsgnind");
  vl = CO;
  m = (IMAT)ARG0(arg);
  ChsgnI(m, &n);
  *rp = n;
}

void ChsgnI(a, c)
IMAT a;
IMAT *c;
{
  /* index type matrix all entry sgn chg */
  IMAT b;
  IMATC ac;
  IENT aent;
  Obj r;
  int ai;

  if( ! a->root ){
    /* zero matrix case */
    *c = a;
  } else {
    NEWIMAT(b);
    b->col = a->col;
    b->row = a->row;
    ai = -1;
    ac = (pointer)a->root;
    GetNextIent(&ac, &aent, &ai);
    while(aent.cr != -1){
      arf_chsgn((Obj)aent.body, (Obj *)&r);
      AppendIent(b, aent.row, aent.col, r);
      GetNextIent(&ac, &aent, &ai);
    }
    *c = b;
  }
}

void Pm2Im(arg, rp)
NODE arg;
IMAT *rp;
{
  /* matrix type convert from MAT to IMAT */
  int i,j, row, col;
  MAT m;
  IMAT Im;
  pointer *a;

  m = (MAT)ARG0(arg);
  row = m->row;
  col = m->col;
  NEWIMAT(Im);
  Im->col = m->col;
  Im->row = m->row;
  for (i = 0; i < row; i++)
    for (j = 0, a = BDY(m)[i]; j < col; j++) 
      if (a[j]) AppendIent(Im, i, j, a[j]);
  *rp = Im;
}

void PIm2m(arg, rp)
NODE arg;
MAT *rp;
{
  /* matrix type convert from IMAT to MAT */
  IMAT Im;
  MAT m;
  int c, i, j;
  IMATC Imc;
  IENT ent;

  Im = (IMAT)ARG0(arg);
  MKMAT(m,Im->row,Im->col);
  if ( Im->root ) {
    /* non zero matrix case */
    Imc = (pointer)Im->root;
    c = -1;
    GetNextIent(&Imc, &ent, &c);
    while (ent.cr != -1) {
      BDY(m)[ent.row][ent.col] = ent.body;
      GetNextIent(&Imc, &ent, &c);
    }
  }
  *rp = m;
}

void AddMatI(vl, a, b, c)
VL vl;
IMAT a, b, *c;
{
  /* add index type matrix */
  int ai, bi;
  IMAT m;
  IMATC ac, bc;
  IENT aent, bent;
  Obj pc;

  if ( a->col != b->col ) {
    error("addmati : colum size mismatch");
    *c = 0;
  } else if ( a->row != b->row ) {
    error("addmati : row size mismatch");
    *c = 0;
  } else if ( ! a->root ) {
    *c = b; /* a : zero matrox */
  } else if ( ! b->root ) {
    *c = a; /* b : zero matrix */
  } else {
    NEWIMAT(m);
    m->col = a->col;
    m->row = a->row;
    ai = -1;
    ac = (pointer)a->root;
    bi = -1;
    bc = (pointer)b->root;
    GetNextIent(&ac, &aent, &ai);
    GetNextIent(&bc, &bent, &bi);
    while( aent.cr != -1 ) {
      if( aent.cr == bent.cr ) {
        arf_add(vl, (Obj)aent.body, (Obj)bent.body, (Obj *)&pc);
        AppendIent( m, aent.row, aent.col, pc);
        GetNextIent( &ac, &aent, &ai );
        GetNextIent( &bc, &bent, &bi );
      } else if ( aent.cr < bent.cr ) {
        AppendIent( m, aent.row, aent.col, (Obj)aent.body);
        GetNextIent( &ac, &aent, &ai);
      } else if ( bent.cr == -1 ) {
        AppendIent( m, aent.row, aent.col, (Obj)aent.body);
        GetNextIent( &ac, &aent, &ai);
      } else {
        AppendIent( m, bent.row, bent.col, (Obj)bent.body);
        GetNextIent( &bc, &bent, &bi);
        if ( bent.cr == -1 ){
          while(aent.cr != -1){
            AppendIent(m, aent.row, aent.col, (Obj)aent.body);
            GetNextIent( &ac, &aent, &ai );
          }
          break;
        }
      }
    }
    while(bent.cr != -1) {
      AppendIent( m, bent.row, bent.col, (Obj)bent.body);
      GetNextIent( &bc, &bent, &bi );
    }
    *c = m;
  }
  return;
}

void SubMatI(vl, a, b, c)
VL vl;
IMAT a, b, *c;
{
  /* subtract index type matrix */
  int ai, bi;
  IMAT m;
  IMATC ac, bc;
  IENT aent, bent;
  Obj obj;

  if ( a->col != b->col ) {
    error("submati : colum size mismatch");
    *c = 0;
  } else if ( a->row != b->row ) {
    error("submati : row size mismatch");
    *c = 0;
  } else if ( ! b->root ) {
    *c = a;
  } else if ( ! a->root ) {
    ChsgnI(b, &m);
    *c = m;
  } else {
    vl = CO;
    NEWIMAT(m);
    m->col = a->col;
    m->row = a->row;
    ai = -1;
    ac = (pointer)a->root;
    bi = -1;
    bc = (pointer)b->root;

    GetNextIent(&ac, &aent, &ai);
    GetNextIent(&bc, &bent, &bi);
    while(aent.cr != -1) {
      if( aent.cr == bent.cr ) {
        arf_sub(vl, (Obj)aent.body, (Obj)bent.body, (Obj *)&obj);
        AppendIent(m, aent.row, aent.col, obj);
        GetNextIent(&ac, &aent, &ai);
        GetNextIent(&bc, &bent, &bi);
      } else if ( aent.cr < bent.cr ) {
        AppendIent( m, aent.row, aent.col, (Obj)aent.body);
        GetNextIent(&ac, &aent, &ai);
      } else if ( bent.cr == -1 ) {
        AppendIent( m, aent.row, aent.col, (Obj)aent.body);
        GetNextIent(&ac, &aent, &ai);
      } else {
        arf_chsgn((Obj)bent.body, (Obj *)&obj);
        AppendIent( m, bent.row, bent.col, (Obj)obj);
        GetNextIent(&bc, &bent, &bi);
        if (bent.cr == -1){
          while(aent.cr != -1){
            AppendIent(m, aent.row, aent.col, (Obj)aent.body);
            GetNextIent(&ac, &aent, &ai);
          }
          break;
        }
      }
    }
    while(bent.cr != -1) {
      arf_chsgn((Obj)bent.body, (Obj *)&obj);
      AppendIent( m, bent.row, bent.col, (Obj)obj);
      GetNextIent(&bc, &bent, &bi);
    }
    *c = m;
  }
}

void MulrMatI(vl, a, b, rp)
VL vl;
Obj a, b, *rp;
{
  /* multiply a expression and a index type matrix */
  IMAT m;
  IENT ent;
  IMATC Im;
  Obj p;
  int ia;

  NEWIMAT(m);
  m->col = ((IMAT)b)->col;
  m->row = ((IMAT)b)->row;
  Im = (IMATC)(((IMAT)b)->root);
  ia = -1;
  GetNextIent(&Im, &ent, &ia);
  while(ent.cr != -1) {
    arf_mul(vl, (Obj)a, (Obj)ent.body, (Obj *)&p);
    AppendIent(m, ent.row, ent.col, (Obj)p);
    GetNextIent(&Im, &ent, &ia);
  }
  *rp = (Obj)m;
}

void MulMatG(vl, a, b, c)
VL vl;
Obj a, b, *c;
{
  /* multiply index type matrix general procedure */
  IMAT m;

  if ( !a ) {
    NEWIMAT(m);
    m->col = ((IMAT)b)->col;
    m->row = ((IMAT)b)->row;
    *c = (Obj)m;
  } else if ( !b ) {
    NEWIMAT(m);
    m->col = ((IMAT)a)->col;
    m->row = ((IMAT)a)->row;
    *c = (Obj)m;
  } else if ( OID(a) <= O_R ) MulrMatI(vl, (Obj)a, (Obj)b, (Obj *)c);
  else if ( OID(b) <= O_R ) MulrMatI(vl, (Obj)b, (Obj)a, (Obj *)c);
  else MulMatS(vl, (IMAT)a, (IMAT)b, (IMAT *)c);
}

void MulMatS(vl, m,n,rp)
VL vl;
IMAT m, n, *rp;
{
  /* multiply index type matrix and index type matrix */
  IMAT r, a11, a12, a21, a22, b11, b12, b21, b22;
  IMAT ans1, ans2, c11, c12, c21, c22, s1, s2, t1, t2, u1, v1, w1;
  pointer u, v, *ab;
  int hmcol, hmrow, hncol, hnrow;
  IENT ent, entn, entm;
  IMATC in,im;
  int *c, ai, bi, cr, row, col, ca11,ca12,ca21,ca22;
  if ( m->col != n->row ) {
    *rp =0; error("mulmati : size mismatch");
  }
  vl = CO;
  NEWIMAT(r);
  r->row = m->row;
  r->col = n->col;
  if( m->clen == 0 || n->clen == 0){
  /* zero matrix case */
    *rp = (pointer)r;
    return;
  } else if(StrassenSize == 0) {
  /* not use Strassen algorithm */
    MulMatI(vl, m, n, &r);
    *rp = (pointer)r;
    return;
#if 0
  } else if(m->row == 1){
    ab = (pointer)MALLOC((n->col +1)*sizeof(Obj));
    ai = -1;
    im = (pointer)m->root;
    GetNextIent(&im, &entm, &ai);
    bi = -1;
    in = (pointer)n->root;
    GetNextIent(&in, &entn, &bi);
    while( entn.cr != -1 ){
      if( entn.row == entm.col ){
        arf_mul(vl,(Obj)entm.body,(Obj)entm.body,(Obj *)&u);
        arf_add(vl,(Obj)u,(Obj)ab[entn.col],(Obj *)&v);
        ab[entn.col] = (pointer)v;
        GetNextIent(&in, &entn, &bi);
      } else if( entn.row < entm.col ) {
        GetNextIent(&in, &entn, &bi);
      } else {
        GetNextIent(&im, &entm, &ai);
      }
    }
    for(ai=0; ai< m->col; ai++){
      if(ab[ai] != 0) AppendIent(r, 1, ai, (Obj)&ab[ai]);
    }
    *rp=r;
    return;
  } else if(n->col == 1){
  /* not yet */
    *rp=0;
    return;
  } else if(m->col == 1){
    ai = -1;
    im = (pointer)m->root;
    GetNextIent(&im, &entm, &ai);
    while( entm.cr != -1 ){
      bi = -1;
      in = (pointer)n->root;
      GetNextIent(&in, &entn, &bi);
      while( entn.cr != -1 ){
        arf_mul(vl,(Obj)entm.body,(Obj)entn.body,(Obj *)&u);
        AppendIent(r, entm.row, entn.col, (Obj)&u);
        GetNextIent(&in, &entn, &bi);
      }
      GetNextIent(&im, &entm, &ai);
    }
    *rp = r;
    return;
#endif
  } else if((m->row<=StrassenSize)||(m->col<=StrassenSize)|| (n->col<=StrassenSize)) {
  /* not use Strassen algorithm */
    MulMatI(vl, m, n, &r);
    *rp = (pointer)r;
    return;
  }
  /* Strassen Algorithm */
  /* calcrate matrix half size */
  hmrow = (m->row + (m->row & 1)) >> 1;
  hmcol = (m->col + (m->col & 1)) >> 1;
  hnrow = (n->row + (n->row & 1)) >> 1;
  hncol = (n->col + (n->col & 1)) >> 1;
  NEWIMAT(a11); a11->col = hmcol; a11->row = hmrow;
  NEWIMAT(a12); a12->col = hmcol; a12->row = hmrow;
  NEWIMAT(a21); a21->col = hmcol; a21->row = hmrow;
  NEWIMAT(a22); a22->col = hmcol; a22->row = hmrow;
  NEWIMAT(b11); b11->col = hncol; b11->row = hnrow;
  NEWIMAT(b12); b12->col = hncol; b12->row = hnrow;
  NEWIMAT(b21); b21->col = hncol; b21->row = hnrow;
  NEWIMAT(b22); b22->col = hncol; b22->row = hnrow;

  /* copy matrix n to 4 small matrix a11,a12,a21,a22 */
  im = (pointer)m->root;
  ai = -1;
  GetNextIent(&im, &ent, &ai);
  while( ent.cr != -1 ){
    if(ent.col < hmcol){
      if(ent.row < hmrow){
        AppendIent(a11,ent.row,ent.col,(Obj)ent.body);
      } else {
        AppendIent(a21,ent.row-hmrow,ent.col,(Obj)ent.body);
      }
    } else {
      if(ent.row < hmrow){
        AppendIent(a12,ent.row,ent.col-hmcol,(Obj)ent.body);
      } else {
        AppendIent(a22,ent.row-hmrow,ent.col-hmcol,(Obj)ent.body);
      }
    }
    GetNextIent(&im, &ent, &ai);
  }
  /* copy matrix m to 4 small matrix b11,b12,b21,b22 */
  im = (pointer)n->root;
  ai = -1;
  GetNextIent(&im, &ent, &ai);
  while( ent.cr != -1 ){
    if(ent.col < hmcol){
      if(ent.row < hnrow){
        AppendIent(b11,ent.row,ent.col,(Obj)ent.body);
      } else {
        AppendIent(b21,ent.row-hnrow,ent.col,(Obj)ent.body);
      }
    } else {
      if(ent.row < hnrow){
        AppendIent(b12,ent.row,ent.col-hncol,(Obj)ent.body);
      } else {
        AppendIent(b22,ent.row-hnrow,ent.col-hncol,(Obj)ent.body);
      }
    }
    GetNextIent(&im, &ent, &ai);
  }
  /* expand matrix by Strassen-Winograd algorithm */
  /* s1 = A21 + A22 */
  AddMatI(vl,a21,a22, &s1);

  /* s2=s1-A11 */
  SubMatI(vl,s1,a11,&s2);

  /* t1=B12-B11 */
  SubMatI(vl,b12,b11,&t1);

  /* t2=B22-t1 */
  SubMatI(vl,b22,t1,&t2);

  /* u=(A11-A21)*(B22-B12) */
  SubMatI(vl,a11,a21,&ans1);
  SubMatI(vl,b22,b12,&ans2);
  MulMatS(vl, ans1,ans2,&u1);

  /* v=s1*t1 */
  MulMatS(vl, s1,t1,&v1);

  /* w=A11*B11+s2*t2 */
  MulMatS(vl, a11,b11,&ans1);
  MulMatS(vl, s2,t2,&ans2);
  AddMatI(vl,ans1,ans2,&w1);

  /* C11 = A11*B11+A12*B21 */
  MulMatS(vl, a12,b21,&ans2);
  AddMatI(vl,ans1,ans2,&c11);

  /* C12 = w1+v1+(A12-s2)*B22 */
  SubMatI(vl,a12,s2,&ans1);
  MulMatS(vl, ans1, b22, &ans2);
  AddMatI(vl, w1, v1, &ans1);
  AddMatI(vl, ans1, ans2, &c12);

  /* C21 = w1+u1+A22*(B21-t2) */
  SubMatI(vl, b21, t2, &ans1);
  MulMatS(vl, a22, ans1, &ans2);
  AddMatI(vl, w1, u1, &ans1);
  AddMatI(vl, ans1, ans2, &c21);

  /* C22 = w1 + u1 + v1 */
  AddMatI(vl, ans1, v1, &c22);

  /* create return i matrix */
  r->col = m->col;
  r->row = n->row;

  /* copy c11 to r */
  ca11 = -1;
  if( c11->root ){
    im = (pointer)c11->root;
    ai = -1;
    GetNextIent(&im, &ent, &ai);
    while( ai != -1 ){
      AppendIent(r,ent.row,ent.col,(Obj)ent.body);
      GetNextIent(&im, &ent, &ai);
    }
  }
  /* copy c21 to r */
  if( c21->root ){
    im = (pointer)c21->root;
    ai = -1;
    GetNextIent(&im, &ent, &ai);
    while( ent.cr != -1 ){
      PutIent(r, ent.row + hmrow, ent.col, (Obj)ent.body);
      GetNextIent(&im, &ent, &ai);
    }
  }
  /* copy c12 to r */
  if( c12->root ){
    im = (pointer)c12->root;
    ai = -1;
    GetNextIent(&im, &ent, &ai);
    while( ent.cr != -1 ){
      PutIent(r, ent.row, ent.col + hncol, (Obj)ent.body);
      GetNextIent(&im, &ent, &ai);
    }
  }
  /* copy c22 to r */
  if( c22->root ){
    im = (pointer)c22->root;
    ai = -1;
    GetNextIent(&im, &ent, &ai);
    while( ent.cr != -1 ){
      PutIent(r, ent.row + hmrow, ent.col + hncol, (Obj)ent.body);
      GetNextIent(&im, &ent, &ai);
    }
  }
  *rp = (pointer)r;
  return;
}

void MulMatI(vl,m,n,r)
VL vl;
IMAT m, n, *r;
{
  /* inner prodocut algorithm for index type multiply procedure */
  IMAT l;
  IMATC im, imt;
  IENT ent, tent;
  int bc, bend, ac, aend;
  int col, row, tcol, trow;
  int i, j, k, ai, bi;
  int bj, ik,kj;
  pointer s, u, v;

  typedef struct oIND {
    int row, col;
    pointer body;
  } *IND;
  IND a,b;
  /* make trans(n) */
  a = (IND)MALLOC(IMATCH * m->clen *sizeof(struct oIND)+1);
  b = (IND)MALLOC(IMATCH * n->clen *sizeof(struct oIND)+1);

  bend = -1;
  row = n->row;
  col = n->col;
  bc = -1;
  im = (pointer)n->root;
  GetNextIent(&im, &ent, &bc);
  b[++bend].body = ent.body;
  b[bend].row = ent.col;
  b[bend].col = ent.row;
  GetNextIent(&im, &ent, &bc);
  while ( ent.cr != -1 ) {
    trow = ent.col;
    tcol = ent.row;
    for( i = bend; i >= 0; i--) {
      if( b[i].row > trow ) {
        b[i + 1] = b[i];
      } else if ( b[i].col > tcol ) {
        b[i + 1] = b[i];
      } else {
        b[i+1].col = tcol;
        b[i+1].row = trow;
        b[i+1].body = ent.body;
        break;
      }
    }
    if ( i == -1 ) {
      b[0].col = tcol;
      b[0].row = trow;
      b[0].body = ent.body;
    }
    bend++;
    GetNextIent(&im, &ent, &bc);
  }

  im = (pointer)m->root;
  ac = -1;
  GetNextIent(&im, &ent, &ac);
  aend = -1;
  while( ent.cr != -1 ){
    a[++aend].col = ent.col;
    a[aend].row = ent.row;
    a[aend].body = ent.body;
    GetNextIent(&im, &ent, &ac);
  }
  /* make return matrix */
  NEWIMAT(l);
  l->row = m->row;
  l->col = n->col;
  ac = -1;
  for( i = 0; i<= aend;) {
    ai = a[i].row;
    bj = b[0].row;
    s = 0;
    ik=i;
    for(j=0; j<= bend;){
      if( a[ik].col == b[j].col) {
        arf_mul(CO,(Obj)a[ik].body, (Obj)b[j].body, (Obj *)&u);
        arf_add(CO,(Obj)s,(Obj)u,(Obj *)&v);
        s = v;
        j++;
        if ( bj != b[j].row ) {
          AppendIent(l, ai, bj, (Obj)s);
          ik = i;
          bj = b[j].row;
          s = 0;
        } else {
          ik++;
          if ( ai != a[ik].row ) {
            AppendIent(l, ai, bj, (Obj)s);
            s = 0;
            while ( bj == b[j].row ) {
              j++;
              if ( j > bend ) break;
            }
            ik = i;
            bj = b[j].row;
          }
        }
      } else if ( a[ik].col > b[j].col ) {
        j++;
        if ( bj != b[j].row ) {
          AppendIent(l, ai, bj, (Obj)s);
          s = 0;
          ik = i;
          bj = b[j].row;
        }
      } else /* if ( a[ik].col < b[j].col ) */ {
        ik++;
        if ( ai != a[ik].row ) {
          AppendIent(l, ai, bj, (Obj)s);
          s = 0;
          while ( bj == b[j].row ) {
            j++;
            if ( j > bend ) break;
          }
          bj = b[j].row;
          ik = i;
        }
      }
    }
    i = ik;
    while ( ai == a[i].row ) {
      i++;
      if( i > aend) break;
    }
  }
  *r = (IMAT)l;
  FREE(a);
  FREE(b);
  return;
}