vendor/gonum.org/v1/gonum/lapack/gonum/dgels.go

   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 (
   8         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/lapack"
  10 )
  11
  12 // Dgels finds a minimum-norm solution based on the matrices A and B using the
  13 // QR or LQ factorization. Dgels returns false if the matrix
  14 // A is singular, and true if this solution was successfully found.
  15 //
  16 // The minimization problem solved depends on the input parameters.
  17 //
  18 //  1. If m >= n and trans == blas.NoTrans, Dgels finds X such that || A*X - B||_2
  19 //     is minimized.
  20 //  2. If m < n and trans == blas.NoTrans, Dgels finds the minimum norm solution of
  21 //     A * X = B.
  22 //  3. If m >= n and trans == blas.Trans, Dgels finds the minimum norm solution of
  23 //     A^T * X = B.
  24 //  4. If m < n and trans == blas.Trans, Dgels finds X such that || A*X - B||_2
  25 //     is minimized.
  26 // Note that the least-squares solutions (cases 1 and 3) perform the minimization
  27 // per column of B. This is not the same as finding the minimum-norm matrix.
  28 //
  29 // The matrix A is a general matrix of size m×n and is modified during this call.
  30 // The input matrix B is of size max(m,n)×nrhs, and serves two purposes. On entry,
  31 // the elements of b specify the input matrix B. B has size m×nrhs if
  32 // trans == blas.NoTrans, and n×nrhs if trans == blas.Trans. On exit, the
  33 // leading submatrix of b contains the solution vectors X. If trans == blas.NoTrans,
  34 // this submatrix is of size n×nrhs, and of size m×nrhs otherwise.
  35 //
  36 // work is temporary storage, and lwork specifies the usable memory length.
  37 // At minimum, lwork >= max(m,n) + max(m,n,nrhs), and this function will panic
  38 // otherwise. A longer work will enable blocked algorithms to be called.
  39 // In the special case that lwork == -1, work[0] will be set to the optimal working
  40 // length.
  41 func (impl Implementation) Dgels(trans blas.Transpose, m, n, nrhs int, a []float64, lda int, b []float64, ldb int, work []float64, lwork int) bool {
  42         notran := trans == blas.NoTrans
  43         checkMatrix(m, n, a, lda)
  44         mn := min(m, n)
  45         checkMatrix(max(m, n), nrhs, b, ldb)
  46
  47         // Find optimal block size.
  48         tpsd := true
  49         if notran {
  50                 tpsd = false
  51         }
  52         var nb int
  53         if m >= n {
  54                 nb = impl.Ilaenv(1, "DGEQRF", " ", m, n, -1, -1)
  55                 if tpsd {
  56                         nb = max(nb, impl.Ilaenv(1, "DORMQR", "LN", m, nrhs, n, -1))
  57                 } else {
  58                         nb = max(nb, impl.Ilaenv(1, "DORMQR", "LT", m, nrhs, n, -1))
  59                 }
  60         } else {
  61                 nb = impl.Ilaenv(1, "DGELQF", " ", m, n, -1, -1)
  62                 if tpsd {
  63                         nb = max(nb, impl.Ilaenv(1, "DORMLQ", "LT", n, nrhs, m, -1))
  64                 } else {
  65                         nb = max(nb, impl.Ilaenv(1, "DORMLQ", "LN", n, nrhs, m, -1))
  66                 }
  67         }
  68         if lwork == -1 {
  69                 work[0] = float64(max(1, mn+max(mn, nrhs)*nb))
  70                 return true
  71         }
  72
  73         if len(work) < lwork {
  74                 panic(shortWork)
  75         }
  76         if lwork < mn+max(mn, nrhs) {
  77                 panic(badWork)
  78         }
  79         if m == 0 || n == 0 || nrhs == 0 {
  80                 impl.Dlaset(blas.All, max(m, n), nrhs, 0, 0, b, ldb)
  81                 return true
  82         }
  83
  84         // Scale the input matrices if they contain extreme values.
  85         smlnum := dlamchS / dlamchP
  86         bignum := 1 / smlnum
  87         anrm := impl.Dlange(lapack.MaxAbs, m, n, a, lda, nil)
  88         var iascl int
  89         if anrm > 0 && anrm < smlnum {
  90                 impl.Dlascl(lapack.General, 0, 0, anrm, smlnum, m, n, a, lda)
  91                 iascl = 1
  92         } else if anrm > bignum {
  93                 impl.Dlascl(lapack.General, 0, 0, anrm, bignum, m, n, a, lda)
  94         } else if anrm == 0 {
  95                 // Matrix is all zeros.
  96                 impl.Dlaset(blas.All, max(m, n), nrhs, 0, 0, b, ldb)
  97                 return true
  98         }
  99         brow := m
 100         if tpsd {
 101                 brow = n
 102         }
 103         bnrm := impl.Dlange(lapack.MaxAbs, brow, nrhs, b, ldb, nil)
 104         ibscl := 0
 105         if bnrm > 0 && bnrm < smlnum {
 106                 impl.Dlascl(lapack.General, 0, 0, bnrm, smlnum, brow, nrhs, b, ldb)
 107                 ibscl = 1
 108         } else if bnrm > bignum {
 109                 impl.Dlascl(lapack.General, 0, 0, bnrm, bignum, brow, nrhs, b, ldb)
 110                 ibscl = 2
 111         }
 112
 113         // Solve the minimization problem using a QR or an LQ decomposition.
 114         var scllen int
 115         if m >= n {
 116                 impl.Dgeqrf(m, n, a, lda, work, work[mn:], lwork-mn)
 117                 if !tpsd {
 118                         impl.Dormqr(blas.Left, blas.Trans, m, nrhs, n,
 119                                 a, lda,
 120                                 work[:n],
 121                                 b, ldb,
 122                                 work[mn:], lwork-mn)
 123                         ok := impl.Dtrtrs(blas.Upper, blas.NoTrans, blas.NonUnit, n, nrhs,
 124                                 a, lda,
 125                                 b, ldb)
 126                         if !ok {
 127                                 return false
 128                         }
 129                         scllen = n
 130                 } else {
 131                         ok := impl.Dtrtrs(blas.Upper, blas.Trans, blas.NonUnit, n, nrhs,
 132                                 a, lda,
 133                                 b, ldb)
 134                         if !ok {
 135                                 return false
 136                         }
 137                         for i := n; i < m; i++ {
 138                                 for j := 0; j < nrhs; j++ {
 139                                         b[i*ldb+j] = 0
 140                                 }
 141                         }
 142                         impl.Dormqr(blas.Left, blas.NoTrans, m, nrhs, n,
 143                                 a, lda,
 144                                 work[:n],
 145                                 b, ldb,
 146                                 work[mn:], lwork-mn)
 147                         scllen = m
 148                 }
 149         } else {
 150                 impl.Dgelqf(m, n, a, lda, work, work[mn:], lwork-mn)
 151                 if !tpsd {
 152                         ok := impl.Dtrtrs(blas.Lower, blas.NoTrans, blas.NonUnit,
 153                                 m, nrhs,
 154                                 a, lda,
 155                                 b, ldb)
 156                         if !ok {
 157                                 return false
 158                         }
 159                         for i := m; i < n; i++ {
 160                                 for j := 0; j < nrhs; j++ {
 161                                         b[i*ldb+j] = 0
 162                                 }
 163                         }
 164                         impl.Dormlq(blas.Left, blas.Trans, n, nrhs, m,
 165                                 a, lda,
 166                                 work,
 167                                 b, ldb,
 168                                 work[mn:], lwork-mn)
 169                         scllen = n
 170                 } else {
 171                         impl.Dormlq(blas.Left, blas.NoTrans, n, nrhs, m,
 172                                 a, lda,
 173                                 work,
 174                                 b, ldb,
 175                                 work[mn:], lwork-mn)
 176                         ok := impl.Dtrtrs(blas.Lower, blas.Trans, blas.NonUnit,
 177                                 m, nrhs,
 178                                 a, lda,
 179                                 b, ldb)
 180                         if !ok {
 181                                 return false
 182                         }
 183                 }
 184         }
 185
 186         // Adjust answer vector based on scaling.
 187         if iascl == 1 {
 188                 impl.Dlascl(lapack.General, 0, 0, anrm, smlnum, scllen, nrhs, b, ldb)
 189         }
 190         if iascl == 2 {
 191                 impl.Dlascl(lapack.General, 0, 0, anrm, bignum, scllen, nrhs, b, ldb)
 192         }
 193         if ibscl == 1 {
 194                 impl.Dlascl(lapack.General, 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb)
 195         }
 196         if ibscl == 2 {
 197                 impl.Dlascl(lapack.General, 0, 0, bignum, bnrm, scllen, nrhs, b, ldb)
 198         }
 199         return true
 200 }