| List of works by author |
Notes:
1. To download matrix routines:- From SWiSHzone.com
- As a swi file
- As a zipped set of files
- Click here or scroll down to the end of this page for a listing of the core of the programs
- copy and paste them directly into your program,
- use them as "include" files,
- or save them in "component" directories for handy use as components.
include "fileName.txt"; to have all the routines inserted into your program before it is compiled.
4. For more detailed directions on their use, read the comment lines in the script of these routines.
5. These routines have been checked for matrices as large as 60x60, requiring a few seconds to run. They should work for larger matrices, given a faster computer and/or more time.
Matrix manipulation package contents:
1. Testing program: lines of code that call on the routines. These mostly output the results via trace statements. They can be helpful as examples for use of the routines.2. Eigenvalue
- Householder routine. This converts a real square symmetric matrix into a tridiagonal matrix, having non-zero elements only along the diagonal and on either side of the diagonal. This prepares a symmetric matrix for the next routine (the QL routine) that extracts the eigenvalues and eigenvectors from the matrix.
- QL routine. This routine extracts the eigenvalues and eigenvectors from a tridiagonal real square symmetric matrix.
- Inverse iteration routine. This uses an eigenvalue to solve for the corresponding eigenvector. It works by repeatedly using the LU decomposition and backsolve routines (also in this package).
- Limited Householder routine. Used with the limited QL routine to solve for the eigenvalues but not the eigenvectors of a real square symmetric tridiagonal matrix.
- Limited QL routine. Used with the limited Householder routine to solve for the eigenvalues but not the eigenvectors of a real square symmetric tridiagonal matrix.
- Eigenvalues routine. Combines the limited Householder routine with the limited QL routine to solve for the eigenvalues (but not eigenvectors) of a real square symmetric tridiagonal matrix.
- Eigenvalue/vector routine. Combines the Householder routine with the QL routine to solve for the eigenvalues and eigenvectors of a real square symmetric tridiagonal matrix.
- LU factoring (decomposition) routine. Factors a real matrix into the product of a lower (having zeros in all elements above the diagonal) and an upper matrix (having zeros in all elements below the diagonal). This is used with the next routine to solve matrix equations. This routine produces a compacted LU matrix, which is a composite of both L and U matrices.
- LU back solve routine. This is used with the LU factoring routine above to solve matrix equations and to invert matrices.
- Separate triangular matrices. This separates the L and U matrices out of the compacted LU form.
- Inverse matrix routine. Uses the above routines to invert a matrix.
- Matrix solve routine. Uses the above routines to solve a matrix equation. In most applications it is faster than the inverse matrix routine.
- Determinant routine. Uses the LU factoring routine to calculate the determinant of a matrix.
- LU factoring (decomposition) routine. Factors a real matrix into the product of a lower (having zeros in all elements above the diagonal) and an upper matrix (having zeros in all elements below the diagonal). This is used with the next routine to solve matrix equations. This routine produces a compacted LU matrix, which is a composite of both L and U matrices.
- LU back solve routine. This is used with the LU factoring routine above to solve matrix equations and to invert matrices.
- Separate triangular matrices. This separates the L and U matrices out of the compacted LU form.
- Inverse matrix routine. Uses the above routines to invert a matrix.
- Matrix solve routine. Uses the above routines to solve a matrix equation. In most applications it is faster than the inverse matrix routine.
- Determinant routine. Uses the LU factoring routine to calculate the determinant of a matrix.
- Ordered matrix produces a matrix with sequential numbers as elements.
- Test matrix produces a 5x5 matrix for testing purposes.
- Householder test matrix produces a given matrix from a Wikipedia example.
- Symmetric test matrix converts the above test matrix into a symmetric matrix.
- Symmetric matrix uses a matrix to produce a symmetric matrix.
- Constant matrix produces a matrix having all its elements equal to a given constant.
- Random matrix produces a matrix having randomly generated elements.
- Pseudo random matrix uses sinusoidal functions to produce a matrix whose elements vary wildly from each other. This function will always produce the same seemingly random matrix. Used for testing procedures.
- Random vector produces a vector (1D array) having randomly generated elements.
- Constant vector produces a vector having all its elements equal to a given constant.
- Ordered vector produces a vector with sequential numbers as elements.
- Print matrix causes a matrix to be printed with trace statements, in a user settable format.
- Print vector causes a vector to be printed with trace statements, in a user settable format.
- Round number rounds a number off to a settable decimal point.
- String number produces a string representing a given number in a user settable format.
- Matrix multiply multiplies two matrices together.
- Matrix vector multiply multiplies a vector by a matrix.
- Identity matrix produces an identity matrix of a specified size.
- Matrix addition adds two matrices together.
- Matrix subtraction subtracts one matrix from another.
- Vector times constant multiplies a vector by a scalar.
- Matrix times constant multiplies a matrix by a scalar.
- Vector magnitude squared produces the square of the magnitude of a vector . You can take the square root of this to get the magnitude.
- Vector dot product creates the dot (or scalar) product of two vectors.
- Vector outer product creates the outer (or matrix) product of two vectors.
- Transpose creates the transpose of a matrix.
- Vector sort re-arranges the elements of a vector (or array) from largest to smallest.
- Reverse vector sort, re-arranges the elements of a vector (or array) from smallest to largest.
- Matrix sort interchanges rows of a matrix so that the specified column has its elements arranged from largest to smallest.
- Reverse matrix sort interchanges rows of a matrix so that the specified column has its elements arranged from smallest to largest.
- Screen print matrix causes a matrix to be displayed on the screen of a swi (Swishmax) or swf (flash) routine in a specified format.
- Screen print horizontal vector causes a vector to be displayed horizontally on the screen of a swi (Swishmax) or swf (flash) routine in a specified format.
- Screen print vertical vector causes a vector to be displayed vertically on the screen of a swi (Swishmax) or swf (flash) routine in a specified format.
- Duplicate matrix duplicates a matrix. Normally setting two matrices (A and B for example) equal (as A=B) merely assigns both variable names to the same matrix. Changing elements in one will change the elements in the other. This routine makes an independent new matrix that initially has the same values as the starting matrix.
- String to 2D array converts a string of values to a 2D array or matrix. The user needs to specify the number of rows and columns of the new matrix. The element separators are assumed to be spaces.
- String to vector converts a string of values to a 1D array or vector. The user needs to specify the number of elements in the new vector. The element separators are assumed to be spaces.
References:
1. A good teaching of the Householder method. http://en.wikipedia.org/wiki/Householder_transformation .2. A very nice website describing the details of the Householder and QL methods http://beige.ucs.indiana.edu/B673/node30.html .
3. Another similar descriptive site: http://math.fullerton.edu/mathews/n2003/HouseholderMod.html .
4. A standard respected reference on most numerical methods including the Householder and QL routine: Numerical Recipes, Press, Flannery, Teukolsky, Vetterling. Web site for this book, http://www.nr.com/1. http://en.wikipedia.org/wiki/LU_decomposition .
5. A good teaching on LU decomposition method http://en.wikipedia.org/wiki/LU_decomposition .
6. A very nice website describing the details of the LU decomposition and similar methods albeit with the examples in C++: http://mymathlib.webtrellis.net/matrices/linearsystems/doolittle.html .
7. A site with all sorts of numerical routines (including the LU decomposition) written for flash actionscript. These routines would need to be changed somewhat to work for Swishmax: http://members.shaw.ca/flashprogramming/wisASLibrary/wis/index.html
8. See my web site on waves and resonances at: http://resonanceswavesandfields.blogspot.com/ .
9. For more on matrix equations and eigenvalues and vectors see: http://en.wikipedia.org/wiki/System_of_linear_equations and http://en.wikipedia.org/wiki/Eigenvalue,_eigenvector_and_eigenspace .
Script:
1. Eigenvalue
/* Eigenvalue/eigenvector package
The following routines are variations of the householder and QL routines,
used to calculate eigenvalues and eigenvectors.
Also included are routines that combine these two to calculate
eigenvalues and eigenvectors more easily.
There is also a function that uses the LU decomposition to
calculate an eigenvector for a specific eigenvalue. This may be
faster if one or only a few eigenvectors of a large array are needed.
References:
1. http://en.wikipedia.org/wiki/Householder_transformation
2. http://beige.ucs.indiana.edu/B673/node30.html --- a very nice
& nbsp; website desribing the details of these methods.
3. Another similar descriptive site:
http://math.fullerton.edu/mathews/n2003/HouseholderMod.html
4. A standard respected reference on most numerical methods
including the Householder and QL routine:
Numerical Recipes, Press, Flannery, Teukolsky, Vetterling
Web site for the above book, http://www.nr.com/
See my web page at http://resonanceswavesandfields.blogspot.com/
*/
function householder(A:Array,EE:Array):Array{
/* This function transforms a symmetric input square matrix A into a
tridiagonal form, suitable for the QL routine to then find
the eigenvalues and eigenvectors for A.
The householder function returns D, the diagonal elements, in
a 1D array.
The off diagonal elements are available as the EE in the argument list.
The 1D EE array needs to be defined in the calling program
before the householder function is called, i.e. var EE=new Array(n); .
In Flash, arrays in the argument list are not protected and
we use this feature to return the EE vector (1D array).
The orthogonal matrix Q is available as the changed input array A.
*/
var n:Number=A.length;
var m:Number=A[0].length;
if(n!=m){trace("This routine only works for square matrices.");return A;}
var l:Number,i:Number,j:Number,k:Number;
var h:Number,scale:Number,f:Number,g:Number,hh:Number,r:Number;
var D=new Array(n);
var sigma:Number,Ail:Number,sqrtTerm:Number,H:Number,sum:Number,K:Number;
var uj:Number,qj:Number;
for(i=n-1;i>0;i--){// step through the columns from right to left (reverse normal order).
l=i-1;// stays as i-1, i.e. as the row number of the off diagonal elements.
h=0;
H=0;
sigma=0;
scale=0;// is the epsilon,
if(l>0){// don't need to do the following for the last (left most) column.
// FIRST: calculate epsilon ... called "scale" in algorithm
for(k=0;k<=l;k++){
scale+=Math.abs(A[i][k]); }// end of k loop
// SECONDLY: if scale is not zero, then scale all the A's with epsilon
// and calculate sigma
if(scale==0){EE[i]=A[i][l]; // skip transformation.
} else{
for(k=0;k<=l;k++){
A[i][k]/=scale; // scales all elements, use scaled A's in transformation.
sigma+=A[i][k]*A[i][k]; // calculate sigma
}//end of k loop
// THIRDLY: use sigma to calculate u AND also H
Ail=A[i][l];// this is the off diagonal element. We need to
// it use a lot here so we make a special variable.
sqrtTerm=-Math.sqrt(sigma)*Math.sign(Ail);// We calculate the square root
// term with the proper sign
// This also almost equals |u|
EE[i]=scale*sqrtTerm;// This is the final result for the single off diagonal
// element ( the EE[i] value ). We calculate
// for this i value.
// Puts other temporary values
// in the yet unused EE terms.
// These EE elements are the off diagonal
// elements of our to-be-created tridiagonal
// matrix. They are calculated
// and are represented as K = uT p/2H.
// This routine also stores the p vectors in the yet
// unused EE elements where p=Au/H.
// These p vectors are calculated in part IV below.
H=sigma-Ail*sqrtTerm;
A[i][l]=Ail-sqrtTerm; // this fixes up the last non-zero u element (other non-zero terms
// are Aik...ok as they are).
// We store u in the i th row of A.
// FOURTH: we work on calculating the vector p = Au/H and K .;
K=0; // K is the off diagonal element ... to be eventually stored in the EE array.
for(j=0;j<=l;j++){
A[j][i]=A[i][j]/H;// store u/h in the i th column of A.
sum=0; // form an element of A.u in g.
for(k=0;k<=j;k++){sum+=A[j][k]*A[i][k];}// this is made more complicated by
// our storage of u in redundant parts of the
// A matrix. The top part of the A
// matrix is not overlaid with u values
// and is used to retrieve real A values.
for(k=j+1;k<=l;k++){sum+=A[k][j]*A[i][k];}
EE[j]=sum/H; // form element of p in temporarily unused element
K+=EE[j]*A[i][j];
}// end of j loop
K/=(H+H);// form k.
// FIFTH: calculate the q matrix and then the A' value
for(j=0;j<=l;j++){
uj=A[i][j]; // these are u[j] values
qj=EE[j]-K*uj;// this is q[j]
EE[j]=qj; // this stores q[j] in the unused EE
//temporarily, over top the p's there.
for(k=0;k<=j;k++){
A[j][k]-=uj*EE[k]+qj*A[i][k];}//<16>
// The above line calculates the new A's, i.e. A'.
}// end of j loop
}// end of if(scale==0)/ }else
}else{EE[i]=A[i][l]} ;// Done if l==0. This statement works
// with the if(l>0){ ... near beginning of routine
// the final } in this line ends this chain.
D[i]=H; // final diagonal values.
}// end of i loop
// SIXTH: do some final things
D[0]=0;
EE[0]=0;// These are the off diagonal elements. Since there
// is one less of these than n, we just zero out
// the unused element (in a vector on length n).
for(i=0;i<n;i++){ // begin accumulation of transformation matrix
// Q, overwriting the remains of the highly
// transformed A.
l=i-1;
if(D[i]!=0){ // this block skipped when i=0
for(j=0;j<=l;j++){// Next we calculate Q
sum=0;
for(k=0;k<=l;k++){sum+=A[i][k]*A[k][j];}// use u and u/h stored in A to form P.Q
for(k=0;k<=l;k++){A[k][j]-=sum*A[k][i];}//
}// end j loop
}// end if
D[i]=A[i][i] ;// Transfer the diagonal elements of Q
// to D. These must be the final diagonal
// elements of our tridiagonal matrix.
A[i][i]=1; //set the diagonal elements of A (now Q) to identity matrix.
for(j=0;j<=l;j++){A[i][j]=0;A[j][i]=0;}
// The above line zeroes out the lower triangle of Q
// to get it in final form.
}// end i loop
return D;
}// end of function householder
// ============================================================
// ============================================================
/* The next routine, the QLroutine, calculates the eigenvalues and eigenvectors
of a symmetric tridiagonal matrix. The above householder routine transforms a symmetric
matrix into a tridiagonal form, suitable for the QLroutine to work on.
The inputs are:
The D array inputs the diagonal elements of a tridiagonal matrix.
The EE array inputs the sub-diagonal elements of a tridaigonal matrix.
The Z matrix is the transform that the householder routine used to change
the matrix into tridiagonal form.
The outputs are:
This routine returns a matrix of the eigenvectors, wherein each column of the
matrix is an eigenvector. These eigenvector are arranged in the same order as
the returned eigenvalues.
The eigenvalues are returned as the transformed input D array (a vector).
In Flash, arrays in the argument list are not protected and
we use this feature to return the eigenvalue list as a transformed D vector (1D array).
*/
function QLroutine(D:Array,EE:Array,Z:Array):Array{
var n:Number=D.length;
var i:Number,l:Number,m:Number,k:Number,iter:Number;
var g:Number,s:Number,c:Number,p:Number,b:Number,r:Number,dd:Number;
var maxIter:Number=30;
for(i=1;i<n;i++){EE[i-1]=EE[i];}//<11>convenient to renumber the elements of EE
// i.e. this shifts the elements down so that
// they range from 0 to n-2 instead of from 1 to n-1.
EE[n-1]=0;
for(l=0;l<n;l++){ // The l (lower case ell) index steps through making smaller and smaller
// matrices to work on. Each iteration of this loop results
// in the top non-zero off diagonal element (the l th element in EE)
// being zeroed. The next iteration will work on the remaining
// non-diagonal matrix (which is one order smaller).
for(iter=0;iter<maxIter;iter++){
for(m=l;m<(n-1);m++){//(12)This small loop looks for a single small off diagonal
// element, at which point the matrix will be split.
// Really, this iteration of the "iter" loop will only
// work on the matrix from l to m where m is the location
// of the first off diagonal zero (past l).
dd=Math.abs(D[m])+Math.abs(D[m+1]);
if((Math.abs(EE[m])+dd)==dd){break} //alter the condition if this is too
//extreme. It will determine when the
//routine considers an off diagonal
//to be effectively zero.
}// end of m loop
if(m==l){break;} // this is true when we have successfully zeroed out the top
// off diagonal element (the l th one). Then it is time to go to
// the very end of the program and start a new l value and work
// on the next smaller matrix.
g=(D[l+1]-D[l])/(2*EE[l]); // form shift
r=Math.sqrt(g*g+1);
g=D[m]-D[l]+EE[l]/(g+Math.abs(r)*Math.sign(g)); // calculates d[m] - k[s].
s=1;
c=1;
p=0;
for(i=m-1;i>=l;i--){// a plane rotation as in the original QL,followed by
// Givens rotations to restore tridiagonal form.
f=s*EE[i];
b=c*EE[i];
if(Math.abs(f)>=Math.abs(g)){
c=g/f;
r=Math.sqrt(c*c+1);
EE[i+1]=f*r;
s=1/r;
c*=s;
}else{ // belongs to if 6 lines above
s=f/g;
r=Math.sqrt(s*s+1);
EE[i+1]=g*r;
c=1/r;
s*=c;
}// end of if/else
g=D[i+1]-p;
r=(D[i]-g)*s+2*c*b;
p=s*r;
D[i+1]=g+p;
g=c*r-b;
for(k=0;k<n;k++){
f=Z[k][i+1];
Z[k][i+1]=s*Z[k][i]+c*f;
Z[k][i]=c*Z[k][i]-s*f;
}// end of k loop
}// end of i loop
D[l]-=p;
EE[l]=g;
EE[m]=0;
}// end of iter loop
// We jump this if we "break" in the line above (close
// to the top) on the condition that m==l.
}// end of l loop
return Z;
}// end of function QLroutine
//==================Beginning of inverseIteration routine
/* This routine uses the LU decomposition routine to iteratively calculate
an eigenvector for a given eigenvalue. It may be faster to use this
routine than the complete QL routine above if only one or a few
eigenvectors are needed. (Use the limited householder and limited
QL routines below to only calculate the eigenvalues when using
the inverse iteration routine.)
*/
function inverseIteration(M:Array,eigenValue:Number,iterations):Array{
// this solves for eigenvectors.
var n:Number=M.length;
var eVector=new Array(n);
var newVector=new Array(n);
var i:Number,j:Number,iter:Number,y:Number,norm:Number,dif:Number;
var dif:Number,magDif:Number;
var I = new Array(n);
var DiffMat = new Array(n);
for(i=0;i<n;i++){I[i]=new Array(n);DiffMat[i]=new Array(n);}
I=_parent.matrixC.identityMatrix(n);
eVector=_parent.matrixC.randomVector(n);// initial guess
norm=Math.sqrt(_parent.matrixC.vectorMagnitudeSquared(eVector));
eVector=_parent.matrixC.vectorTimesConstant(1/norm,eVector);
DiffMat=_parent.matrixC.matrixSubtraction(M,_parent.matrixC.matrixTimesConstant(eigenValue,I));
DiffMat=_parent.LUdecomposition.LUfactoring(DiffMat);// DiffMat is now in the compact LU form
// and we need to use _global.indx which is also created by "LUfactoring"
for(iter=0;iter<iterations;iter++){
newVector=_parent.LUdecomposition.LUbackSolve(DiffMat,_global.indx,eVector);
norm=Math.sqrt(_parent.matrixC.vectorMagnitudeSquared(newVector));
newVector=_parent.matrixC.vectorTimesConstant(1/norm,newVector);
if(newVector[0]<0){newVector=_parent.matrixC.vectorTimesConstant(-1,newVector);}
magDif=0;
for(i=0;i<n;i++){
dif=newVector[i]-eVector[i];
magDif+=dif*dif;
eVector[i]=newVector[i];}// end of i loop
if((magDif>0.3)and(iter>3)){iter=0;eVector=_parent.matrixC.randomVector(n);
trace("try new random starting eigenvector"); };
}// end of iter loop
_global.eigenChange=Math.sqrt(magDif);
if(eVector[0]<0){eVector=_parent.matrixC.vectorTimesConstant(-1,eVector);}
return eVector;
}// end of function inverse iteration
// ==================== Beginning of "limitedHouseholder" function ================
function limitedHouseholder(A:Array):Array{
// Only useful for finding eigenvalues, not eigenvectors...
// does not form transform array Q
// This returns D (the diagonal elements), the off diagonal element, EE.
// A dummy vector needs to be input as EE of length n.
// The orginal matrix A is destroyed.
// references: Numerical methods and
// http://en.wikipedia.org/wiki/Householder_transformation
var n:Number=A.length;
var m:Number=A[0].length;
if(n!=m){trace("This routine only works for square matrices.");return A;}
var l:Number,i:Number,j:Number,k:Number;
var h:Number,scale:Number,f:Number,g:Number,hh:Number,r:Number;
var EE=new Array(n);
var D=new Array(n);
var sigma:Number,Ail:Number,sqrtTerm:Number,H:Number,sum:Number,K:Number;
var uj:Number,qj:Number;
for(i=n-1;i>0;i--){// Step through the columns from right to left (reverse normal order).
l=i-1;// stays as i-1 as the row of the off diagonal elements.
h=0;
H=0;
sigma=0;
scale=0;// is the epsilon,
if(l>0){// don't need to do the following for the last (left most) column... remember
// we are working right to left.
// FIRST: calculate epsilon, called "scale" here.
for(k=0;k<=l;k++){
scale+=Math.abs(A[i][k]); }// end of k loop
// SECOND: if scale is not zero, then scale all the A's with epsilon
// and calculate sigma.
if(scale==0){EE[i]=A[i][l]; // skip transformation.
} else{
for(k=0;k<=l;k++){
A[i][k]/=scale; // Scales all elements. Use scaled A's in transformation.
sigma+=A[i][k]*A[i][k]; // Calculate sigma.
}//end of k loop
// THIRD: use sigma to calculate u .
Ail=A[i][l];// This is the off diagonal element. We need to
// it use a lot here so we make a special variable.
sqrtTerm=-Math.sqrt(sigma)*Math.sign(Ail);// We calculate the square root
// term in eq.11.2.14 with the proper sign.
// This also almost equals |u|.
EE[i]=scale*sqrtTerm; // This is the final result for the single off diagonal
// element ( the EE[i] value ). We calculate
// for this i value.
// Puts other temporary values
// in the yet unused EE terms.
// These EE elements are the off diagonal
// elements of our to-be-created tridiagonal matrix. They are calculated
// using Eq. 11.2.11 and are represented as K = uT p/2H
// in the algebra in the book, Numerical Recipes.
// The factor of two in the denominator may cancel the missing two above (for H),
// letting the program just ignore the 2 in both.
// This routine also stores the p vectors in the yet
// unused EE elements where p=Au/H as per 11.2.10 .
// These p vectors are calculated in part IV below.
H=sigma-Ail*sqrtTerm;
A[i][l]=Ail-sqrtTerm; // this fixes up the last non-zero u element (other non-zero terms
// are Aik...ok as they are).
// We store u in the i th row of A.
// FOURTH: we work on calculating the vector p as in Eq. 11.2.10 p = Au/H and K in 11.2.11 .
K=0; // K is the off diagonal element ... to be eventually stored in the EE array.
for(j=0;j<=l;j++){
sum=0; // form an element of A.u in g.
for(k=0;k<=j;k++){sum+=A[j][k]*A[i][k];}// This is made more complicated by
// our storage of u in redundant parts of the
// A matrix. The top part of the A
// matrix is not overlaid with u values
// and is used to retrieve real A values.
for(k=j+1;k<=l;k++){sum+=A[k][j]*A[i][k];}//
EE[j]=sum/H; // form element of p in temporarily unused element
K+=EE[j]*A[i][j];
}// end of j loop
K/=(H+H);// form k, eqn 11.2.11 . Note: this has the factor of 2 in the denominator.
// FIFTH: calculate the q matrix 11.2.12 and then the A' value 11.2.13 .
for(j=0;j<=l;j++){
uj=A[i][j]; // these are u[j] values
qj=EE[j]-K*uj;// this is q[j] as per Eq. 11.2.12 .
EE[j]=qj; // this stores q[j] in the unused EE
//temporarily, over top the p's there.
for(k=0;k<=j;k++){
A[j][k]-=uj*EE[k]+qj*A[i][k];}
// The above line calculates the new A's, i.e. A'
// using 11.2.13
}// end of j loop
}// end of if(scale==0)/ }else
}else{EE[i]=A[i][l]} ;// Done if l==0. This statement works
// with the if(l>0){ ... near beginning of routine
// the final } in this line ends this chain.
D[i]=H; // final diagonal values.
}// end of i loop
// SIXTH: do some final things
EE[0]=0;// These are the off diagonal elements. Since there
// is one less of these than n, we just zero out
// the unused element (in a vector of length n).
for(i=0;i<n;i++){// Transfer the diagonal elements of Q
D[i]=A[i][i]; // to D. These must be the final diagonal
// elements of our tridiagonal matrix.
}// end i loop
for(i=0;i<n;i++){A[i]=EE[i];}// Overwrite the EE's (off diagonal elements)
// onto the A array, which is now a vector instead
// of a 2D array.
return D;
}// end of function limited householder
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function limitedQLroutine(D:Array,EE:Array):Array{
// a QL routine without eigenvector capabilities
var n:Number=D.length;
var i:Number,l:Number,m:Number,k:Number,iter:Number;
var g:Number,s:Number,c:Number,p:Number,b:Number,r:Number,dd:Number;
for(i=1;i<n;i++){EE[i-1]=EE[i];}//convenient to renumber the elements of EE
// i.e. this shifts the elements down on so that
// they range from 0 to n-2 instead of from 1 to n-1.
EE[n-1]=0;
for(l=0;l<n;l++){ // The l index steps through making smaller and smaller
// matrices to work on. Each iteration of this loop results
// in the top non-zero off diagonal element (the l th element in EE)
// being zeroed. The next iteration will work on the remaining
// non-diagonal matrix (which is one order smaller).
for(iter=0;iter<6;iter++){
for(m=l;m<(n-1);m++){//This small loop looks for a single small off diagonal
// element (in the range l to n-1), at which point the
// matrix will be split.
// The rest of the iteration of the iter loop will only
// work on the matrix from l to m where m is the location
// of the first off diagonal zero (past l) and ignore
// the rest for now.
dd=Math.abs(D[m])+Math.abs(D[m+1]);
if((Math.abs(EE[m])+dd)==dd){break}// alter the condition if this is too
//extreme. It will determine when the
//routine considers an off diagonal
//to be effectively zero.
}// end of m loop
if(m==l){break} // this is true when we have successfully zeroed out the top
// off diagonal element (the l th one). Then it is time to go to
// the very end of the program and start a new l value and work
// on the next smaller matrix.
g=(D[l+1]-D[l])/(2*EE[l]); // form shift
r=Math.sqrt(g*g+1);
g=D[m]-D[l]+EE[l]/(g+Math.abs(r)*Math.sign(g)); // this is d[m] - k[s].
s=1;
c=1;
p=0;
for(i=m-1;i>=l;i--){// a plane rotation as in the original QL, followed by
// Givens rotations to restore tridiagonal form.
f=s*EE[i];
b=c*EE[i];
if(Math.abs(f)>=Math.abs(g)){
c=g/f;
r=Math.sqrt(c*c+1);
EE[i+1]=f*r;
s=1/r;
c*=s;
}else{ // belongs to the "if", 6 lines above.
s=f/g;
r=Math.sqrt(s*s+1);
EE[i+1]=g*r;
c=1/r;
s*=c;
}// end of if/else
g=D[i+1]-p;
r=(D[i]-g)*s+2*c*b;
p=s*r;
D[i+1]=g+p;
g=c*r-b;
}// end of i loop
D[l]-=p;
EE[l]=g;
EE[m]=0;
}//end of iter loop.
}// end of l loop
return D;
}// end of function QLroutine
// ============================================================
// The following routine uses the routines above to calculate the eigenvalues of
// matrix, M. It returns a vector of these eigenvalues. It does not protect
// the input matrix.
function eigenvalues(M:Array):Array{
var n:Number=M.length;
var vector=new Array(n);
vector=limitedHouseholder(M);//vector now becomes the list of diagonal elements
vector=limitedQLroutine(vector,M);// M is transformed into the off
// diagonal elements by the previous step.
return vector// this is the list of eigenvalues
}// end of function eigenValue
//=================================================
// The following routine uses the above routines to calculate both the
// eigenvalues and the eigenvectors for matrix M. It does not
// protect the input matrix M.
function eigenValueVectors(M:Array):Array{
var n:Number=M.length;
var arr=new Array(n);
var EE=new Array(n);
arr=householder(M,EE);// "arr" now becomes the list of diagonal elements.
// "EE" is initially undefined, and becomes the
// off diagonal elements as calculated by
// the householder routine.
// "M" is initially the matrix for which we want the
// eigenvalues and vectors. It is changed into the
// Q's, the 2D transform generated by the
// householder routine.
M=QLroutine(arr,EE,M);
return arr// This returns the eigenvalues (the vector "arr").
// The argument, M, is changed into a column-by-column array
// of eigenvectors.
// If only a few eigenvectors are needed, the function
// inverseIteration (above) may be more efficient to use.
} // end of eigenValueVectors routine
2. LUdecomposition
LU decomposition package/*
The following routines are used to solve matrix equations such as
M x = b where M is a known square matrix, x is an unknown vector
(one dimensional array), and b is a known vector. These methods
solve for x.
These methods have been successfully run in Swishmax3
for matrices up to 60x60 in size, in a few seconds. Larger matrices
should be able to be accomodated with faster computers and/or more
time.
References:
1. http://en.wikipedia.org/wiki/LU_decomposition
2. http://mymathlib.webtrellis.net/matrices/linearsystems/doolittle.html
--- a very nice website describing the details of these methods albeit
with the examples in C++ instead of swishmax script.
3. A site with all sorts of numerical routines (including the
LU decomposition) written for flash actionscript. They would need to
be changed somewhat to work for Swishmax:
http://members.shaw.ca/flashprogramming/wisASLibrary/wis/index.html
4. A standard respected reference on most numerical methods
including the LU decomposition:
Numerical Recipes, Press, Flannery, Teukolsky, Vetterling
Web site for the above book, http://www.nr.com/
See my web page at http://resonanceswavesandfields.blogspot.com/
*/
/* The next function factors a square matrix into two triangular matrices,
one a lower triangular matrix (the L matrix, having 1's along its
diagonal and zeros above the diagonal) and a upper triangular matrix (the U one
having zeros below the diagonal): M = L U .
This method is more commonly called LU decomposition.
The input argument is the matrix to be decomposed or factored. The output
is a compacted merging of the L and U matrices into a single matrix.
In this compacted form, the 1's along the diagonal of the L matrix
are deleted (replaced by zeros) and the two matrices are written on top
o f each other. The non-zero elements of the U matrix replace the
zeros of the L matrix. The row swapping done by LUfactoring is made
available as the vector _global.indx and the even/odd number of swaps
in the sign of _global.d .
This version of LUfactoring does not modify the input matrix.
This factoring algorithm is normally used in conjunction with
a back solving routine, called LUbackSolve here. After the factoring,
solving the equation M x = L U x = b is very straightforward (by
back solving).
*/
function LUfactoring(InputArray:Array):Array {
var n:Number=InputArray.length;
var m:Number=InputArray[0].length;
if(m != n){trace("Array is not square.");return InputArray};
var i:Number, j:Number, k:Number;
_global.indx=new Array(n);
_global.d=1;
var norm:Array=new Array(n);
var temp:Number, imax:Number, max:Number, sum:Number;
var tiny:Number=0.00001;
var A:Array=new Array(n);// the purpose of this first block is to make a working copy of the input matrix
// so that the original argument is not changed
for(i=0;i<n;i++){
A[i]=new Array(n);
for(j=0;j<n;j++){ A[i][j]=InputArray[i][j];}
}//end of i loop
for (i=0;i<n;i++){
max=0;
for (j=0;j<n;j++){
temp=Math.abs(A[i][j]);
if(temp>max){max=temp};
}// end of j loop
if(max==0) {trace("singular matrix"); return A;}
norm[i] = 1/max;
} // end of i loop
for (j = 0; j < n; j++) { // very large loop... involves rest of function
for (i = 0; i<j; i++) {
sum = A[i][j];
for (k=0; k<i; k++){
sum -= A[i][k]*A[k][j];
} // end of k loop
A[i][j] = sum
} // end of i loop
max = 0;
for (i=j; i<n; i++)
sum=A[i][j];
for (k=0; k<j; k++) {
sum -= A[i][k]*A[k][j];
} // end of k loop
A[i][j] = sum;
temp = norm[i]*Math.abs(sum);
if(temp>=max) {
imax = i;
max = temp }
} // end of i loop
if(j != imax){
for (k=0; k<n; k++){
temp = A[imax][k];
A[imax][k] = A[j][k];
A[j][k] = temp;
} // end k loop
_global.d = -_global.d;
norm[imax] = norm[j];
} //end of if
_global.indx[j] = imax;
if(A[j][j] == 0) {A[j][j] = tiny;}
if(j != n-1) {
temp = 1/A[j][j];
for(i = j+1; i<n; i++) {A[i][j] *= temp}
}// end of if
}// end of large j loop
return A;
} // end of function LUfactoring
/* ===================== Second function ... uses factored matrix to solve for x ==============================
To be used with function LUfactoring above (which factors a matrix into
lower, L, and upper, U, triangular matrices) to solve the the matrix equation
A x = L U x = L x' = b, where A and b are known and x'=U x. The matrix, A, needs to
be in the form that the output of the above function generates, as compacted lower/upper
triangular matrices. The result, the vector x, overwrites B in the output.
The b vector is shown as B in the program.
The intermediate result x' = _global.xPrime is available outside the routine
(delete two lines at the end of the program if you do not want this feature.)
*/
function LUbackSolve(A:Array, indx:Array, B:Array):Array {
var isZero:Boolean, i:Number, j:Number, k:Number, sum:Number;
var n:Number=A.length;
var xPrime=new Array(n);
var out=new Array(n);// this is x, the solution (a vector)
isZero=true;
for(i=0;i<n;i++){xPrime[i]=B[i];}// copies B to protect this input.
for (i=0; i<n; i++){// this loop solves L x' = b for x', where x' = Ux
k = indx[i];
sum = xPrime[k]; // the missing identity matrix times b, with permutations.
xPrime[k] = xPrime[i]; // rearranges the b's after using the original B[k] in the line above.
if(!isZero){
for(j=0; j<i; j++) {
sum -= A[i][j]*xPrime[j];}// steps through each element of the lower matrix
} else if(sum != 0) {isZero = false} // end if and else if
xPrime[i] = sum;
} // end of i loop
for(i=n-1; i>=0; i--) { // this loop solves U x = x' for x
sum = xPrime[i];
for(j=i+1; j<n; j++) {
sum -= A[i][j]*out[j];}// steps through each element of the upper matrix
out[i] = sum/A[i][i];
} // end of i loop
_global.xPrime=new Array(n);
f or(i=0;i<n;i++){_global.xPrime[i]=xPrime[i];}
return out
}// end of function LUbackSolve
/* The next routine separates out the L and U matrices from their
compacted form. Its input is the compacted output of
LUfactoring. It returns the L (lower) matrix and makes the U(upper)
matrix globally available as _global.U .
*/
function separateTriangleMatrices(A:Array):Array{
var n:Number=A.length;
var L:Array=new Array(n);
var U:Array=new Array(n);
var i:Number,j:Number;
for (i=0;i<n;i++){
L[i]=new Array(n);
U[i]=new Array(n);
for (j=0;j<n;j++){
if(j<i){
L[i][j]=A[i][j];
U[i][j]=0;
}else if (i==j){
L[i][j]=1;
U[i][j]=A[i][j];
}else {
L[i][j]=0;
U[i][j]=A[i][j];
}// end of if/elseif/else chain
}// end of j loop
}// end of i loop
_global.U=U; // _global.U must be externally defined
return L;
}// end of function separateTriangleMatrices
/* The next function returns the inverse matrix of a square
matrix. If the user only needs to solve one particular
equation, they should use matrixSolve instead, because it is much faster,
particularly for large matrices.
*/
function inverseMatrix(A:Array):Array{
var n:Number=A.length;
var m:Number=A[0].length;
var I=new Array(n);
var LU=new Array(n);
var InverseMat=new Array(n);
var partial=new Array(n);
if(m!=n){trace("This needs to be a square matrix to allow inversion.");}
for(var i:Number=0;i<n;i++){LU[i]=new Array(n);InverseMat[i]=new Array(n);}
LU=LUfactoring(A)
for(var j:Number=0;j<n;j++){
for(i=0;i<n;i++){I[i]=0;if(i==j){I[i]=1}}
partial=LUbackSolve(LU,_global.indx,I);
for(i=0;i<n;i++){InverseMat[i][j]=partial[i];}
}// end of j loop;
return InverseMat;
} // end of function inverseMatrix
/* The next routine solves the matrix equation A x = b.
Its input arguments are A and b. It returns the
vector x.
*/
function matrixSolve(A:Array,b:Array):Array{
var n:Number=A.length;
var m:Number=A[0].length;
var LU=new Array(n);
var out=new Array(n);
if(m!=n){trace("This needs to be a square matrix to allow inversion.");}
for(var i:Number=0;i<n;i++){LU[i]=new Array(n);}
LU=LUfactoring(A);
out=LUbackSolve(LU,_global.indx,b);
return out;
}// end of function matrix solve
/* The following routine uses the LU factoring routine
(LU decomposition) to calculate the determinant of
input square matrix A.
*/
function determinant(A:Array):Number{
var n:Number=A.length;
var m:Number=A[0].length;
var LU=new Array(n);
if(m!=n){trace("This needs to be a square matrix.")}
for(var i:Number=0;i<n;i++){LU[i]=new Array(n);}
LU=LUfactoring(A);
var product:Number=_global.d;
for(var i:Number=0;i<n;i++){product*=LU[i][i];}
return product;
}// end of function determinant
3. matrixC
General matrix manipulation package// See my web page at http://resonanceswavesandfields.blogspot.com/
function orderedMatrix(rows:Number,columns:Number):Array{
// creates and fills array with sequential numbers
var value:Number=0;
var element:Array = new Array(rows);
for (var i:Number=0;i<rows;i++){
element[i]=new Array(columns);
for (var j:Number=0;j<columns;j++){
element[i][j]=value;
value++}
}// end of loops
return element
}// end of function
function testMatrix():Array {
// creates a specific 5x5 matrix used to test some of the routines
var n:Number=5;
var A:Array=new Array(n);
A[0]=new Array(2.2,1,0,-3,1.5);
A[1]=new Array(-0.9,2.6,-0.1,-3,0.99);
A[2]=new Array(4,-1,0.33,-4,-1.55);
A[3]=new Array(1,3.7,-0.33,-4.3,1.4);
A[4]=new Array(0,5,-2,-3,9);
return A;
}// end of function
function householderTestMatrix():Array{
// an 4x4 example matrix used in
// http://en.wikipedia.org/wiki/Householder_transformation
var n:Number=4;
var A:Array=new Array(n);
A[0]=new Array(4,1,-2,2);
A[1]=new Array(1,2,0,1);
A[2]=new Array(-2,0,3,-2);
A[3]=new Array(2,1,-2,-1);
return A;
}// end of function householderTestMatrix
function symmetricTestMatrix():Array{
// converts the above test matrix into a symmetric matrix
// to allow testing some routines that work on symmetric matrices.
var out=new Array;
out=testMatrix();
out=symmetricMatrix(out);
return out;
}// end of function symmetric test matrix
function symmetricMatrix(A:Array):Array{
// will convert any matrix A into a symmetric matrix.
var n:Number=A.length;// number of rows
var m:Number=A[0].length;//number of columns
if(m!=n){trace("must be a square matrix to make symmetric");return A;}
var out=new Array(n);
for(var i:Number=0;i<n;i++){
out[i]=new Array(n);
for(var j:Number=0;j<n;j++){
out[i][j]=(A[i][j]+A[j][i])/2;
}// end of j loop
}// end of i loop
return out;
} // end of function symmetric metrix
function constantMatrix(rows:Number,columns:Number,value:Number):Array{
// creates a matrix in which all elements have the same value.
var element:Array = new Array(rows);
for (var i:Number=0;i<rows;i++){
element[i]=new Array(columns);
for (var j:Number=0;j<columns;j++){
element[i][j]=value; } }// end of loops
return element;
} // end of function
function randomMatrix(rows:Number,columns:Number):Array{
// creates and fills an array with random numbers
// in the range of -50 to +50 .
var element:Array = new Array(rows);
for (var i:Number=0;i<rows;i++){
element[i]=new Array(columns);
for (var j:Number=0;j<columns;j++){
element[i][j]=Math.randomInt(100)-50; } }// end of loops
return element;
}// end of function
function pseudoRandomMatrix(rows:Number,columns:Number):Array{
// creates and fills an array with varying numbers in the range
// of -2 to +2. This routine will always return the same
// matrix (if the number of rows and columns are the same)
// each time it is called (the reason for the "pseudo" title.)
var element:Array = new Array(rows);
var k1:Number,k2:Number;
for (var i:Number=0;i<rows;i++){
element[i]=new Array(columns);
}
k1=1;
k2=Math.sqrt(3)*k1;
element=orderedMatrix(rows,columns);
for (i=0;i<rows;i++){
for (var j:Number=0;j<columns;j++){
element[i][j]=Math.sin(k1*i*j)+Math.cos(k2*(j+i));}}
return element
}// end of function
function randomVector(rows:Number):Array{
// creates and fills a vector (a 1D array) with random numbers in the
// range of -50 to +50.
var element:Array = new Array(rows);
for (var j:Number=0;j<rows;j++){
element[j]=Math.randomInt(100)-50; } // end of loop
return element
}// end of function
function constantVector(rows:Number,value:Number):Array{
// creates and fills a vector (a one dimensional
// array) with a constant value.
var element:Array = new Array(rows);
for (var j:Number=0;j<rows;j++){
element[j]=value; } // end of loop
return element
}// end of function
function orderedVector(rows:Number,value:Number):Array{
// creates and fills a vector with an increasing set of numbers
var element:Array = new Array(rows);
var value:Number=5;
for (var j:Number=0;j<rows;j++){
element[j]=value;
value++ } // end of loop
return element
}// end of function
function printMatrix(Mat:Array,field:Number,decPosition:Number,separatorStg:String):Void{
// causes the trace function to print out a matrix in a structured format.
// "field" is the length in character spaces of each number or matrix element.
// "decPosition" is the number of digits to be displayed to the right of the decimal place.
// "separatorStg" is the string to be used between numbers. Spaces are a good separator.
// A typical use of this routine might be printMatrix(M,6,2," ");
// This would cause the Matrix M to be printed out with 6 character spaces reserved for
// each number, 2 digits to the right of the decimal, and two space between adjacent numbers.
if(Mat[0][0]==undefined){trace("toVectorPrint");printVector(Mat,field,decPosition,separatorStg);return}
var rows:Number, columns:Number;
rows=Mat.length;
columns=Mat[0].length;
var out:String;
for (var i:Number=0;i<rows;i++){
out="";
for (var j:Number=0;j<columns;j++){
out+=stringNumber(field,decPosition,Mat[i][j])+separatorStg; }// end of j loop
trace(out)}// end of i loop
}// end of function
function printVector(vec:Array,field:Number,decPosition:Number,separatorStg:String):Void{
// causes the trace function to print out a vector (1D array) in a structured format.
// See the instruction for print Matrix (above) for more details.
if(vec[0][0] != undefined){trace("toMatrixPrint");printMatrix(vec,field,decPosition,separatorStg);return}
var rows:Number;
rows=vec.length;
var out:String="";
for (var i:Number=0;i<rows;i++){
out+=stringNumber(field,decPosition,vec[i])+separatorStg; }// end of j loop
trace(out)
}// end of function
function roundNum(d:Number,value:Number):Number {
// returns a number rounded off to a specified number of digits.
// where n is an integer number of decimal digits
// n can be negative for rounding to the left of the decimal place
var factor:Number=Math.pow(10,d);
return Math.round(value*factor)/factor;
} // end of function
function stringNumber(field:Number,decPlace:Number,value:Number):String{
// converts a number to a formated string for use in several of the above output
// routines. For more specifics see the instructions for printMatrix above.
// This routine as it stands now could be expanded to gracefully handle scientific
// notation.
var strg:String;
var targetDecimalPosition:Number=field-decPlace-1;
strg=String(roundNum(decPlace,value));
// check length, if trying to clip off numbers or scientific notation
// (make a special function for this)...Do this later
// see where the decimal point is and add zeroes on the right as necessary
var decimalPosition:Number
decimalPosition=strg.indexOf(".");
if((decimalPosition==-1) && ((decPlace>0) && (strg.length<field))) {strg+=".";}
decimalPosition=strg.indexOf(".");
if(strg.length<field and decimalPosition<targetDecimalPosition){
for(var i=0;i<targetDecimalPosition-decimalPosition;i++){strg=" "+strg}}
if(strg.length<field){
for(var i=0;i<=field-strg.length;i++){strg+="0"}}
while(strg.length>field and strg.charAt(0)==" "){strg=strg.substring(1)}
return strg
}// end of function
function matrixMultiply(A:Array,B:Array):Array {
// returns the matrix that is the product of A times B.
var rowsA:Number, colnsA:Number, rowsB:Number, colnsB:Number;
rowsA=A.length;
colnsA=A[0].length;
rowsB=B.length;
colnsB=B[0].length;
if(colnsA != rowsB){ trace("ERROR: The number of columns of A does not equal the number of rows of B!");}
if((B[0][0]==undefined)and(colnsB=1)){return matrixVectorMultiply(A,B)}
var C:Array=new Array(rowsA);
for(var i:Number=0;i<rowsA;i++){ //for all rows of output matrix
C[i]=new Array(colnsB);
for(var j:Number=0;j<colnsB;j++) { // for all columns of output matrix
C[i][j]=0;
for(var k:Number=0;k<colnsA;k++){ // for all summed-over indices
C[i][j]+=A[i][k]*B[k][j];
}// end of k loop
}//end of j loop
}// end of i loop
return C;
}// end of function
function matrixVectorMultiply(A:Array,V:Array):Array {
// returns the products of matrix A times vector V.
var rowsA:Number, colnsA:Number, rowsV:Number;
rowsV=V.length;
var vector:Array=new Array(rowsV);
if(V[0][0]==undefined){
for(var j:Number=0;j<rowsV;j++){vector[j]=V[j]}
}else{for(j=0;j<rowsV;j++){vector[j]=V[j][0]}}//converts a 2D array of nx1 elements to a 1D array
rowsA=A.length;
colnsA=A[0].length;
rowsV=V.length;
if(colnsA != rowsV){ trace("ERROR: The number of columns of A does not equal the number of rows of V!");}
var W:Array=new Array(rowsA);
for(var i:Number=0;i<rowsA;i++){ //for all rows of output matrix
W[i]=0;
for(var j:Number=0;j<colnsA;j++) { // for all columns of output matrix
W[i]+=A[i][j]*vector[j];
}//end of j loop
}// end of i loop
return W;
}// end of function
function identityMatrix(order:Number):Array {
// returns the identity matrix of size "order" x "order".
// An identity matrix has all its elements zero, except those along the
// diagonal are ones. When a vector or matrix is multiplied by the
// identity matrix, the product is the same vector or matrix.
var I:Array=new Array(order);
for(var i:Number=0;i<order;i++){ //for all rows of output matrix
I[i]=new Array(order);
for(var j:Number=0;j<order;j++) { // for all columns of output matrix
if(i!=j)I[i][j]=0;
else I[i][j]=1
}//end of j loop
}// end of i loop
return I;
} // end of function
function matrixAddition(A:Array,B:Array):Array {
// returns the matrix sum of A + B . Both A and B need to be of the same dimensions.
var rowsA:Number, colnsA:Number, rowsB:Number, colnsB:Number;
rowsA=A.length;
colnsA=A[0].length;
rowsB=B.length;
colnsB=B[0].length;
if(colnsA != colnsB){ trace("ERROR: The number of columns of A does not equal the number of columns of B!");}
if(rowsA != rowsB){ trace("ERROR: The number of rows of A does not equal the number of rows of B!");}
trace("line 119 "+rowsA+" "+colnsA+" "+rowsB+" "+colnsB);
var C:Array=new Array(rowsA);
for(var i:Number=0;i<rowsA;i++){ //for all rows of output matrix
C[i]=new Array(colnsA);
for(var j:Number=0;j<colnsA;j++) { // for all columns of output matrix
C[i][j]=A[i][j]+B[i][j];
}//end of j loop
}// end of i loop
return C;
} // end of function
function matrixSubtraction(A:Array,B:Array):Array {
// returns the matrix difference of A - B . Both A and B need to be of the same dimensions.
var rowsA:Number, colnsA:Number, rowsB:Number, colnsB:Number;
rowsA=A.length;
colnsA=A[0].length;
rowsB=B.length;
colnsB=B[0].length;
if(colnsA != colnsB){ trace("ERROR: The number of columns of A does not equal the number of columns of B!");}
if(rowsA != rowsB){ trace("ERROR: The number of rows of A does not equal the number of rows of B!");}
trace("line 119 "+rowsA+" "+colnsA+" "+rowsB+" "+colnsB);
var C:Array=new Array(rowsA);
for(var i:Number=0;i<rowsA;i++){ //for all rows of output matrix
C[i]=new Array(colnsA);
for(var j:Number=0;j<colnsA;j++) { // for all columns of output matrix
C[i][j]=A[i][j]-B[i][j];
}//end of j loop
}// end of i loop
return C;
} // end of function
function vectorTimesConstant(s:Number,V:Array):Array{
// returns the vector formed by multiplying a vector V by the scalar s.
var rows:Number=V.length;
var C:Array=new Array(rows);
for(var j:Number=0;j<rows;j++) { // for all columns of output matrix
C[j]=s*V[j];
}//end of j loop
return C;
} // end of function
function matrixTimesConstant(s:Number,A:Array):Array{
// returns the matrix formed by multiplying a matrix A by the scalar s.
var rows:Number=A.length;
var colns:Number=A[0].length;
var C:Array=new Array(rows);
for(var i:Number=0;i<rows;i++){ //for all rows of output matrix
C[i]=new Array(colns);
for(var j:Number=0;j<colns;j++) { // for all columns of output matrix
C[i][j]=s*A[i][j];
}//end of j loop
}// end of i loop
return C;
} // end of function
function vectorMagnitudeSquared(V:Array):Number{
// returns the magnitude squared of a vector, i.e. the sum of all
// its elements squared.
var rows:Number=V.length;
var sum:Number=0;
for(var j:Number=0;j<rows;j++) { // for all columns of output matrix
sum+=V[j]*V[j];
}//end of j loop
return sum;
} // end of function
function vectorDotProduct(V:Array,W:Array):Number{
// returns the dot (inner or scalar) product of vector V times vector W.
// This product is a simple number.
var rowsV:Number=V.length;
var rowsW:Number=W.length;
if(rowsW != rowsV){trace("The number of elements in the two vectors are unequal");return -1;}
var sum:Number=0;
for(var j:Number=0;j<rowsV;j++) { // for all columns of output matrix
sum+=V[j]*W[j];
}//end of j loop
return sum;
} // end of function
function vectorOuterProduct(V:Array,W:Array):Array{
// returns the outer or matrix product of vector V and vector W. The result
// is a matrix.
var rowsV:Number=V.length;
var rowsW:Number=W.length;
if(rowsW != rowsV){trace("The number of elements in the two vectors are unequal");return V;}
var C:Array=new Array(rowsV);
for(var i:Number=0;i<rowsV;i++) {
C[i]=new Array(rowsV);
for(var j:Number=0;j<rowsV;j++) { // for all columns of output matrix
C[i][j]=V[i]*W[j];
}}//end of loops
return C;
} // end of function
function transpose(A:array):Array{
// returns the transpose of matrix A, formed by interchanging the rows and columns of A.
var n:Number=A.length;// number of rows of A
var m:Number=A[0].length;// number of columns of A
var out=new Array(n);
for(var i:Number=0;i<n;i++){
out[i]=new Array(m);
for(var j:Number=0;j<m;j++){
out[j][i]=A[i][j];
}// end of j loop
}// end of i loop
return out;
}// end of function transpose
function vectorSort(V:Array):Array{
// sorts vector V (a 1D array) to return a vector in which the the elements are
// arranged from largest to smallest.
var n:Number=V.length;
var temp:Number=0;
var AA=new Array(n);
for(var i:Number=0;i<n;i++){AA[i]=V[i];}
for(var iter:Number=0;iter<n;iter++){
for(i=0;i<n-1;i++){
if(AA[i+1]>AA[i]){
temp=AA[i];
AA[i]=AA[i+1];
AA[i+1]=temp;}// end if
}}// end of loops
return AA;
}// end of function sort
function reverseVectorSort(V:Array):Array{
// sorts vector A (a 1D array) to return a vector in which the the elements are
// arranged from smallest to largest.
var n:Number=V.length;
var out=new Array(n);
for(var i:Number=0;i<n;i++){out[i]=-V[i]};
out = vectorSort(out);
for(i=0;i<n;i++){out[i]=-out[i]};
return out;
}// end of function reverse sort
function matrixSort(M:Array,sortColumn:Number):Array{
// sorts matrix by interchanging rows so that the elements in the matrix
// column "sortColumn" are arranged from largest(at the top) to smallest.
var n:Number=M.length;
var m:Number=M[0].length;
var temp:Number=0;
var AA=new Array(n);
for(var i:Number=0;i<n;i++){AA[i]=new Array(m);
for(var j:Number=0;j<m;j++){
AA[i][j]=M[i][j];
}}// end of loops
for(var iter:Number=0;iter<n;iter++){
for(i=0;i<n-1;i++){
if(AA[i+1][sortColumn]>AA[i][sortColumn]){
for(j=0;j<m;j++){
temp=AA[i][j];
AA[i][j]=AA[i+1][j];
AA[i+1][j]=temp;}// end if
}}}// end of loops
return AA;
}// end of function sort
function reverseMatrixSort(M:Array,sortColumn:Number):Array{
// sorts matrix by interchanging rows so that the elements in the matrix
// column "sortColumn" are arranged from smallest(at the top) to largest.
var n:Number=M.length;
var m:Number=M[0].length;
var temp:Number=0;
var AA=new Array(n);
for(var i:Number=0;i<n;i++){AA[i]=new Array(m);
for(var j:Number=0;j<m;j++){
AA[i][j]=M[i][j];
}}// end of loops
for(var iter:Number=0;iter<n;iter++){
for(i=0;i<n-1;i++){
if(AA[i+1][sortColumn]<AA[i][sortColumn]){
for(j=0;j<m;j++){
temp=AA[i][j];
AA[i][j]=AA[i+1][j];
AA[i+1][j]=temp;}// end if
}}}// end of loops
return AA;
}// end of function sort
function screenPrintMatrix(Mat:Array,field:Number,decPosition:Number,separatorStg:String,
fontSize:Number,x:Number,y:Number,id:Number):Void{
// causes a matrix to appear on the screen in a structured format.
// "field" is the length in character spaces of each number or matrix element.
// "decPosition" is the number of digits to be displayed to the right of the decimal place.
// "separatorStg" is the string to be used between numbers. Spaces are a good separator.
// fontSize is the size of font to be used... try 12 to start with.
// x and y are the top left corner of the matrix to be created.
// "id" is an integer used to uniquely name each text field. Set this equal to 0, 1, or 2 or some
// positive integer that is unique for each text field.
// A typical use of this routine might be screenPrintMatrix(M,6,2," ",12,10,20,3);
// This would cause the Matrix M to be printed out with 6 character spaces reserved for
// each number, 2 digits to the right of the decimal, and two spaces between adjacent numbers.
// The top left corner of the matrix would be located at screen position x,y = 10,20 .
// The "id" of this text field would be "txtField3".
if(Mat[0][0]==undefined){trace("toVectorPrint");printVector(Mat,field,decPosition,separatorStg);return}
var rows:Number, columns:Number;
rows=Mat.length;
columns=Mat[0].length;
var out:String;
var depthN:Number=this.getDepth()+id+1;
// trace("461 depth = "+depthN);
var widthN:Number=fontSize*(field+separatorStg.length)*columns;
var heightN:Number=fontSize*rows*2.5;
this.createTextField("txtField"+String(id), depthN, x, y, widthN, heightN);
txt=eval("txtField"+String(id));
txt.type="input"; // allows the user/viewer to modify this text field,
// but it can still be set by the program via myTextField.text=newString;
// The default setting for "type" is "dynamic" that does allow program
// resetting, but not user resetting.
var myTxtFormat=new TextFormat();
myTxtFormat=txt.getTextFormat();
myTxtFormat.size=fontSize;
out="";
for (var i:Number=0;i<rows;i++){
for (var j:Number=0;j<columns;j++){
out+=stringNumber(field,decPosition,Mat[i][j])+separatorStg; }// end of j loop
out+=chr(13)+chr(10);
}// end of i loop
txt.text=out;
txt.setTextFormat(myTxtFormat);//this must be after adding the contents
}// end of function
function screenPrintHorizontalVector(V:Array,field:Number,decPosition:Number,separatorStg:String,
fontSize:Number,x:Number,y:Number,id:Number):Void{
// causes a horizontal vector (1D array) to appear on the screen in a structured format.
// "field" is the length in character spaces of each number or matrix element.
// "decPosition" is the number of digits to be displayed to the right of the decimal place.
// "separatorStg" is the string to be used between numbers. Spaces are a good separator.
// fontSize is the size of font to be used... try 12 to start with.
// x and y are the top left corner of the matrix to be created.
// "id" is an integer used to uniquely name each text field. Set this equal to 0, 1, or 2 or some
// positive integer that is unique for each text field.
// A typical use of this routine might be screenPrintHorizontalVector(M,6,2," ",12,10,20,2);
// This would cause the vector V to be printed out with 6 character spaces reserved for
// each number, 2 digits to the right of the decimal, and two space between adjacent numbers.
// The top left corner of the vector would be located at screen position x,y = 10,20 .
// The "id" of this text field would be "txtField2".
if(V[0]==undefined){trace("toVectorPrint");return}
var columns:Number;
columns=V.length;
var out:String;
var depthN:Number=this.getDepth()+id+1;
// trace("498 depth = "+depthN);
var widthN:Number=fontSize*(field+separatorStg.length)*columns;
var heightN:Number=fontSize*3;
this.createTextField("txtField"+String(id), depthN, x, y, widthN, heightN);
txt=eval("txtField"+String(id));
var myTxtFormat=new TextFormat();
myTxtFormat=txt.getTextFormat();
myTxtFormat.size=fontSize;
out="";
for (var j:Number=0;j<columns;j++){
out+=stringNumber(field,decPosition,V[j])+separatorStg; }// end of j loop
txt.text=out;
txt.setTextFormat(myTxtFormat);
}// end of function
function screenPrintVerticalVector(V:Array,field:Number,decPosition:Number,separatorStg:String,
fontSize:Number,x:Number,y:Number,id:Number):Void{
// causes a vertically oriented vector (array) to appear on the screen in a structured format.
// "field" is the length in character spaces of each number or matrix element.
// "decPosition" is the number of digits to be displayed to the right of the decimal place.
// "separatorStg" is the string to be used between numbers. Spaces are a good separator.
// fontSize is the size of font to be used... try 12 to start with.
// x and y are the top left corner of the matrix to be created.
// "id" is an integer used to uniquely name each text field. Set this equal to 0, 1, or 2 or some
// positive integer that is unique for each text field.
// A typical use of this routine might be screenPrintHorizontalVector(M,6,2," ",12,10,20,4);
// This would cause the vector V to be printed out with 6 character spaces reserved for
// each number, 2 digits to the right of the decimal, and two space between adjacent numbers.
// The top left corner of the matrix would be located at screen position x,y = 10,20 .
// The "id" of this text field would be "txtField4".
if(V[0]==undefined){trace("toVectorPrint");return}
var rows:Number;
rows=V.length;
var out:String;
var depthN:Number=this.getDepth()+id+1;
var widthN:Number=fontSize*(field+separatorStg.length);
var heightN:Number=fontSize*rows*5;
this.createTextField("txtField"+String(id), depthN, x, y, widthN, heightN);
txt=eval("txtField"+String(id));
txt.type="input";// allows the user/viewer to modify this text field.
var myTxtFormat=new TextFormat();
myTxtFormat=txt.getTextFormat();
myTxtFormat.size=fontSize;
out="";
for (var i:Number=0;i<rows;i++){
out+=stringNumber(field,decPosition,V[i])+separatorStg+chr(13)+chr(10); }// end of i loop
txt.text=out;
txt.setTextFormat(myTxtFormat);
}// end of function
function duplicateMatrix(M:Array):Array{
var n:Number=M.length;
var m:Number=M[0].length;
var out=new Array(n);
for(var i:Number=0;i<n;i++){
out[i]=new Array(m);
for(var j:Number=0;j<m;j++){
out[i][j]=M[i][j];
}}// end of loops
return out;
}// end of function duplicate matrix
function stringTo2dArray(inString:String,nRows:Number,mColns:Number):Array{
//assumes elements are followed by spaces and rows by carriage returns
var strg:String=inString; // copy the inString to protect it while changing the working copy
var out:Array=new Array(nRows);
for(var i:Number=0;i<nRows;i++){
out[i]=new Array(mColns);
strg=strg.trimLeft();
for(var j:Number=0;j<mColns;j++){
var indx:Number=strg.indexOf(" ");
var num:Number=Number(strg.slice(0,indx));
out[i][j]=num;
strg=strg.slice(indx);
strg=strg.trimLeft();
}// end of j loop
}// end of i loop
return out
} // end of function stringTo2dArray
function stringToVector(inString:String,n:Number):Array{
var tempM:Array=new Array(n);
var out:Array=new Array(n);
for(var i:Number=0;i<n;i++){
tempM[i]=new Array[1];}
tempM=stringTo2dArray(inString,1,n);
for(j=0;j<n;j++){
out[j]=tempM[0][j];}
return out;
}// end of function stringToVector
