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

Example 6 with SparseMatrixNd

use of maspack.matrix.SparseMatrixNd 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 7 with SparseMatrixNd

use of maspack.matrix.SparseMatrixNd 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 8 with SparseMatrixNd

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

the class IncompleteLUDecompositionTest method main.

public static void main(String[] args) {
    SparseMatrixNd A = new SparseMatrixNd(1, 1);
    Random random = new Random(51);
    try {
        A.scan(new ReaderTokenizer(new FileReader("/ubc/ece/home/hct/other/elliote/A.matrix")));
        System.out.println("read matrix " + A.rowSize() + " " + A.colSize());
    } catch (IOException e) {
        System.out.println(e.getMessage());
        System.exit(0);
    }
    int n = A.rowSize();
    System.out.println("loaded matrix");
    // int n = 400;
    // A.setSize(n, n);
    // SVDecomposition svd;
    // do
    // {
    // A.setRandom(1, 2, n*2, random);
    // for(int i = 0; i < n; i++)
    // A.set(i, i, random.nextDouble() + 1);
    // 
    // svd = new SVDecomposition(A);
    // }
    // while(svd.determinant() < 1e-10);
    // A.mulTranspose(A);
    // System.out.println("created spd matrix");
    // waitforuser();
    IncompleteLUDecomposition ilud = new IncompleteLUDecomposition();
    ilud.factor(A, 0.00001);
    if (ilud.L.containsNaN() || ilud.U.containsNaN()) {
        System.out.println("LU contains NaN");
    }
    System.out.println("factored matrix");
    // waitforuser();
    SparseMatrixNd tmp = new SparseMatrixNd(n, n);
    tmp.mul(ilud.L, ilud.U);
    // System.out.println(new MatrixNd(tmp));
    tmp.sub(A);
    System.out.println("residual matrix one norm " + tmp.oneNorm());
    VectorNd x = new VectorNd(n);
    VectorNd b = new VectorNd(n);
    VectorNd xc = new VectorNd(n);
    try {
        System.out.println("writing solve matrix and incomplete cholesky decomposition");
        PrintWriter fwr;
        fwr = new PrintWriter("/ubc/ece/home/hct/other/elliote/LUA.matrix");
        // fwr.append("[ ");
        A.write(fwr, new NumberFormat(), WriteFormat.Dense);
        // fwr.append(" ]");
        fwr.close();
        fwr = new PrintWriter("/ubc/ece/home/hct/other/elliote/L.matrix");
        // fwr.append("[ ");
        ilud.L.write(fwr, new NumberFormat(), WriteFormat.Dense);
        // fwr.append(" ]");
        fwr.close();
        fwr = new PrintWriter("/ubc/ece/home/hct/other/elliote/U.matrix");
        // fwr.append("[ ");
        ilud.U.write(fwr, new NumberFormat(), WriteFormat.Dense);
        // fwr.append(" ]");
        fwr.close();
        System.out.println("finished writing");
    } catch (IOException e) {
        System.out.println(e.getMessage());
        System.exit(0);
    }
    // test backsolves
    System.out.println();
    CGSolver isolver = new CGSolver();
    DirectSolver dsolver = new UmfpackSolver();
    // dsolver.factor(A);
    b.setRandom(-1, 1, random);
    System.out.println("solving L * x = b");
    ilud.solveL(x, b);
    dsolver.analyzeAndFactor(ilud.L);
    dsolver.solve(xc, b);
    System.out.println("b " + b);
    System.out.println("x " + x);
    System.out.println("xc " + xc);
    if (!x.epsilonEquals(xc, 1e-6)) {
        System.out.println("backsolve failed");
    }
    System.out.println();
    System.out.println("solving U * x = b");
    ilud.solveU(x, b);
    dsolver.analyzeAndFactor(ilud.U);
    dsolver.solve(xc, b);
    System.out.println("b " + b);
    System.out.println("x " + x);
    System.out.println("xc " + xc);
    if (!x.epsilonEquals(xc, 1e-6)) {
        System.out.println("backsolve failed");
    }
    // test upcg solver
    System.out.println();
    System.out.println("solving A * x = b");
    double tol = 1e-2;
    int maxit = 1500;
    System.out.println("preconditioned solve untransformed");
    x.setZero();
    isolver.solve(x, A, b, tol, maxit, ilud);
    System.out.println("iterations " + isolver.getNumIterations());
    System.out.println("b " + b);
    System.out.println("x " + x);
    System.out.println();
    // System.out.println("preconditioned solve transformed");
    // x.setZero();
    // isolver.solveTransformed(x, A, b, tol, maxit, icd);
    // System.out.println("iterations " + isolver.getNumIterations());
    // System.out.println("b " + b);
    // System.out.println("x " + x);
    // 
    // System.out.println();
    System.out.println("unpreconditioned solve");
    x.setZero();
    isolver.solve(x, A, b, tol, maxit);
    System.out.println("iterations " + isolver.getNumIterations());
    System.out.println("b " + b);
    System.out.println("x " + x);
    System.out.println();
    System.out.println("direct solve");
    x.setZero();
    dsolver.analyzeAndFactor(A);
    dsolver.solve(x, b);
    System.out.println("b " + b);
    System.out.println("x " + x);
}

Example 9 with SparseMatrixNd

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

the class MFreeModel3d method updateNodeMasses.

public void updateNodeMasses(double totalMass, VectorNd massMatrixDiag) {
    if (totalMass <= 0) {
        totalMass = integrateMass();
    }
    if (massMatrixDiag == null) {
        SparseMatrixNd massMatrix = computeConsistentMassMatrix();
        massMatrixDiag = new VectorNd(massMatrix.rowSize());
        for (int i = 0; i < massMatrix.rowSize(); i++) {
            double rowSum = 0;
            for (int j = 0; j < massMatrix.colSize(); j++) {
                rowSum += massMatrix.get(i, j);
            }
            // rowSum += massMatrix.get(i, i);
            massMatrixDiag.set(i, rowSum);
        }
    }
    double mTrace = massMatrixDiag.sum();
    for (FemNode3d node : myNodes) {
        int i = node.getNumber();
        double m = totalMass / mTrace * massMatrixDiag.get(i);
        node.setMass(m);
        node.setMassExplicit(true);
    }
}

Example 10 with SparseMatrixNd

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