gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/gonum/dptts2.go (about)

     1  // Copyright ©2023 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 "gonum.org/v1/gonum/blas/blas64"
     8  
     9  // dptts2 solves a tridiagonal system of the form
    10  //
    11  //	A * X = B
    12  //
    13  // using the L*D*Lᵀ factorization of A computed by Dpttrf. D is a diagonal
    14  // matrix specified in d, L is a unit bidiagonal matrix whose subdiagonal is
    15  // specified in e, and X and B are n×nrhs matrices.
    16  func (impl Implementation) dptts2(n, nrhs int, d, e []float64, b []float64, ldb int) {
    17  	// Quick return if possible.
    18  	if n <= 1 {
    19  		if n == 1 {
    20  			bi := blas64.Implementation()
    21  			bi.Dscal(nrhs, 1/d[0], b, 1)
    22  		}
    23  		return
    24  	}
    25  
    26  	// Solve A * X = B using the factorization A = L*D*Lᵀ, overwriting each
    27  	// right hand side vector with its solution.
    28  	for j := 0; j < nrhs; j++ {
    29  		// Solve L * x = b.
    30  		for i := 1; i < n; i++ {
    31  			b[i*ldb+j] -= b[(i-1)*ldb+j] * e[i-1]
    32  		}
    33  		// Solve D * Lᵀ * x = b.
    34  		b[(n-1)*ldb+j] /= d[n-1]
    35  		for i := n - 2; i >= 0; i-- {
    36  			b[i*ldb+j] = b[i*ldb+j]/d[i] - b[(i+1)*ldb+j]*e[i]
    37  		}
    38  	}
    39  }