github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dorgbr.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 native
     6  
     7  import "github.com/gonum/lapack"
     8  
     9  // Dorgbr generates one of the matrices Q or P^T computed by Dgebrd
    10  // computed from the decomposition Dgebrd. See Dgebd2 for the description of
    11  // Q and P^T.
    12  //
    13  // If vect == lapack.ApplyQ, then a is assumed to have been an m×k matrix and
    14  // Q is of order m. If m >= k, then Dorgbr returns the first n columns of Q
    15  // where m >= n >= k. If m < k, then Dorgbr returns Q as an m×m matrix.
    16  //
    17  // If vect == lapack.ApplyP, then A is assumed to have been a k×n matrix, and
    18  // P^T is of order n. If k < n, then Dorgbr returns the first m rows of P^T,
    19  // where n >= m >= k. If k >= n, then Dorgbr returns P^T as an n×n matrix.
    20  //
    21  // Dorgbr is an internal routine. It is exported for testing purposes.
    22  func (impl Implementation) Dorgbr(vect lapack.DecompUpdate, m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
    23  	mn := min(m, n)
    24  	wantq := vect == lapack.ApplyQ
    25  	if wantq {
    26  		if m < n || n < min(m, k) || m < min(m, k) {
    27  			panic(badDims)
    28  		}
    29  	} else {
    30  		if n < m || m < min(n, k) || n < min(n, k) {
    31  			panic(badDims)
    32  		}
    33  	}
    34  	if wantq {
    35  		if m >= k {
    36  			checkMatrix(m, k, a, lda)
    37  		} else {
    38  			checkMatrix(m, m, a, lda)
    39  		}
    40  	} else {
    41  		if n >= k {
    42  			checkMatrix(k, n, a, lda)
    43  		} else {
    44  			checkMatrix(n, n, a, lda)
    45  		}
    46  	}
    47  	work[0] = 1
    48  	if wantq {
    49  		if m >= k {
    50  			impl.Dorgqr(m, n, k, a, lda, tau, work, -1)
    51  		} else if m > 1 {
    52  			impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, -1)
    53  		}
    54  	} else {
    55  		if k < n {
    56  			impl.Dorglq(m, n, k, a, lda, tau, work, -1)
    57  		} else if n > 1 {
    58  			impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, -1)
    59  		}
    60  	}
    61  	lworkopt := int(work[0])
    62  	lworkopt = max(lworkopt, mn)
    63  	if lwork == -1 {
    64  		work[0] = float64(lworkopt)
    65  		return
    66  	}
    67  	if len(work) < lwork {
    68  		panic(badWork)
    69  	}
    70  	if lwork < mn {
    71  		panic(badWork)
    72  	}
    73  	if m == 0 || n == 0 {
    74  		work[0] = 1
    75  		return
    76  	}
    77  	if wantq {
    78  		// Form Q, determined by a call to Dgebrd to reduce an m×k matrix.
    79  		if m >= k {
    80  			impl.Dorgqr(m, n, k, a, lda, tau, work, lwork)
    81  		} else {
    82  			// Shift the vectors which define the elementary reflectors one
    83  			// column to the right, and set the first row and column of Q to
    84  			// those of the unit matrix.
    85  			for j := m - 1; j >= 1; j-- {
    86  				a[j] = 0
    87  				for i := j + 1; i < m; i++ {
    88  					a[i*lda+j] = a[i*lda+j-1]
    89  				}
    90  			}
    91  			a[0] = 1
    92  			for i := 1; i < m; i++ {
    93  				a[i*lda] = 0
    94  			}
    95  			if m > 1 {
    96  				// Form Q[1:m-1, 1:m-1]
    97  				impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, lwork)
    98  			}
    99  		}
   100  	} else {
   101  		// Form P^T, determined by a call to Dgebrd to reduce a k×n matrix.
   102  		if k < n {
   103  			impl.Dorglq(m, n, k, a, lda, tau, work, lwork)
   104  		} else {
   105  			// Shift the vectors which define the elementary reflectors one
   106  			// row downward, and set the first row and column of P^T to
   107  			// those of the unit matrix.
   108  			a[0] = 1
   109  			for i := 1; i < n; i++ {
   110  				a[i*lda] = 0
   111  			}
   112  			for j := 1; j < n; j++ {
   113  				for i := j - 1; i >= 1; i-- {
   114  					a[i*lda+j] = a[(i-1)*lda+j]
   115  				}
   116  				a[j] = 0
   117  			}
   118  			if n > 1 {
   119  				impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork)
   120  			}
   121  		}
   122  	}
   123  	work[0] = float64(lworkopt)
   124  }