import javax.xml.ws.soap.MTOM;
/**
* Created by natalia on 11.09.14.
*/
public class Holesskiy {
double alpha = 0.6;
int n = 25;
Matrix a;
int k = 0;
int mkmin1 = 0;
int mk;
int nk;
int iMu1;
int jMu1;
Matrix Akmin1NOSwap = new Matrix(n, n);
Matrix L = new Matrix(n, n);//current mul of Lk-s and Pk-s ~= L1*P1*L2*P2*...
Matrix E;
boolean printResults = false;
public void findMu1() {
if (n - mkmin1 < 2)
return;
double max1 = a.GetElement(mkmin1,mkmin1);
iMu1 = mkmin1;
jMu1 = mkmin1;
for(int i = mkmin1; i < n; i++) {
int j = i;
if(Math.abs(a.GetElement(i, j)) > max1) {
max1 = Math.abs(a.GetElement(i, j));
iMu1 = i;
jMu1 = j;
}
}
nk = 1;
}
private boolean findPk() {
if (printResults)
System.out.println("step " + k);
k++;
//int mkmin1 = 0;
findMu1();
//=====for testing==
Akmin1NOSwap = new Matrix(a);
//============
a.swapRaws(mkmin1, iMu1);
a.swapColumns(mkmin1, iMu1);
return true;
}
private void findAk() {
mk = mkmin1 + nk;
Matrix ak = new Matrix(n, n);
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
ak.setElement(i, j, 0);
for(int i = 0; i < mk; i++)
for (int j = 0; j < mk; j++)
ak.setElement(i, j, a.GetElement(i, j));
//ak.show();
//cut Bk
Matrix Bk = new Matrix(n-mk, nk);
for(int i = 0; i < n-mk; i++)
for(int j = 0; j < nk; j++)
Bk.setElement(i, j, a.GetElement(mk+i, j + mkmin1));
if (printResults) {
System.out.print("\nBk:\n");
Bk.show();
}
//cut Ck
Matrix Ck = new Matrix(n-mk, n-mk);
for(int i = 0; i < n-mk; i++)
for(int j = 0; j < n-mk; j++)
Ck.setElement(i, j, a.GetElement(mk+i, mk+j));
if (printResults) {
System.out.print("\nCk:\n");
Ck.show();
}
//cut Tk
Matrix Tk = new Matrix(nk, nk);
for(int i = 0; i < nk; i++)
for(int j = 0; j < nk; j++)
Tk.setElement(i, j, a.GetElement(i+mkmin1, j+mkmin1));
if (printResults) {
System.out.print("\nTk:\n");
Tk.show();
}
//Tk^-1
Matrix Tkmin1 = new Matrix(Tk);
Tkmin1 = Tkmin1.degMin1();
if (printResults) {
System.out.println("Tk^-1: ");
Tkmin1.show();
System.out.println("\nmust be E:");
Matrix check = Tkmin1.times(Tk);
check.show();
}
//count A~(k)
//count Ck - Bk * Tk^-1 * BkSopryazhonnaya
Matrix BktimesTkmin1 = Bk.times(Tkmin1);
Matrix BkSopr = Bk.transpose();
Matrix forMul = BktimesTkmin1.times(BkSopr);
Matrix podAk = Ck.minus(forMul);
if (printResults) {
System.out.println("\npodmatrix of Ak: ");
podAk.show();
}
//======================count=Lk==and==L==============
Matrix minBktimesTkmin1 = BktimesTkmin1.unMinus();//----------------------
Matrix Lk = new Matrix(n - mkmin1, n - mkmin1);
//count Lk
for (int i = 0; i < n - mkmin1; i++) {
for(int j = 0; j < n - mkmin1; j++) {
if(i == j)
Lk.setElement(i, i, 1);
else if(i >= nk && j < nk)
Lk.setElement(i, j, minBktimesTkmin1.GetElement(i - nk, j));
else
Lk.setElement(i, j, 0);
}
}
if (printResults) {
System.out.println("\n-Bk*Tk^-1: ");
minBktimesTkmin1.show();
System.out.println("\nLk: ");
Lk.show();
}
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
Matrix LkNoSwap = new Matrix(Lk);
Matrix wholeLkNoSwap = new Matrix(n, n);
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
if(i == j)
wholeLkNoSwap.setElement(i, i, 1);
else
wholeLkNoSwap.setElement(i, j, 0);
for(int i = 0; i < n - mkmin1; i++)
for(int j = 0; j < n - mkmin1; j++)
wholeLkNoSwap.setElement(i + mkmin1 , j + mkmin1 , LkNoSwap.GetElement(i, j));
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Lk.swapColumns(0, iMu1 - mkmin1);
Matrix wholeLk = new Matrix(n, n);
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
if(i == j)
wholeLk.setElement(i, i, 1);
else
wholeLk.setElement(i, j, 0);
for(int i = 0; i < n - mkmin1; i++)
for(int j = 0; j < n - mkmin1; j++)
wholeLk.setElement(i + mkmin1 , j + mkmin1 , Lk.GetElement(i, j));
L = wholeLk.times(L);
if (printResults) {
System.out.println("bigLk: \n");
wholeLk.show();
System.out.println("L :");
L.show();
}
//=============end=of=count=L===============
for (int i = 0; i < n - mk; i++)
for (int j = 0; j < n - mk; j++)
ak.setElement(i+mk, j+mk, podAk.GetElement(i, j));
if (printResults) {
System.out.println("\nresult Ak: ");
ak.show();
}
//check:
Matrix whLkTranspose = new Matrix(wholeLkNoSwap);
whLkTranspose = whLkTranspose.transpose();
Matrix b = wholeLkNoSwap.times(a);
Matrix checkKResult = b.times(whLkTranspose);
if (printResults) {
System.out.println("must be 0");
checkKResult.minus(ak).show();
}
//change values for the next step:
a = ak;
mkmin1 = mk;
}
public Matrix solveAxEqualB(Matrix b) {
Matrix Ltrans = L.transpose();
Matrix result = Ltrans.times(a.degMin1().times(L));//A^-1 = L^t * Ak^-1 * L
Matrix X = result.times(b);
return X;
}
public void testH() {
//a = Matrix.Hilbert(n, n);
a = Matrix.matrix25();
//a = new Matrix(new double[][] {{1,5,4}, {5,3,2}, {4,2,8}});
//a = randomSymmetricMatrix();
System.out.println("matrix a (initial):");
for (int i = 0; i < a.M; i++) {
for (int j = 0; j < a.N; j++)
System.out.print(" "+ a.data[i][j]);
System.out.println();
}
System.out.println();
Matrix initial = new Matrix(a);
mk = 1;
mkmin1 = 0;
a.show();
for(int i = 0;i < n; i++)
for(int j = 0; j < n; j++)
L.setElement(i, j, 0);
for (int i = 0 ; i < n; i++)
L.setElement(i, i, 1);
E = new Matrix(L);
while (mk < n -1) {
if(findPk())
findAk();
else
break;
}
Matrix Ltrans = L.transpose();
Matrix result = Ltrans.times(a.degMin1().times(L));
Matrix check = initial.degMin1().minus(result);
if (printResults) {
System.out.println("check : A^-1 - ( L^t * Ak^-1 * L )");
check.show();
}
Matrix B = new Matrix(n, 1);
for (int i = 0; i < n; i++) {
B.setElement(i, 0, 1);
}
Matrix X = solveAxEqualB(B);
if (printResults) {
System.out.println("x = A^-1 * b: ");
X.show();
}
System.out.println("check: b - A*x: ");
check = B.minus(initial.times(X));
System.out.println("diff");
check.show();
float normx = (float) 0.0;
for(int i = 0; i < B.M; i++){
normx += check.data[i][0] * check.data[i][0];
}
normx = (float)Math.sqrt(normx);
float score = (float) ((normx - Math.sqrt(n)) / Math.sqrt(n));
System.out.println(score);
System.out.println("hi, we've gone throw " + k + " steps;)");
k = 0;
}
public static void main(String args[]) {
Holesskiy b = new Holesskiy();
b.testH();
}
}