fcn make_array(n,m,v){ (m).pump(List.createLong(m).write,v)*n }
fcn eye(n){ // Creates a nxn identity matrix.
   I:=make_array(n,n,0.0);
   foreach j in (n){ I[j][j]=1.0 }
   I
}

// Creates the pivoting matrix for A.
fcn pivotize(A){
   n:=A.len();	// rows
   P:=eye(n);
   foreach i in (n){
      max,row:=A[i][i],i;
      foreach j in ([i..n-1]){
         if(A[j][i]>max) max,row=A[j][i],j;
      }
      if(i!=row) P.swap(i,row);
   }
   // Return P.
   P
}

// Decomposes a square matrix A by PA=LU and returns L, U and P.
fcn lu(A){
   n:=A.len();
   L:=eye(n);
   U:=make_array(n,n,0.0);
   P:=pivotize(A);
   A=matMult(P,A);

   foreach j in (n){
      foreach i in (j+1){
         U[i][j]=A[i][j] - (i).reduce('wrap(s,k){ s + U[k][j]*L[i][k] },0.0);
      }
      foreach i in ([j..n-1]){
         L[i][j]=( A[i][j] -
		   (j).reduce('wrap(s,k){ s + U[k][j]*L[i][k] },0.0) ) /
		 U[j][j];
      }
   }
   // Return L, U and P.
   return(L,U,P);
}

fcn matMult(a,b){
   n,m,p:=a[0].len(),a.len(),b[0].len();
   ans:=make_array(n,m,0.0);
   foreach i,j,k in (m,p,n){ ans[i][j]+=a[i][k]*b[k][j]; }
   ans
}