/*
Copyright (C) 2008 Richard Gomes

This source code is release under the BSD License.

This file is part of JQuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://jquantlib.org/

JQuantLib is free software: you can redistribute it and/or modify it
under the terms of the JQuantLib license.  You should have received a
copy of the license along with this program; if not, please email
. The license is also available online at
.

This program 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 license for more details.

JQuantLib is based on QuantLib. http://quantlib.org/
When applicable, the original copyright notice follows this notice.
 */

/*
 Copyright (C) 2008 Klaus Spanderen

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 . The license is also available online at
 .

 This program 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 license for more details.
*/

package org.jquantlib.math.matrixutilities;

import java.util.Arrays;

import org.jquantlib.QL;
import org.jquantlib.math.optimization.Minpack;


/**
 * QR decompoisition
 * 

 * This implementation is based on MINPACK/J
 * 

 * This subroutine uses householder transformations with optional column pivoting to compute a QR factorization of
 * the {@latex$ m} by {@latex$ n} matrix {@latex$ A}.
 * 

 * That is, Minpack_f77#qrfac determines an orthogonal matrix {@latex$ Q}, a permutation matrix {@latex$ P}, and
 * an upper trapezoidal matrix {@latex$ R} with diagonal elements of nonincreasing magnitude, such that {@latex$ A*P = Q*R}.
 * 

 * Return value ipvt is an integer array of length {@latex$ n}, which defines the permutation matrix {@latex$ P}
 * such that {@latex$ A*P = Q*R}.
 * 

 * Column {@latex$ j} of {@latex$ P} is column ipvt[j] of the identity matrix.
 * 

 * See lmdiff.cpp for further details.
 *
 * @see MathWorld
 * @see Wikipedia
 * @see MINPACK/J
 * @see MINPACK
 *
 * @author Richard Gomes
 */
public class QRDecomposition {

    private final int m;
    private final int n;
    private final Matrix A;
    private final Matrix Q;
    private final Matrix R;
    private final Matrix P;
    private final int ipvt[];
    private final boolean isNonSingular;


    /**
     * Calculates the QR-decomposition of the given matrix.
     * @param A The matrix to decompose.
     */
    public QRDecomposition(final Matrix matrix) {
        this(matrix, false);
    }


    /**
     * Calculates the QR-decomposition of the given matrix.
     * @param M The matrix to decompose.
     */
    public QRDecomposition(final Matrix A, final boolean pivot) {
        this.A = A;
        this.m = A.rows();
        this.n = A.cols();
        this.ipvt = new int[n];

        final Matrix mT = A.clone().toJava().transpose();
        final Array rdiag = new Array(n);
        final Array wa = new Array(n);

        System.out.println("mT (BEFORE) = " + mT);
        Minpack.qrfac(
                m, n,
                mT,                      // input/output parameter (sorry for that :~ )
                pivot,
                ipvt, rdiag, rdiag, wa  // output parameters (sorry for that :~ )
                );

        System.out.println("mT (AFTER)  = " + mT);
        System.out.println("Array ipvt = " + Arrays.toString(ipvt));
        System.out.println("Array rdiag = " + rdiag);
        System.out.println("Array wa = " + wa);

        // obtain R
        final double[][] r = new double[n][n];
        for (int i=0; i < n; i++) {
            // r[i][i] = rdiag[i];
            r[i][i] = rdiag.get(i);
            if (i < m) {
                // std::copy(mT.column_begin(i)+i+1, mT.column_end(i), r.row_begin(i)+i+1);
                for (int k = i + 1; k < n; k++) {
                    // r[i][k] = mT[k][i];
                    r[i][k] = mT.get(k, i);
                }
            }
        }

        // obtain Q
        final double[][] q = new double[m][n];
        final double w[] = new double[m];
        for (int k=0; k < m; k++) {
            Arrays.fill(w, 0.0);
            w[k] = 1.0;
            for (int j=0; j < Math.min(n, m); j++) {
                final double t3 = mT.get(j, j);
                if (t3 != 0.0) {
                    // final double t = std::inner_product(mT.row_begin(j)+j, mT.row_end(j), w.begin()+j, 0.0)/t3;
                    double t = 0.0;
                    for (int p = j; p < m; p++) {
                        t += mT.get(j, p) * w[p];
                    }
                    t /= t3;

                    for (int i=j; i