diff -u soy.i.orig soy.i
--- soy.i.orig
+++ soy.i
@@ -269,10 +269,13 @@ func ruoadd(a,b,ur=,un=)
   /* DOCUMENT func ruoadd(a,b)
      Addition of two square symmetric RUO matrices "a" and "b".
      SEE ALSO: ruoadd_float.c, ruoadd_double.c (soy.c)
+     Modified by Marcos van Dam, April 2011.
+     Now adds diagonal matrices
   */
 {
   if (typeof(*a.xn) != typeof(*b.xn)) error,"Mixed Data Types";
   if (a.r != b.r) error,"Matrices have incompatible dimensions!";
+
   argc = 5;
   if (ur == []) ur = MR;
   if (un == []) un = MN;
@@ -284,10 +287,14 @@ func ruoadd(a,b,ur=,un=)
     c.ix = &(array(long,ur));
     c.jx = &(array(long,un));
     c.xn = &(array(float,un));
-    c.xd = &(array(float,ur));
-    tt = array(float,a.n+b.n);
-    t = [&a,&b,&c,&tt,&s];
-    tmp = ruoadd_float(argc, &t);
+    if (a.n+b.n == 0){ // fix bug when adding two diagonal matrices
+      c.xd = &(*a.xd + *b.xd);
+    } else {
+      c.xd = &(array(float,ur));
+      tt = array(float,a.n+b.n);
+      t = [&a,&b,&c,&tt,&s];
+      tmp = ruoadd_float(argc, &t);
+    }
     return c;}
   else if (typeof(*a.xn) == "double") {
     c = ruo_d();
@@ -296,10 +303,14 @@ func ruoadd(a,b,ur=,un=)
     c.ix = &(array(long,ur));
     c.jx = &(array(long,un));
     c.xn = &(array(double,un));
-    c.xd = &(array(double,ur));
-    tt = array(double,a.n+b.n);
-    t = [&a,&b,&c,&tt,&s];
-    tmp = ruoadd_double(argc, &t);
+    if (a.n+b.n == 0){ // fix bug when adding two diagonal matrices
+      c.xd = &(*a.xd + *b.xd);
+    } else {
+      c.xd = &(array(double,ur));
+      tt = array(double,a.n+b.n);
+      t = [&a,&b,&c,&tt,&s];
+      tmp = ruoadd_double(argc, &t);
+    }
     return c;}
   else error,"Unsupported Data Type";
 }
@@ -433,7 +444,7 @@ func rcotr(a)
   at.t = a.t;
   at.ix = &(array(long,ur));
   at.jx = &(array(long,un));
-  if (a.n > 0){    
+  if (a.n > 0){
     sjx = long(sort((*a.jx)(1:a.n)));
     hjx = (*a.jx)(sjx);
     ax = array(long,a.c);
@@ -523,7 +534,7 @@ func ruoinf(a)
 {
   if (typeof(*a.xn) == "float") x = array(float,a.r,a.r);
   else if (typeof(*a.xn) == "double") x = array(double,a.r,a.r);
-  else error,"Unsupported Data Type";  
+  else error,"Unsupported Data Type";
   for (i=1; i<=a.r; i++) x(i,i) = (*a.xd)(i);
   for (i=1; i<a.r; i++) {
     if ((*a.ix)(i+1) > (*a.ix)(i)) {
@@ -763,6 +774,8 @@ func save_ruo(a,fn)
      Saves an RUO structure a to the binary file fn by converting
      all of its elements to float (double) and putting them into a
      single vector.
+     Modified by Marcos van Dam, April 2011.
+     Now saves diagonal matrices
   */
 {
   r = a.r;
@@ -770,16 +783,20 @@ func save_ruo(a,fn)
   if (typeof(*a.xn) == "float") {
     v = array(float,n*2+r*2+4);
     v(1:2) = float([n,r]);
-    v(3:n+2) = (*a.xn)(1:n);
-    v(n+3:2*n+2) = float((*a.jx)(1:n));
+    if (n > 0){ // fixed bug if n == 0
+      v(3:n+2) = (*a.xn)(1:n);
+      v(n+3:2*n+2) = float((*a.jx)(1:n));
+    }
     v(2*n+3:2*n+r+4) = float((*a.ix)(1:r+2));
     v(2*n+r+5:2*n+2*r+4) = (*a.xd)(1:r);
   }
   else if (typeof(*a.xn) == "double") {
     v = array(double,n*2+r*2+4);
     v(1:2) = double([n,r]);
-    v(3:n+2) = (*a.xn)(1:n);
-    v(n+3:2*n+2) = double((*a.jx)(1:n));
+    if (n > 0){ // fixed bug if n == 0
+      v(3:n+2) = (*a.xn)(1:n);
+      v(n+3:2*n+2) = double((*a.jx)(1:n));
+    }
     v(2*n+3:2*n+r+4) = double((*a.ix)(1:r+2));
     v(2*n+r+5:2*n+2*r+4) = (*a.xd)(1:r);
   }
@@ -812,11 +829,11 @@ func restore_rco(fn, ur=, un=)
   xn = v(4:a.n+3);
   jx = long(v(a.n+4:2*a.n+3));
   ix = long(v(2*a.n+4:2*a.n+a.r+5));
-  
+
   (*a.xn)(1:a.n) = xn;
   (*a.jx)(1:numberof(jx)) = jx;
   (*a.ix)(1:numberof(ix)) = ix;
-  
+
   return a;
 }
 
@@ -826,6 +843,8 @@ func restore_ruo(fn, ur=, un=)
      Returns the RUO structure saved in the file fn by save_rco.
      Modified by Marcos van Dam, August 2010.
      Now pads out the pointers, so that the matrices can be manipulated
+     Modified by Marcos van Dam, April 2011.
+     Now restores diagonal matrices
   */
 {
   if (ur == []) ur = MR;
@@ -842,13 +861,14 @@ func restore_ruo(fn, ur=, un=)
   a.xn = &(array(float,un));
   a.xd = &(array(float,ur));
 
-  xn = v(3:a.n+2);
-  jx = long(v(a.n+3:2*a.n+2));
+  if (a.n > 0){ // fixed bug if n == 0
+    xn = v(3:a.n+2);
+    jx = long(v(a.n+3:2*a.n+2));
+    (*a.xn)(1:numberof(xn)) = xn;
+    (*a.jx)(1:numberof(jx)) = jx;
+  }
   ix = long(v(2*a.n+3:2*a.n+a.r+4));
   xd = v(2*a.n+a.r+5:2*a.n+2*a.r+4);
-
-  (*a.xn)(1:numberof(xn)) = xn;
-  (*a.jx)(1:numberof(jx)) = jx;
   (*a.ix)(1:numberof(ix)) = ix;
   (*a.xd)(1:numberof(xd)) = xd;
 
@@ -873,7 +893,7 @@ func rcodr(&a,r)
     if (r == a.r) {
       (*a.jx)(a.n-nel+1:a.n) *= 0;
       (*a.xn)(a.n-nel+1:a.n) *= 0.0f;
-      (*a.ix)(a.r+1) = 0; 
+      (*a.ix)(a.r+1) = 0;
     } else if (r == 1) {
       (*a.jx)(1:a.n-nel) = (*a.jx)(nel+1:a.n);
       (*a.xn)(1:a.n-nel) = (*a.xn)(nel+1:a.n);
@@ -885,7 +905,7 @@ func rcodr(&a,r)
       //(*a.xn)((*a.ix)(r):a.n-nel-1) = (*a.xn)((*a.ix)(r+1):a.n-1); //orig
       (*a.jx)((*a.ix)(r)+1:a.n) = (*a.jx)((*a.ix)(r)+1+nel:a.n+nel); //rev
       (*a.xn)((*a.ix)(r)+1:a.n) = (*a.xn)((*a.ix)(r)+1+nel:a.n+nel); //rev
-      
+
       (*a.ix)(r+1:a.r) = (*a.ix)(r+2:a.r+1)-nel;
       (*a.ix)(a.r+1) = 0;
     }
@@ -981,7 +1001,7 @@ func ruo2rco(a)
   u.xn = &(dxn);
   u.n += u.r;
   b = rcoadd(u,l);
-  
+
   return b;
 }
 //==================================================================
@@ -1101,14 +1121,14 @@ func spunit(n, precision=)
     spidentity = spruo(double(unit(1)));
   }
   else {
-    spidentity = spruo(float(unit(1)));    
+    spidentity = spruo(float(unit(1)));
   }
 
   spidentity.r = n;
   (*spidentity.xd)(1:n) = 1;
-  
+
   return spidentity;
-}              
+}
 
 
 //==================================================================
@@ -1125,7 +1145,7 @@ func spzeros(m, n, precision=)
   else {
     zero_row = sprco(array(float,[2,n,1]));
   }
-  
+
   for (row_counter=1;row_counter<=m;row_counter++){
     if (row_counter == 1){
       zero_matrix = zero_row;
@@ -1135,4 +1155,4 @@ func spzeros(m, n, precision=)
     }
   }
   return zero_matrix;
-}  
+}