vendor/gonum.org/v1/gonum/mat/dense_arithmetic.go

   1 // Copyright ©2013 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 mat
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/blas/blas64"
  12         "gonum.org/v1/gonum/lapack/lapack64"
  13 )
  14
  15 // Add adds a and b element-wise, placing the result in the receiver. Add
  16 // will panic if the two matrices do not have the same shape.
  17 func (m *Dense) Add(a, b Matrix) {
  18         ar, ac := a.Dims()
  19         br, bc := b.Dims()
  20         if ar != br || ac != bc {
  21                 panic(ErrShape)
  22         }
  23
  24         aU, _ := untranspose(a)
  25         bU, _ := untranspose(b)
  26         m.reuseAs(ar, ac)
  27
  28         if arm, ok := a.(RawMatrixer); ok {
  29                 if brm, ok := b.(RawMatrixer); ok {
  30                         amat, bmat := arm.RawMatrix(), brm.RawMatrix()
  31                         if m != aU {
  32                                 m.checkOverlap(amat)
  33                         }
  34                         if m != bU {
  35                                 m.checkOverlap(bmat)
  36                         }
  37                         for ja, jb, jm := 0, 0, 0; ja < ar*amat.Stride; ja, jb, jm = ja+amat.Stride, jb+bmat.Stride, jm+m.mat.Stride {
  38                                 for i, v := range amat.Data[ja : ja+ac] {
  39                                         m.mat.Data[i+jm] = v + bmat.Data[i+jb]
  40                                 }
  41                         }
  42                         return
  43                 }
  44         }
  45
  46         var restore func()
  47         if m == aU {
  48                 m, restore = m.isolatedWorkspace(aU)
  49                 defer restore()
  50         } else if m == bU {
  51                 m, restore = m.isolatedWorkspace(bU)
  52                 defer restore()
  53         }
  54
  55         for r := 0; r < ar; r++ {
  56                 for c := 0; c < ac; c++ {
  57                         m.set(r, c, a.At(r, c)+b.At(r, c))
  58                 }
  59         }
  60 }
  61
  62 // Sub subtracts the matrix b from a, placing the result in the receiver. Sub
  63 // will panic if the two matrices do not have the same shape.
  64 func (m *Dense) Sub(a, b Matrix) {
  65         ar, ac := a.Dims()
  66         br, bc := b.Dims()
  67         if ar != br || ac != bc {
  68                 panic(ErrShape)
  69         }
  70
  71         aU, _ := untranspose(a)
  72         bU, _ := untranspose(b)
  73         m.reuseAs(ar, ac)
  74
  75         if arm, ok := a.(RawMatrixer); ok {
  76                 if brm, ok := b.(RawMatrixer); ok {
  77                         amat, bmat := arm.RawMatrix(), brm.RawMatrix()
  78                         if m != aU {
  79                                 m.checkOverlap(amat)
  80                         }
  81                         if m != bU {
  82                                 m.checkOverlap(bmat)
  83                         }
  84                         for ja, jb, jm := 0, 0, 0; ja < ar*amat.Stride; ja, jb, jm = ja+amat.Stride, jb+bmat.Stride, jm+m.mat.Stride {
  85                                 for i, v := range amat.Data[ja : ja+ac] {
  86                                         m.mat.Data[i+jm] = v - bmat.Data[i+jb]
  87                                 }
  88                         }
  89                         return
  90                 }
  91         }
  92
  93         var restore func()
  94         if m == aU {
  95                 m, restore = m.isolatedWorkspace(aU)
  96                 defer restore()
  97         } else if m == bU {
  98                 m, restore = m.isolatedWorkspace(bU)
  99                 defer restore()
 100         }
 101
 102         for r := 0; r < ar; r++ {
 103                 for c := 0; c < ac; c++ {
 104                         m.set(r, c, a.At(r, c)-b.At(r, c))
 105                 }
 106         }
 107 }
 108
 109 // MulElem performs element-wise multiplication of a and b, placing the result
 110 // in the receiver. MulElem will panic if the two matrices do not have the same
 111 // shape.
 112 func (m *Dense) MulElem(a, b Matrix) {
 113         ar, ac := a.Dims()
 114         br, bc := b.Dims()
 115         if ar != br || ac != bc {
 116                 panic(ErrShape)
 117         }
 118
 119         aU, _ := untranspose(a)
 120         bU, _ := untranspose(b)
 121         m.reuseAs(ar, ac)
 122
 123         if arm, ok := a.(RawMatrixer); ok {
 124                 if brm, ok := b.(RawMatrixer); ok {
 125                         amat, bmat := arm.RawMatrix(), brm.RawMatrix()
 126                         if m != aU {
 127                                 m.checkOverlap(amat)
 128                         }
 129                         if m != bU {
 130                                 m.checkOverlap(bmat)
 131                         }
 132                         for ja, jb, jm := 0, 0, 0; ja < ar*amat.Stride; ja, jb, jm = ja+amat.Stride, jb+bmat.Stride, jm+m.mat.Stride {
 133                                 for i, v := range amat.Data[ja : ja+ac] {
 134                                         m.mat.Data[i+jm] = v * bmat.Data[i+jb]
 135                                 }
 136                         }
 137                         return
 138                 }
 139         }
 140
 141         var restore func()
 142         if m == aU {
 143                 m, restore = m.isolatedWorkspace(aU)
 144                 defer restore()
 145         } else if m == bU {
 146                 m, restore = m.isolatedWorkspace(bU)
 147                 defer restore()
 148         }
 149
 150         for r := 0; r < ar; r++ {
 151                 for c := 0; c < ac; c++ {
 152                         m.set(r, c, a.At(r, c)*b.At(r, c))
 153                 }
 154         }
 155 }
 156
 157 // DivElem performs element-wise division of a by b, placing the result
 158 // in the receiver. DivElem will panic if the two matrices do not have the same
 159 // shape.
 160 func (m *Dense) DivElem(a, b Matrix) {
 161         ar, ac := a.Dims()
 162         br, bc := b.Dims()
 163         if ar != br || ac != bc {
 164                 panic(ErrShape)
 165         }
 166
 167         aU, _ := untranspose(a)
 168         bU, _ := untranspose(b)
 169         m.reuseAs(ar, ac)
 170
 171         if arm, ok := a.(RawMatrixer); ok {
 172                 if brm, ok := b.(RawMatrixer); ok {
 173                         amat, bmat := arm.RawMatrix(), brm.RawMatrix()
 174                         if m != aU {
 175                                 m.checkOverlap(amat)
 176                         }
 177                         if m != bU {
 178                                 m.checkOverlap(bmat)
 179                         }
 180                         for ja, jb, jm := 0, 0, 0; ja < ar*amat.Stride; ja, jb, jm = ja+amat.Stride, jb+bmat.Stride, jm+m.mat.Stride {
 181                                 for i, v := range amat.Data[ja : ja+ac] {
 182                                         m.mat.Data[i+jm] = v / bmat.Data[i+jb]
 183                                 }
 184                         }
 185                         return
 186                 }
 187         }
 188
 189         var restore func()
 190         if m == aU {
 191                 m, restore = m.isolatedWorkspace(aU)
 192                 defer restore()
 193         } else if m == bU {
 194                 m, restore = m.isolatedWorkspace(bU)
 195                 defer restore()
 196         }
 197
 198         for r := 0; r < ar; r++ {
 199                 for c := 0; c < ac; c++ {
 200                         m.set(r, c, a.At(r, c)/b.At(r, c))
 201                 }
 202         }
 203 }
 204
 205 // Inverse computes the inverse of the matrix a, storing the result into the
 206 // receiver. If a is ill-conditioned, a Condition error will be returned.
 207 // Note that matrix inversion is numerically unstable, and should generally
 208 // be avoided where possible, for example by using the Solve routines.
 209 func (m *Dense) Inverse(a Matrix) error {
 210         // TODO(btracey): Special case for RawTriangular, etc.
 211         r, c := a.Dims()
 212         if r != c {
 213                 panic(ErrSquare)
 214         }
 215         m.reuseAs(a.Dims())
 216         aU, aTrans := untranspose(a)
 217         switch rm := aU.(type) {
 218         case RawMatrixer:
 219                 if m != aU || aTrans {
 220                         if m == aU || m.checkOverlap(rm.RawMatrix()) {
 221                                 tmp := getWorkspace(r, c, false)
 222                                 tmp.Copy(a)
 223                                 m.Copy(tmp)
 224                                 putWorkspace(tmp)
 225                                 break
 226                         }
 227                         m.Copy(a)
 228                 }
 229         default:
 230                 m.Copy(a)
 231         }
 232         ipiv := getInts(r, false)
 233         defer putInts(ipiv)
 234         ok := lapack64.Getrf(m.mat, ipiv)
 235         if !ok {
 236                 return Condition(math.Inf(1))
 237         }
 238         work := getFloats(4*r, false) // must be at least 4*r for cond.
 239         lapack64.Getri(m.mat, ipiv, work, -1)
 240         if int(work[0]) > 4*r {
 241                 l := int(work[0])
 242                 putFloats(work)
 243                 work = getFloats(l, false)
 244         } else {
 245                 work = work[:4*r]
 246         }
 247         defer putFloats(work)
 248         lapack64.Getri(m.mat, ipiv, work, len(work))
 249         norm := lapack64.Lange(CondNorm, m.mat, work)
 250         rcond := lapack64.Gecon(CondNorm, m.mat, norm, work, ipiv) // reuse ipiv
 251         if rcond == 0 {
 252                 return Condition(math.Inf(1))
 253         }
 254         cond := 1 / rcond
 255         if cond > ConditionTolerance {
 256                 return Condition(cond)
 257         }
 258         return nil
 259 }
 260
 261 // Mul takes the matrix product of a and b, placing the result in the receiver.
 262 // If the number of columns in a does not equal the number of rows in b, Mul will panic.
 263 func (m *Dense) Mul(a, b Matrix) {
 264         ar, ac := a.Dims()
 265         br, bc := b.Dims()
 266
 267         if ac != br {
 268                 panic(ErrShape)
 269         }
 270
 271         aU, aTrans := untranspose(a)
 272         bU, bTrans := untranspose(b)
 273         m.reuseAs(ar, bc)
 274         var restore func()
 275         if m == aU {
 276                 m, restore = m.isolatedWorkspace(aU)
 277                 defer restore()
 278         } else if m == bU {
 279                 m, restore = m.isolatedWorkspace(bU)
 280                 defer restore()
 281         }
 282         aT := blas.NoTrans
 283         if aTrans {
 284                 aT = blas.Trans
 285         }
 286         bT := blas.NoTrans
 287         if bTrans {
 288                 bT = blas.Trans
 289         }
 290
 291         // Some of the cases do not have a transpose option, so create
 292         // temporary memory.
 293         // C = A^T * B = (B^T * A)^T
 294         // C^T = B^T * A.
 295         if aUrm, ok := aU.(RawMatrixer); ok {
 296                 amat := aUrm.RawMatrix()
 297                 if restore == nil {
 298                         m.checkOverlap(amat)
 299                 }
 300                 if bUrm, ok := bU.(RawMatrixer); ok {
 301                         bmat := bUrm.RawMatrix()
 302                         if restore == nil {
 303                                 m.checkOverlap(bmat)
 304                         }
 305                         blas64.Gemm(aT, bT, 1, amat, bmat, 0, m.mat)
 306                         return
 307                 }
 308                 if bU, ok := bU.(RawSymmetricer); ok {
 309                         bmat := bU.RawSymmetric()
 310                         if aTrans {
 311                                 c := getWorkspace(ac, ar, false)
 312                                 blas64.Symm(blas.Left, 1, bmat, amat, 0, c.mat)
 313                                 strictCopy(m, c.T())
 314                                 putWorkspace(c)
 315                                 return
 316                         }
 317                         blas64.Symm(blas.Right, 1, bmat, amat, 0, m.mat)
 318                         return
 319                 }
 320                 if bU, ok := bU.(RawTriangular); ok {
 321                         // Trmm updates in place, so copy aU first.
 322                         bmat := bU.RawTriangular()
 323                         if aTrans {
 324                                 c := getWorkspace(ac, ar, false)
 325                                 var tmp Dense
 326                                 tmp.SetRawMatrix(amat)
 327                                 c.Copy(&tmp)
 328                                 bT := blas.Trans
 329                                 if bTrans {
 330                                         bT = blas.NoTrans
 331                                 }
 332                                 blas64.Trmm(blas.Left, bT, 1, bmat, c.mat)
 333                                 strictCopy(m, c.T())
 334                                 putWorkspace(c)
 335                                 return
 336                         }
 337                         m.Copy(a)
 338                         blas64.Trmm(blas.Right, bT, 1, bmat, m.mat)
 339                         return
 340                 }
 341                 if bU, ok := bU.(*VecDense); ok {
 342                         m.checkOverlap(bU.asGeneral())
 343                         bvec := bU.RawVector()
 344                         if bTrans {
 345                                 // {ar,1} x {1,bc}, which is not a vector.
 346                                 // Instead, construct B as a General.
 347                                 bmat := blas64.General{
 348                                         Rows:   bc,
 349                                         Cols:   1,
 350                                         Stride: bvec.Inc,
 351                                         Data:   bvec.Data,
 352                                 }
 353                                 blas64.Gemm(aT, bT, 1, amat, bmat, 0, m.mat)
 354                                 return
 355                         }
 356                         cvec := blas64.Vector{
 357                                 Inc:  m.mat.Stride,
 358                                 Data: m.mat.Data,
 359                         }
 360                         blas64.Gemv(aT, 1, amat, bvec, 0, cvec)
 361                         return
 362                 }
 363         }
 364         if bUrm, ok := bU.(RawMatrixer); ok {
 365                 bmat := bUrm.RawMatrix()
 366                 if restore == nil {
 367                         m.checkOverlap(bmat)
 368                 }
 369                 if aU, ok := aU.(RawSymmetricer); ok {
 370                         amat := aU.RawSymmetric()
 371                         if bTrans {
 372                                 c := getWorkspace(bc, br, false)
 373                                 blas64.Symm(blas.Right, 1, amat, bmat, 0, c.mat)
 374                                 strictCopy(m, c.T())
 375                                 putWorkspace(c)
 376                                 return
 377                         }
 378                         blas64.Symm(blas.Left, 1, amat, bmat, 0, m.mat)
 379                         return
 380                 }
 381                 if aU, ok := aU.(RawTriangular); ok {
 382                         // Trmm updates in place, so copy bU first.
 383                         amat := aU.RawTriangular()
 384                         if bTrans {
 385                                 c := getWorkspace(bc, br, false)
 386                                 var tmp Dense
 387                                 tmp.SetRawMatrix(bmat)
 388                                 c.Copy(&tmp)
 389                                 aT := blas.Trans
 390                                 if aTrans {
 391                                         aT = blas.NoTrans
 392                                 }
 393                                 blas64.Trmm(blas.Right, aT, 1, amat, c.mat)
 394                                 strictCopy(m, c.T())
 395                                 putWorkspace(c)
 396                                 return
 397                         }
 398                         m.Copy(b)
 399                         blas64.Trmm(blas.Left, aT, 1, amat, m.mat)
 400                         return
 401                 }
 402                 if aU, ok := aU.(*VecDense); ok {
 403                         m.checkOverlap(aU.asGeneral())
 404                         avec := aU.RawVector()
 405                         if aTrans {
 406                                 // {1,ac} x {ac, bc}
 407                                 // Transpose B so that the vector is on the right.
 408                                 cvec := blas64.Vector{
 409                                         Inc:  1,
 410                                         Data: m.mat.Data,
 411                                 }
 412                                 bT := blas.Trans
 413                                 if bTrans {
 414                                         bT = blas.NoTrans
 415                                 }
 416                                 blas64.Gemv(bT, 1, bmat, avec, 0, cvec)
 417                                 return
 418                         }
 419                         // {ar,1} x {1,bc} which is not a vector result.
 420                         // Instead, construct A as a General.
 421                         amat := blas64.General{
 422                                 Rows:   ar,
 423                                 Cols:   1,
 424                                 Stride: avec.Inc,
 425                                 Data:   avec.Data,
 426                         }
 427                         blas64.Gemm(aT, bT, 1, amat, bmat, 0, m.mat)
 428                         return
 429                 }
 430         }
 431
 432         row := getFloats(ac, false)
 433         defer putFloats(row)
 434         for r := 0; r < ar; r++ {
 435                 for i := range row {
 436                         row[i] = a.At(r, i)
 437                 }
 438                 for c := 0; c < bc; c++ {
 439                         var v float64
 440                         for i, e := range row {
 441                                 v += e * b.At(i, c)
 442                         }
 443                         m.mat.Data[r*m.mat.Stride+c] = v
 444                 }
 445         }
 446 }
 447
 448 // strictCopy copies a into m panicking if the shape of a and m differ.
 449 func strictCopy(m *Dense, a Matrix) {
 450         r, c := m.Copy(a)
 451         if r != m.mat.Rows || c != m.mat.Cols {
 452                 // Panic with a string since this
 453                 // is not a user-facing panic.
 454                 panic(ErrShape.Error())
 455         }
 456 }
 457
 458 // Exp calculates the exponential of the matrix a, e^a, placing the result
 459 // in the receiver. Exp will panic with matrix.ErrShape if a is not square.
 460 //
 461 // Exp uses the scaling and squaring method described in section 3 of
 462 // http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf.
 463 func (m *Dense) Exp(a Matrix) {
 464         r, c := a.Dims()
 465         if r != c {
 466                 panic(ErrShape)
 467         }
 468
 469         var w *Dense
 470         if m.IsZero() {
 471                 m.reuseAsZeroed(r, r)
 472                 w = m
 473         } else {
 474                 w = getWorkspace(r, r, true)
 475         }
 476         for i := 0; i < r*r; i += r + 1 {
 477                 w.mat.Data[i] = 1
 478         }
 479
 480         const (
 481                 terms   = 10
 482                 scaling = 4
 483         )
 484
 485         small := getWorkspace(r, r, false)
 486         small.Scale(math.Pow(2, -scaling), a)
 487         power := getWorkspace(r, r, false)
 488         power.Copy(small)
 489
 490         var (
 491                 tmp   = getWorkspace(r, r, false)
 492                 factI = 1.
 493         )
 494         for i := 1.; i < terms; i++ {
 495                 factI *= i
 496
 497                 // This is OK to do because power and tmp are
 498                 // new Dense values so all rows are contiguous.
 499                 // TODO(kortschak) Make this explicit in the NewDense doc comment.
 500                 for j, v := range power.mat.Data {
 501                         tmp.mat.Data[j] = v / factI
 502                 }
 503
 504                 w.Add(w, tmp)
 505                 if i < terms-1 {
 506                         tmp.Mul(power, small)
 507                         tmp, power = power, tmp
 508                 }
 509         }
 510         putWorkspace(small)
 511         putWorkspace(power)
 512         for i := 0; i < scaling; i++ {
 513                 tmp.Mul(w, w)
 514                 tmp, w = w, tmp
 515         }
 516         putWorkspace(tmp)
 517
 518         if w != m {
 519                 m.Copy(w)
 520                 putWorkspace(w)
 521         }
 522 }
 523
 524 // Pow calculates the integral power of the matrix a to n, placing the result
 525 // in the receiver. Pow will panic if n is negative or if a is not square.
 526 func (m *Dense) Pow(a Matrix, n int) {
 527         if n < 0 {
 528                 panic("matrix: illegal power")
 529         }
 530         r, c := a.Dims()
 531         if r != c {
 532                 panic(ErrShape)
 533         }
 534
 535         m.reuseAs(r, c)
 536
 537         // Take possible fast paths.
 538         switch n {
 539         case 0:
 540                 for i := 0; i < r; i++ {
 541                         zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c])
 542                         m.mat.Data[i*m.mat.Stride+i] = 1
 543                 }
 544                 return
 545         case 1:
 546                 m.Copy(a)
 547                 return
 548         case 2:
 549                 m.Mul(a, a)
 550                 return
 551         }
 552
 553         // Perform iterative exponentiation by squaring in work space.
 554         w := getWorkspace(r, r, false)
 555         w.Copy(a)
 556         s := getWorkspace(r, r, false)
 557         s.Copy(a)
 558         x := getWorkspace(r, r, false)
 559         for n--; n > 0; n >>= 1 {
 560                 if n&1 != 0 {
 561                         x.Mul(w, s)
 562                         w, x = x, w
 563                 }
 564                 if n != 1 {
 565                         x.Mul(s, s)
 566                         s, x = x, s
 567                 }
 568         }
 569         m.Copy(w)
 570         putWorkspace(w)
 571         putWorkspace(s)
 572         putWorkspace(x)
 573 }
 574
 575 // Scale multiplies the elements of a by f, placing the result in the receiver.
 576 //
 577 // See the Scaler interface for more information.
 578 func (m *Dense) Scale(f float64, a Matrix) {
 579         ar, ac := a.Dims()
 580
 581         m.reuseAs(ar, ac)
 582
 583         aU, aTrans := untranspose(a)
 584         if rm, ok := aU.(RawMatrixer); ok {
 585                 amat := rm.RawMatrix()
 586                 if m == aU || m.checkOverlap(amat) {
 587                         var restore func()
 588                         m, restore = m.isolatedWorkspace(a)
 589                         defer restore()
 590                 }
 591                 if !aTrans {
 592                         for ja, jm := 0, 0; ja < ar*amat.Stride; ja, jm = ja+amat.Stride, jm+m.mat.Stride {
 593                                 for i, v := range amat.Data[ja : ja+ac] {
 594                                         m.mat.Data[i+jm] = v * f
 595                                 }
 596                         }
 597                 } else {
 598                         for ja, jm := 0, 0; ja < ac*amat.Stride; ja, jm = ja+amat.Stride, jm+1 {
 599                                 for i, v := range amat.Data[ja : ja+ar] {
 600                                         m.mat.Data[i*m.mat.Stride+jm] = v * f
 601                                 }
 602                         }
 603                 }
 604                 return
 605         }
 606
 607         for r := 0; r < ar; r++ {
 608                 for c := 0; c < ac; c++ {
 609                         m.set(r, c, f*a.At(r, c))
 610                 }
 611         }
 612 }
 613
 614 // Apply applies the function fn to each of the elements of a, placing the
 615 // resulting matrix in the receiver. The function fn takes a row/column
 616 // index and element value and returns some function of that tuple.
 617 func (m *Dense) Apply(fn func(i, j int, v float64) float64, a Matrix) {
 618         ar, ac := a.Dims()
 619
 620         m.reuseAs(ar, ac)
 621
 622         aU, aTrans := untranspose(a)
 623         if rm, ok := aU.(RawMatrixer); ok {
 624                 amat := rm.RawMatrix()
 625                 if m == aU || m.checkOverlap(amat) {
 626                         var restore func()
 627                         m, restore = m.isolatedWorkspace(a)
 628                         defer restore()
 629                 }
 630                 if !aTrans {
 631                         for j, ja, jm := 0, 0, 0; ja < ar*amat.Stride; j, ja, jm = j+1, ja+amat.Stride, jm+m.mat.Stride {
 632                                 for i, v := range amat.Data[ja : ja+ac] {
 633                                         m.mat.Data[i+jm] = fn(j, i, v)
 634                                 }
 635                         }
 636                 } else {
 637                         for j, ja, jm := 0, 0, 0; ja < ac*amat.Stride; j, ja, jm = j+1, ja+amat.Stride, jm+1 {
 638                                 for i, v := range amat.Data[ja : ja+ar] {
 639                                         m.mat.Data[i*m.mat.Stride+jm] = fn(i, j, v)
 640                                 }
 641                         }
 642                 }
 643                 return
 644         }
 645
 646         for r := 0; r < ar; r++ {
 647                 for c := 0; c < ac; c++ {
 648                         m.set(r, c, fn(r, c, a.At(r, c)))
 649                 }
 650         }
 651 }
 652
 653 // RankOne performs a rank-one update to the matrix a and stores the result
 654 // in the receiver. If a is zero, see Outer.
 655 //  m = a + alpha * x * y'
 656 func (m *Dense) RankOne(a Matrix, alpha float64, x, y Vector) {
 657         ar, ac := a.Dims()
 658         xr, xc := x.Dims()
 659         if xr != ar || xc != 1 {
 660                 panic(ErrShape)
 661         }
 662         yr, yc := y.Dims()
 663         if yr != ac || yc != 1 {
 664                 panic(ErrShape)
 665         }
 666
 667         if a != m {
 668                 aU, _ := untranspose(a)
 669                 if rm, ok := aU.(RawMatrixer); ok {
 670                         m.checkOverlap(rm.RawMatrix())
 671                 }
 672         }
 673
 674         var xmat, ymat blas64.Vector
 675         fast := true
 676         xU, _ := untranspose(x)
 677         if rv, ok := xU.(RawVectorer); ok {
 678                 xmat = rv.RawVector()
 679                 m.checkOverlap((&VecDense{mat: xmat, n: x.Len()}).asGeneral())
 680         } else {
 681                 fast = false
 682         }
 683         yU, _ := untranspose(y)
 684         if rv, ok := yU.(RawVectorer); ok {
 685                 ymat = rv.RawVector()
 686                 m.checkOverlap((&VecDense{mat: ymat, n: y.Len()}).asGeneral())
 687         } else {
 688                 fast = false
 689         }
 690
 691         if fast {
 692                 if m != a {
 693                         m.reuseAs(ar, ac)
 694                         m.Copy(a)
 695                 }
 696                 blas64.Ger(alpha, xmat, ymat, m.mat)
 697                 return
 698         }
 699
 700         m.reuseAs(ar, ac)
 701         for i := 0; i < ar; i++ {
 702                 for j := 0; j < ac; j++ {
 703                         m.set(i, j, a.At(i, j)+alpha*x.AtVec(i)*y.AtVec(j))
 704                 }
 705         }
 706 }
 707
 708 // Outer calculates the outer product of the column vectors x and y,
 709 // and stores the result in the receiver.
 710 //  m = alpha * x * y'
 711 // In order to update an existing matrix, see RankOne.
 712 func (m *Dense) Outer(alpha float64, x, y Vector) {
 713         xr, xc := x.Dims()
 714         if xc != 1 {
 715                 panic(ErrShape)
 716         }
 717         yr, yc := y.Dims()
 718         if yc != 1 {
 719                 panic(ErrShape)
 720         }
 721
 722         r := xr
 723         c := yr
 724
 725         // Copied from reuseAs with use replaced by useZeroed
 726         // and a final zero of the matrix elements if we pass
 727         // the shape checks.
 728         // TODO(kortschak): Factor out into reuseZeroedAs if
 729         // we find another case that needs it.
 730         if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
 731                 // Panic as a string, not a mat.Error.
 732                 panic("mat: caps not correctly set")
 733         }
 734         if m.IsZero() {
 735                 m.mat = blas64.General{
 736                         Rows:   r,
 737                         Cols:   c,
 738                         Stride: c,
 739                         Data:   useZeroed(m.mat.Data, r*c),
 740                 }
 741                 m.capRows = r
 742                 m.capCols = c
 743         } else if r != m.mat.Rows || c != m.mat.Cols {
 744                 panic(ErrShape)
 745         }
 746
 747         var xmat, ymat blas64.Vector
 748         fast := true
 749         xU, _ := untranspose(x)
 750         if rv, ok := xU.(RawVectorer); ok {
 751                 xmat = rv.RawVector()
 752                 m.checkOverlap((&VecDense{mat: xmat, n: x.Len()}).asGeneral())
 753         } else {
 754                 fast = false
 755         }
 756         yU, _ := untranspose(y)
 757         if rv, ok := yU.(RawVectorer); ok {
 758                 ymat = rv.RawVector()
 759                 m.checkOverlap((&VecDense{mat: ymat, n: y.Len()}).asGeneral())
 760         } else {
 761                 fast = false
 762         }
 763
 764         if fast {
 765                 for i := 0; i < r; i++ {
 766                         zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c])
 767                 }
 768                 blas64.Ger(alpha, xmat, ymat, m.mat)
 769                 return
 770         }
 771
 772         for i := 0; i < r; i++ {
 773                 for j := 0; j < c; j++ {
 774                         m.set(i, j, alpha*x.AtVec(i)*y.AtVec(j))
 775                 }
 776         }
 777 }