<!--Define JavaScript functions. -->

function sasum(N, vec){//Returns the norm of the input vector
 var dummy = 0.0;
 for (var i = 0; i < N; i++)  dummy += Math.abs(vec[i]);
 return dummy;
}//End sasum

function isamax(N, Matrix, col){
// Considers elements of column, col, of the input matrix, Matrix, on or below the diagonal.
// Returns the row number in which the largest element found.

var dummy = col, xmag, smax = Math.abs(Matrix[col][col]);

for (var i = col + 1; i < N; i++){
 xmag = Math.abs(Matrix[i][col]);
 if (xmag > smax){
  dummy = i;
  smax = xmag;
 } //End if (xmag > smax)
}//End for i
return dummy;
}//End isamax

function saxbr(N, A, b, iw, oPar){

var anorm = 0.0, dummy = 0.0, ynorm = 1.0, s, wk, wkm, sm, t;
var w = new Array(N);
var l = 1;

for (var i = 0; i < N; i++){
 w[i] = dummy = 0.0;
 for (var j = 0; j < N; j++)  dummy += Math.abs(A[j][i]);
 anorm = ((dummy > anorm) ? dummy : anorm);
}//End for i

if (anorm == 0.0){
 oPar.rcond = 0;
 oPar.outEr = -1;
 return;
}//End if (anorm == 0)

for (var i = 0; i < (N - 1); i++){
 l = isamax(N, A, i);
 // l becomes pivot element. A zero pivot implies this column is already triangularized.
 iw[i] = l;

 if (A[l][i] != 0){
  if (l != i){        //Interchange if necessary
   dummy = A[l][i];
   A[l][i] = A[i][i];
   A[i][i] = dummy;
  }//End if (l != i)
  dummy = -1.0/A[i][i];   //Compute multipliers

  for (var j = (i + 1); j < N; j++)  A[j][i] *= dummy;

  for (var j = (i + 1); j < N; j++){     //Row elimination with column indexing
   dummy = A[l][j];
   if (l != i){
    A[l][j] = A[i][j];
    A[i][j] = dummy;
   }//End if (l != i)
   for (var k = (i + 1); k < N; k++)  A[k][j] += dummy*A[k][i];
  }///End for j.
 }//End if A[l][i] != 0)
}//End for i

iw[N - 1] = N - 1;

// rcond = 1/Norm(A)*estimate of Norm(Inverse(A)).
// Estimate = Norm(Z)/Norm(Y) where A*Z = Y and Trans(A)*Y = E.
// Trans(A) is the transpose of A. The components of E are chosen to cause minimum local growth in the elements of W, where Trans(U)*W = E. The vectors are frequently rescaled to avoid overflow.

dummy = 1.0;
for (var i = 0; i < N; i++){
 if (w[i] != 0.0)  dummy = ((w[i] > 0.0) ? -dummy : dummy);
 if (Math.abs(-w[i] + dummy) > Math.abs(A[i][i])){
  s = Math.abs(A[i][i])/Math.abs(-w[i] + dummy);
  for (var j = 0; j < N; j++)  w[j] *= s;
  dummy *= s;
 }//End if abs
 wk = -w[i] + dummy;
 wkm = -(dummy + w[i]);
 s = Math.abs(wk);
 sm = Math.abs(wkm);
 if (A[i][i] != 0.0){
  wk /= A[i][i];
  wkm /= A[i][i];
 }//End if (A[i][i] != 0)
 else
  wkm = wk = 1.0;
 if (i < (N - 1)){
  for (var j = (i + 1); j < N; j++){
   sm += Math.abs(w[j] + wkm*A[i][j]);
   w[j] += wk*A[i][j];
   s += Math.abs(w[j]);
  }//End for j
  if (s < sm){
   t = -wk + wkm;
   wk = wkm;
   for (var j = (i + 1); j < N; j++)  w[j] += t*A[i][j];
  }//End if (s < sm)
 }//End if (i < (N - 1))
 w[i] = wk;
}//End for i

s = 1.0/sasum(N, w);

for (var i = 0; i < N; i++) w[i] *= s;   //Rescale w

//Solve Trans(L)*Y = W

for (var i = N - 1; i >= 0; i--){
 if (i < (N - 1)){
  dummy = 0.0;
  for (var k = (i + 1); k < N; k++)  dummy += A[k][i]*w[k];
  w[i] += dummy;
 }//End if (i < (N - 1))

 if (Math.abs(w[i]) > 1.0){
  s = 1.0/Math.abs(w[i]);
  for (var k = 0; k < N; k++)  w[k] *= s;
 }//End if abs(w[i] > 1)
 l = iw[i];
 t = w[l];
 w[l] = w[i];
 w[i] = t;
}//End for i

s = 1.0/sasum(N, w);

for (var i = 0; i < N; i++) w[i] *= s;

//Solve L*V = Y

for (var i = 0; i < N; i++){
 l = iw[i];
 t = w[l];
 w[l] = w[i];
 w[i] = t;
 for (var j = (i + 1); j < N; j++)  w[j] += t*A[j][i];
 if (Math.abs(w[i]) > 1.0){
  s = 1.0/Math.abs(w[i]);
  for (var j = 0; j < N; j++)  w[j] *=s;
  ynorm *= s;
 }//End if (abs(w[i]) > 1)
}//End for i

s = 1.0/sasum(N, w);

for (var i = 0; i < N; i++) w[i] *= s;

ynorm *= s;

//Solve U*Z = V

for (var i = N - 1; i >= 0; i--){
 if (Math.abs(w[i]) > Math.abs(A[i][i])){
  s = Math.abs(A[i][i])/Math.abs(w[i]);
  for (var k = 0; k < N; k++)  w[k] *= s;
  ynorm *= s;
 }//End if abs(w[i] > abs(A[i][i])
 w[i] = ((A[i][i] != 0.0) ? w[i]/A[i][i] : 1.0);
 t = -w[i];
 for (var k = 0; k < i; k++)  w[k] += t*A[k][i];
}//End for i

ynorm /= sasum(N, w);
dummy = ynorm/anorm;
oPar.rcond = dummy;

// Now find machine epsilon by brute force; the algorithm should be general enough
// to port to other machines, operating systems, etc.
// epsmch is the smallest positive number such that 1.0 + epsmch is not equal to 1.0
// In the MSDN Library, DBL_EPSILON is 2.2204460492503131e-16
// accessed by using #include <float.h>

sm = 1.0;
do {
 ynorm = sm;
 sm /= 2.0;
 s = 1.0 + sm;
} while (s > 1.0); // End do-while loop
//At this point, ynorm should equal the machine precision, epsmch

dummy = ynorm/dummy;
dummy = -Math.log(dummy)/Math.LN10;

oPar.outEr = ((dummy > 0) ? Math.floor(dummy) : -2);

//Solve Ax = b by solving L*Y = b

for (var i = 0; i < (N - 1); i++){
 l = iw[i];
 t = b[l];
 if (l != i){
  b[l] = b[i];
  b[i] = t;
 }//End if (l != i) 
 for (var j = (i + 1); j < N; j++)  b[j] += t*A[j][i];
}//End for i

for (var j = (N - 1); j >= 0; j--){
 b[j] /= A[j][j];
 t = -b[j];
 for (var k = 0; k < j; k++)  b[k] += t*A[k][j];
}//End for j

return;
}  //End of saxbr

// end of JavaScript-->