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

Example 6 with VectorNd

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

the class CRSolver method solve.

public boolean solve(VectorNd x, LinearTransformNd A, VectorNd b, double tol, int maxIter, LinearTransformNd P) {
    int n = A.colSize();
    VectorNd xkp = new VectorNd(n);
    VectorNd pkp = new VectorNd(n);
    VectorNd rkp = new VectorNd(n);
    VectorNd wkp = new VectorNd(n);
    VectorNd zkp = new VectorNd(n);
    VectorNd xk = new VectorNd(n);
    VectorNd pk = new VectorNd(n);
    VectorNd rk = new VectorNd(n);
    VectorNd wk = new VectorNd(n);
    VectorNd zk = new VectorNd(n);
    xkp.setZero();
    pkp.set(b);
    rkp.set(b);
    A.mul(zkp, rkp);
    A.mul(wkp, pkp);
    double phikp = b.norm();
    double mukp = rkp.dot(zkp);
    double phik = phikp;
    double muk = mukp;
    k = 1;
    while (k < maxIter && phik > tol) {
        // System.out.println("cr iteration " + k + " " + maxIter + " " +
        // phik);
        double alphak = mukp / wkp.normSquared();
        xk.scaledAdd(alphak, pkp, xkp);
        rk.scaledAdd(-alphak, wkp, rk);
        phik = rk.norm();
        A.mul(zk, rk);
        muk = rk.dot(zk);
        double betak = muk / mukp;
        pk.scaledAdd(betak, pkp, rk);
        wk.scaledAdd(betak, wkp, zk);
        k++;
        xkp.set(xk);
        pkp.set(pk);
        rkp.set(rk);
        zkp.set(zk);
        wkp.set(wk);
        phikp = phik;
        mukp = muk;
    }
    x.set(xk);
    phi = phik;
    return true;
}

Example 7 with VectorNd

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

the class CRSolverTest method test.

public void test() {
    // use a small matrix size for now
    int size = 20;
    double eps = 1e-9;
    CRSolver solver = new CRSolver();
    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 8 with VectorNd

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

the class CRSolverTest method timing.

public void timing() {
    CGSolver solver = new CGSolver();
    FunctionTimer timer = new FunctionTimer();
    for (int nbod = 10; nbod <= 50; nbod += 10) {
        SparseMatrixNd Sc = createBlockDiagonal(nbod);
        VectorNd xc = new VectorNd(Sc.rowSize());
        VectorNd bc = new VectorNd(Sc.rowSize());
        bc.setRandom(-0.5, 0.5, randGen);
        int loopCnt = 10;
        boolean solved = true;
        timer.start();
        for (int i = 0; i < loopCnt; i++) {
            xc.setZero();
            if (!solver.solve(xc, Sc, bc, 1e-8, 10 * xc.size())) {
                solved = false;
                break;
            }
        }
        timer.stop();
        if (!solved) {
            System.out.println("" + nbod + " bodies: No convergence");
        } else {
            System.out.println("" + nbod + " bodies: " + timer.resultMsec(loopCnt));
            System.out.println("         iterations=" + solver.getNumIterations() + " error=" + solver.getRelativeResidual());
        }
    }
}

Example 9 with VectorNd

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

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

the class SetHandler method initializeWidget.

public static boolean initializeWidget(LabeledComponentBase widget, Property prop) {
    String name = prop.getName();
    PropertyInfo info = prop.getInfo();
    Class<?> type = info.getValueClass();
    if (widget instanceof LabeledWidget) {
        LabeledWidget lwidget = (LabeledWidget) widget;
        if (textIsEmpty(lwidget.getLabelText())) {
            lwidget.setLabelText(name);
        }
        if (textIsEmpty(lwidget.getToolTipText())) {
            lwidget.setToolTipText(info.getDescription());
        }
    }
    if (widget instanceof LabeledControl) {
        removeOldListeners((LabeledControl) widget);
    }
    try {
        if (String.class.isAssignableFrom(type)) {
            // if String has a range, then use StringSelector
            // otherwise, use a simple StringField
            Range stringRange = prop.getRange();
            if (info.isInheritable() || stringRange == null || !(stringRange instanceof StringRange) || ((StringRange) stringRange).isWildcard()) {
                StringField stringField = (StringField) widget;
                stringField.setColumns(20);
                stringField.addValueChangeListener(new PropChangeListener(prop) {

                    public void valueChange(ValueChangeEvent e) {
                        if (e.getValue() == null || e.getValue().equals("")) {
                            super.valueChange(new ValueChangeEvent(e.getSource(), ""));
                        } else {
                            super.valueChange(e);
                        }
                    }
                });
                stringField.setStretchable(true);
            } else {
                String[] constants = ((StringRange) stringRange).getValidStrings();
                StringSelector selector = (StringSelector) widget;
                selector.setSelections(constants, null);
                selector.addValueChangeListener(new PropChangeListener(prop));
            }
        } else if (type == double.class || type == float.class || type == Double.class || type == Float.class) {
            DoubleField doubleField = (DoubleField) widget;
            Range range = prop.getRange();
            if (range instanceof NumericInterval) {
                doubleField.setRange((NumericInterval) range);
            }
            // }
            if (info.getPrintFormat() != null && formatIsDefault(doubleField)) {
                doubleField.setFormat(info.getPrintFormat());
            }
            GuiUtils.setFixedWidth(doubleField.getTextField(), 100);
            doubleField.addValueChangeListener(new PropChangeListener(prop));
        } else if (type == int.class || type == Integer.class) {
            IntegerField intField = (IntegerField) widget;
            GuiUtils.setFixedWidth(intField.getTextField(), 100);
            intField.addValueChangeListener(new PropChangeListener(prop));
        } else if (type == boolean.class || type == Boolean.class) {
            if (info.isReadOnly()) {
                StringField stringField = (StringField) widget;
                stringField.setColumns(5);
            } else {
                BooleanSelector selector = (BooleanSelector) widget;
                selector.addValueChangeListener(new PropChangeListener(prop));
            }
        } else if (VectorBase.class.isAssignableFrom(type) && info.getDimension() != -1) {
            VectorBase resultVec;
            try {
                resultVec = (VectorBase) type.newInstance();
            } catch (Exception e) {
                throw new InternalErrorException("Error creating no-args instance of " + type);
            }
            if (resultVec instanceof VectorNd) {
                ((VectorNd) resultVec).setSize(info.getDimension());
            }
            VectorField vectorField = (VectorField) widget;
            // from scratch)
            if (vectorField.getVectorSize() != info.getDimension()) {
                vectorField.setVectorSize(info.getDimension());
            } else {
                VectorNd existingValue = vectorField.getVectorValue();
                if (existingValue != null) {
                    resultVec.set(existingValue);
                }
            }
            vectorField.setResultHolder(resultVec);
            vectorField.addValueChangeListener(new PropChangeListener(prop));
            if (info.getPrintFormat() != null && formatIsDefault(vectorField)) {
                vectorField.setFormat(info.getPrintFormat());
            }
            vectorField.setStretchable(true);
        } else if (VectoriBase.class.isAssignableFrom(type) && info.getDimension() != -1) {
            VectoriBase resultVec;
            try {
                resultVec = (VectoriBase) type.newInstance();
            } catch (Exception e) {
                throw new InternalErrorException("Error creating no-args instance of " + type);
            }
            if (resultVec instanceof VectorNi) {
                ((VectorNi) resultVec).setSize(info.getDimension());
            }
            VectoriField vectorField = (VectoriField) widget;
            // from scratch)
            if (vectorField.getVectorSize() != info.getDimension()) {
                vectorField.setVectorSize(info.getDimension());
            } else {
                VectorNi existingValue = vectorField.getVectorValue();
                if (existingValue != null) {
                    resultVec.set(existingValue);
                }
            }
            vectorField.setResultHolder(resultVec);
            vectorField.addValueChangeListener(new PropChangeListener(prop));
            if (info.getPrintFormat() != null && formatIsDefault(vectorField)) {
                vectorField.setFormat(info.getPrintFormat());
            }
            vectorField.setStretchable(true);
        } else if (VectorBase.class.isAssignableFrom(type) && info.getDimension() == -1) {
            VectorBase resultVec;
            try {
                resultVec = (VectorBase) type.newInstance();
            } catch (Exception e) {
                throw new InternalErrorException("Error creating no-args instance of " + type);
            }
            VariableVectorField vectorField = (VariableVectorField) widget;
            VectorNd existingValue = vectorField.getVectorValue();
            if (vectorField.getVectorSize() != resultVec.size()) {
                resultVec.setSize(vectorField.getVectorSize());
            }
            if (existingValue != null) {
                resultVec.set(existingValue);
            }
            vectorField.setResultHolder(resultVec);
            vectorField.addValueChangeListener(new PropChangeListener(prop));
            if (info.getPrintFormat() != null && formatIsDefault(vectorField)) {
                vectorField.setFormat(info.getPrintFormat());
            }
            vectorField.setStretchable(true);
        } else if (SymmetricMatrix3d.class.isAssignableFrom(type)) {
            SymmetricMatrix3dField matrixField = (SymmetricMatrix3dField) widget;
            matrixField.addValueChangeListener(new PropChangeListener(prop));
            if (info.getPrintFormat() != null && formatIsDefault(matrixField)) {
                matrixField.setFormat(info.getPrintFormat());
            }
            matrixField.setStretchable(true);
        } else if (RigidTransform3d.class.isAssignableFrom(type)) {
            RigidTransformWidget transformField = (RigidTransformWidget) widget;
            transformField.addValueChangeListener(new PropChangeListener(prop));
            transformField.setStretchable(true);
        } else if (AffineTransform3d.class.isAssignableFrom(type)) {
            AffineTransformWidget transformField = (AffineTransformWidget) widget;
            transformField.addValueChangeListener(new PropChangeListener(prop));
            transformField.setStretchable(true);
        } else if (Rectangle2d.class.isAssignableFrom(type)) {
            RectangleField rectField = (RectangleField) widget;
            rectField.addValueChangeListener(new PropChangeListener(prop));
            if (info.getPrintFormat() != null && formatIsDefault(rectField)) {
                rectField.setFormat(info.getPrintFormat());
            }
            rectField.setStretchable(true);
        } else // }
        if (AxisAngle.class.isAssignableFrom(type)) {
            AxisAngleField orientationField = (AxisAngleField) widget;
            orientationField.addValueChangeListener(new PropChangeListener(prop));
            orientationField.setStretchable(true);
        } else if (Enum.class.isAssignableFrom(type)) {
            Enum<?>[] constants = null;
            Range range = prop.getRange();
            if (range != null && range instanceof EnumRange) {
                constants = ((EnumRange<?>) range).getValidEnums();
            } else {
                constants = (Enum[]) type.getEnumConstants();
            }
            if (info.isReadOnly()) {
                StringField stringField = (StringField) widget;
                int ncols = 0;
                for (int i = 0; i < constants.length; i++) {
                    int len = constants[i].toString().length();
                    if (len > ncols) {
                        ncols = len;
                    }
                }
                stringField.setColumns(ncols);
            } else {
                EnumSelector selector = (EnumSelector) widget;
                selector.setSelections(constants, null);
                selector.addValueChangeListener(new PropChangeListener(prop));
            }
        } else if (Color.class.isAssignableFrom(type)) {
            ColorSelector selector = (ColorSelector) widget;
            if (info.getNullValueOK()) {
                selector.enableNullColors();
            }
            selector.addValueChangeListener(new PropChangeListener(prop));
        } else if (IntegerInterval.class.isAssignableFrom(type)) {
            IntegerIntervalField rangeField = (IntegerIntervalField) widget;
            rangeField.addValueChangeListener(new PropChangeListener(prop));
            rangeField.setStretchable(true);
        } else if (NumericInterval.class.isAssignableFrom(type)) {
            DoubleIntervalField rangeField = (DoubleIntervalField) widget;
            rangeField.addValueChangeListener(new PropChangeListener(prop));
            rangeField.setStretchable(true);
        } else if (GLGridResolution.class.isAssignableFrom(type)) {
            GridResolutionField resField = (GridResolutionField) widget;
            resField.addValueChangeListener(new PropChangeListener(prop));
        } else if (Font.class.isAssignableFrom(type)) {
            FontField fontField = (FontField) widget;
            fontField.addValueChangeListener(new PropChangeListener(prop));
        } else if (CompositeProperty.class.isAssignableFrom(type)) {
            if (widget instanceof CompositePropertyWidget) {
                CompositePropertyWidget compProp = (CompositePropertyWidget) widget;
                compProp.setProperty(prop);
            } else {
                CompositePropertyPanel compProp = (CompositePropertyPanel) widget;
                compProp.setExpandState(info.getWidgetExpandState());
                compProp.initializeSelection(prop);
            }
        } else {
            return false;
        }
    } catch (ClassCastException e) {
        throw new IllegalArgumentException("widget type " + widget.getClass() + " inappropriate for property type " + type);
    }
    // finishWidget (widget, prop);
    return true;
}