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

Example 11 with SparseMatrixNd

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

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

the class IncompleteCholeskyDecompositionTest method main.

public static void main(String[] args) {
    SparseMatrixNd A = new SparseMatrixNd(1, 1);
    Random random = new Random(123);
    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 = 70;
    // 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();
    IncompleteCholeskyDecomposition icd = new IncompleteCholeskyDecomposition();
    icd.factor(A, 0.01);
    if (icd.C.containsNaN()) {
        System.out.println("L contains NaN");
    }
    System.out.println("factored matrix");
    // waitforuser();
    SparseMatrixNd tmp = new SparseMatrixNd(n, n);
    tmp.mulTransposeRight(icd.C, icd.C);
    // 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 = new
    // PrintWriter("/ubc/ece/home/hct/other/elliote/A.matrix");
    // fwr.append("[ ");
    // A.write(fwr, new NumberFormat(), WriteFormat.Sparse);
    // fwr.append(" ]");
    // fwr.close();
    // 
    // fwr = new PrintWriter("/ubc/ece/home/hct/other/elliote/L.matrix");
    // fwr.append("[ ");
    // icd.L.write(fwr, new NumberFormat(), WriteFormat.Sparse);
    // 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");
    icd.solve(x, b);
    dsolver.analyzeAndFactor(icd.C);
    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 L' * x = b");
    icd.solveTranspose(x, b);
    tmp.transpose(icd.C);
    dsolver.analyzeAndFactor(tmp);
    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-20;
    int maxit = 500;
    System.out.println("preconditioned solve untransformed");
    x.setZero();
    isolver.solve(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("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 13 with SparseMatrixNd

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

the class CGSolverTest method test.

public void test() {
    // use a small matrix size for now
    int size = 20;
    double eps = 1e-9;
    CGSolver solver = new CGSolver();
    LUDecomposition LU = new LUDecomposition();
    MatrixNd M = new MatrixNd(size, size);
    VectorNd x = new VectorNd(size);
    VectorNd b = new VectorNd(size);
    VectorNd xcheck = new VectorNd(size);
    // identity matrix to try as trivial pre-conditioner
    MatrixNd I = new MatrixNd(size, size);
    I.setIdentity();
    // inverse matrix to try as pre-conditioner
    MatrixNd Minv = new MatrixNd(size, size);
    for (int i = 0; i < 10; i++) {
        int numIter;
        M.setRandom(-0.5, 0.5, randGen);
        // make sure matrix is SPD
        M.mulTranspose(M);
        b.setRandom(-0.5, 0.5, randGen);
        x.setZero();
        if (!solver.solve(x, M, b, eps, size * size)) {
            throw new TestException("No convergence: error=" + solver.getRelativeResidual());
        }
        numIter = solver.getNumIterations();
        LU.factor(M);
        LU.solve(xcheck, b);
        if (!xcheck.epsilonEquals(x, x.norm() * eps)) {
            throw new TestException("Solver gave wrong answer. Expected\n" + xcheck.toString("%8.3f") + "\nGot\n" + x.toString("%8.3f"));
        }
        // System.out.println ("condEst=" + LU.conditionEstimate (M));
        SparseMatrixNd Sc = createBlockDiagonal(10);
        VectorNd bc = new VectorNd(Sc.colSize());
        VectorNd xc = new VectorNd(Sc.colSize());
        bc.setRandom(-0.5, 0.5, randGen);
        if (!solver.solve(xc, Sc, bc, eps, 10 * xc.size())) {
            throw new TestException("No convergence, constraint problem: error=" + solver.getRelativeResidual());
        }
        x.setZero();
        if (!solver.solve(x, M, b, eps, size * size, I)) {
            throw new TestException("No convergence, identity preconditioner: error=" + solver.getRelativeResidual());
        }
        if (numIter != solver.getNumIterations()) {
            throw new TestException("Different iteration count with identity preconditioner: " + solver.getNumIterations() + " vs. " + numIter);
        }
        LU.inverse(Minv);
        x.setZero();
        if (!solver.solve(x, M, b, eps, size * size, Minv)) {
            throw new TestException("No convergence, inverse preconditioner: error=" + solver.getRelativeResidual());
        }
        if (solver.getNumIterations() > 2) {
            throw new TestException("Num iterations > 2 with inverse preconditioner");
        }
    // System.out.println ("num iter=" + solver.getNumIterations());
    // try
    // { PrintWriter pw =
    // new PrintWriter (new FileWriter ("mat.m"));
    // pw.println ("M=[\n" + Sc.toString ("%12.9f") + "]");
    // pw.println ("b=[\n" + bc.toString ("%12.9f") + "]");
    // pw.println ("x=[\n" + xc.toString ("%12.9f") + "]");
    // pw.close();
    // }
    // catch (Exception e)
    // {
    // }
    }
}

Example 14 with SparseMatrixNd

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

the class CGSolverTest method createBlockDiagonal.

public SparseMatrixNd createBlockDiagonal(int nblocks) {
    SparseMatrixNd M = new SparseMatrixNd(5 * nblocks, 5 * nblocks);
    MatrixNd block = new MatrixNd(5, 5);
    for (int k = 0; k < nblocks; k++) {
        block.setRandom(-0.5, 0.5, randGen);
        block.mulTranspose(block);
        for (int i = 0; i < 5; i++) {
            for (int j = 0; j < 5; j++) {
                M.set(5 * k + i, 5 * k + j, block.get(i, j));
            }
        }
    }
    return M;
}

Example 15 with SparseMatrixNd

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

the class CRSolverTest method createBlockDiagonal.

public SparseMatrixNd createBlockDiagonal(int nblocks) {
    SparseMatrixNd M = new SparseMatrixNd(5 * nblocks, 5 * nblocks);
    MatrixNd block = new MatrixNd(5, 5);
    for (int k = 0; k < nblocks; k++) {
        block.setRandom(-0.5, 0.5, randGen);
        block.mulTranspose(block);
        for (int i = 0; i < 5; i++) {
            for (int j = 0; j < 5; j++) {
                M.set(5 * k + i, 5 * k + j, block.get(i, j));
            }
        }
    }
    return M;
}