| // Copyright ©2015 The Gonum Authors. All rights reserved. |
| // Use of this source code is governed by a BSD-style |
| // license that can be found in the LICENSE file. |
| |
| package gonum |
| |
| import ( |
| "gonum.org/v1/gonum/blas" |
| "gonum.org/v1/gonum/blas/blas64" |
| "gonum.org/v1/gonum/lapack" |
| ) |
| |
| // Dlarfb applies a block reflector to a matrix. |
| // |
| // In the call to Dlarfb, the mxn c is multiplied by the implicitly defined matrix h as follows: |
| // c = h * c if side == Left and trans == NoTrans |
| // c = c * h if side == Right and trans == NoTrans |
| // c = hᵀ * c if side == Left and trans == Trans |
| // c = c * hᵀ if side == Right and trans == Trans |
| // h is a product of elementary reflectors. direct sets the direction of multiplication |
| // h = h_1 * h_2 * ... * h_k if direct == Forward |
| // h = h_k * h_k-1 * ... * h_1 if direct == Backward |
| // The combination of direct and store defines the orientation of the elementary |
| // reflectors. In all cases the ones on the diagonal are implicitly represented. |
| // |
| // If direct == lapack.Forward and store == lapack.ColumnWise |
| // V = [ 1 ] |
| // [v1 1 ] |
| // [v1 v2 1] |
| // [v1 v2 v3] |
| // [v1 v2 v3] |
| // If direct == lapack.Forward and store == lapack.RowWise |
| // V = [ 1 v1 v1 v1 v1] |
| // [ 1 v2 v2 v2] |
| // [ 1 v3 v3] |
| // If direct == lapack.Backward and store == lapack.ColumnWise |
| // V = [v1 v2 v3] |
| // [v1 v2 v3] |
| // [ 1 v2 v3] |
| // [ 1 v3] |
| // [ 1] |
| // If direct == lapack.Backward and store == lapack.RowWise |
| // V = [v1 v1 1 ] |
| // [v2 v2 v2 1 ] |
| // [v3 v3 v3 v3 1] |
| // An elementary reflector can be explicitly constructed by extracting the |
| // corresponding elements of v, placing a 1 where the diagonal would be, and |
| // placing zeros in the remaining elements. |
| // |
| // t is a k×k matrix containing the block reflector, and this function will panic |
| // if t is not of sufficient size. See Dlarft for more information. |
| // |
| // work is a temporary storage matrix with stride ldwork. |
| // work must be of size at least n×k side == Left and m×k if side == Right, and |
| // this function will panic if this size is not met. |
| // |
| // Dlarfb is an internal routine. It is exported for testing purposes. |
| func (Implementation) Dlarfb(side blas.Side, trans blas.Transpose, direct lapack.Direct, store lapack.StoreV, m, n, k int, v []float64, ldv int, t []float64, ldt int, c []float64, ldc int, work []float64, ldwork int) { |
| nv := m |
| if side == blas.Right { |
| nv = n |
| } |
| switch { |
| case side != blas.Left && side != blas.Right: |
| panic(badSide) |
| case trans != blas.Trans && trans != blas.NoTrans: |
| panic(badTrans) |
| case direct != lapack.Forward && direct != lapack.Backward: |
| panic(badDirect) |
| case store != lapack.ColumnWise && store != lapack.RowWise: |
| panic(badStoreV) |
| case m < 0: |
| panic(mLT0) |
| case n < 0: |
| panic(nLT0) |
| case k < 0: |
| panic(kLT0) |
| case store == lapack.ColumnWise && ldv < max(1, k): |
| panic(badLdV) |
| case store == lapack.RowWise && ldv < max(1, nv): |
| panic(badLdV) |
| case ldt < max(1, k): |
| panic(badLdT) |
| case ldc < max(1, n): |
| panic(badLdC) |
| case ldwork < max(1, k): |
| panic(badLdWork) |
| } |
| |
| if m == 0 || n == 0 { |
| return |
| } |
| |
| nw := n |
| if side == blas.Right { |
| nw = m |
| } |
| switch { |
| case store == lapack.ColumnWise && len(v) < (nv-1)*ldv+k: |
| panic(shortV) |
| case store == lapack.RowWise && len(v) < (k-1)*ldv+nv: |
| panic(shortV) |
| case len(t) < (k-1)*ldt+k: |
| panic(shortT) |
| case len(c) < (m-1)*ldc+n: |
| panic(shortC) |
| case len(work) < (nw-1)*ldwork+k: |
| panic(shortWork) |
| } |
| |
| bi := blas64.Implementation() |
| |
| transt := blas.Trans |
| if trans == blas.Trans { |
| transt = blas.NoTrans |
| } |
| // TODO(btracey): This follows the original Lapack code where the |
| // elements are copied into the columns of the working array. The |
| // loops should go in the other direction so the data is written |
| // into the rows of work so the copy is not strided. A bigger change |
| // would be to replace work with workᵀ, but benchmarks would be |
| // needed to see if the change is merited. |
| if store == lapack.ColumnWise { |
| if direct == lapack.Forward { |
| // V1 is the first k rows of C. V2 is the remaining rows. |
| if side == blas.Left { |
| // W = Cᵀ V = C1ᵀ V1 + C2ᵀ V2 (stored in work). |
| |
| // W = C1. |
| for j := 0; j < k; j++ { |
| bi.Dcopy(n, c[j*ldc:], 1, work[j:], ldwork) |
| } |
| // W = W * V1. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, |
| n, k, 1, |
| v, ldv, |
| work, ldwork) |
| if m > k { |
| // W = W + C2ᵀ V2. |
| bi.Dgemm(blas.Trans, blas.NoTrans, n, k, m-k, |
| 1, c[k*ldc:], ldc, v[k*ldv:], ldv, |
| 1, work, ldwork) |
| } |
| // W = W * Tᵀ or W * T. |
| bi.Dtrmm(blas.Right, blas.Upper, transt, blas.NonUnit, n, k, |
| 1, t, ldt, |
| work, ldwork) |
| // C -= V * Wᵀ. |
| if m > k { |
| // C2 -= V2 * Wᵀ. |
| bi.Dgemm(blas.NoTrans, blas.Trans, m-k, n, k, |
| -1, v[k*ldv:], ldv, work, ldwork, |
| 1, c[k*ldc:], ldc) |
| } |
| // W *= V1ᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, n, k, |
| 1, v, ldv, |
| work, ldwork) |
| // C1 -= Wᵀ. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < n; i++ { |
| for j := 0; j < k; j++ { |
| c[j*ldc+i] -= work[i*ldwork+j] |
| } |
| } |
| return |
| } |
| // Form C = C * H or C * Hᵀ, where C = (C1 C2). |
| |
| // W = C1. |
| for i := 0; i < k; i++ { |
| bi.Dcopy(m, c[i:], ldc, work[i:], ldwork) |
| } |
| // W *= V1. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, m, k, |
| 1, v, ldv, |
| work, ldwork) |
| if n > k { |
| bi.Dgemm(blas.NoTrans, blas.NoTrans, m, k, n-k, |
| 1, c[k:], ldc, v[k*ldv:], ldv, |
| 1, work, ldwork) |
| } |
| // W *= T or Tᵀ. |
| bi.Dtrmm(blas.Right, blas.Upper, trans, blas.NonUnit, m, k, |
| 1, t, ldt, |
| work, ldwork) |
| if n > k { |
| bi.Dgemm(blas.NoTrans, blas.Trans, m, n-k, k, |
| -1, work, ldwork, v[k*ldv:], ldv, |
| 1, c[k:], ldc) |
| } |
| // C -= W * Vᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, m, k, |
| 1, v, ldv, |
| work, ldwork) |
| // C -= W. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < m; i++ { |
| for j := 0; j < k; j++ { |
| c[i*ldc+j] -= work[i*ldwork+j] |
| } |
| } |
| return |
| } |
| // V = (V1) |
| // = (V2) (last k rows) |
| // Where V2 is unit upper triangular. |
| if side == blas.Left { |
| // Form H * C or |
| // W = Cᵀ V. |
| |
| // W = C2ᵀ. |
| for j := 0; j < k; j++ { |
| bi.Dcopy(n, c[(m-k+j)*ldc:], 1, work[j:], ldwork) |
| } |
| // W *= V2. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, n, k, |
| 1, v[(m-k)*ldv:], ldv, |
| work, ldwork) |
| if m > k { |
| // W += C1ᵀ * V1. |
| bi.Dgemm(blas.Trans, blas.NoTrans, n, k, m-k, |
| 1, c, ldc, v, ldv, |
| 1, work, ldwork) |
| } |
| // W *= T or Tᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, transt, blas.NonUnit, n, k, |
| 1, t, ldt, |
| work, ldwork) |
| // C -= V * Wᵀ. |
| if m > k { |
| bi.Dgemm(blas.NoTrans, blas.Trans, m-k, n, k, |
| -1, v, ldv, work, ldwork, |
| 1, c, ldc) |
| } |
| // W *= V2ᵀ. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, n, k, |
| 1, v[(m-k)*ldv:], ldv, |
| work, ldwork) |
| // C2 -= Wᵀ. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < n; i++ { |
| for j := 0; j < k; j++ { |
| c[(m-k+j)*ldc+i] -= work[i*ldwork+j] |
| } |
| } |
| return |
| } |
| // Form C * H or C * Hᵀ where C = (C1 C2). |
| // W = C * V. |
| |
| // W = C2. |
| for j := 0; j < k; j++ { |
| bi.Dcopy(m, c[n-k+j:], ldc, work[j:], ldwork) |
| } |
| |
| // W = W * V2. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, m, k, |
| 1, v[(n-k)*ldv:], ldv, |
| work, ldwork) |
| if n > k { |
| bi.Dgemm(blas.NoTrans, blas.NoTrans, m, k, n-k, |
| 1, c, ldc, v, ldv, |
| 1, work, ldwork) |
| } |
| // W *= T or Tᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, trans, blas.NonUnit, m, k, |
| 1, t, ldt, |
| work, ldwork) |
| // C -= W * Vᵀ. |
| if n > k { |
| // C1 -= W * V1ᵀ. |
| bi.Dgemm(blas.NoTrans, blas.Trans, m, n-k, k, |
| -1, work, ldwork, v, ldv, |
| 1, c, ldc) |
| } |
| // W *= V2ᵀ. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, m, k, |
| 1, v[(n-k)*ldv:], ldv, |
| work, ldwork) |
| // C2 -= W. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < m; i++ { |
| for j := 0; j < k; j++ { |
| c[i*ldc+n-k+j] -= work[i*ldwork+j] |
| } |
| } |
| return |
| } |
| // Store = Rowwise. |
| if direct == lapack.Forward { |
| // V = (V1 V2) where v1 is unit upper triangular. |
| if side == blas.Left { |
| // Form H * C or Hᵀ * C where C = (C1; C2). |
| // W = Cᵀ * Vᵀ. |
| |
| // W = C1ᵀ. |
| for j := 0; j < k; j++ { |
| bi.Dcopy(n, c[j*ldc:], 1, work[j:], ldwork) |
| } |
| // W *= V1ᵀ. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, n, k, |
| 1, v, ldv, |
| work, ldwork) |
| if m > k { |
| bi.Dgemm(blas.Trans, blas.Trans, n, k, m-k, |
| 1, c[k*ldc:], ldc, v[k:], ldv, |
| 1, work, ldwork) |
| } |
| // W *= T or Tᵀ. |
| bi.Dtrmm(blas.Right, blas.Upper, transt, blas.NonUnit, n, k, |
| 1, t, ldt, |
| work, ldwork) |
| // C -= Vᵀ * Wᵀ. |
| if m > k { |
| bi.Dgemm(blas.Trans, blas.Trans, m-k, n, k, |
| -1, v[k:], ldv, work, ldwork, |
| 1, c[k*ldc:], ldc) |
| } |
| // W *= V1. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, n, k, |
| 1, v, ldv, |
| work, ldwork) |
| // C1 -= Wᵀ. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < n; i++ { |
| for j := 0; j < k; j++ { |
| c[j*ldc+i] -= work[i*ldwork+j] |
| } |
| } |
| return |
| } |
| // Form C * H or C * Hᵀ where C = (C1 C2). |
| // W = C * Vᵀ. |
| |
| // W = C1. |
| for j := 0; j < k; j++ { |
| bi.Dcopy(m, c[j:], ldc, work[j:], ldwork) |
| } |
| // W *= V1ᵀ. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, m, k, |
| 1, v, ldv, |
| work, ldwork) |
| if n > k { |
| bi.Dgemm(blas.NoTrans, blas.Trans, m, k, n-k, |
| 1, c[k:], ldc, v[k:], ldv, |
| 1, work, ldwork) |
| } |
| // W *= T or Tᵀ. |
| bi.Dtrmm(blas.Right, blas.Upper, trans, blas.NonUnit, m, k, |
| 1, t, ldt, |
| work, ldwork) |
| // C -= W * V. |
| if n > k { |
| bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n-k, k, |
| -1, work, ldwork, v[k:], ldv, |
| 1, c[k:], ldc) |
| } |
| // W *= V1. |
| bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, m, k, |
| 1, v, ldv, |
| work, ldwork) |
| // C1 -= W. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < m; i++ { |
| for j := 0; j < k; j++ { |
| c[i*ldc+j] -= work[i*ldwork+j] |
| } |
| } |
| return |
| } |
| // V = (V1 V2) where V2 is the last k columns and is lower unit triangular. |
| if side == blas.Left { |
| // Form H * C or Hᵀ C where C = (C1 ; C2). |
| // W = Cᵀ * Vᵀ. |
| |
| // W = C2ᵀ. |
| for j := 0; j < k; j++ { |
| bi.Dcopy(n, c[(m-k+j)*ldc:], 1, work[j:], ldwork) |
| } |
| // W *= V2ᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, n, k, |
| 1, v[m-k:], ldv, |
| work, ldwork) |
| if m > k { |
| bi.Dgemm(blas.Trans, blas.Trans, n, k, m-k, |
| 1, c, ldc, v, ldv, |
| 1, work, ldwork) |
| } |
| // W *= T or Tᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, transt, blas.NonUnit, n, k, |
| 1, t, ldt, |
| work, ldwork) |
| // C -= Vᵀ * Wᵀ. |
| if m > k { |
| bi.Dgemm(blas.Trans, blas.Trans, m-k, n, k, |
| -1, v, ldv, work, ldwork, |
| 1, c, ldc) |
| } |
| // W *= V2. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, n, k, |
| 1, v[m-k:], ldv, |
| work, ldwork) |
| // C2 -= Wᵀ. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < n; i++ { |
| for j := 0; j < k; j++ { |
| c[(m-k+j)*ldc+i] -= work[i*ldwork+j] |
| } |
| } |
| return |
| } |
| // Form C * H or C * Hᵀ where C = (C1 C2). |
| // W = C * Vᵀ. |
| // W = C2. |
| for j := 0; j < k; j++ { |
| bi.Dcopy(m, c[n-k+j:], ldc, work[j:], ldwork) |
| } |
| // W *= V2ᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, m, k, |
| 1, v[n-k:], ldv, |
| work, ldwork) |
| if n > k { |
| bi.Dgemm(blas.NoTrans, blas.Trans, m, k, n-k, |
| 1, c, ldc, v, ldv, |
| 1, work, ldwork) |
| } |
| // W *= T or Tᵀ. |
| bi.Dtrmm(blas.Right, blas.Lower, trans, blas.NonUnit, m, k, |
| 1, t, ldt, |
| work, ldwork) |
| // C -= W * V. |
| if n > k { |
| bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n-k, k, |
| -1, work, ldwork, v, ldv, |
| 1, c, ldc) |
| } |
| // W *= V2. |
| bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, m, k, |
| 1, v[n-k:], ldv, |
| work, ldwork) |
| // C1 -= W. |
| // TODO(btracey): This should use blas.Axpy. |
| for i := 0; i < m; i++ { |
| for j := 0; j < k; j++ { |
| c[i*ldc+n-k+j] -= work[i*ldwork+j] |
| } |
| } |
| } |
