Logo Search packages:      
Sourcecode: quantlib version File versions  Download package

Disposable< std::vector< Size > > QuantLib::qrDecomposition ( const Matrix &  A,
Matrix &  q,
Matrix &  r,
bool  pivot = true 
)

QR decompoisition.

This implementation is based on MINPACK (<http://www.netlib.org/minpack>, <http://www.netlib.org/cephes/linalg.tgz>)

This subroutine uses householder transformations with column pivoting (optional) to compute a qr factorization of the m by n matrix A. That is, qrfac determines an orthogonal matrix q, a permutation matrix p, and an upper trapezoidal matrix r with diagonal elements of nonincreasing magnitude, such that A*p = q*r.

Return value ipvt is an integer array of length n, which defines the permutation matrix p such that A*p = q*r. Column j of p is column ipvt(j) of the identity matrix.

See lmdiff.cpp for further details.

Definition at line 29 of file qrdecomposition.cpp.

References QuantLib::Array::begin(), QuantLib::Matrix::begin(), QuantLib::Matrix::column_begin(), QuantLib::Matrix::column_end(), QuantLib::Matrix::columns(), QuantLib::Array::end(), QuantLib::Matrix::row_begin(), QuantLib::Matrix::row_end(), and QuantLib::Matrix::rows().

Referenced by qrSolve().

                                                               {
        Matrix mT = transpose(M);
        const Size m = M.rows();
        const Size n = M.columns();

        boost::scoped_array<int> lipvt(new int[n]);
        boost::scoped_array<double> rdiag(new double[n]);
        boost::scoped_array<double> wa(new double[n]);

        MINPACK::qrfac(m, n, mT.begin(), 0, (pivot)?1:0,
                       lipvt.get(), n, rdiag.get(), rdiag.get(), wa.get());

        if (r.columns() != n || r.rows() !=n)
            r = Matrix(n, n);

        for (Size i=0; i < n; ++i) {
            std::fill(r.row_begin(i), r.row_begin(i)+i, 0.0);
            r[i][i] = rdiag[i];
            if (i < m) {
                std::copy(mT.column_begin(i)+i+1, mT.column_end(i),
                          r.row_begin(i)+i+1);
            }
            else {
                std::fill(r.row_begin(i)+i+1, r.row_end(i), 0.0);
            }
        }

        if (q.rows() != m || q.columns() != n)
            q = Matrix(m, n);

        Array w(m);
        for (Size k=0; k < m; ++k) {
            std::fill(w.begin(), w.end(), 0.0);
            w[k] = 1.0;

            for (Size j=0; j < std::min(n, m); ++j) {
                const Real t3 = mT[j][j];
                if (t3 != 0.0) {
                    const Real t
                        = std::inner_product(mT.row_begin(j)+j, mT.row_end(j),
                                             w.begin()+j, 0.0)/t3;
                    for (Size i=j; i<m; ++i) {
                        w[i]-=mT[j][i]*t;
                    }
                }
                q[k][j] = w[j];
            }
            std::fill(q.row_begin(k) + std::min(n, m), q.row_end(k), 0.0);
        }

        std::vector<Size> ipvt(n);
        if (pivot) {
            std::copy(lipvt.get(), lipvt.get()+n, ipvt.begin());
        }
        else {
            for (Size i=0; i < n; ++i)
                ipvt[i] = i;
        }

        return ipvt;
    }


Generated by  Doxygen 1.6.0   Back to index