no.uib.cipr.matrix.QR Maven / Gradle / Ivy

Go to download
/*
 * 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 com.github.fommil.netlib.LAPACK;
import org.netlib.util.intW;

/**
 * Computers QR decompositions
 */
public class QR extends OrthogonalComputer {

    /**
     * Constructs an empty QR decomposition
     * 
     * @param m
     *            Number of rows. Must be larger than or equal the number of
     *            columns
     * @param n
     *            Number of columns
     */
    public QR(int m, int n) {
        super(m, n, true);

        if (n > m)
            throw new IllegalArgumentException("n > m");

        int lwork;

        // Query optimal workspace. First for computing the factorization
        {
            work = new double[1];
            intW info = new intW(0);
            LAPACK.getInstance().dgeqrf(m, n, new double[0], Matrices.ld(m),
                    new double[0], work, -1, info);

            if (info.val != 0)
                lwork = n;
            else
                lwork = (int) work[0];
            lwork = Math.max(1, lwork);
            work = new double[lwork];
        }

        // Workspace needed for generating an explicit orthogonal matrix
        {
            workGen = new double[1];
            intW info = new intW(0);
            LAPACK.getInstance().dorgqr(m, n, k, new double[0], Matrices.ld(m),
                    new double[0], workGen, -1, info);

            if (info.val != 0)
                lwork = n;
            else
                lwork = (int) workGen[0];
            lwork = Math.max(1, lwork);
            workGen = new double[lwork];
        }

    }

    /**
     * Convenience method to compute a QR decomposition
     * 
     * @param A
     *            Matrix to decompose. Not modified
     * @return Newly allocated decomposition
     */
    public static QR factorize(Matrix A) {
        return new QR(A.numRows(), A.numColumns()).factor(new DenseMatrix(A));
    }

    @Override
    public QR factor(DenseMatrix A) {

        if (Q.numRows() != A.numRows())
            throw new IllegalArgumentException("Q.numRows() != A.numRows()");
        else if (Q.numColumns() != A.numColumns())
            throw new IllegalArgumentException(
                    "Q.numColumns() != A.numColumns()");
        else if (R == null)
            throw new IllegalArgumentException("R == null");

        /*
         * Calculate factorisation, and extract the triangular factor
         */
        intW info = new intW(0);
        LAPACK.getInstance().dgeqrf(m, n, A.getData(), Matrices.ld(m), tau,
                work, work.length, info);

        if (info.val < 0)
            throw new IllegalArgumentException();

        R.zero();
        for (MatrixEntry e : A)
            if (e.row() <= e.column())
                R.set(e.row(), e.column(), e.get());

        /*
         * Generate the orthogonal matrix
         */
        info.val = 0;
        LAPACK.getInstance().dorgqr(m, n, k, A.getData(), Matrices.ld(m), tau,
                workGen, workGen.length, info);

        if (info.val < 0)
            throw new IllegalArgumentException();

        Q.set(A);

        return this;
    }

    /**
     * Returns the upper triangular factor
     */
    public UpperTriangDenseMatrix getR() {
        return R;
    }
}