qrwyk.go

// Copyright (c) Harri Rautila, 2013

// This file is part of github.com/hrautila/matops package. It is free software,
// distributed under the terms of GNU Lesser General Public License Version 3, or
// any later version. See the COPYING tile included in this archive.

package matops

import (
    "github.com/hrautila/matrix"
    "errors"
    //"fmt"
)

// Build Q in place by applying elementary reflectors in reverse order to
// an implied identity matrix.  This forms Q = H(1)H(2) ... H(k)
// 
// this is compatibe with lapack.DORG2R
func unblockedBuildQ(A, tau, w *matrix.FloatMatrix, kb int) error {
    var err error = nil
    var ATL, ATR, ABL, ABR matrix.FloatMatrix
    var A00, a01, A02, a10t, a11, a12t, A20, a21, A22 matrix.FloatMatrix
    var tT, tB matrix.FloatMatrix
    var t0, tau1, t2, w1 matrix.FloatMatrix
    var mb int
    var rowvec bool

    mb  = A.Rows() - A.Cols()
    rowvec =  tau.Rows() == 1

    partition2x2(
        &ATL, &ATR,
        &ABL, &ABR,  A, mb, 0, pBOTTOMRIGHT)

    if rowvec {
        partition1x2(
            &tT, &tB,    tau, 0, pRIGHT)
    } else {
        partition2x1(
            &tT,
            &tB,    tau, 0, pBOTTOM)
    }

    // clearing of the columns of the right and setting ABR to unit diagonal
    // (only if not applying all reflectors, kb > 0)

    for ATL.Rows() > 0 && ATL.Cols() > 0 {
        repartition2x2to3x3(&ATL,
            &A00,  &a01, &A02,
            &a10t, &a11, &a12t,
            &A20,  &a21, &A22,   A, 1, pTOPLEFT)
        if rowvec {
            repartition1x2to1x3(&tT,
                &t0,  &tau1, &t2,     tau, 1, pLEFT)
        } else {
            repartition2x1to3x1(&tT,
                &t0,
                &tau1,
                &t2,     tau, 1, pTOP)
        }

        // --------------------------------------------------------
        
        // adjust workspace to correct size
        w.SubMatrix(&w1, 0, 0, 1, a12t.Cols())  
        // apply Householder reflection from left
        applyHHTo2x1(&tau1, &a21, &a12t, &A22, &w1, LEFT)

        // apply (in-place) current elementary reflector to unit vector
        a21.Scale(-tau1.Float())
        a11.SetAt(0, 0, 1.0 - tau1.Float())

        // zero the upper part
        a01.SetIndexes(0.0)
        
        // --------------------------------------------------------
        continue3x3to2x2(
            &ATL, &ATR,
            &ABL, &ABR,   &A00, &a11, &A22,   A, pTOPLEFT)
        if rowvec {
            continue1x3to1x2(
                &tT, &tB,    &t0, &tau1,   tau, pLEFT)
        } else {
            continue3x1to2x1(
                &tT,
                &tB,    &t0, &tau1,   tau, pTOP)
        }
    }
    return err
}

func blockedBuildQ(A, tau, W *matrix.FloatMatrix, nb int) error {
    var err error = nil
    var ATL, ATR, ABL, ABR, AL matrix.FloatMatrix
    var A00, A01, A02, A10, A11, A12, A20, A21, A22 matrix.FloatMatrix
    var tT, tB matrix.FloatMatrix
    var t0, tau1, t2, Tw, Wrk matrix.FloatMatrix
    var mb int

    mb  = A.Rows() - A.Cols()
    Twork := matrix.FloatZeros(nb, nb)

    partition2x2(
        &ATL, &ATR,
        &ABL, &ABR,  A, mb, 0, pBOTTOMRIGHT)
    partition2x1(
        &tT,
        &tB,    tau, 0, pBOTTOM)

    // clearing of the columns of the right and setting ABR to unit diagonal
    // (only if not applying all reflectors, kb > 0)

    for ATL.Rows() > 0 && ATL.Cols() > 0 {
        repartition2x2to3x3(&ATL,
            &A00,  &A01, &A02,
            &A10,  &A11, &A12,
            &A20,  &A21, &A22,   A, nb, pTOPLEFT)
        repartition2x1to3x1(&tT,
            &t0,
            &tau1,
            &t2,     tau, nb, pTOP)

        // --------------------------------------------------------
        
        // build block reflector from current block
        merge2x1(&AL, &A11, &A21)
        Twork.SubMatrix(&Tw, 0, 0, A11.Cols(), A11.Cols())
        unblkQRBlockReflector(&Tw, &AL, &tau1)
        
        // update with current block reflector (I - Y*T*Y.T)*Atrailing
        W.SubMatrix(&Wrk, 0, 0, A12.Cols(), A11.Cols())
        updateWithQT(&A12, &A22, &A11, &A21, &Tw, &Wrk, nb, false)
        
        // use unblocked version to compute current block
        W.SubMatrix(&Wrk, 0, 0, 1, A11.Cols())
        unblockedBuildQ(&AL, &tau1, &Wrk, 0)
        
        // zero upper part
        A01.SetIndexes(0.0)

        // --------------------------------------------------------
        continue3x3to2x2(
            &ATL, &ATR,
            &ABL, &ABR,   &A00, &A11, &A22,   A, pTOPLEFT)
        continue3x1to2x1(
            &tT,
            &tB,    &t0, &tau1,   tau, pTOP)
    }
    return err
}

func blockedBuildQT(A, T, W *matrix.FloatMatrix, nb int) error {
    var err error = nil
    var ATL, ATR, ABL, ABR, AL matrix.FloatMatrix
    var A00, A01, A11, A12, A21, A22 matrix.FloatMatrix
    var TTL, TTR, TBL, TBR matrix.FloatMatrix
    var T00, T01, T02, T11, T12, T22 matrix.FloatMatrix
    var tau1, Wrk matrix.FloatMatrix
    var mb int

    mb  = A.Rows() - A.Cols()

    partition2x2(
        &ATL, &ATR,
        &ABL, &ABR,  A, mb, 0, pBOTTOMRIGHT)
    partition2x2(
        &TTL, &TTR,
        &TBL, &TBR,  T, 0, 0, pBOTTOMRIGHT)

    // clearing of the columns of the right and setting ABR to unit diagonal
    // (only if not applying all reflectors, kb > 0)

    for ATL.Rows() > 0 && ATL.Cols() > 0 {
        repartition2x2to3x3(&ATL,
            &A00, &A01, nil,
            nil,  &A11, &A12,
            nil,  &A21, &A22,   A, nb, pTOPLEFT)
        repartition2x2to3x3(&TTL,
            &T00, &T01, &T02,
            nil,  &T11, &T12,
            nil,  nil,  &T22,   T, nb, pTOPLEFT)

        // --------------------------------------------------------
        
        // update with current block reflector (I - Y*T*Y.T)*Atrailing
        W.SubMatrix(&Wrk, 0, 0, A12.Cols(), A11.Cols())
        updateWithQT(&A12, &A22, &A11, &A21, &T11, &Wrk, nb, false)
        
        // use unblocked version to compute current block
        W.SubMatrix(&Wrk, 0, 0, 1, A11.Cols())
        // elementary scalar coefficients on the diagonal, column vector
        T11.Diag(&tau1)
        merge2x1(&AL, &A11, &A21)
        // do an unblocked update to current block
        unblockedBuildQ(&AL, &tau1, &Wrk, 0)
        
        // zero upper part
        A01.SetIndexes(0.0)
        // --------------------------------------------------------
        continue3x3to2x2(
            &ATL, &ATR,
            &ABL, &ABR,   &A00, &A11, &A22,   A, pTOPLEFT)
        continue3x3to2x2(
            &TTL, &TTR,
            &TBL, &TBR,   &T00, &T11, &T22,   T, pTOPLEFT)
    }
    return err
}

/*
 * Generate an M-by-N real matrix Q with ortonormal columns, which is
 * defined as the product of k elementary reflectors and block reflector T
 *
 *    Q = H(1) H(2) . . . H(k)
 *
 * as returned by DecomposeQRT().
 *
 * Arguments:
 *  A     On entry, QR factorization as returned by DecomposeQRT() where the lower
 *        trapezoidal  part holds the elementary reflectors. On exit, the M-by-N
 *        matrix Q.
 *
 *  tau   The scalar factors of elementary reflectors as returned by   DecomposeQR()
 *
 *  W     Workspace, size A.Cols()-by-nb.
 *
 *  nb    Blocksize for blocked invocations. If nb == 0 unblocked algorith is used
 *
 * Compatible with lapack.DORGQR
 */
func BuildQ(A, tau, W *matrix.FloatMatrix, nb int) (*matrix.FloatMatrix, error) {
    var err error = nil
    if nb != 0 && W == nil {
        return nil, errors.New("workspace not defined")
    }
    // default is from LEFT
    if nb != 0 && (W.Cols() < nb || W.Rows() < A.Cols()) {
        return nil, errors.New("workspace too small")
    }

    if nb == 0 {
        w := matrix.FloatZeros(1, A.Cols())
        err = unblockedBuildQ(A, tau, w, 0)
    } else {
        var Wrk matrix.FloatMatrix
        Wrk.SubMatrixOf(W, 0, 0, A.Cols(), nb)
        err = blockedBuildQ(A, tau, &Wrk, nb)
    }
    return A, err
}


/*
 * Generate an M-by-N real matrix Q with ortonormal columns, which is
 * defined as the product of k elementary reflectors and block reflector T
 *
 *    Q = H(1) H(2) . . . H(k)
 *
 * generated using the compact WY representaion as returned by DecomposeQRT().
 *
 * Arguments:
 *  A     On entry, QR factorization as returned by DecomposeQRT() where the lower
 *        trapezoidal  part holds the elementary reflectors. On exit, the M-by-N
 *        matrix Q.
 *
 *  T     The block reflector computed from elementary reflectors as returned by
 *        DecomposeQRT() or computed from elementary reflectors and scalar coefficients
 *        by BuildT()
 *
 *  W     Workspace, size A.Cols()-by-nb.
 *
 *  nb    Blocksize for blocked invocations. If nb == 0 default value T.Cols() 
 *        is used.
 *
 * Compatible with lapack.DORGQR
 */
func BuildQT(A, T, W *matrix.FloatMatrix, nb int) (*matrix.FloatMatrix, error) {
    var err error = nil
    if nb == 0 {
        nb = A.Cols()
    }
    // default is from LEFT
    if nb != 0 && (W.Cols() < nb || W.Rows() < A.Cols()) {
        return nil, errors.New("workspace too small")
    }

    var Wrk matrix.FloatMatrix
    Wrk.SubMatrixOf(W, 0, 0, A.Cols(), nb)
    err = blockedBuildQT(A, T, &Wrk, nb)
    return A, err
}


// Local Variables:
// tab-width: 4
// indent-tabs-mode: nil
// End: