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

Example 1 with SparseVectorNd

use of maspack.matrix.SparseVectorNd 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);
    }
}