github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/testlapack/matgen.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  	"math/rand"
    10  
    11  	"github.com/gonum/blas"
    12  	"github.com/gonum/blas/blas64"
    13  )
    14  
    15  // Dlatm1 computes the entries of dst as specified by mode, cond and rsign.
    16  //
    17  // mode describes how dst will be computed:
    18  //  |mode| == 1: dst[0] = 1 and dst[1:n] = 1/cond
    19  //  |mode| == 2: dst[:n-1] = 1/cond and dst[n-1] = 1
    20  //  |mode| == 3: dst[i] = cond^{-i/(n-1)}, i=0,...,n-1
    21  //  |mode| == 4: dst[i] = 1 - i*(1-1/cond)/(n-1)
    22  //  |mode| == 5: dst[i] = random number in the range (1/cond, 1) such that
    23  //                    their logarithms are uniformly distributed
    24  //  |mode| == 6: dst[i] = random number from the distribution given by dist
    25  // If mode is negative, the order of the elements of dst will be reversed.
    26  // For other values of mode Dlatm1 will panic.
    27  //
    28  // If rsign is true and mode is not ±6, each entry of dst will be multiplied by 1
    29  // or -1 with probability 0.5
    30  //
    31  // dist specifies the type of distribution to be used when mode == ±6:
    32  //  dist == 1: Uniform[0,1)
    33  //  dist == 2: Uniform[-1,1)
    34  //  dist == 3: Normal(0,1)
    35  // For other values of dist Dlatm1 will panic.
    36  //
    37  // rnd is used as a source of random numbers.
    38  func Dlatm1(dst []float64, mode int, cond float64, rsign bool, dist int, rnd *rand.Rand) {
    39  	amode := mode
    40  	if amode < 0 {
    41  		amode = -amode
    42  	}
    43  	if amode < 1 || 6 < amode {
    44  		panic("testlapack: invalid mode")
    45  	}
    46  	if cond < 1 {
    47  		panic("testlapack: cond < 1")
    48  	}
    49  	if amode == 6 && (dist < 1 || 3 < dist) {
    50  		panic("testlapack: invalid dist")
    51  	}
    52  
    53  	n := len(dst)
    54  	if n == 0 {
    55  		return
    56  	}
    57  
    58  	switch amode {
    59  	case 1:
    60  		dst[0] = 1
    61  		for i := 1; i < n; i++ {
    62  			dst[i] = 1 / cond
    63  		}
    64  	case 2:
    65  		for i := 0; i < n-1; i++ {
    66  			dst[i] = 1
    67  		}
    68  		dst[n-1] = 1 / cond
    69  	case 3:
    70  		dst[0] = 1
    71  		if n > 1 {
    72  			alpha := math.Pow(cond, -1/float64(n-1))
    73  			for i := 1; i < n; i++ {
    74  				dst[i] = math.Pow(alpha, float64(i))
    75  			}
    76  		}
    77  	case 4:
    78  		dst[0] = 1
    79  		if n > 1 {
    80  			condInv := 1 / cond
    81  			alpha := (1 - condInv) / float64(n-1)
    82  			for i := 1; i < n; i++ {
    83  				dst[i] = float64(n-i-1)*alpha + condInv
    84  			}
    85  		}
    86  	case 5:
    87  		alpha := math.Log(1 / cond)
    88  		for i := range dst {
    89  			dst[i] = math.Exp(alpha * rnd.Float64())
    90  		}
    91  	case 6:
    92  		switch dist {
    93  		case 1:
    94  			for i := range dst {
    95  				dst[i] = rnd.Float64()
    96  			}
    97  		case 2:
    98  			for i := range dst {
    99  				dst[i] = 2*rnd.Float64() - 1
   100  			}
   101  		case 3:
   102  			for i := range dst {
   103  				dst[i] = rnd.NormFloat64()
   104  			}
   105  		}
   106  	}
   107  
   108  	if rsign && amode != 6 {
   109  		for i, v := range dst {
   110  			if rnd.Float64() < 0.5 {
   111  				dst[i] = -v
   112  			}
   113  		}
   114  	}
   115  
   116  	if mode < 0 {
   117  		for i := 0; i < n/2; i++ {
   118  			dst[i], dst[n-i-1] = dst[n-i-1], dst[i]
   119  		}
   120  	}
   121  }
   122  
   123  // Dlagsy generates an n×n symmetric matrix A, by pre- and post- multiplying a
   124  // real diagonal matrix D with a random orthogonal matrix:
   125  //  A = U * D * U^T.
   126  //
   127  // work must have length at least 2*n, otherwise Dlagsy will panic.
   128  //
   129  // The parameter k is unused but it must satisfy
   130  //  0 <= k <= n-1.
   131  func Dlagsy(n, k int, d []float64, a []float64, lda int, rnd *rand.Rand, work []float64) {
   132  	checkMatrix(n, n, a, lda)
   133  	if k < 0 || max(0, n-1) < k {
   134  		panic("testlapack: invalid value of k")
   135  	}
   136  	if len(d) != n {
   137  		panic("testlapack: bad length of d")
   138  	}
   139  	if len(work) < 2*n {
   140  		panic("testlapack: insufficient work length")
   141  	}
   142  
   143  	// Initialize lower triangle of A to diagonal matrix.
   144  	for i := 1; i < n; i++ {
   145  		for j := 0; j < i; j++ {
   146  			a[i*lda+j] = 0
   147  		}
   148  	}
   149  	for i := 0; i < n; i++ {
   150  		a[i*lda+i] = d[i]
   151  	}
   152  
   153  	bi := blas64.Implementation()
   154  
   155  	// Generate lower triangle of symmetric matrix.
   156  	for i := n - 2; i >= 0; i-- {
   157  		for j := 0; j < n-i; j++ {
   158  			work[j] = rnd.NormFloat64()
   159  		}
   160  		wn := bi.Dnrm2(n-i, work[:n-i], 1)
   161  		wa := math.Copysign(wn, work[0])
   162  		var tau float64
   163  		if wn != 0 {
   164  			wb := work[0] + wa
   165  			bi.Dscal(n-i-1, 1/wb, work[1:n-i], 1)
   166  			work[0] = 1
   167  			tau = wb / wa
   168  		}
   169  
   170  		// Apply random reflection to A[i:n,i:n] from the left and the
   171  		// right.
   172  		//
   173  		// Compute y := tau * A * u.
   174  		bi.Dsymv(blas.Lower, n-i, tau, a[i*lda+i:], lda, work[:n-i], 1, 0, work[n:2*n-i], 1)
   175  
   176  		// Compute v := y - 1/2 * tau * ( y, u ) * u.
   177  		alpha := -0.5 * tau * bi.Ddot(n-i, work[n:2*n-i], 1, work[:n-i], 1)
   178  		bi.Daxpy(n-i, alpha, work[:n-i], 1, work[n:2*n-i], 1)
   179  
   180  		// Apply the transformation as a rank-2 update to A[i:n,i:n].
   181  		bi.Dsyr2(blas.Lower, n-i, -1, work[:n-i], 1, work[n:2*n-i], 1, a[i*lda+i:], lda)
   182  	}
   183  
   184  	// Store full symmetric matrix.
   185  	for i := 1; i < n; i++ {
   186  		for j := 0; j < i; j++ {
   187  			a[j*lda+i] = a[i*lda+j]
   188  		}
   189  	}
   190  }
   191  
   192  // Dlagge generates a real general m×n matrix A, by pre- and post-multiplying
   193  // a real diagonal matrix D with random orthogonal matrices:
   194  //  A = U*D*V.
   195  //
   196  // d must have length min(m,n), and work must have length m+n, otherwise Dlagge
   197  // will panic.
   198  //
   199  // The parameters ku and kl are unused but they must satisfy
   200  //  0 <= kl <= m-1,
   201  //  0 <= ku <= n-1.
   202  func Dlagge(m, n, kl, ku int, d []float64, a []float64, lda int, rnd *rand.Rand, work []float64) {
   203  	checkMatrix(m, n, a, lda)
   204  	if kl < 0 || max(0, m-1) < kl {
   205  		panic("testlapack: invalid value of kl")
   206  	}
   207  	if ku < 0 || max(0, n-1) < ku {
   208  		panic("testlapack: invalid value of ku")
   209  	}
   210  	if len(d) != min(m, n) {
   211  		panic("testlapack: bad length of d")
   212  	}
   213  	if len(work) < m+n {
   214  		panic("testlapack: insufficient work length")
   215  	}
   216  
   217  	// Initialize A to diagonal matrix.
   218  	for i := 0; i < m; i++ {
   219  		for j := 0; j < n; j++ {
   220  			a[i*lda+j] = 0
   221  		}
   222  	}
   223  	for i := 0; i < min(m, n); i++ {
   224  		a[i*lda+i] = d[i]
   225  	}
   226  
   227  	// Quick exit if the user wants a diagonal matrix.
   228  	// if kl == 0 && ku == 0 {
   229  	// 	return
   230  	// }
   231  
   232  	bi := blas64.Implementation()
   233  
   234  	// Pre- and post-multiply A by random orthogonal matrices.
   235  	for i := min(m, n) - 1; i >= 0; i-- {
   236  		if i < m-1 {
   237  			for j := 0; j < m-i; j++ {
   238  				work[j] = rnd.NormFloat64()
   239  			}
   240  			wn := bi.Dnrm2(m-i, work[:m-i], 1)
   241  			wa := math.Copysign(wn, work[0])
   242  			var tau float64
   243  			if wn != 0 {
   244  				wb := work[0] + wa
   245  				bi.Dscal(m-i-1, 1/wb, work[1:m-i], 1)
   246  				work[0] = 1
   247  				tau = wb / wa
   248  			}
   249  
   250  			// Multiply A[i:m,i:n] by random reflection from the left.
   251  			bi.Dgemv(blas.Trans, m-i, n-i,
   252  				1, a[i*lda+i:], lda, work[:m-i], 1,
   253  				0, work[m:m+n-i], 1)
   254  			bi.Dger(m-i, n-i,
   255  				-tau, work[:m-i], 1, work[m:m+n-i], 1,
   256  				a[i*lda+i:], lda)
   257  		}
   258  		if i < n-1 {
   259  			for j := 0; j < n-i; j++ {
   260  				work[j] = rnd.NormFloat64()
   261  			}
   262  			wn := bi.Dnrm2(n-i, work[:n-i], 1)
   263  			wa := math.Copysign(wn, work[0])
   264  			var tau float64
   265  			if wn != 0 {
   266  				wb := work[0] + wa
   267  				bi.Dscal(n-i-1, 1/wb, work[1:n-i], 1)
   268  				work[0] = 1
   269  				tau = wb / wa
   270  			}
   271  
   272  			// Multiply A[i:m,i:n] by random reflection from the right.
   273  			bi.Dgemv(blas.NoTrans, m-i, n-i,
   274  				1, a[i*lda+i:], lda, work[:n-i], 1,
   275  				0, work[n:n+m-i], 1)
   276  			bi.Dger(m-i, n-i,
   277  				-tau, work[n:n+m-i], 1, work[:n-i], 1,
   278  				a[i*lda+i:], lda)
   279  		}
   280  	}
   281  
   282  	// TODO(vladimir-ch): Reduce number of subdiagonals to kl and number of
   283  	// superdiagonals to ku.
   284  }
   285  
   286  // dlarnv fills dst with random numbers from a uniform or normal distribution
   287  // specified by dist:
   288  //  dist=1: uniform(0,1),
   289  //  dist=2: uniform(-1,1),
   290  //  dist=3: normal(0,1).
   291  // For other values of dist dlarnv will panic.
   292  func dlarnv(dst []float64, dist int, rnd *rand.Rand) {
   293  	switch dist {
   294  	default:
   295  		panic("testlapack: invalid dist")
   296  	case 1:
   297  		for i := range dst {
   298  			dst[i] = rnd.Float64()
   299  		}
   300  	case 2:
   301  		for i := range dst {
   302  			dst[i] = 2*rnd.Float64() - 1
   303  		}
   304  	case 3:
   305  		for i := range dst {
   306  			dst[i] = rnd.NormFloat64()
   307  		}
   308  	}
   309  }
   310  
   311  // dlattr generates an n×n triangular test matrix A with its properties uniquely
   312  // determined by imat and uplo, and returns whether A has unit diagonal. If diag
   313  // is blas.Unit, the diagonal elements are set so that A[k,k]=k.
   314  //
   315  // trans specifies whether the matrix A or its transpose will be used.
   316  //
   317  // If imat is greater than 10, dlattr also generates the right hand side of the
   318  // linear system A*x=b, or A^T*x=b. Valid values of imat are 7, and all between 11
   319  // and 19, inclusive.
   320  //
   321  // b mush have length n, and work must have length 3*n, and dlattr will panic
   322  // otherwise.
   323  func dlattr(imat int, uplo blas.Uplo, trans blas.Transpose, n int, a []float64, lda int, b, work []float64, rnd *rand.Rand) (diag blas.Diag) {
   324  	checkMatrix(n, n, a, lda)
   325  	if len(b) != n {
   326  		panic("testlapack: bad length of b")
   327  	}
   328  	if len(work) < 3*n {
   329  		panic("testlapack: insufficient length of work")
   330  	}
   331  	if uplo != blas.Upper && uplo != blas.Lower {
   332  		panic("testlapack: bad uplo")
   333  	}
   334  	if trans != blas.Trans && trans != blas.NoTrans {
   335  		panic("testlapack: bad trans")
   336  	}
   337  
   338  	if n == 0 {
   339  		return blas.NonUnit
   340  	}
   341  
   342  	ulp := dlamchE * dlamchB
   343  	smlnum := dlamchS
   344  	bignum := (1 - ulp) / smlnum
   345  
   346  	bi := blas64.Implementation()
   347  
   348  	switch imat {
   349  	default:
   350  		// TODO(vladimir-ch): Implement the remaining cases.
   351  		panic("testlapack: invalid or unimplemented imat")
   352  	case 7:
   353  		// Identity matrix. The diagonal is set to NaN.
   354  		diag = blas.Unit
   355  		switch uplo {
   356  		case blas.Upper:
   357  			for i := 0; i < n; i++ {
   358  				a[i*lda+i] = math.NaN()
   359  				for j := i + 1; j < n; j++ {
   360  					a[i*lda+j] = 0
   361  				}
   362  			}
   363  		case blas.Lower:
   364  			for i := 0; i < n; i++ {
   365  				for j := 0; j < i; j++ {
   366  					a[i*lda+j] = 0
   367  				}
   368  				a[i*lda+i] = math.NaN()
   369  			}
   370  		}
   371  	case 11:
   372  		// Generate a triangular matrix with elements between -1 and 1,
   373  		// give the diagonal norm 2 to make it well-conditioned, and
   374  		// make the right hand side large so that it requires scaling.
   375  		diag = blas.NonUnit
   376  		switch uplo {
   377  		case blas.Upper:
   378  			for i := 0; i < n-1; i++ {
   379  				dlarnv(a[i*lda+i:i*lda+n], 2, rnd)
   380  			}
   381  		case blas.Lower:
   382  			for i := 1; i < n; i++ {
   383  				dlarnv(a[i*lda:i*lda+i+1], 2, rnd)
   384  			}
   385  		}
   386  		for i := 0; i < n; i++ {
   387  			a[i*lda+i] = math.Copysign(2, a[i*lda+i])
   388  		}
   389  		// Set the right hand side so that the largest value is bignum.
   390  		dlarnv(b, 2, rnd)
   391  		imax := bi.Idamax(n, b, 1)
   392  		bscal := bignum / math.Max(1, b[imax])
   393  		bi.Dscal(n, bscal, b, 1)
   394  	case 12:
   395  		// Make the first diagonal element in the solve small to cause
   396  		// immediate overflow when dividing by T[j,j]. The off-diagonal
   397  		// elements are small (cnorm[j] < 1).
   398  		diag = blas.NonUnit
   399  		tscal := 1 / math.Max(1, float64(n-1))
   400  		switch uplo {
   401  		case blas.Upper:
   402  			for i := 0; i < n; i++ {
   403  				dlarnv(a[i*lda+i:i*lda+n], 2, rnd)
   404  				bi.Dscal(n-i-1, tscal, a[i*lda+i+1:], 1)
   405  				a[i*lda+i] = math.Copysign(1, a[i*lda+i])
   406  			}
   407  			a[(n-1)*lda+n-1] *= smlnum
   408  		case blas.Lower:
   409  			for i := 0; i < n; i++ {
   410  				dlarnv(a[i*lda:i*lda+i+1], 2, rnd)
   411  				bi.Dscal(i, tscal, a[i*lda:], 1)
   412  				a[i*lda+i] = math.Copysign(1, a[i*lda+i])
   413  			}
   414  			a[0] *= smlnum
   415  		}
   416  		dlarnv(b, 2, rnd)
   417  	case 13:
   418  		// Make the first diagonal element in the solve small to cause
   419  		// immediate overflow when dividing by T[j,j]. The off-diagonal
   420  		// elements are O(1) (cnorm[j] > 1).
   421  		diag = blas.NonUnit
   422  		switch uplo {
   423  		case blas.Upper:
   424  			for i := 0; i < n; i++ {
   425  				dlarnv(a[i*lda+i:i*lda+n], 2, rnd)
   426  				a[i*lda+i] = math.Copysign(1, a[i*lda+i])
   427  			}
   428  			a[(n-1)*lda+n-1] *= smlnum
   429  		case blas.Lower:
   430  			for i := 0; i < n; i++ {
   431  				dlarnv(a[i*lda:i*lda+i+1], 2, rnd)
   432  				a[i*lda+i] = math.Copysign(1, a[i*lda+i])
   433  			}
   434  			a[0] *= smlnum
   435  		}
   436  		dlarnv(b, 2, rnd)
   437  	case 14:
   438  		// T is diagonal with small numbers on the diagonal to
   439  		// make the growth factor underflow, but a small right hand side
   440  		// chosen so that the solution does not overflow.
   441  		diag = blas.NonUnit
   442  		switch uplo {
   443  		case blas.Upper:
   444  			for i := 0; i < n; i++ {
   445  				for j := i + 1; j < n; j++ {
   446  					a[i*lda+j] = 0
   447  				}
   448  				if (n-1-i)&0x2 == 0 {
   449  					a[i*lda+i] = smlnum
   450  				} else {
   451  					a[i*lda+i] = 1
   452  				}
   453  			}
   454  		case blas.Lower:
   455  			for i := 0; i < n; i++ {
   456  				for j := 0; j < i; j++ {
   457  					a[i*lda+j] = 0
   458  				}
   459  				if i&0x2 == 0 {
   460  					a[i*lda+i] = smlnum
   461  				} else {
   462  					a[i*lda+i] = 1
   463  				}
   464  			}
   465  		}
   466  		// Set the right hand side alternately zero and small.
   467  		switch uplo {
   468  		case blas.Upper:
   469  			b[0] = 0
   470  			for i := n - 1; i > 0; i -= 2 {
   471  				b[i] = 0
   472  				b[i-1] = smlnum
   473  			}
   474  		case blas.Lower:
   475  			for i := 0; i < n-1; i += 2 {
   476  				b[i] = 0
   477  				b[i+1] = smlnum
   478  			}
   479  			b[n-1] = 0
   480  		}
   481  	case 15:
   482  		// Make the diagonal elements small to cause gradual overflow
   483  		// when dividing by T[j,j]. To control the amount of scaling
   484  		// needed, the matrix is bidiagonal.
   485  		diag = blas.NonUnit
   486  		texp := 1 / math.Max(1, float64(n-1))
   487  		tscal := math.Pow(smlnum, texp)
   488  		switch uplo {
   489  		case blas.Upper:
   490  			for i := 0; i < n; i++ {
   491  				a[i*lda+i] = tscal
   492  				if i < n-1 {
   493  					a[i*lda+i+1] = -1
   494  				}
   495  				for j := i + 2; j < n; j++ {
   496  					a[i*lda+j] = 0
   497  				}
   498  			}
   499  		case blas.Lower:
   500  			for i := 0; i < n; i++ {
   501  				for j := 0; j < i-1; j++ {
   502  					a[i*lda+j] = 0
   503  				}
   504  				if i > 0 {
   505  					a[i*lda+i-1] = -1
   506  				}
   507  				a[i*lda+i] = tscal
   508  			}
   509  		}
   510  		dlarnv(b, 2, rnd)
   511  	case 16:
   512  		// One zero diagonal element.
   513  		diag = blas.NonUnit
   514  		switch uplo {
   515  		case blas.Upper:
   516  			for i := 0; i < n; i++ {
   517  				dlarnv(a[i*lda+i:i*lda+n], 2, rnd)
   518  				a[i*lda+i] = math.Copysign(2, a[i*lda+i])
   519  			}
   520  		case blas.Lower:
   521  			for i := 0; i < n; i++ {
   522  				dlarnv(a[i*lda:i*lda+i+1], 2, rnd)
   523  				a[i*lda+i] = math.Copysign(2, a[i*lda+i])
   524  			}
   525  		}
   526  		iy := n / 2
   527  		a[iy*lda+iy] = 0
   528  		dlarnv(b, 2, rnd)
   529  		bi.Dscal(n, 2, b, 1)
   530  	case 17:
   531  		// Make the offdiagonal elements large to cause overflow when
   532  		// adding a column of T. In the non-transposed case, the matrix
   533  		// is constructed to cause overflow when adding a column in
   534  		// every other step.
   535  		diag = blas.NonUnit
   536  		tscal := (1 - ulp) / (dlamchS / ulp)
   537  		texp := 1.0
   538  		switch uplo {
   539  		case blas.Upper:
   540  			for i := 0; i < n; i++ {
   541  				for j := i; j < n; j++ {
   542  					a[i*lda+j] = 0
   543  				}
   544  			}
   545  			for j := n - 1; j >= 1; j -= 2 {
   546  				a[j] = -tscal / float64(n+1)
   547  				a[j*lda+j] = 1
   548  				b[j] = texp * (1 - ulp)
   549  				a[j-1] = -tscal / float64(n+1) / float64(n+2)
   550  				a[(j-1)*lda+j-1] = 1
   551  				b[j-1] = texp * float64(n*n+n-1)
   552  				texp *= 2
   553  			}
   554  			b[0] = float64(n+1) / float64(n+2) * tscal
   555  		case blas.Lower:
   556  			for i := 0; i < n; i++ {
   557  				for j := 0; j <= i; j++ {
   558  					a[i*lda+j] = 0
   559  				}
   560  			}
   561  			for j := 0; j < n-1; j += 2 {
   562  				a[(n-1)*lda+j] = -tscal / float64(n+1)
   563  				a[j*lda+j] = 1
   564  				b[j] = texp * (1 - ulp)
   565  				a[(n-1)*lda+j+1] = -tscal / float64(n+1) / float64(n+2)
   566  				a[(j+1)*lda+j+1] = 1
   567  				b[j+1] = texp * float64(n*n+n-1)
   568  				texp *= 2
   569  			}
   570  			b[n-1] = float64(n+1) / float64(n+2) * tscal
   571  		}
   572  	case 18:
   573  		// Generate a unit triangular matrix with elements between -1
   574  		// and 1, and make the right hand side large so that it requires
   575  		// scaling. The diagonal is set to NaN.
   576  		diag = blas.Unit
   577  		switch uplo {
   578  		case blas.Upper:
   579  			for i := 0; i < n; i++ {
   580  				a[i*lda+i] = math.NaN()
   581  				dlarnv(a[i*lda+i+1:i*lda+n], 2, rnd)
   582  			}
   583  		case blas.Lower:
   584  			for i := 0; i < n; i++ {
   585  				dlarnv(a[i*lda:i*lda+i], 2, rnd)
   586  				a[i*lda+i] = math.NaN()
   587  			}
   588  		}
   589  		// Set the right hand side so that the largest value is bignum.
   590  		dlarnv(b, 2, rnd)
   591  		iy := bi.Idamax(n, b, 1)
   592  		bnorm := math.Abs(b[iy])
   593  		bscal := bignum / math.Max(1, bnorm)
   594  		bi.Dscal(n, bscal, b, 1)
   595  	case 19:
   596  		// Generate a triangular matrix with elements between
   597  		// bignum/(n-1) and bignum so that at least one of the column
   598  		// norms will exceed bignum.
   599  		// Dlatrs cannot handle this case for (typically) n>5.
   600  		diag = blas.NonUnit
   601  		tleft := bignum / math.Max(1, float64(n-1))
   602  		tscal := bignum * (float64(n-1) / math.Max(1, float64(n)))
   603  		switch uplo {
   604  		case blas.Upper:
   605  			for i := 0; i < n; i++ {
   606  				dlarnv(a[i*lda+i:i*lda+n], 2, rnd)
   607  				for j := i; j < n; j++ {
   608  					aij := a[i*lda+j]
   609  					a[i*lda+j] = math.Copysign(tleft, aij) + tscal*aij
   610  				}
   611  			}
   612  		case blas.Lower:
   613  			for i := 0; i < n; i++ {
   614  				dlarnv(a[i*lda:i*lda+i+1], 2, rnd)
   615  				for j := 0; j <= i; j++ {
   616  					aij := a[i*lda+j]
   617  					a[i*lda+j] = math.Copysign(tleft, aij) + tscal*aij
   618  				}
   619  			}
   620  		}
   621  		dlarnv(b, 2, rnd)
   622  		bi.Dscal(n, 2, b, 1)
   623  	}
   624  
   625  	// Flip the matrix if the transpose will be used.
   626  	if trans == blas.Trans {
   627  		switch uplo {
   628  		case blas.Upper:
   629  			for j := 0; j < n/2; j++ {
   630  				bi.Dswap(n-2*j-1, a[j*lda+j:], 1, a[(j+1)*lda+n-j-1:], -lda)
   631  			}
   632  		case blas.Lower:
   633  			for j := 0; j < n/2; j++ {
   634  				bi.Dswap(n-2*j-1, a[j*lda+j:], lda, a[(n-j-1)*lda+j+1:], -1)
   635  			}
   636  		}
   637  	}
   638  
   639  	return diag
   640  }
   641  
   642  func checkMatrix(m, n int, a []float64, lda int) {
   643  	if m < 0 {
   644  		panic("testlapack: m < 0")
   645  	}
   646  	if n < 0 {
   647  		panic("testlapack: n < 0")
   648  	}
   649  	if lda < max(1, n) {
   650  		panic("testlapack: lda < max(1, n)")
   651  	}
   652  	if len(a) < (m-1)*lda+n {
   653  		panic("testlapack: insufficient matrix slice length")
   654  	}
   655  }