<!--Define JavaScript functions. -->

var N = 5; //Global variable, dimension of matrix.

function sasum(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(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 SimEq5Solve(dataForm){

var A = new Array(N);
for (var i = 0; i < N; i++)
 A[i] = new Array(N);

A[0][0] = parseFloat(dataForm.a11.value);
A[0][1] = parseFloat(dataForm.a12.value);
A[0][2] = parseFloat(dataForm.a13.value);
A[0][3] = parseFloat(dataForm.a14.value);
A[0][4] = parseFloat(dataForm.a15.value);

A[1][0] = parseFloat(dataForm.a21.value);
A[1][1] = parseFloat(dataForm.a22.value);
A[1][2] = parseFloat(dataForm.a23.value);
A[1][3] = parseFloat(dataForm.a24.value);
A[1][4] = parseFloat(dataForm.a25.value);

A[2][0] = parseFloat(dataForm.a31.value);
A[2][1] = parseFloat(dataForm.a32.value);
A[2][2] = parseFloat(dataForm.a33.value);
A[2][3] = parseFloat(dataForm.a34.value);
A[2][4] = parseFloat(dataForm.a35.value);

A[3][0] = parseFloat(dataForm.a41.value);
A[3][1] = parseFloat(dataForm.a42.value);
A[3][2] = parseFloat(dataForm.a43.value);
A[3][3] = parseFloat(dataForm.a44.value);
A[3][4] = parseFloat(dataForm.a45.value);

A[4][0] = parseFloat(dataForm.a51.value);
A[4][1] = parseFloat(dataForm.a52.value);
A[4][2] = parseFloat(dataForm.a53.value);
A[4][3] = parseFloat(dataForm.a54.value);
A[4][4] = parseFloat(dataForm.a55.value);

var b = new Array(N);

b[0] = parseFloat(dataForm.b1.value);
b[1] = parseFloat(dataForm.b2.value);
b[2] = parseFloat(dataForm.b3.value);
b[3] = parseFloat(dataForm.b4.value);
b[4] = parseFloat(dataForm.b5.value);

var anorm = 0.0, dummy;

for (var i = 0; i < N; 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){
 dataForm.rcond.value = 0;
 dataForm.erCode.value = -1;
 return;
}//End if (anorm == 0)

var iw = new Array(N), l;

for (var i = 0; i < (N - 1); i++){
 l = isamax(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.
*/

var s, wk, wkm, sm, t;
var w = new Array(N);

for (var i = 0; i < N; i++)
 w[i] = 0.0;

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(w);

for (var i = 0; i < N; i++)   //Rescale w
 w[i] *= s;

//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(w);

for (var i = 0; i < N; i++)
 w[i] *= s;

var ynorm = 1.0;

//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(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(w);

dummy = ynorm/anorm;
dataForm.rcond.value = 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;

dataForm.erCode.value = ((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

dataForm.x1.value = b[0];
dataForm.x2.value = b[1];
dataForm.x3.value = b[2];
dataForm.x4.value = b[3];
dataForm.x5.value = b[4];

return;
}  //End of SimEq5Solve

// end of JavaScript-->