Linear algebra libraries







Numerical Methods

Linear algebra libraries

It is convenient to bundle functions into libraries. The functions that were described in the previous sections have been put in a library called matrix.js. The functions are:

copy(A) - Returns a copy of matrix A.
det(A) - Returns the determinant of square matrix A.
dot(A,B) - Returns a matrix that corresponds to the matrix product of matrices A and B.
identity(N) - Returns an NxN identity matrix.
inverse(A) - Returns the inverse of a square matrix.
LU_matrix(A) - Returns the LU matrix of square matrix A. This function overwrites A.
function LUP_decomp(A) - LU decomposition with partial pivoting. Returns the LU matrix, the permutation matrix P, and sgn=(-1) times the number of row swaps. This function overwrites A.
LU_solve(LU,b) - LU is a square NxN LU-matrix, and b is NxM matrix of M column vectors. Returns x, a NxM matrix of M column vectors. Each column of x is a solution of Ax=b.
matrix2tex(A) - Formats a matrix as a tex bmatrix.
new_matrix(N,M) - Declares an NxM matrix.
permutation(pv) - Returns a permutation matrix. pv is a vector of integers. pv[i] specifies the column of the 1 in row i.
pivot(A) - exchanges rows of A to put elements with large absolute values along the diagonal. Returns the permutation matrix P and PA = the matrix product of P times A. This function overwrites A.
transpose(A) - Returns the transpose of A.
write_matrix(A) - Formats a matrix as text.

The contents of this library are shown in the text area below.

function det(A) { //returns the determinant of square matrix A
  LUP = LUP_decomp(A);
  D = 1;
  for (i=0; i<A.length; i++) {//product the diagonal elements of the LU matrix
    D = D*LUP.LU[i][i]; 
  }
  D = D*LUP.sgn;
  return D;
}

function inverse(A) { // returns the inverse of a square matrix
  LUP = LUP_decomp(A);
  Pb = dot(LUP.P,identity(A.length));
  Ainv = LU_solve(LUP.LU,Pb);
  return Ainv;
}

function LUP_decomp(A) { //LU decomposition with partial pivoting
  //returns the LU matrix, the permutation matrix of matrix A, and (-1) times the number of row swaps 
  //This function overwrites A
  N = A.length; 
  P = identity(N); 
  sgn = 1; //even number of row exchanges sign = 1, odd number of exchanges sgn=-1

for (i=0; i<N; i++) {

//start pivot section
    Umax = 0;
    for (r=i; r<N; r++) {
      Uii=A[r][i];
      q = 0;
      while (q<i) {
        Uii -= A[r][q]*A[q][r];
        q++;
      } 
      if (Math.abs(Uii)>Umax) {
        Umax = Math.abs(Uii);
        row = r;
      }
    }
    if (i!=row) {//swap rows
      sgn = -sgn;
      for (q=0; q<N; q++) {
        tmp = P[i][q];
        P[i][q]=P[row][q];
        P[row][q]=tmp;
        tmp = A[i][q];
        A[i][q]=A[row][q];
        A[row][q]=tmp;
      } 
    }
   //end pivot section

j = i;
    while (j<N) { //Determine U across row i
      q = 0;
      while (q<i) {
        A[i][j] -= A[i][q]*A[q][j];
        q++;
      } 
      j++;
    }
    j = i+1;
    while (j<N) { //Determine L down column i
      q = 0;
      while (q<i) {
        A[j][i] -= A[j][q]*A[q][i];
        q++;
      } 
      A[j][i] = A[j][i]/A[i][i];
      j++;
    } 
  }
  return {LU:A,P:P,sgn:sgn};
}

function LU_solve(LU,b) {
  //LU is a square NxN LU-matrix, the LU decomposition of matrix A
  //b is NxM matrix of M column vectors
  //LU_solve(LU,b) returns x, a NxM matrix of M column vectors
  //each column of x is a solution of Ax=b

//initialize x and y to the same size as b
  rowb = b.length
  colb = b[0].length;
  x = new_matrix(rowb,colb);
  y = new_matrix(rowb,colb);
  N = LU.length;
  for (q=0; q<colb; q++) {
    for (i=0; i<N; i++) { //forward substitution
      y[i][q]=b[i][q];
      j = 0;
      while (j<i) {
        y[i][q] -= LU[i][j]*y[j][q];
        j++;
      }
    }
    for (i=N-1; i>-1; i--) { //back substitution
      x[i][q]=y[i][q];
      j = i+1;
      while (j<N) {
        x[i][q] -= LU[i][j]*x[j][q];
        j++;
      }
     x[i][q]=x[i][q]/LU[i][i];
    }
  }
  return x;
}

function new_matrix(N,M) { // declares an NxM 
  var newA = new Array(N);
  for (i=0; i<N; i++) {
    newA[i] = new Array(M);
  }
  return newA;
}

function pivot(A) { //determines pivot matrix for A 
  N = A.length; 
  P = identity(N);

for (j = 0; j<N; j++) {
    row = j;
    Amax = 0;
    for (i = j; i <N; i++) {//look down column j for largest element
      absA = Math.abs(A[i][j]);
      if (absA>Amax) {
        Amax = absA;
        row = i;
      }
    }
    if (j!=row) {//swap rows
      for (q=0; q<N; q++) {
        tmp = P[j][q];
        P[j][q]=P[row][q];
        P[row][q]=tmp;
        tmp = A[j][q];
        A[j][q]=A[row][q];
        A[row][q]=tmp;
      }
    }
  }
  return {P:P,PA:A};//returns permutation matrix P and PA
}

function transpose(A) {
  N = A.length;
  M = A[0].length;
  T = new_matrix(M,N);
  for (i=0; i<N; i++) { 
    for (j=0; j<M; j++) {
      T[j][i] = A[i][j];
    }
  }
  return T;
}

//***** display functions

The library can be imported with <script> tags. It is possible to import it as a local file or over the internet by specifying the URL of the file. Examples of both methods of importing the library are given below.

<p>A permutation matrix has the property that one element is equal to 1 in every row and every column. All other elements are zero. A permuation matrix is declared with the function permutation([3,1,4,2]). This function looks at the number $N$ of elements in the vector that it has been passed and creates a $N\times N$ matrix with a 1 in the third column of the first row, a 1 in the first column of the second row, a 1 in the fouth column of the third row, and a 1 in the second column of the fourth row. </p>

P = permutation([3,1,4,2]);

str = "\\[ \\begin{equation}";
str += "P ="+matrix2tex(P);
str += "\\end{equation} \\]";
document.write(str);
</script>

<p>The function transpose(P) is called to calculate the transpose of $P$. </p>

<p>The function inverse(P) is called to calculate the inverse of $P$. </p>

<p>Permutation matrices have the property that inverse equals the transpose.</p>

<p>Multiplying a matrix by a permutation matrix from the left reorders the rows of the matrix.</p>

<p>Multiplying a matrix by a permutation matrix from the right reorders the columns of the matrix.</p>

C/C++ libraries

All major computer languages have scientific libraries. Below the Linear Algebra section of the GNU Scientific library manual is embedded.