github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dlags2.go (about)

     1  // Copyright ©2017 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 testlapack
     6  
     7  import (
     8  	"math"
     9  	"testing"
    10  
    11  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/blas/blas64"
    15  	"github.com/gopherd/gonum/floats/scalar"
    16  )
    17  
    18  type Dlags2er interface {
    19  	Dlags2(upper bool, a1, a2, a3, b1, b2, b3 float64) (csu, snu, csv, snv, csq, snq float64)
    20  }
    21  
    22  func Dlags2Test(t *testing.T, impl Dlags2er) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	for _, upper := range []bool{true, false} {
    25  		for i := 0; i < 100; i++ {
    26  			// Generate randomly the elements of a 2×2 matrix A
    27  			//  [ a1 a2 ] or [ a1 0  ]
    28  			//  [ 0  a3 ]    [ a2 a3 ]
    29  			a1 := rnd.Float64()
    30  			a2 := rnd.Float64()
    31  			a3 := rnd.Float64()
    32  			// Generate randomly the elements of a 2×2 matrix B.
    33  			//  [ b1 b2 ] or [ b1 0  ]
    34  			//  [ 0  b3 ]    [ b2 b3 ]
    35  			b1 := rnd.Float64()
    36  			b2 := rnd.Float64()
    37  			b3 := rnd.Float64()
    38  
    39  			// Compute orthogal matrices U, V, Q
    40  			//  U = [  csu  snu ], V = [  csv snv ], Q = [  csq   snq ]
    41  			//      [ -snu  csu ]      [ -snv csv ]      [ -snq   csq ]
    42  			// that transform A and B.
    43  			csu, snu, csv, snv, csq, snq := impl.Dlags2(upper, a1, a2, a3, b1, b2, b3)
    44  
    45  			// Check that U, V, Q are orthogonal matrices (their
    46  			// determinant is equal to 1).
    47  			detU := det2x2(csu, snu, -snu, csu)
    48  			if !scalar.EqualWithinAbsOrRel(math.Abs(detU), 1, 1e-14, 1e-14) {
    49  				t.Errorf("U not orthogonal: det(U)=%v", detU)
    50  			}
    51  			detV := det2x2(csv, snv, -snv, csv)
    52  			if !scalar.EqualWithinAbsOrRel(math.Abs(detV), 1, 1e-14, 1e-14) {
    53  				t.Errorf("V not orthogonal: det(V)=%v", detV)
    54  			}
    55  			detQ := det2x2(csq, snq, -snq, csq)
    56  			if !scalar.EqualWithinAbsOrRel(math.Abs(detQ), 1, 1e-14, 1e-14) {
    57  				t.Errorf("Q not orthogonal: det(Q)=%v", detQ)
    58  			}
    59  
    60  			// Create U, V, Q explicitly as dense matrices.
    61  			u := blas64.General{
    62  				Rows:   2,
    63  				Cols:   2,
    64  				Stride: 2,
    65  				Data:   []float64{csu, snu, -snu, csu},
    66  			}
    67  			v := blas64.General{
    68  				Rows:   2,
    69  				Cols:   2,
    70  				Stride: 2,
    71  				Data:   []float64{csv, snv, -snv, csv},
    72  			}
    73  			q := blas64.General{
    74  				Rows:   2,
    75  				Cols:   2,
    76  				Stride: 2,
    77  				Data:   []float64{csq, snq, -snq, csq},
    78  			}
    79  
    80  			// Create A and B explicitly as dense matrices.
    81  			a := blas64.General{Rows: 2, Cols: 2, Stride: 2}
    82  			b := blas64.General{Rows: 2, Cols: 2, Stride: 2}
    83  			if upper {
    84  				a.Data = []float64{a1, a2, 0, a3}
    85  				b.Data = []float64{b1, b2, 0, b3}
    86  			} else {
    87  				a.Data = []float64{a1, 0, a2, a3}
    88  				b.Data = []float64{b1, 0, b2, b3}
    89  			}
    90  
    91  			tmp := blas64.General{Rows: 2, Cols: 2, Stride: 2, Data: make([]float64, 4)}
    92  			// Transform A as Uᵀ*A*Q.
    93  			blas64.Gemm(blas.Trans, blas.NoTrans, 1, u, a, 0, tmp)
    94  			blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, q, 0, a)
    95  			// Transform B as Vᵀ*A*Q.
    96  			blas64.Gemm(blas.Trans, blas.NoTrans, 1, v, b, 0, tmp)
    97  			blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, q, 0, b)
    98  
    99  			// Extract elements of transformed A and B that should be equal to zero.
   100  			var gotA, gotB float64
   101  			if upper {
   102  				gotA = a.Data[1]
   103  				gotB = b.Data[1]
   104  			} else {
   105  				gotA = a.Data[2]
   106  				gotB = b.Data[2]
   107  			}
   108  			// Check that they are indeed zero.
   109  			if !scalar.EqualWithinAbsOrRel(gotA, 0, 1e-14, 1e-14) {
   110  				t.Errorf("unexpected non-zero value for zero triangle of Uᵀ*A*Q: %v", gotA)
   111  			}
   112  			if !scalar.EqualWithinAbsOrRel(gotB, 0, 1e-14, 1e-14) {
   113  				t.Errorf("unexpected non-zero value for zero triangle of Vᵀ*B*Q: %v", gotB)
   114  			}
   115  		}
   116  	}
   117  }
   118  
   119  func det2x2(a, b, c, d float64) float64 { return a*d - b*c }