Ticket #33904: patch-soy.i.diff

File patch-soy.i.diff, 6.0 KB (added by paumard, 12 years ago)

  • soy.i 

    diff -u soy.i.orig soy.i
    old new func ruoadd(a,b,ur=,un=) 
    269269  /* DOCUMENT func ruoadd(a,b)
    270270     Addition of two square symmetric RUO matrices "a" and "b".
    271271     SEE ALSO: ruoadd_float.c, ruoadd_double.c (soy.c)
     272     Modified by Marcos van Dam, April 2011.
     273     Now adds diagonal matrices
    272274  */
    273275{
    274276  if (typeof(*a.xn) != typeof(*b.xn)) error,"Mixed Data Types";
    275277  if (a.r != b.r) error,"Matrices have incompatible dimensions!";
     278
    276279  argc = 5;
    277280  if (ur == []) ur = MR;
    278281  if (un == []) un = MN;
    … … func ruoadd(a,b,ur=,un=) 
    284287    c.ix = &(array(long,ur));
    285288    c.jx = &(array(long,un));
    286289    c.xn = &(array(float,un));
    287     c.xd = &(array(float,ur));
    288     tt = array(float,a.n+b.n);
    289     t = [&a,&b,&c,&tt,&s];
    290     tmp = ruoadd_float(argc, &t);
     290    if (a.n+b.n == 0){ // fix bug when adding two diagonal matrices
     291      c.xd = &(*a.xd + *b.xd);
     292    } else {
     293      c.xd = &(array(float,ur));
     294      tt = array(float,a.n+b.n);
     295      t = [&a,&b,&c,&tt,&s];
     296      tmp = ruoadd_float(argc, &t);
     297    }
    291298    return c;}
    292299  else if (typeof(*a.xn) == "double") {
    293300    c = ruo_d();
    … … func ruoadd(a,b,ur=,un=) 
    296303    c.ix = &(array(long,ur));
    297304    c.jx = &(array(long,un));
    298305    c.xn = &(array(double,un));
    299     c.xd = &(array(double,ur));
    300     tt = array(double,a.n+b.n);
    301     t = [&a,&b,&c,&tt,&s];
    302     tmp = ruoadd_double(argc, &t);
     306    if (a.n+b.n == 0){ // fix bug when adding two diagonal matrices
     307      c.xd = &(*a.xd + *b.xd);
     308    } else {
     309      c.xd = &(array(double,ur));
     310      tt = array(double,a.n+b.n);
     311      t = [&a,&b,&c,&tt,&s];
     312      tmp = ruoadd_double(argc, &t);
     313    }
    303314    return c;}
    304315  else error,"Unsupported Data Type";
    305316}
    … … func rcotr(a) 
    433444  at.t = a.t;
    434445  at.ix = &(array(long,ur));
    435446  at.jx = &(array(long,un));
    436   if (a.n > 0){   
     447  if (a.n > 0){
    437448    sjx = long(sort((*a.jx)(1:a.n)));
    438449    hjx = (*a.jx)(sjx);
    439450    ax = array(long,a.c);
    … … func ruoinf(a) 
    523534{
    524535  if (typeof(*a.xn) == "float") x = array(float,a.r,a.r);
    525536  else if (typeof(*a.xn) == "double") x = array(double,a.r,a.r);
    526   else error,"Unsupported Data Type"; 
     537  else error,"Unsupported Data Type";
    527538  for (i=1; i<=a.r; i++) x(i,i) = (*a.xd)(i);
    528539  for (i=1; i<a.r; i++) {
    529540    if ((*a.ix)(i+1) > (*a.ix)(i)) {
    … … func save_ruo(a,fn) 
    763774     Saves an RUO structure a to the binary file fn by converting
    764775     all of its elements to float (double) and putting them into a
    765776     single vector.
     777     Modified by Marcos van Dam, April 2011.
     778     Now saves diagonal matrices
    766779  */
    767780{
    768781  r = a.r;
    … … func save_ruo(a,fn) 
    770783  if (typeof(*a.xn) == "float") {
    771784    v = array(float,n*2+r*2+4);
    772785    v(1:2) = float([n,r]);
    773     v(3:n+2) = (*a.xn)(1:n);
    774     v(n+3:2*n+2) = float((*a.jx)(1:n));
     786    if (n > 0){ // fixed bug if n == 0
     787      v(3:n+2) = (*a.xn)(1:n);
     788      v(n+3:2*n+2) = float((*a.jx)(1:n));
     789    }
    775790    v(2*n+3:2*n+r+4) = float((*a.ix)(1:r+2));
    776791    v(2*n+r+5:2*n+2*r+4) = (*a.xd)(1:r);
    777792  }
    778793  else if (typeof(*a.xn) == "double") {
    779794    v = array(double,n*2+r*2+4);
    780795    v(1:2) = double([n,r]);
    781     v(3:n+2) = (*a.xn)(1:n);
    782     v(n+3:2*n+2) = double((*a.jx)(1:n));
     796    if (n > 0){ // fixed bug if n == 0
     797      v(3:n+2) = (*a.xn)(1:n);
     798      v(n+3:2*n+2) = double((*a.jx)(1:n));
     799    }
    783800    v(2*n+3:2*n+r+4) = double((*a.ix)(1:r+2));
    784801    v(2*n+r+5:2*n+2*r+4) = (*a.xd)(1:r);
    785802  }
    … … func restore_rco(fn, ur=, un=) 
    812829  xn = v(4:a.n+3);
    813830  jx = long(v(a.n+4:2*a.n+3));
    814831  ix = long(v(2*a.n+4:2*a.n+a.r+5));
    815  
     832
    816833  (*a.xn)(1:a.n) = xn;
    817834  (*a.jx)(1:numberof(jx)) = jx;
    818835  (*a.ix)(1:numberof(ix)) = ix;
    819  
     836
    820837  return a;
    821838}
    822839
    … … func restore_ruo(fn, ur=, un=) 
    826843     Returns the RUO structure saved in the file fn by save_rco.
    827844     Modified by Marcos van Dam, August 2010.
    828845     Now pads out the pointers, so that the matrices can be manipulated
     846     Modified by Marcos van Dam, April 2011.
     847     Now restores diagonal matrices
    829848  */
    830849{
    831850  if (ur == []) ur = MR;
    … … func restore_ruo(fn, ur=, un=) 
    842861  a.xn = &(array(float,un));
    843862  a.xd = &(array(float,ur));
    844863
    845   xn = v(3:a.n+2);
    846   jx = long(v(a.n+3:2*a.n+2));
     864  if (a.n > 0){ // fixed bug if n == 0
     865    xn = v(3:a.n+2);
     866    jx = long(v(a.n+3:2*a.n+2));
     867    (*a.xn)(1:numberof(xn)) = xn;
     868    (*a.jx)(1:numberof(jx)) = jx;
     869  }
    847870  ix = long(v(2*a.n+3:2*a.n+a.r+4));
    848871  xd = v(2*a.n+a.r+5:2*a.n+2*a.r+4);
    849 
    850   (*a.xn)(1:numberof(xn)) = xn;
    851   (*a.jx)(1:numberof(jx)) = jx;
    852872  (*a.ix)(1:numberof(ix)) = ix;
    853873  (*a.xd)(1:numberof(xd)) = xd;
    854874
    … … func rcodr(&a,r) 
    873893    if (r == a.r) {
    874894      (*a.jx)(a.n-nel+1:a.n) *= 0;
    875895      (*a.xn)(a.n-nel+1:a.n) *= 0.0f;
    876       (*a.ix)(a.r+1) = 0; 
     896      (*a.ix)(a.r+1) = 0;
    877897    } else if (r == 1) {
    878898      (*a.jx)(1:a.n-nel) = (*a.jx)(nel+1:a.n);
    879899      (*a.xn)(1:a.n-nel) = (*a.xn)(nel+1:a.n);
    … … func rcodr(&a,r) 
    885905      //(*a.xn)((*a.ix)(r):a.n-nel-1) = (*a.xn)((*a.ix)(r+1):a.n-1); //orig
    886906      (*a.jx)((*a.ix)(r)+1:a.n) = (*a.jx)((*a.ix)(r)+1+nel:a.n+nel); //rev
    887907      (*a.xn)((*a.ix)(r)+1:a.n) = (*a.xn)((*a.ix)(r)+1+nel:a.n+nel); //rev
    888      
     908
    889909      (*a.ix)(r+1:a.r) = (*a.ix)(r+2:a.r+1)-nel;
    890910      (*a.ix)(a.r+1) = 0;
    891911    }
    … … func ruo2rco(a) 
    9811001  u.xn = &(dxn);
    9821002  u.n += u.r;
    9831003  b = rcoadd(u,l);
    984  
     1004
    9851005  return b;
    9861006}
    9871007//==================================================================
    … … func spunit(n, precision=) 
    11011121    spidentity = spruo(double(unit(1)));
    11021122  }
    11031123  else {
    1104     spidentity = spruo(float(unit(1)));   
     1124    spidentity = spruo(float(unit(1)));
    11051125  }
    11061126
    11071127  spidentity.r = n;
    11081128  (*spidentity.xd)(1:n) = 1;
    1109  
     1129
    11101130  return spidentity;
    1111 }             
     1131}
    11121132
    11131133
    11141134//==================================================================
    … … func spzeros(m, n, precision=) 
    11251145  else {
    11261146    zero_row = sprco(array(float,[2,n,1]));
    11271147  }
    1128  
     1148
    11291149  for (row_counter=1;row_counter<=m;row_counter++){
    11301150    if (row_counter == 1){
    11311151      zero_matrix = zero_row;
    … … func spzeros(m, n, precision=) 
    11351155    }
    11361156  }
    11371157  return zero_matrix;
    1138 } 
     1158}