github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dorm2r.go (about)

     1  // Copyright ©2015 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package gonum
     6  
     7  import "github.com/jingcheng-WU/gonum/blas"
     8  
     9  // Dorm2r multiplies a general matrix C by an orthogonal matrix from a QR factorization
    10  // determined by Dgeqrf.
    11  //  C = Q * C   if side == blas.Left and trans == blas.NoTrans
    12  //  C = Qᵀ * C  if side == blas.Left and trans == blas.Trans
    13  //  C = C * Q   if side == blas.Right and trans == blas.NoTrans
    14  //  C = C * Qᵀ  if side == blas.Right and trans == blas.Trans
    15  // If side == blas.Left, a is a matrix of size m×k, and if side == blas.Right
    16  // a is of size n×k.
    17  //
    18  // tau contains the Householder factors and is of length at least k and this function
    19  // will panic otherwise.
    20  //
    21  // work is temporary storage of length at least n if side == blas.Left
    22  // and at least m if side == blas.Right and this function will panic otherwise.
    23  //
    24  // Dorm2r is an internal routine. It is exported for testing purposes.
    25  func (impl Implementation) Dorm2r(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64) {
    26  	left := side == blas.Left
    27  	switch {
    28  	case !left && side != blas.Right:
    29  		panic(badSide)
    30  	case trans != blas.Trans && trans != blas.NoTrans:
    31  		panic(badTrans)
    32  	case m < 0:
    33  		panic(mLT0)
    34  	case n < 0:
    35  		panic(nLT0)
    36  	case k < 0:
    37  		panic(kLT0)
    38  	case left && k > m:
    39  		panic(kGTM)
    40  	case !left && k > n:
    41  		panic(kGTN)
    42  	case lda < max(1, k):
    43  		panic(badLdA)
    44  	case ldc < max(1, n):
    45  		panic(badLdC)
    46  	}
    47  
    48  	// Quick return if possible.
    49  	if m == 0 || n == 0 || k == 0 {
    50  		return
    51  	}
    52  
    53  	switch {
    54  	case left && len(a) < (m-1)*lda+k:
    55  		panic(shortA)
    56  	case !left && len(a) < (n-1)*lda+k:
    57  		panic(shortA)
    58  	case len(c) < (m-1)*ldc+n:
    59  		panic(shortC)
    60  	case len(tau) < k:
    61  		panic(shortTau)
    62  	case left && len(work) < n:
    63  		panic(shortWork)
    64  	case !left && len(work) < m:
    65  		panic(shortWork)
    66  	}
    67  
    68  	if left {
    69  		if trans == blas.NoTrans {
    70  			for i := k - 1; i >= 0; i-- {
    71  				aii := a[i*lda+i]
    72  				a[i*lda+i] = 1
    73  				impl.Dlarf(side, m-i, n, a[i*lda+i:], lda, tau[i], c[i*ldc:], ldc, work)
    74  				a[i*lda+i] = aii
    75  			}
    76  			return
    77  		}
    78  		for i := 0; i < k; i++ {
    79  			aii := a[i*lda+i]
    80  			a[i*lda+i] = 1
    81  			impl.Dlarf(side, m-i, n, a[i*lda+i:], lda, tau[i], c[i*ldc:], ldc, work)
    82  			a[i*lda+i] = aii
    83  		}
    84  		return
    85  	}
    86  	if trans == blas.NoTrans {
    87  		for i := 0; i < k; i++ {
    88  			aii := a[i*lda+i]
    89  			a[i*lda+i] = 1
    90  			impl.Dlarf(side, m, n-i, a[i*lda+i:], lda, tau[i], c[i:], ldc, work)
    91  			a[i*lda+i] = aii
    92  		}
    93  		return
    94  	}
    95  	for i := k - 1; i >= 0; i-- {
    96  		aii := a[i*lda+i]
    97  		a[i*lda+i] = 1
    98  		impl.Dlarf(side, m, n-i, a[i*lda+i:], lda, tau[i], c[i:], ldc, work)
    99  		a[i*lda+i] = aii
   100  	}
   101  }