svd.java

另一个功能更强大的矩阵运算软件开源代码

/*
 * Copyright (C) 2003-2006 Bjørn-Ove Heimsund
 * 
 * This file is part of MTJ.
 * 
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
 * 
 * This library is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
 * for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */

package no.uib.cipr.matrix;

import org.netlib.lapack.LAPACK;
import org.netlib.util.intW;

/**
 * Computes singular value decompositions
 */
public class SVD {

    /**
     * Work array
     */
    private final double[] work;

    /**
     * Work array
     */
    private final int[] iwork;

    /**
     * Matrix dimension
     */
    private final int m, n;

    /**
     * Compute the singular vectors fully?
     */
    private final boolean vectors;

    /**
     * Job to do
     */
    private final JobSVD job;

    /**
     * The singular values
     */
    private final double[] S;

    /**
     * Singular vectors
     */
    private final DenseMatrix U, Vt;

    /**
     * Creates an empty SVD which will compute all singular values and vectors
     * 
     * @param m
     *            Number of rows
     * @param n
     *            Number of columns
     */
    public SVD(int m, int n) {
        this(m, n, true);
    }

    /**
     * Creates an empty SVD
     * 
     * @param m
     *            Number of rows
     * @param n
     *            Number of columns
     * @param vectors
     *            True to compute the singular vectors, false for just the
     *            singular values
     */
    public SVD(int m, int n, boolean vectors) {
        this.m = m;
        this.n = n;
        this.vectors = vectors;

        // Allocate space for the decomposition
        S = new double[Math.min(m, n)];
        if (vectors) {
            U = new DenseMatrix(m, m);
            Vt = new DenseMatrix(n, n);
        } else
            U = Vt = null;

        job = vectors ? JobSVD.All : JobSVD.None;

        // Find workspace requirements
        iwork = new int[8 * Math.min(m, n)];

        // Query optimal workspace
        double[] worksize = new double[1];
        intW info = new intW(0);
        LAPACK.getInstance().dgesdd(job.netlib(), m, n, new double[0],
                Matrices.ld(m), new double[0], new double[0], Matrices.ld(m),
                new double[0], Matrices.ld(n), worksize, -1,
                iwork, info);

        // Allocate workspace
        int lwork = -1;
        if (info.val != 0) {
            if (vectors)
                lwork = 3
                        * Math.min(m, n)
                        * Math.min(m, n)
                        + Math.max(Math.max(m, n), 4 * Math.min(m, n)
                                * Math.min(m, n) + 4 * Math.min(m, n));
            else
                lwork = 3
                        * Math.min(m, n)
                        * Math.min(m, n)
                        + Math.max(Math.max(m, n), 5 * Math.min(m, n)
                                * Math.min(m, n) + 4 * Math.min(m, n));
        } else
            lwork = (int) worksize[0];

        lwork = Math.max(lwork, 1);
        work = new double[lwork];
    }

    /**
     * Convenience method for computing a full SVD
     * 
     * @param A
     *            Matrix to decompose, not modified
     * @return Newly allocated factorization
     * @throws NotConvergedException
     */
    public static SVD factorize(Matrix A) throws NotConvergedException {
        return new SVD(A.numRows(), A.numColumns()).factor(new DenseMatrix(A));
    }

    /**
     * Computes an SVD
     * 
     * @param A
     *            Matrix to decompose. Size must conform, and it will be
     *            overwritten on return. Pass a copy to avoid this
     * @return The current decomposition
     * @throws NotConvergedException
     */
    public SVD factor(DenseMatrix A) throws NotConvergedException {
        if (A.numRows() != m)
            throw new IllegalArgumentException("A.numRows() != m");
        else if (A.numColumns() != n)
            throw new IllegalArgumentException("A.numColumns() != n");

        intW info = new intW(0);
        LAPACK.getInstance().dgesdd(job.netlib(), m, n, A.getData(), Matrices.ld(m), S,
                vectors ? U.getData() : new double[0], Matrices.ld(m),
                vectors ? Vt.getData() : new double[0], Matrices.ld(n), work, work.length,
                iwork, info);

        if (info.val > 0)
            throw new NotConvergedException(
                    NotConvergedException.Reason.Iterations);
        else if (info.val < 0)
            throw new IllegalArgumentException();

        return this;
    }

    /**
     * True if singular vectors are stored
     */
    public boolean hasSingularVectors() {
        return U != null;
    }

    /**
     * Returns the left singular vectors, column-wise. Not available for partial
     * decompositions
     * 
     * @return Matrix of size m*m
     */
    public DenseMatrix getU() {
        return U;
    }

    /**
     * Returns the right singular vectors, row-wise. Not available for partial
     * decompositions
     * 
     * @return Matrix of size n*n
     */
    public DenseMatrix getVt() {
        return Vt;
    }

    /**
     * Returns the singular values (stored in descending order)
     * 
     * @return Array of size min(m,n)
     */
    public double[] getS() {
        return S;
    }

}