# gaussian elimination
# revised 1/23/08 bic

# test matrices
mat1 <- rbind(c(9,3,4,7),
             c(4,3,4,8),
             c(1,1,1,3))

mat2 <- rbind(c(1,-1,2,-1,-8),
             c(2,-2,3,-3,-20),
             c(1,1,1,0,-2),
             c(1,-1,4,3,4))

mat3 <- rbind(c(1,1,1,8),
              c(1,0,1,7),
              c(0,0,0,1))

mat4 <- rbind(c(0.003000,59.14,59.17),
              c(5.291,-6.130,46.78))

mat5 <- rbind(c(.00001,1,1),
              c(1,1,2))

mat6 <- rbind(
              c(1,1,2),
              c(.00001,1,1))

mat7 <- rbind(c(.0001,1,1),
              c(1,1,2))

mat9 <- rbind(c(10^-20,1,1), 
              c(1,1,2))

mat10 <- rbind(c(.3000*10^-3,.1566*10^1,.1569*10^1),
               c(.3454*10^0,-.2436*10^1,.1018*10^1))

mat10a <- rbind(c(.3454*10^0,-.2436*10^1,.1018*10^1),
                c(.3000*10^-3,.1566*10^1,.1569*10^1))


# code
swapRows <- function(mat, to, from)
  {
    new = mat;
    new[to,] = mat[from,];
    new[from,] = mat[to,];
    return(new);
  }

addRows <- function(mat, to, from, scale=1)   
  {
    mat[to,] = scale*mat[from,]+mat[to,];
    
    return(mat)
  }


backsub <- function(mat)
  {
    n = NROW(mat);
    x = array(0, dim=n);
    x[n] = mat[n,n+1]/mat[n,n];

    for(i in (n-1):1) {
      s = 0;
      for(j in (i+1):n) {
        s = s + mat[i,j]*x[j];
      }
      x[i] = (mat[i,n+1] - s)/mat[i,i];
    }

    return(x);
  }
        

rowReduce <- function(mat, debug=F)
  {
    n = NROW(mat);
    for (i in 1:(n-1)) {
      p = i;
      while (p < n+1) {
        if (mat[p, i] != 0) {
          break;
        }
        else {
          p = p+1;
        }
      }
      if (p > n) {
        return(-1);  # can't find a row to work
      }
      if (p != i) {
        if (debug) {
          cat(sprintf("swapRows(mat,i=%d p=%d)\n", i, p));
        }
        mat = swapRows(mat,i,p);
      }
      for (j in (i+1):n) {
        scale = -1*mat[j,i]/mat[i,i];
        mat = addRows(mat,j,i,scale);
        if (debug) {
          cat(sprintf(
             "addRows(mat,to=%d from=%d, scale=%f)\n",
                               j,i, scale));
        }
      }
    }
    return(mat);
  }
        
        
ge <- function(mat, debug=F) {
  m = rowReduce(mat, debug);
  if (length(attributes(m)) == 0) {
  # check for a -1 by looking at attribute list.
    #  m == -1 doesn't work
    cat("row reduce doesn't work\n");
  }
  else backsub(m);
}
     

        
sos <- function(){
  source("gauss-elim-start-2008.r")
}