Examples with SparseMatrixCell maspack.matrix.SparseMatrixCell used on opensource projects

Example 1 with SparseMatrixCell

use of maspack.matrix.SparseMatrixCell in project artisynth_core by artisynth.

the class IncompleteLUDecomposition method factor.

/**
 * Peforms an Cholesky decomposition on the Matrix M.
 *
 * @param A
 * matrix on which to perform the Cholesky decomposition
 */
public void factor(SparseMatrixNd A) {
    n = A.rowSize();
    LU = new SparseMatrixNd(A);
    for (int k = 0; k < n; k++) {
        SparseMatrixCell cellik = LU.getCol(k);
        while (cellik != null && cellik.i < k) cellik = cellik.down;
        double akk = cellik.value;
        cellik = cellik.down;
        while (cellik != null) {
            int i = cellik.i;
            cellik.value /= akk;
            double aik = cellik.value;
            SparseMatrixCell cellij = LU.getRow(i);
            SparseMatrixCell cellkj = LU.getRow(k);
            while (cellij != null && cellij.j <= k) cellij = cellij.next;
            while (cellkj != null && cellkj.j <= k) cellkj = cellkj.next;
            while (cellij != null && cellkj != null) {
                if (cellij.j == cellkj.j)
                    cellij.value -= aik * cellkj.value;
                if (cellij.j < cellkj.j)
                    cellij = cellij.next;
                else
                    cellkj = cellkj.next;
            }
            cellik = cellik.down;
        }
    }
}

Example 2 with SparseMatrixCell

use of maspack.matrix.SparseMatrixCell in project artisynth_core by artisynth.

the class IncompleteLUDecomposition method factor.

public void factor(SparseMatrixNd A, double droptol) {
    n = A.rowSize();
    L = new SparseMatrixNd(n, n);
    U = new SparseMatrixNd(n, n);
    for (int i = 0; i < n; i++) {
        SparseVectorNd w = new SparseVectorNd(n);
        SparseVectorCell cellprev = null;
        SparseMatrixCell cellik = A.getRow(i);
        double itol = 0;
        while (cellik != null) {
            SparseVectorCell wk = new SparseVectorCell(cellik.j, cellik.value);
            w.addEntry(wk, cellprev);
            cellprev = wk;
            // w.set(cellik.j, cellik.value);
            itol += cellik.value * cellik.value;
            cellik = cellik.next;
        }
        itol = Math.sqrt(itol) * droptol;
        // cellik = A.rows[i];
        // while(cellik != null)
        // {
        // System.out.println(cellik.j + " " + cellik.value + " " +
        // w.get(cellik.j));
        // cellik = cellik.next;
        // }
        cellprev = null;
        SparseVectorCell wk = w.getCells();
        while (wk != null && wk.i < i) {
            int k = wk.i;
            wk.value = wk.value / U.get(k, k);
            if (Math.abs(wk.value) > itol) {
                SparseMatrixCell uk = U.getRow(k);
                while (uk != null && uk.j <= k) uk = uk.next;
                while (uk != null) {
                    // while(wj != null && wj.i < uk.j && wj.next != null &&
                    // wj.next.i <= uk.j)
                    // wj = wj.next;
                    // 
                    // if(wj.i == uk.j)
                    // wj.value = wj.value - wk.value * uk.value;
                    // else
                    // {
                    // w.set(uk.j, -wk.value * uk.value);
                    // w.addEntry(new SparseVectorCell(uk.j, -wk.value *
                    // uk.value), wj);
                    // }
                    double newval = w.get(uk.j) - wk.value * uk.value;
                    w.set(uk.j, newval);
                    // System.out.println(newval);
                    uk = uk.next;
                }
                cellprev = wk;
                wk = wk.next();
            } else {
                // System.out.println("dropping " + itol);
                w.set(wk.i, 0);
                // w.removeEntry(wk, cellprev);
                if (cellprev != null)
                    wk = cellprev.next();
                else
                    wk = w.getCells();
            }
        }
        cellprev = null;
        SparseMatrixCell mcellprev = null;
        wk = w.getCells();
        while (wk != null) {
            if (Math.abs(wk.value) > itol || wk.i == i) {
                SparseMatrixCell newcell = new SparseMatrixCell(i, wk.i, wk.value);
                if (wk.i == i)
                    mcellprev = null;
                if (wk.i < i) {
                    // L.set(i, wk.i, wk.value);
                    L.addEntry(newcell, mcellprev, L.prevColEntry(i, wk.i));
                } else {
                    // U.set(i, wk.i, wk.value);
                    U.addEntry(newcell, mcellprev, U.prevColEntry(i, wk.i));
                }
                mcellprev = newcell;
            }
            wk = wk.next();
        }
        L.set(i, i, 1.0);
    }
}

Example 3 with SparseMatrixCell

use of maspack.matrix.SparseMatrixCell in project artisynth_core by artisynth.

the class IncompleteCholeskyDecomposition method factor.

public void factor(SparseMatrixNd A, double tol) {
    n = A.rowSize();
    // modified ILUT as per matlab manual
    // L.set(A);
    // 
    // //slow but functional version
    // VectorNd vec = new VectorNd(n);
    // for(int i = 0; i < n; i++)
    // {
    // L.getRow(i, vec);
    // double toli = vec.norm() * tol;
    // for(int k = 0; k < i; k++)
    // {
    // double aik = vec.get(k) / L.get(k,k);
    // vec.set(k, 0);
    // if(Math.abs(aik) > toli)
    // {
    // for(int j = k+1; j < n; j++)
    // {
    // double aij = vec.get(j) - aik * L.get(k, j);
    // vec.set(j, aij);
    // }
    // }
    // }
    // 
    // for(int j = 0; j < n; j++)
    // {
    // if(Math.abs(vec.get(j)) < toli)
    // {
    // vec.set(j, 0);
    // }
    // }
    // 
    // L.setRow(i, vec);
    // }
    // 
    // SparseMatrixNd goodL = new SparseMatrixNd(L);
    // //update L to get incomplete cholesky
    // VectorNd vec = new VectorNd(n);
    // for(int i = 0; i < n; i++)
    // {
    // L.getRow(i, vec);
    // double diag = vec.get(i);
    // System.out.println(diag);
    // vec.scale(1.0/Math.sqrt(diag));
    // L.setRow(i, vec);
    // }
    C = new SparseMatrixNd(A);
    for (int i = 0; i < n; i++) {
        SparseMatrixCell cellik = C.getRow(i);
        double itol = 0;
        while (cellik != null) {
            itol += cellik.value * cellik.value;
            cellik = cellik.next;
        }
        itol = Math.sqrt(itol) * tol;
        // System.out.println("first itol " + itol);
        cellik = C.getRow(i);
        while (cellik != null && cellik.j < i) {
            int k = cellik.j;
            double aik = cellik.value / A.get(k, k);
            if (Math.abs(aik) > itol) {
                SparseMatrixCell cellkj = C.getRow(k);
                while (cellkj != null && cellkj.j <= k) cellkj = cellkj.next;
                SparseMatrixCell cellij = cellik;
                while (cellkj != null) {
                    int j = cellkj.j;
                    while (cellij != null && cellij.j < j && cellij.next != null && cellij.next.j <= j) cellij = cellij.next;
                    double diff = aik * cellkj.value;
                    if (cellij != null && cellij.j == j) {
                        cellij.value -= diff;
                    } else {
                        SparseMatrixCell cellnew = new SparseMatrixCell(i, j, -diff);
                        C.addEntry(cellnew, cellij, C.prevColEntry(i, j));
                    }
                    cellkj = cellkj.next;
                }
            }
            C.removeEntry(cellik, null, i == 0 ? null : C.getRow(k));
            cellik = cellik.next;
        // L.set(i, k, 0);
        }
        int maxj = i;
        double max = 0;
        cellik = C.getRow(i);
        itol = 0;
        while (cellik != null) {
            double tmp = Math.abs(cellik.value);
            if (tmp > max) {
                maxj = cellik.j;
                max = tmp;
            }
            itol += cellik.value * cellik.value;
            cellik = cellik.next;
        }
        itol = Math.sqrt(itol) * tol;
        // System.out.println(maxj);
        // System.out.println("second itol " + itol);
        SparseMatrixCell lastcellik = C.getRow(i);
        // if(lastcellik.value <= 0)
        // lastcellik.value = Math.sqrt(itol);
        cellik = lastcellik.next;
        while (cellik != null) {
            if (Math.abs(cellik.value) < itol && cellik.j != maxj) {
                C.removeEntry(cellik, lastcellik, C.prevColEntry(i, cellik.j));
            } else {
                lastcellik = cellik;
            }
            cellik = cellik.next;
        }
    }
    for (int i = 0; i < n; i++) {
        SparseMatrixCell cellij = C.getRow(i);
        double diagdiv = 1.0 / Math.sqrt(cellij.value);
        while (cellij != null) {
            cellij.value *= diagdiv;
            cellij = cellij.next;
        }
    }
    C.transpose();
// System.out.println("Cholesky L");
// System.out.println(new MatrixNd(this.L));
}

Example 4 with SparseMatrixCell

use of maspack.matrix.SparseMatrixCell in project artisynth_core by artisynth.

the class IncompleteCholeskyDecomposition method factor.

/**
 * Peforms an Cholesky decomposition on the Matrix M.
 *
 * @param A
 * matrix on which to perform the Cholesky decomposition
 * @throws ImproperSizeException
 * if M is not square
 * @throws IllegalArgumentException
 * if M is detected to be not symmetric positive definite
 */
public void factor(SparseMatrixNd A) {
    // int n = A.rowSize();
    // L = new SparseMatrixNd(n, n);
    // 
    // for(int i = 0; i < n; i++)
    // {
    // double v;
    // 
    // SparseMatrixCell icell = A.rows[i];
    // 
    // SparseMatrixCell tolcell = icell;
    // double inorm = 0;
    // while(tolcell != null)
    // {
    // inorm += tolcell.value * tolcell.value;
    // tolcell = tolcell.next;
    // }
    // inorm = Math.sqrt(inorm);
    // 
    // //compute the diagonal offset parts of the row
    // while(icell.j < i || icell != null)
    // {
    // v = icell.value;
    // // System.out.println("initial value " + v);
    // 
    // SparseMatrixCell ikcell = L.rows[i];
    // SparseMatrixCell jkcell = L.rows[icell.j];
    // 
    // while(ikcell != null && ikcell.j < icell.j && jkcell.j < icell.j)
    // {
    // 
    // if(ikcell.j == jkcell.j)
    // {
    // v -= ikcell.value * jkcell.value;
    // // System.out.println("multiplying values " + ikcell.value + " " +
    // jkcell.value);
    // }
    // if(ikcell.j < jkcell.j)
    // {
    // ikcell = ikcell.next;
    // }
    // else
    // {
    // jkcell = jkcell.next;
    // }
    // }
    // 
    // while(jkcell.j != icell.j)
    // jkcell = jkcell.next;
    // 
    // v /= jkcell.value;
    // 
    // L.set(i, icell.j, v);
    // icell = icell.next;
    // }
    // 
    // //compute diagonal
    // v = icell.value;
    // SparseMatrixCell ikcell = L.rows[i];
    // while(ikcell != null && ikcell.j < i)
    // {
    // v -= ikcell.value * ikcell.value;
    // ikcell = ikcell.next;
    // }
    // 
    // if(v <= 0)
    // v = inorm * 0.001;//A.get(i, i);
    // 
    // v = Math.sqrt(v);
    // 
    // if(Double.isNaN(v) || v == 0)
    // {
    // System.out.println("NaN created");
    // System.out.println(v);
    // }
    // 
    // L.set(i, i, v);
    // }
    n = A.rowSize();
    C = new SparseMatrixNd(A);
    for (int i = 0; i < n; i++) {
        SparseMatrixCell cellij = C.getRow(i);
        while (cellij != null && cellij.j < i) {
            double vij = cellij.value;
            int j = cellij.j;
            SparseMatrixCell cellik = C.getRow(i);
            SparseMatrixCell celljk = C.getRow(j);
            while (cellik != null && cellik.j < cellij.j && celljk != null && celljk.j < cellij.j) {
                if (cellik.j == celljk.j) {
                    vij -= cellik.value * celljk.value;
                }
                if (cellik.j < celljk.j) {
                    cellik = cellik.next;
                } else {
                    celljk = celljk.next;
                }
            }
            while (celljk.j != j) {
                celljk = celljk.next;
            }
            vij /= celljk.value;
            if (Double.isNaN(vij)) {
                System.out.println("NaN " + celljk.value);
            }
            cellij.value = vij;
            cellij = cellij.next;
        }
        double vii = cellij.value;
        SparseMatrixCell cellik = C.getRow(i);
        while (cellik != null && cellik.j < i) {
            vii -= cellik.value * cellik.value;
            cellik = cellik.next;
        }
        if (vii <= 0) {
            vii = cellij.value;
        }
        System.out.println(vii);
        vii = Math.sqrt(vii);
        cellij.value = vii;
        SparseMatrixCell cellii = cellij;
        cellij = cellij.next;
        while (cellij != null) {
            C.removeEntry(cellij, cellii, null);
            cellij = cellij.next;
        }
    }
}

Example 5 with SparseMatrixCell

use of maspack.matrix.SparseMatrixCell in project artisynth_core by artisynth.

the class IncompleteCholeskyDecomposition method solveTranspose.

// public boolean solveDense(Vector x, Vector b) throws
// ImproperStateException, ImproperSizeException
// {
// for(int i = 0; i < n; i++)
// {
// double v = b.get(i);
// 
// for(int j = 0; j < i; j++)
// {
// v -= x.get(j) * L.get(i, j);
// }
// 
// v /= L.get(i, i);
// 
// x.set(i, v);
// }
// 
// return true;
// }
/**
 * Solves L' * x = b
 */
public boolean solveTranspose(VectorNd x, VectorNd b) throws ImproperStateException, ImproperSizeException {
    double[] xbuf = x.getBuffer();
    double[] bbuf = b.getBuffer();
    for (int i = (n - 1); i >= 0; i--) {
        double v = bbuf[i];
        SparseMatrixCell jicell = C.getCol(i);
        jicell = jicell.down;
        while (jicell != null) {
            v -= xbuf[jicell.i] * jicell.value;
            jicell = jicell.down;
        }
        v /= C.get(i, i);
        xbuf[i] = v;
    }
    return true;
}