Mark Bishop

Analytical Chemistry and

Signal Acquisition/Analysis



/**
 * Copyright 2009 Mark Bishop This program is free software: you can
 * redistribute it and/or modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
 * details: http://www.gnu.org/licenses
 * 
 * The author makes no warranty for the accuracy, completeness, safety, or
 * usefulness of any information provided and does not represent that its use
 * would not infringe privately owned right.
 */
 
public class CpuMatrixOps {
 
	/**
	 * 
	 * @param A
	 *            A[m,n]
	 * @param B
	 *            B[m,n]
	 * @return C[m,n] = A + B
	 */
	public static float[][] Add(float[][] A, float[][] B) {
		int m = A.length;
		int n = B[0].length;
		float[][] C = new float[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j] + B[i][j];
			}
		}
		return C;
	}
 
	/**
	 * 
	 * @param a
	 *            a[n]
	 * @param b
	 *            b[n]
	 * @return v[n] = a + b
	 */
	public static float[] Add(float[] a, float[] b) {
		int n = a.length;
		float[] v = new float[n];
		for (int i = 0; i < n; i++) {
			v[i] = a[i] + b[i];
		}
		return v;
	}
 
	/**
	 * 
	 * @param A
	 *            A[m.n]
	 * @return A Transpose [n,m]
	 */
	public static float[][] Transpose(float[][] A) {
		int m = A.length;
		int n = A[0].length;
		float[][] aT = new float[n][m];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				aT[j][i] = A[i][j];
			}
		}
		return aT;
	}
 
	/**
	 * Loop order k i j
	 * 
	 * @param A
	 *            A[m,p]
	 * @param B
	 *            B[p,n]
	 * @return C[m,n] = AB
	 */
	public static float[][] Multiply(float[][] A, float[][] B) {
		int m = A.length;
		int n = B[0].length;
		int p = A[0].length;
		float[][] C = new float[m][n];
		for (int k = 0; k < p; k++) {
			for (int i = 0; i < m; i++) {
				for (int j = 0; j < n; j++) {
					C[i][j] = A[i][k] * B[k][j] + C[i][j];
				}
			}
		}
		return C;
	}
 
	/**
	 * 
	 * @param A
	 *            A[m,n]
	 * @param b
	 *            b[n]
	 * @return v[n] = Ab
	 */
	public static float[] Multiply(float[][] A, float[] b) {
		int m = A.length;
		int n = A[0].length;
		float[] v = new float[m];
		for (int k = 0; k < n; k++) {
			for (int i = 0; i < m; i++) {
				v[i] = A[i][k] * b[k] + v[i];
			}
		}
		return v;
	}
 
	/**
	 * 
	 * @param a
	 *            a[n]
	 * @param b
	 *            a[n]
	 * @return ab[n] = a dot b
	 */
	public static float VectorInnerProduct(float[] a, float[] b) {
		int n = a.length;
		float ab = 0;
		for (int i = 0; i < n; i++) {
			ab += a[i] * b[i];
		}
		return ab;
	}
 
	/**
	 * Tensor product
	 * 
	 * @param vC
	 *            a column vector of length m
	 * @param vR
	 *            a row matrix of length n
	 * @return C[m,n] = vC X vR
	 */
	public static float[][] VectorOuterProduct(float[] vC, float[] vR) {
		int m = vC.length;
		int n = vR.length;
		float[][] C = new float[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = vC[i] * vR[j] + C[i][j];
			}
		}
		return C;
	}
 
	/**
	 * 
	 * @param v
	 * @return vector 2 norm of v
	 */
	public static float Vector2Norm(float[] v) {
		int n = v.length;
		float norm = 0;
		for (int i = 0; i < n; i++) {
			norm = (float) (norm + Math.pow(v[i], 2));
		}
		norm = (float) Math.pow(norm, 0.5);
		return norm;
	}
 
	/**
	 * 
	 * @param v
	 *            v[n]
	 * @return InfinityNorm(v) ... absolute value of highest magnitude element
	 */
	public static float VectorInfinityNorm(float[] v) {
		float lFinity = 0;
		int n = v.length;
		float[] temp = new float[n];
		for (int i = 0; i < n; i++) {
			temp[i] = Math.abs(v[i]);
		}
		lFinity = temp[0];
		for (int i = 1; i < n; i++) {
			if (temp[i] > lFinity) {
				lFinity = temp[i];
			}
		}
		return lFinity;
	}
 
	/**
	 * 
	 * @param alpha
	 *            a scalar
	 * @param v
	 *            v[n]
	 * @return alpha*v[n]
	 */
	public static float[] Scale(float alpha, float[] v) {
		int n = v.length;
		float[] alphaV = new float[n];
		for (int i = 0; i < n; i++) {
			alphaV[i] = alpha * v[i];
		}
		return alphaV;
	}
 
	/**
	 * 
	 * @param alpha
	 *            a scalar
	 * @param A
	 *            A[m,n]
	 * @return alpha*A[m,n]
	 */
	public static float[][] Scale(float alpha, float[][] A) {
		int m = A.length;
		int n = A[0].length;
		float[][] alphaA = new float[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				alphaA[i][j] = alpha * A[i][j];
			}
		}
		return alphaA;
	}
 
	/**
	 * 
	 * @param A
	 *            A[m.n]
	 * @return Frobenius Norm(A)
	 */
	public static float FNorm(float[][] A) {
		int m = A.length;
		int n = A[0].length;
		float fNorm = 0;
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				fNorm = fNorm + (float) Math.pow(A[i][j], 2);
			}
		}
		fNorm = (float) Math.pow(fNorm, 0.5);
		return fNorm;
	}
 
	/**
	 * 
	 * @return gives an upper bound on the relative error due to rounding in
	 *         floating point arithmetic
	 */
	public static float MachineEpsilonfloat() {
		float eps = 1.0f;
		do
			eps /= 2.0;
		while ((float) (1.0 + (eps / 2.0)) != 1.0);
		return eps;
	}
 
	/**
	 * 
	 * @param n
	 * @return I[n,n]
	 */
	public static float[][] Eye(int n) {
		float[][] eye = new float[n][n];
		for (int i = 0; i < n; i++) {
			eye[i][i] = 1.0f;
		}
		return eye;
	}
 
	/**
	 * 
	 * @param A
	 *            A[m,n] where A is upper triangular
	 * @param b
	 *            b[m]
	 * @return x[n] where Ax = b
	 */
	public static float[] BackCalculateX(float[][] upperTriangularMatrix,
			float[] b) {
		int n = upperTriangularMatrix[0].length;
		float temp = 0;
		b[n - 1] = b[n - 1] / upperTriangularMatrix[n - 1][n - 1];
		for (int i = n - 2; i >= 0; i += -1) {
			for (int j = n - 1; j >= i + 1; j += -1) {
				temp = temp - upperTriangularMatrix[i][j] * b[j];
			}
			temp = temp + b[i];
			b[i] = temp / upperTriangularMatrix[i][i];
			temp = 0;
		}
		return b;
	}
 
	/**
	 * Makes a rotator matrix. G(i,k,theta)
	 * 
	 * @param n
	 * @param i
	 * @param k
	 * @param radians
	 * @return
	 */
	public static float[][] Gikt(int n, int i, int k, double theta) {
		double sin = Math.sin(theta);
		double cos = Math.cos(theta);
		float[][] G = CpuMatrixOps.Eye(n);
		G[i][i] = (float) cos;
		G[i][k] = (float) sin;
		G[k][i] = (float) (-1 * sin);
		G[k][k] = (float) cos;
		return G;
	}
}
 
 
ID may be exposed