<!--Define JavaScript functions. -->

var N = 3; //Global variable, dimension of matrix.

function fObj(A, cvec, xvec){ // Calculates value of objective function, F
 // Since the POSITION of the maximum does not depend upon the scalar, alpha, do not bother
 // passing it in to the function or including it in calculations. The MAGNITUDE of the
 // maximum depends upon alpha, so alpha will be added at very end of program

var dummy = 0.0;            // A dummy variable
var dumvec = new Array(N);  // A dummy vector

for (var i = 0; i < N; i++) {
 dumvec[i] = 0.0;
 for (var j = 0; j < N; j++)
  dumvec[i] += A[i][j]*xvec[j];
} // End for i

for (var i = 0; i < N; i++)
 dummy += xvec[i]*dumvec[i];

dummy/= 2;

for (var i = 0; i < N; i++)
 dummy += cvec[i]*xvec[i];

return dummy;
} // End fObj

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 SimEq3Solve(A, b, x, codVec){

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){
 codVec[0] = 0; // rcond
 codVec[1] = -1; //erCode
 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)). = codVec[0]
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;
codVec[0] = 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;

codVec[1] = ((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

x[0] = b[0];
x[1] = b[1];
x[2] = b[2];

return;
}  //End of SimEq3Solve

function qp3Solve(dataForm){

 var H = new Array(N);       // Define the Hessian matrix, H
 var Hwork = new Array(N);   // Create a working copy of the H matrix

 for (var i = 0; i < N; i++) {
  H[i] = new Array(N);
  Hwork[i] = new Array(N);
 }

 H[0][0] = parseFloat(dataForm.h11.value);

 H[1][0] = parseFloat(dataForm.h21.value);
 H[1][1] = parseFloat(dataForm.h22.value);

 H[2][0] = parseFloat(dataForm.h31.value);
 H[2][1] = parseFloat(dataForm.h32.value);
 H[2][2] = parseFloat(dataForm.h33.value);

 H[0][1] = H[1][0];
 H[0][2] = H[2][0];

 H[1][2] = H[2][1];

 for (var i = 0; i < N; i++)  // Copy H into Hwork
  for (var j = 0; j < N; j++)
   Hwork[i][j] = H[i][j];

 var c = new Array(N);
 var cwork = new Array(N);   // Create a working copy of the c vector
 var x = new Array(N);

 c[0] = parseFloat(dataForm.c1.value);
 c[1] = parseFloat(dataForm.c2.value);
 c[2] = parseFloat(dataForm.c3.value);

 for (var i = 0; i < N; i++)  // Copy c vector into cwork vector
  cwork[i] = c[i];

 var alpha;
 var pVec = new Array(2);  // Parameter vector
 alpha = parseFloat(dataForm.alpha.value);

 if (dataForm.constrFlag[0].checked){ // Constraints present
  alert("This utility presently does not accept constraints. This capability will be added in the future. No further action          taken.");
 } // End if constraint Flag true
 else {                           // No constraints; simply solve [H](x) = -(c)
  alert("No constraints present.");

  for (var i = 0; i < N; i++)
   cwork[i] = -cwork[i];

  SimEq3Solve(Hwork, cwork, x, pVec);

  dataForm.x1.value = x[0];
  dataForm.x2.value = x[1];
  dataForm.x3.value = x[2];

  dataForm.fObjval.value = alpha + fObj(H, c, x);

  dataForm.rcond.value = pVec[0];
  dataForm.erCode.value = pVec[1];

 } // End else no constraints

 return;
}  //End of qp3Solve

// end of JavaScript-->