Skip to content

I
imagej-elphel
Commit 99291b6d authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

Simplified, direct decomposition

parent f616c5d7
Hide whitespace changes
Inline Side-by-side
Showing with 119 additions and 100 deletions
src/main/java/com/elphel/imagej/common/CholeskyBlock.java
View file @ 99291b6d
......@@ -170,93 +170,8 @@ public class CholeskyBlock {
return j * (m * np) + ((j >= nf) ? (np-nf*m): m) * m * i;
}
/**
* Get a new diagonal square submatrix
* @param arr A or L packed array (line-scan order for each tile column)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagSquare(
double [] arr,
int i){
int l = (i >= nf) ? (np-nf*m): m;
double [][] a_diag = new double[l][l];
return getDiagSquare (arr, a_diag, i);
}
/**
* Get a new diagonal square submatrix
* @param arr A or L packed array (line-scan order for each tile column)
* @param a_diag - preallocated array (should be smaller if needed)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagSquare(
double [] arr,
double [][] a_diag,
int i){
int indx = indx_IJ(i,i);
int l = a_diag.length;
for (int k = 0; k < l; k++) {
int kl = k*l;
System.arraycopy( arr, indx+kl, a_diag[k], 0, l);
}
return a_diag;
}
/**
* Get a new diagonal square submatrix, copy only lower triangle including diagonal
* @param arr A or L packed array (line-scan order for each tile column)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagTriangle(
double [] arr,
int i){
int l = (i >= nf) ? (np-nf*m): m;
double [][] a_diag = new double[l][l];
return getDiagLTriangle (arr, a_diag, i);
}
/**
* Get a new diagonal square submatrix, copy only lower triangle including diagonal
* @param arr A or L packed array (line-scan order for each tile column)
* @param a_diag - preallocated array (should be smaller if needed)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagLTriangle(
double [] arr,
double [][] a_diag,
int i){
int indx = indx_IJ(i,i);
int l = a_diag.length;
for (int k = 0; k < l; k++) {
int kl = k*l;
System.arraycopy( arr, indx+kl, a_diag[k], 0, k+1);
}
return a_diag;
}
/**
* Save calculated tile L lower diagonal matrix to a packed array
* @param arr A or L packed array (line-scan order for each tile column)
* @param l_diag lower triangular array with Cholesky decomposition
* @param i tile index on a diagonal
*/
public void putDiagLTriangle(
double [] arr,
double [][] l_diag,
int i) {
int indx = indx_IJ(i,i);
int l = l_diag.length;
for (int k = 0; k < l; k++) {
int kl = k*l;
System.arraycopy(l_diag[k], 0, arr, indx+kl, k+1);
}
}
public void setL21(
int i, // i > j,
......@@ -310,6 +225,8 @@ public class CholeskyBlock {
public void choleskyBlockMulti() {
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
cholesky_single(0);
/*
final double [][] Am = new double[m][m];
final double [][] Ah = (n > nf) ? new double[(np-nf*m)][(np-nf*m)] : null;
final double [][] A1 = (np < m) ? Ah : Am; // smaller than a tile
......@@ -330,6 +247,7 @@ public class CholeskyBlock {
L, // double [] arr,
L1, // double [][] l_diag,
0); // int i)
*/
// Calculate first column under diagonal (L21) - maybe use m
ai.set(1); // start from the second tile row
for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
......@@ -356,15 +274,13 @@ public class CholeskyBlock {
threads[ithread] = new Thread() {
public void run() {
for (int nRow= ai.getAndIncrement(); nRow < n; nRow = ai.getAndIncrement()) {
// if ((nRow == (n-1)) && (ftile_diag == (n-1))) {
// System.out.println("calling setA22() first: diag="+(ftile_diag-1)+", nRow="+nRow);
// }
setA22(
ftile_diag-1, // int diag,// < col, < row
nRow, // int row, // >= col
ftile_diag); // int col)
if (nRow == ftile_diag) {
cholesky_single(ftile_diag);
/*
double [][] At = (ftile_diag < nf) ? Am: Ah; // full or reduced array
double [][] Lt = (ftile_diag < nf) ? Lm: Lh; // full or reduced array
// Extract top-left tile (only lower triangle)
......@@ -381,11 +297,6 @@ public class CholeskyBlock {
L, // double [] arr,
Lt, // double [][] l_diag,
ftile_diag); // int i)
/*
// Calculate first column under diagonal (L21) - maybe use m
for (int tr = ftile_diag; tr < n; tr++) {
setL21(tr, ftile_diag);
}
*/
}
}
......@@ -417,9 +328,6 @@ public class CholeskyBlock {
int ncol = (ntile-1) - (nrow * (nrow + 1) /2);
int row = ftile_diag + nrow + 1;
int col = ftile_diag + ncol + 1;
// if ((row == (n-1)) && (col == (n-1))) {
// System.out.println("calling setA22(): diag="+(ftile_diag-1)+", row="+row);
// }
setA22(
ftile_diag-1, // int diag,// < col, < row
row, // int row, // >= col
......@@ -437,11 +345,37 @@ public class CholeskyBlock {
// Cholesky-Banachiewicz Algorithm ?
// Single-threaded
public void cholesky_single(int diag) {
int h = (diag < nf) ? m : (np-nf*m);
int indx = indx_IJ(diag, diag);
Arrays.fill(L, indx, indx+h*h-1, 0);
for (int j = 0; j < h; j++) {
int indx_j = indx+ j * h;
int indx_jj = indx_j + j;
double d = 0.0;
for (int k = 0; k < j; k++) {
int indx_k = indx+ k * h;
double s = 0.0;
for (int i = 0; i < k; i++) {
s += L[indx_k + i] * L[indx_j + i];
}
s = (A[indx_j + k]-s)/L[indx_k+k];
L[indx_j + k] = s;
d = d + s*s;
}
d = A[indx_jj] - d;
L[indx_jj] = Math.sqrt(Math.max(d,0.0));
}
}
// ======================== unused =====================
public static double [][] cholesky_single (
double [][] a,
double [][] l) {
int n = a.length;
// double [][] l = new double[n][n];
for (int j = 0; j < n; j++) {
Arrays.fill(l[j], 0.0);
}
......@@ -463,9 +397,94 @@ public class CholeskyBlock {
}
return l;
}
/**
* Get a new diagonal square submatrix, copy only lower triangle including diagonal
* @param arr A or L packed array (line-scan order for each tile column)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagTriangle(
double [] arr,
int i){
int l = (i >= nf) ? (np-nf*m): m;
double [][] a_diag = new double[l][l];
return getDiagLTriangle (arr, a_diag, i);
}
/**
* Get a new diagonal square submatrix, copy only lower triangle including diagonal
* @param arr A or L packed array (line-scan order for each tile column)
* @param a_diag - preallocated array (should be smaller if needed)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagLTriangle(
double [] arr,
double [][] a_diag,
int i){
int indx = indx_IJ(i,i);
int l = a_diag.length;
for (int k = 0; k < l; k++) {
int kl = k*l;
System.arraycopy( arr, indx+kl, a_diag[k], 0, k+1);
}
return a_diag;
}
/**
* Save calculated tile L lower diagonal matrix to a packed array
* @param arr A or L packed array (line-scan order for each tile column)
* @param l_diag lower triangular array with Cholesky decomposition
* @param i tile index on a diagonal
*/
public void putDiagLTriangle(
double [] arr,
double [][] l_diag,
int i) {
int indx = indx_IJ(i,i);
int l = l_diag.length;
for (int k = 0; k < l; k++) {
int kl = k*l;
System.arraycopy(l_diag[k], 0, arr, indx+kl, k+1);
}
}
/**
* Get a new diagonal square submatrix
* @param arr A or L packed array (line-scan order for each tile column)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagSquare(
double [] arr,
int i){
int l = (i >= nf) ? (np-nf*m): m;
double [][] a_diag = new double[l][l];
return getDiagSquare (arr, a_diag, i);
}
/**
* Get a new diagonal square submatrix
* @param arr A or L packed array (line-scan order for each tile column)
* @param a_diag - preallocated array (should be smaller if needed)
* @param i tile index on the diagonal
* @return new array with tile data (m x m, or smaller for the bottom right)
*/
public double [][] getDiagSquare(
double [] arr,
double [][] a_diag,
int i){
int indx = indx_IJ(i,i);
int l = a_diag.length;
for (int k = 0; k < l; k++) {
int kl = k*l;
System.arraycopy( arr, indx+kl, a_diag[k], 0, l);
}
return a_diag;
}
}
src/main/java/com/elphel/imagej/common/CholeskyLLTMulti.java
View file @ 99291b6d
......@@ -268,7 +268,7 @@ public class CholeskyLLTMulti {
Matrix wjtjlambda_in,
Matrix jty_in,
String title) {
int block_size = 70; // 100; // 70;//70~best // 64; // 60; // 70; // 80; // 120; // 150; // 200; // 100; // 4; // 10; // 100; // 28;
int block_size = 100; // 70;//70~best // 64; // 60; // 70; // 80; // 120; // 150; // 200; // 100; // 4; // 10; // 100; // 28;
boolean truncate = false;
int trunc_size = 199; // 0 to use full size
Matrix wjtjlambda,jty;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment