vendor/gonum.org/v1/gonum/lapack/testlapack/matgen.go

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