/*
* Open Source Physics software is free software as described near the bottom of this code file.
*
* For additional information and documentation on Open Source Physics please see:
* <http://www.opensourcephysics.org/>
*/
/**
Adapted to OSP by Javier E. Hasbun, 2009.
<P>
@Copyright (c) 2009
This software is to support the Open Source Physics library
http://www.opensourcephysics.org under the terms of the GNU General Public
License (GPL) as published by the Free Software Foundation.
**/
/**
ComplexMatrix.java
solve, invert, etc. last argument is output
void solve(Complex A[][], Complex Y[], Complex X[]); X=A^-1*Y
void invert(Complex A[][]); A=A^-1
Complex determinant(Complex A[]); d=det A
void eigenvalues(Complex A[][], ComplexV[][], Complex Y[]); V,Y=eigen A
void eigenCheck(Complex A[][], V[][], Complex Y[]); printout
void multiply(Complex A[][], Complex B[][], Complex C[][]); C=A*B
void add(Complex A[][], Complex B[][], Complex C[][]); C=A+B
void subtract(Complex C[][], Complex A[][], complex B[][]); C=A-B
double norm1(Complex A[][]); d=norm1 A
double norm2(Complex A[][]); sqrt largest eigenvalue A^T*A d=norm2 A
double normFro(Complex A[][]); Frobenius d=normFro A
double normInf(Complex A[][]); d=normInf A
void copy(Complex A[][], Complex B[][]); B=A
void boolean equals(Complex A[][], Complex B[][]); A==B
void fromDouble(double A[], double b[][], Complex C[][]); C=(A,b)
void fromDouble(double A[], Complex B[][]); B=A
void identity(A[][]); A=I
void zero(A[][]); A=0
void print(Complex A[][]); A
void multiply(Complex A[][], Complex X[], Complex Y[]); Y=A*X
void add(Complex X[], Complex Y[], Complex Z[]); Z=X+Y
void subtract(Complex X[], Complex Y[], Complex Z[]); Z=X-Y
double norm1(Complex X[]); d=norm1 X
double norm2(Complex X[]); d=norm2 X
double normInf(Complex X[]); d=normInf X
void copy(Complex X[], Complex Y[]); Y=X
void boolean equals(Complex X[], Complex Y[]); X==Y
void fromDouble(double X[], double Y[], Complex Z[]); Z=(X,Y)
void fromDouble(double X[], Complex Y[]); Y=(X,0)
void unitVector(X[],I); X[I]=1, else 0
void zero(X[]); X=0
void print(Complex X[]); X
void readSize(String, rowCol[])
void read(String, Complex[][])
void closeInput()
void write(String, Complex[][])
void closeOutput()
**/
package org.opensourcephysics.numerics;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.PrintWriter;
/**
* Class description
*
*/
public strictfp class ComplexMatrix {
public static void solve(final Complex A[][], final Complex Y[], Complex X[]) {
// solve complex linear equations for X where Y = A * X
// method: Gauss-Jordan elimination using maximum pivot
// usage: ComplexMatrix.solve(A,Y,X);
// Translated to java by : Jon Squire , 4 April 2003
// First written by Jon Squire December 1959 for IBM 650, translated to
// other languages e.g. Fortran converted to Ada converted to C
// converted to java
int n = A.length;
int m = n+1;
Complex B[][] = new Complex[n][m]; // working matrix
int row[] = new int[n]; // row interchange indicies
int hold, I_pivot; // pivot indicies
Complex pivot; // pivot element value
double abs_pivot;
if((A[0].length!=n)||(Y.length!=n)||(X.length!=n)) {
System.out.println("Error in ComplexMatrix.solve inconsistent array sizes."); //$NON-NLS-1$
}
// build working data structure
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
B[i][j] = A[i][j];
}
B[i][n] = Y[i];
}
// set up row interchange vectors
for(int k = 0; k<n; k++) {
row[k] = k;
}
// begin main reduction loop
for(int k = 0; k<n; k++) {
// find largest element for pivot
pivot = B[row[k]][k];
abs_pivot = pivot.abs();
I_pivot = k;
for(int i = k+1; i<n; i++) {
if(B[row[i]][k].abs()>abs_pivot) {
I_pivot = i;
pivot = B[row[i]][k];
abs_pivot = pivot.abs();
}
}
// have pivot, interchange row indicies
hold = row[k];
row[k] = row[I_pivot];
row[I_pivot] = hold;
// check for near singular
if(abs_pivot<1.0E-10) {
for(int j = k+1; j<n+1; j++) {
B[row[k]][j] = new Complex(0.0, 0.0);
}
System.out.println("redundant row (singular) "+row[k]); //$NON-NLS-1$
} // singular, delete row
else {
// reduce about pivot
for(int j = k+1; j<n+1; j++) {
B[row[k]][j] = B[row[k]][j].div(B[row[k]][k]);
}
// inner reduction loop
for(int i = 0; i<n; i++) {
if(i!=k) {
for(int j = k+1; j<n+1; j++) {
B[row[i]][j] = B[row[i]][j].subtract(B[row[i]][k].mul(B[row[k]][j]));
}
}
}
}
// finished inner reduction
}
// end main reduction loop
// build X for return, unscrambling rows
for(int i = 0; i<n; i++) {
X[i] = B[row[i]][n];
}
} // end solve
public static final void invert(Complex A[][]) {
int n = A.length;
int row[] = new int[n];
int col[] = new int[n];
Complex temp[] = new Complex[n];
int hold, I_pivot, J_pivot;
Complex pivot;
double abs_pivot;
if(A[0].length!=n) {
System.out.println("Error in Complex.Matrix.invert,"+ //$NON-NLS-1$
" matrix not square."); //$NON-NLS-1$
}
// set up row and column interchange vectors
for(int k = 0; k<n; k++) {
row[k] = k;
col[k] = k;
}
// begin main reduction loop
for(int k = 0; k<n; k++) {
// find largest element for pivot
pivot = A[row[k]][col[k]];
I_pivot = k;
J_pivot = k;
for(int i = k; i<n; i++) {
for(int j = k; j<n; j++) {
abs_pivot = pivot.abs();
if(A[row[i]][col[j]].abs()>abs_pivot) {
I_pivot = i;
J_pivot = j;
pivot = A[row[i]][col[j]];
}
}
}
if(pivot.abs()<1.0E-10) {
System.out.println("ComplexMatrix is singular !"); //$NON-NLS-1$
return;
}
hold = row[k];
row[k] = row[I_pivot];
row[I_pivot] = hold;
hold = col[k];
col[k] = col[J_pivot];
col[J_pivot] = hold;
// reduce about pivot
A[row[k]][col[k]] = (new Complex(1.0, 0.0)).div(pivot);
for(int j = 0; j<n; j++) {
if(j!=k) {
A[row[k]][col[j]] = A[row[k]][col[j]].mul(A[row[k]][col[k]]);
}
}
// inner reduction loop
for(int i = 0; i<n; i++) {
if(k!=i) {
for(int j = 0; j<n; j++) {
if(k!=j) {
A[row[i]][col[j]] = A[row[i]][col[j]].subtract(A[row[i]][col[k]].mul(A[row[k]][col[j]]));
}
}
A[row[i]][col[k]] = A[row[i]][col[k]].mul(A[row[k]][col[k]]).neg();
}
}
}
// end main reduction loop
// unscramble rows
for(int j = 0; j<n; j++) {
for(int i = 0; i<n; i++) {
temp[col[i]] = A[row[i]][j];
}
for(int i = 0; i<n; i++) {
A[i][j] = temp[i];
}
}
// unscramble columns
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
temp[row[j]] = A[i][col[j]];
}
for(int j = 0; j<n; j++) {
A[i][j] = temp[j];
}
}
} // end invert
public static final Complex determinant(final Complex A[][]) {
int n = A.length;
Complex D = new Complex(1.0, 0.0); // determinant
Complex B[][] = new Complex[n][n]; // working matrix
int row[] = new int[n]; // row interchange indicies
int hold, I_pivot; // pivot indicies
Complex pivot; // pivot element value
double abs_pivot;
if(A[0].length!=n) {
System.out.println("Error in ComplexMatrix.determinant,"+ //$NON-NLS-1$
" inconsistent array sizes."); //$NON-NLS-1$
}
// build working matrix
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
B[i][j] = A[i][j];
}
}
// set up row interchange vectors
for(int k = 0; k<n; k++) {
row[k] = k;
}
// begin main reduction loop
for(int k = 0; k<n-1; k++) {
// find largest element for pivot
pivot = B[row[k]][k];
abs_pivot = pivot.abs();
I_pivot = k;
for(int i = k; i<n; i++) {
if(B[row[i]][k].abs()>abs_pivot) {
I_pivot = i;
pivot = B[row[i]][k];
abs_pivot = pivot.abs();
}
}
// have pivot, interchange row indicies
if(I_pivot!=k) {
hold = row[k];
row[k] = row[I_pivot];
row[I_pivot] = hold;
D = D.neg();
}
// check for near singular
if(abs_pivot<1.0E-10) {
return new Complex(0.0, 0.0);
}
D = D.mul(pivot);
// reduce about pivot
for(int j = k+1; j<n; j++) {
B[row[k]][j] = B[row[k]][j].div(B[row[k]][k]);
}
// inner reduction loop
for(int i = 0; i<n; i++) {
if(i!=k) {
for(int j = k+1; j<n; j++) {
B[row[i]][j] = B[row[i]][j].subtract(B[row[i]][k].mul(B[row[k]][j]));
}
}
}
// finished inner reduction
}
// end of main reduction loop
return D.mul(B[row[n-1]][n-1]);
} // end determinant
public static final void eigenvalues(final Complex A[][], Complex V[][], Complex Y[]) {
// cyclic Jacobi iterative method of finding eigenvalues
// advertized for symmetric real
int n = A.length;
Complex AA[][] = new Complex[n][n];
// double norm; // removed by W. Christian
Complex c = new Complex(1.0, 0.0);
Complex s = new Complex(0.0, 0.0);
if((A[0].length!=n)||(V.length!=n)||(V[0].length!=n)||(Y.length!=n)) {
System.out.println("Error in ComplexMatrix.eigenvalues,"+ //$NON-NLS-1$
" inconsistent array sizes."); //$NON-NLS-1$
}
identity(V); // start V as identity matrix
copy(A, AA);
for(int k = 0; k<n; k++) {
// norm=norm4(AA); // removed by W. Christian
for(int i = 0; i<n-1; i++) {
for(int j = i+1; j<n; j++) {
schur2(AA, i, j, c, s);
mat44(i, j, c, s, AA, V);
}
} // end one iteration
}
// norm = norm4(AA); // final quality check if desired // removed by W. Christian
for(int i = 0; i<n; i++) {
// copy eigenvalues back to caller
Y[i] = AA[i][i];
}
} // end eigenvalues
public static final void eigenCheck(final Complex A[][], final Complex V[][], final Complex Y[]) {
if((A==null)||(V==null)||(Y==null)) {
return;
}
// check A * X = lambda X lambda=Y[i] X=V[i]
// check determinant(A- lambda I) = 0
int n = A.length;
Complex B[][] = new Complex[n][n];
Complex C[][] = new Complex[n][n];
Complex X[] = new Complex[n];
Complex Z[] = new Complex[n];
Complex T[] = new Complex[n];
double norm = 0.0;
if((A[0].length!=n)||(V.length!=n)||(V[0].length!=n)||(Y.length!=n)) {
System.out.println("Error in ComplexMatrix.eigenCheck,"+ //$NON-NLS-1$
" inconsistent array sizes."); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
X[j] = V[j][i];
}
mul(A, X, T);
for(int j = 0; j<n; j++) {
Z[j] = T[j].subtract(Y[i].mul(X[j]));
}
System.out.println("check for near zero norm of Z["+i+"]="+Z[i]); //$NON-NLS-1$ //$NON-NLS-2$
}
norm = norm2(Z);
System.out.println("norm ="+norm+" is eigen vector error indication 1."); //$NON-NLS-1$ //$NON-NLS-2$
System.out.println("det V = "+ComplexMatrix.determinant(V)); //$NON-NLS-1$
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
Z[j] = V[j][i];
}
System.out.println("check for 1.0 = "+norm2(Z)); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
identity(B);
mul(B, Y[i], C);
subtract(A, C, B);
Z[i] = determinant(B);
}
norm = norm2(Z);
System.out.println("norm ="+norm+" is eigen value error indication."); //$NON-NLS-1$ //$NON-NLS-2$
} // end eigenCheck
static void schur2(final Complex A[][], final int p, final int q, Complex c, Complex s) {
Complex tau;
Complex tau_tau_1;
Complex t;
if(A[0].length!=A.length) {
System.out.println("Error in schur2 of Complex jacobi,"+ //$NON-NLS-1$
" inconsistent array sizes."); //$NON-NLS-1$
}
if(A[p][q].abs()!=0.0) {
tau = (A[q][q].subtract(A[p][p])).div((A[p][q].mul(2.0)));
tau_tau_1 = (tau.mul(tau).add(1.0)).sqrt();
if(tau.abs()>=0.0) {
t = tau.add(tau_tau_1).invert();
} else {
t = (tau_tau_1.subtract(tau)).invert().neg();
}
c = (t.mul(t)).add(1.0).sqrt().invert();
s = t.mul(c);
} else {
c = new Complex(1.0, 0.0);
s = new Complex(0.0, 0.0);
}
} // end schur2
static void mat22(final Complex c, final Complex s, final Complex A[][], Complex B[][]) {
if((A.length!=2)||(A[0].length!=2)||(B.length!=2)||(B[0].length!=2)) {
System.out.println("Error in mat22 of Jacobi, not both 2 by 2"); //$NON-NLS-1$
}
Complex T[][] = new Complex[2][2];
T[0][0] = c.mul(A[0][0]).subtract(s.mul(A[0][1]));
T[0][1] = s.mul(A[0][0]).add(c.mul(A[0][1]));
T[1][0] = c.mul(A[1][0]).subtract(s.mul(A[1][1]));
T[1][1] = s.mul(A[1][0]).add(c.mul(A[1][1]));
B[0][0] = c.mul(T[0][0]).subtract(s.mul(T[1][0]));
B[0][1] = c.mul(T[0][1]).subtract(s.mul(T[1][1]));
B[1][0] = s.mul(T[0][0]).add(c.mul(T[1][0]));
B[1][1] = s.mul(T[0][1]).add(c.mul(T[1][1]));
} // end mat2
static void mat44(final int p, final int q, final Complex c, final Complex s, final Complex A[][], Complex V[][]) {
int n = A.length;
Complex B[][] = new Complex[n][n];
Complex J[][] = new Complex[n][n];
if((A[0].length!=n)||(V.length!=n)||(V[0].length!=n)) {
System.out.println("Error in mat44 of Complex Jacobi,"+ //$NON-NLS-1$
" A or V not same and square"); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
J[i][j] = new Complex(0.0, 0.0);
}
J[i][i] = new Complex(1.0, 0.0);
}
J[p][p] = c; /* J transpose */
J[p][q] = s.neg();
J[q][q] = c;
J[q][p] = s;
mul(J, A, B);
J[p][q] = s;
J[q][p] = s.neg();
mul(B, J, A);
mul(V, J, B);
copy(B, V);
} // end mat44
static double norm4(final Complex A[][]) // for Jacobi
{
int n = A.length;
int nr = A[0].length;
double nrm = 0.0;
if(n!=nr) {
System.out.println("Error in Complex norm4, non square A["+ //$NON-NLS-1$
n+"]["+nr+"]"); //$NON-NLS-1$ //$NON-NLS-2$
}
for(int i = 0; i<n-1; i++) {
for(int j = i+1; j<n; j++) {
nrm = nrm+A[i][j].abs()+A[j][i].abs();
}
}
return nrm/(n*n-n);
} // end norm4
public static final void mul(final Complex A[][], final Complex B[][], Complex C[][]) {
int ni = A.length;
int nk = A[0].length;
int nj = B[0].length;
if((B.length!=nk)||(C.length!=ni)||(C[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.mul,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = new Complex(0.0, 0.0);
for(int k = 0; k<nk; k++) {
C[i][j] = C[i][j].add(A[i][k].mul(B[k][j]));
}
}
}
} // end mul
public static final void mul(final Complex A[][], final Complex B, Complex C[][]) {
int ni = A.length;
int nj = A[0].length;
if((C.length!=ni)||(C[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.mul,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = A[i][j].mul(B);
}
}
} // end mul
public static final void add(final Complex A[][], final Complex B[][], Complex C[][]) {
int ni = A.length;
int nj = A[0].length;
if((B.length!=ni)||(C.length!=ni)||(B[0].length!=nj)||(C[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.add,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = A[i][j].add(B[i][j]);
}
}
} // end add
public static final void add(final Complex A[][], final Complex B, Complex C[][]) {
int ni = A.length;
int nj = A[0].length;
if((C.length!=ni)||(C[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.add,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = A[i][j].add(B);
}
}
} // end add
public static final void subtract(final Complex A[][], final Complex B[][], Complex C[][]) {
int ni = A.length;
int nj = A[0].length;
if((B.length!=ni)||(C.length!=ni)||(B[0].length!=nj)||(C[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.subtract,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = A[i][j].subtract(B[i][j]);
}
}
} // end subtract
public static final void subtract(final Complex A[][], final Complex B, Complex C[][]) {
int ni = A.length;
int nj = A[0].length;
if((C.length!=ni)||(C[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.subtract,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = A[i][j].subtract(B);
}
}
} // end subtract
public static final double norm1(final Complex A[][]) {
double norm = 0.0;
double colSum;
int ni = A.length;
int nj = A[0].length;
for(int j = 0; j<nj; j++) {
colSum = 0.0;
for(int i = 0; i<ni; i++) {
colSum = colSum+A[i][j].abs();
}
norm = Math.max(norm, colSum);
}
return norm;
} // end norm1
public static final double normInf(final Complex A[][]) {
double norm = 0.0;
double rowSum;
int ni = A.length;
int nj = A[0].length;
for(int i = 0; i<ni; i++) {
rowSum = 0.0;
for(int j = 0; j<nj; j++) {
rowSum = rowSum+A[i][j].abs();
}
norm = Math.max(norm, rowSum);
}
return norm;
} // end normInf
public static final double normFro(final Complex A[][]) {
double norm = 0.0;
int n = A.length;
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
norm = norm+A[i][j].abs()*A[i][j].abs();
}
}
return Math.sqrt(norm);
} // end normFro
public static final double norm2(final Complex A[][]) {
double r = 0.0; // largest eigenvalue
int n = A.length;
Complex B[][] = new Complex[n][n];
Complex V[][] = new Complex[n][n];
Complex BI[] = new Complex[n];
if(A[0].length!=n) {
System.out.println("Error in ComplexMatrix.norm2,"+ //$NON-NLS-1$
" matrix not square."); //$NON-NLS-1$
}
for(int i = 0; i<n;
i++) // B = A^T * A
{
for(int j = 0; j<n; j++) {
B[i][j] = new Complex(0.0, 0.0);
for(int k = 0; k<n; k++) {
B[i][j] = B[i][j].add(A[k][i].mul(A[k][j]));
}
}
}
eigenvalues(B, V, BI);
for(int i = 0; i<n; i++) {
r = Math.max(r, BI[i].abs());
}
return Math.sqrt(r);
} // end subtract
public static final void copy(final Complex A[][], Complex B[][]) {
int ni = A.length;
int nj = A[0].length;
if((B.length!=ni)||(B[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.copy,"+ //$NON-NLS-1$
" inconsistent sizes."); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
B[i][j] = A[i][j];
}
}
} // end copy
public static final boolean equals(final Complex A[][], final Complex B[][]) {
int ni = A.length;
int nj = A[0].length;
boolean same = true;
if((B.length!=ni)||(B[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.equals,"+ //$NON-NLS-1$
" inconsistent sizes."); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
same = same&&A[i][j].equals(B[i][j]);
}
}
return same;
} // end equals
public static final void fromDouble(final double A[][], final double B[][], Complex C[][]) {
int ni = A.length;
int nj = A[0].length;
if((C.length!=ni)||(C[0].length!=nj)||(B.length!=ni)||(B[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.fromDouble,"+ //$NON-NLS-1$
" inconsistent sizes."); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = new Complex(A[i][j], B[i][j]);
}
}
} // end fromDouble
public static final void fromDouble(final double A[][], Complex C[][]) {
int ni = A.length;
int nj = A[0].length;
if((C.length!=ni)||(C[0].length!=nj)) {
System.out.println("Error in ComplexMatrix.fromDouble,"+ //$NON-NLS-1$
" inconsistent sizes."); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
C[i][j] = new Complex(A[i][j]);
}
}
} // end fromDouble
public static final void identity(Complex A[][]) {
int n = A.length;
if(n!=A[0].length) {
System.out.println("Error in ComplexMatrix.identity,"+ //$NON-NLS-1$
" inconsistent sizes."); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
for(int j = 0; j<n; j++) {
A[i][j] = new Complex(0.0);
}
A[i][i] = new Complex(1.0);
}
} // end identity
public static final void zero(Complex A[][]) {
int ni = A.length;
int nj = A[0].length;
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
A[i][j] = new Complex(0.0);
}
}
} // end zero
public static final void print(final Complex A[][]) {
int ni = A.length;
int nj = A[0].length;
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
System.out.println("A["+i+"]["+j+"]="+A[i][j]); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
}
}
} // end print
public static final void mul(final Complex A[][], final Complex B[], Complex C[]) {
int ni = A.length;
int nj = A[0].length;
if((B.length!=nj)||(C.length!=ni)) {
System.out.println("Error in ComplexMatrix.mul,"+ //$NON-NLS-1$
" incompatible sizes."); //$NON-NLS-1$
}
for(int i = 0; i<ni; i++) {
C[i] = new Complex(0.0, 0.0);
for(int j = 0; j<nj; j++) {
C[i] = C[i].add(A[i][j].mul(B[j]));
}
}
} // end mul
public static final void add(final Complex X[], final Complex Y[], Complex Z[]) {
int n = X.length;
if((Y.length!=n)||(Z.length!=n)) {
System.out.println("Error in ComplexMatrix.add,"+ //$NON-NLS-1$
" incompatible sizes."); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
Z[i] = X[i].add(Y[i]);
}
} // end add
public static final void subtract(final Complex X[], final Complex Y[], Complex Z[]) {
int n = X.length;
if((Y.length!=n)||(Z.length!=n)) {
System.out.println("Error in ComplexMatrix.subtract,"+ //$NON-NLS-1$
" incompatible sizes."); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
Z[i] = X[i].subtract(Y[i]);
}
} // end subtract
public static final double norm1(final Complex X[]) {
double norm = 0.0;
int n = X.length;
for(int i = 0; i<n; i++) {
norm = norm+X[i].abs();
}
return norm;
} // end norm1
public static final double norm2(final Complex X[]) {
double norm = 0.0;
int n = X.length;
for(int i = 0; i<n; i++) {
norm = norm+X[i].abs()*X[i].abs();
}
return StrictMath.sqrt(norm);
} // end norm2
public static final double normInf(final Complex X[]) {
double norm = 0.0;
int n = X.length;
for(int i = 0; i<n; i++) {
norm = Math.max(norm, X[i].abs());
}
return norm;
} // end normInf
public static final void copy(final Complex X[], Complex Y[]) {
int n = X.length;
if(Y.length!=n) {
System.out.println("Error in ComplexMatrix.copy,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
Y[i] = X[i];
}
} // end copy
public static final boolean equals(final Complex X[], final Complex Y[]) {
int n = X.length;
boolean same = true;
if(Y.length!=n) {
System.out.println("Error in ComplexMatrix.equals,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
same = same&&X[i].equals(Y[i]);
}
return same;
} // end equals
public static final void fromDouble(final double X[], Complex Z[]) {
int n = X.length;
if(Z.length!=n) {
System.out.println("Error in ComplexMatrix.fromDouble,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
Z[i] = new Complex(X[i]);
}
} // end fromDouble
public static final void fromDouble(final double X[], final double Y[], Complex Z[]) {
int n = X.length;
if((Z.length!=n)||(Y.length!=n)) {
System.out.println("Error in ComplexMatrix.fromDouble,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
for(int i = 0; i<n; i++) {
Z[i] = new Complex(X[i], Y[i]);
}
} // end fromDouble
public static void fromRoots(Complex X[], Complex Y[]) {
int n = X.length;
if(Y.length!=n+1) {
System.out.println("Error in ComplexMatrix.fromRoots,"+ //$NON-NLS-1$
" incompatible sizes"); //$NON-NLS-1$
}
Y[0] = X[0].neg();
Y[1] = new Complex(1.0);
if(n==1) {
return;
}
for(int i = 1; i<n; i++) {
Y[i+1] = new Complex(0.0);
for(int j = 0; j<=i; j++) {
Y[i+1-j] = Y[i-j].subtract(Y[i+1-j].mul(X[i]));
}
Y[0] = Y[0].mul(X[i]).neg();
}
}
public static final void unitVector(Complex X[], int j) {
int n = X.length;
for(int i = 0; i<n; i++) {
X[i] = new Complex(0.0);
}
X[j] = new Complex(1.0);
} // end unitVector
public static final void zero(Complex X[]) {
int n = X.length;
for(int i = 0; i<n; i++) {
X[i] = new Complex(0.0);
}
} // end zero
public static final void print(final Complex X[]) {
int n = X.length;
for(int i = 0; i<n; i++) {
System.out.println("X["+i+"]="+X[i]); //$NON-NLS-1$ //$NON-NLS-2$
}
} // end print
private static String input_file_name; // for read
private static BufferedReader in; // for read
private static String output_file_name; // for write
private static BufferedWriter out; // for write
private static PrintWriter file_out; // for write
public static final void readSize(String file_name, int rowCol[]) {
String input_line = new String("@"); // unread //$NON-NLS-1$
int len; // input_line length
int index; // start of field
int last; // end of field
String intStr; // for parseInt (does not ignore leading and
// trailing white space)
int ni; // number of rows
int nj; // number of columns
if((input_file_name==null)||!file_name.equals(input_file_name)) {
input_file_name = file_name;
try {
in = new BufferedReader(new FileReader(file_name));
} catch(Exception e) {
System.out.println("ComplexMatrix.read unable to open file "+file_name); //$NON-NLS-1$
return;
}
}
ni = 0;
nj = 0;
try {
input_line = in.readLine(); // first read before 'while'
while(input_line!=null) // allow leading blank lines
{
input_line = input_line.trim();
len = input_line.length();
if(len==0) {
input_line = in.readLine();
continue; // skip blank lines
}
if(input_line.charAt(0)=='(') {
System.out.println("ComplexMatrix.readSize unable to get size "+file_name); //$NON-NLS-1$
break;
}
index = 0;
last = input_line.indexOf(' '); // first blank after number
if(last==-1) {
last = len;
}
intStr = input_line.substring(index, last);
ni = Integer.parseInt(intStr);
input_line = input_line.substring(last, len);
input_line = input_line.trim();
len = input_line.length();
if(len==0) {
nj = ni;
break;
}
index = 0;
last = input_line.indexOf(' '); // first blank after number
if(last==-1) {
last = len;
}
intStr = input_line.substring(index, last);
nj = Integer.parseInt(intStr);
break;
} // end while
} catch(Exception e) {
System.out.println("ComplexMatrix.readSize unable to get size "+file_name); //$NON-NLS-1$
}
rowCol[0] = ni;
rowCol[1] = nj;
} // end readSize
public static final void read(String file_name, Complex A[][]) {
String input_line = new String("@"); // unread //$NON-NLS-1$
int len; // input_line length
int index; // start of field
int last; // end of field
String intStr; // for parseInt (does not ignore leading and
// trailing white space)
int i, ni; // 0..ni-1
int j, nj; // 0..nj-1
boolean have_line = false;
if((input_file_name==null)||!file_name.equals(input_file_name)) {
input_file_name = file_name;
try {
in = new BufferedReader(new FileReader(file_name));
} catch(Exception e) {
System.out.println("ComplexMatrix.read unable to open file "+file_name); //$NON-NLS-1$
return;
}
}
ni = 0;
nj = 0;
try {
input_line = in.readLine(); // first read before 'while'
while(input_line!=null) // allow leading blank lines
{
input_line = input_line.trim();
len = input_line.length();
if(len==0) {
input_line = in.readLine();
continue; // skip blank lines
}
if(input_line.charAt(0)=='(') {
ni = A.length; // no size, just complex data
nj = A[0].length;
have_line = true;
break;
}
index = 0;
last = input_line.indexOf(' '); // first blank after number
if(last==-1) {
last = len;
}
intStr = input_line.substring(index, last);
ni = Integer.parseInt(intStr);
input_line = input_line.substring(last, len);
input_line = input_line.trim();
len = input_line.length();
if(len==0) {
nj = ni;
break;
}
index = 0;
last = input_line.indexOf(' '); // first blank after number
if(last==-1) {
last = len;
}
intStr = input_line.substring(index, last);
nj = Integer.parseInt(intStr);
break;
} // end while
} catch(Exception e) {
System.out.println("ComplexMatrix.read unable to get size "+file_name); //$NON-NLS-1$
}
// now read complex numbers
i = 0;
j = 0;
if((A.length!=ni)||(A[0].length!=nj)) {
System.out.println("incompatible size in ComplexMatrix.read"); //$NON-NLS-1$
return;
}
try {
if(!have_line) {
input_line = in.readLine(); // first read before 'while'
}
have_line = false;
while(input_line!=null) {
input_line = input_line.trim();
len = input_line.length();
if(len==0) {
input_line = in.readLine();
continue; // skip blank lines
}
index = 0;
last = input_line.indexOf(')'); // closing )
if(last==-1) {
input_line = in.readLine();
continue;
}
intStr = input_line.substring(index, last+1);
A[i][j] = Complex.parseComplex(intStr);
j++;
if(j==nj) {
j = 0;
i++;
}
if(i==ni) {
break;
}
input_line = input_line.substring(last+1);
} // end while
} catch(Exception e) {
System.out.println("ComplexMatrix.read unable to read data "+file_name); //$NON-NLS-1$
}
} // end read
public static final void closeInput() {
try {
in.close();
} catch(Exception e) {
System.out.println("ComplexMatrix.closeInput not closed"); //$NON-NLS-1$
}
input_file_name = null;
} // end closeInput
public static final void write(String file_name, Complex A[][]) {
int ni = A.length;
int nj = A[0].length;
if((output_file_name==null)||!file_name.equals(output_file_name)) {
output_file_name = file_name;
try {
out = new BufferedWriter(new FileWriter(file_name));
file_out = new PrintWriter(out);
} catch(Exception e) {
System.out.println("ComplexMatrix.write unable to open file "+file_name); //$NON-NLS-1$
return;
}
}
// write size
if(ni==nj) {
file_out.println(ni);
} else {
file_out.println(ni+" "+nj); //$NON-NLS-1$
}
// now write complex numbers
try {
for(int i = 0; i<ni; i++) {
for(int j = 0; j<nj; j++) {
file_out.println(A[i][j].toString());
}
}
file_out.println();
} catch(Exception e) {
System.out.println("ComplexMatrix.write unable to write data "+file_name); //$NON-NLS-1$
}
} // end write
public static final void closeOutput() {
file_out.close();
output_file_name = null;
} // end closeOutput
} // end class ComplexMatrix
/*
* Open Source Physics software is free software; you can redistribute
* it and/or modify it under the terms of the GNU General Public License (GPL) as
* published by the Free Software Foundation; either version 2 of the License,
* or(at your option) any later version.
* Code that uses any portion of the code in the org.opensourcephysics package
* or any subpackage (subdirectory) of this package must must also be be released
* under the GNU GPL license.
*
* This software 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.
*
* You should have received a copy of the GNU General Public License
* along with this; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
* or view the license online at http://www.gnu.org/copyleft/gpl.html
*
* Copyright (c) 2007 The Open Source Physics project
* http://www.opensourcephysics.org
*/