github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dgetf2.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 testlapack
     6  
     7  import (
     8  	"testing"
     9  
    10  	"math/rand"
    11  
    12  	"github.com/gopherd/gonum/blas"
    13  	"github.com/gopherd/gonum/blas/blas64"
    14  	"github.com/gopherd/gonum/floats"
    15  )
    16  
    17  type Dgetf2er interface {
    18  	Dgetf2(m, n int, a []float64, lda int, ipiv []int) bool
    19  }
    20  
    21  func Dgetf2Test(t *testing.T, impl Dgetf2er) {
    22  	rnd := rand.New(rand.NewSource(1))
    23  	for _, test := range []struct {
    24  		m, n, lda int
    25  	}{
    26  		{10, 10, 0},
    27  		{10, 5, 0},
    28  		{10, 5, 0},
    29  
    30  		{10, 10, 20},
    31  		{5, 10, 20},
    32  		{10, 5, 20},
    33  	} {
    34  		m := test.m
    35  		n := test.n
    36  		lda := test.lda
    37  		if lda == 0 {
    38  			lda = n
    39  		}
    40  		a := make([]float64, m*lda)
    41  		for i := range a {
    42  			a[i] = rnd.Float64()
    43  		}
    44  		aCopy := make([]float64, len(a))
    45  		copy(aCopy, a)
    46  
    47  		mn := min(m, n)
    48  		ipiv := make([]int, mn)
    49  		for i := range ipiv {
    50  			ipiv[i] = rnd.Int()
    51  		}
    52  		ok := impl.Dgetf2(m, n, a, lda, ipiv)
    53  		checkPLU(t, ok, m, n, lda, ipiv, a, aCopy, 1e-14, true)
    54  	}
    55  
    56  	// Test with singular matrices (random matrices are almost surely non-singular).
    57  	for _, test := range []struct {
    58  		m, n, lda int
    59  		a         []float64
    60  	}{
    61  		{
    62  			m:   2,
    63  			n:   2,
    64  			lda: 2,
    65  			a: []float64{
    66  				1, 0,
    67  				0, 0,
    68  			},
    69  		},
    70  		{
    71  			m:   2,
    72  			n:   2,
    73  			lda: 2,
    74  			a: []float64{
    75  				1, 5,
    76  				2, 10,
    77  			},
    78  		},
    79  		{
    80  			m:   3,
    81  			n:   3,
    82  			lda: 3,
    83  			// row 3 = row1 + 2 * row2
    84  			a: []float64{
    85  				1, 5, 7,
    86  				2, 10, -3,
    87  				5, 25, 1,
    88  			},
    89  		},
    90  		{
    91  			m:   3,
    92  			n:   4,
    93  			lda: 4,
    94  			// row 3 = row1 + 2 * row2
    95  			a: []float64{
    96  				1, 5, 7, 9,
    97  				2, 10, -3, 11,
    98  				5, 25, 1, 31,
    99  			},
   100  		},
   101  	} {
   102  		if impl.Dgetf2(test.m, test.n, test.a, test.lda, make([]int, min(test.m, test.n))) {
   103  			t.Log("Returned ok with singular matrix.")
   104  		}
   105  	}
   106  }
   107  
   108  // checkPLU checks that the PLU factorization contained in factorize matches
   109  // the original matrix contained in original.
   110  func checkPLU(t *testing.T, ok bool, m, n, lda int, ipiv []int, factorized, original []float64, tol float64, print bool) {
   111  	var hasZeroDiagonal bool
   112  	for i := 0; i < min(m, n); i++ {
   113  		if factorized[i*lda+i] == 0 {
   114  			hasZeroDiagonal = true
   115  			break
   116  		}
   117  	}
   118  	if hasZeroDiagonal && ok {
   119  		t.Error("Has a zero diagonal but returned ok")
   120  	}
   121  	if !hasZeroDiagonal && !ok {
   122  		t.Error("Non-zero diagonal but returned !ok")
   123  	}
   124  
   125  	// Check that the LU decomposition is correct.
   126  	mn := min(m, n)
   127  	l := make([]float64, m*mn)
   128  	ldl := mn
   129  	u := make([]float64, mn*n)
   130  	ldu := n
   131  	for i := 0; i < m; i++ {
   132  		for j := 0; j < n; j++ {
   133  			v := factorized[i*lda+j]
   134  			switch {
   135  			case i == j:
   136  				l[i*ldl+i] = 1
   137  				u[i*ldu+i] = v
   138  			case i > j:
   139  				l[i*ldl+j] = v
   140  			case i < j:
   141  				u[i*ldu+j] = v
   142  			}
   143  		}
   144  	}
   145  
   146  	LU := blas64.General{
   147  		Rows:   m,
   148  		Cols:   n,
   149  		Stride: n,
   150  		Data:   make([]float64, m*n),
   151  	}
   152  	U := blas64.General{
   153  		Rows:   mn,
   154  		Cols:   n,
   155  		Stride: ldu,
   156  		Data:   u,
   157  	}
   158  	L := blas64.General{
   159  		Rows:   m,
   160  		Cols:   mn,
   161  		Stride: ldl,
   162  		Data:   l,
   163  	}
   164  	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, L, U, 0, LU)
   165  
   166  	p := make([]float64, m*m)
   167  	ldp := m
   168  	for i := 0; i < m; i++ {
   169  		p[i*ldp+i] = 1
   170  	}
   171  	for i := len(ipiv) - 1; i >= 0; i-- {
   172  		v := ipiv[i]
   173  		blas64.Swap(blas64.Vector{N: m, Inc: 1, Data: p[i*ldp:]},
   174  			blas64.Vector{N: m, Inc: 1, Data: p[v*ldp:]})
   175  	}
   176  	P := blas64.General{
   177  		Rows:   m,
   178  		Cols:   m,
   179  		Stride: m,
   180  		Data:   p,
   181  	}
   182  	aComp := blas64.General{
   183  		Rows:   m,
   184  		Cols:   n,
   185  		Stride: lda,
   186  		Data:   make([]float64, m*lda),
   187  	}
   188  	copy(aComp.Data, factorized)
   189  	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, P, LU, 0, aComp)
   190  	if !floats.EqualApprox(aComp.Data, original, tol) {
   191  		if print {
   192  			t.Errorf("PLU multiplication does not match original matrix.\nWant: %v\nGot: %v", original, aComp.Data)
   193  			return
   194  		}
   195  		t.Error("PLU multiplication does not match original matrix.")
   196  	}
   197  }