vendor/gonum.org/v1/gonum/lapack/gonum/dlarfx.go

   1 // Copyright ©2016 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 gonum
   6
   7 import "gonum.org/v1/gonum/blas"
   8
   9 // Dlarfx applies an elementary reflector H to a real m×n matrix C, from either
  10 // the left or the right, with loop unrolling when the reflector has order less
  11 // than 11.
  12 //
  13 // H is represented in the form
  14 //  H = I - tau * v * v^T,
  15 // where tau is a real scalar and v is a real vector. If tau = 0, then H is
  16 // taken to be the identity matrix.
  17 //
  18 // v must have length equal to m if side == blas.Left, and equal to n if side ==
  19 // blas.Right, otherwise Dlarfx will panic.
  20 //
  21 // c and ldc represent the m×n matrix C. On return, C is overwritten by the
  22 // matrix H * C if side == blas.Left, or C * H if side == blas.Right.
  23 //
  24 // work must have length at least n if side == blas.Left, and at least m if side
  25 // == blas.Right, otherwise Dlarfx will panic. work is not referenced if H has
  26 // order < 11.
  27 //
  28 // Dlarfx is an internal routine. It is exported for testing purposes.
  29 func (impl Implementation) Dlarfx(side blas.Side, m, n int, v []float64, tau float64, c []float64, ldc int, work []float64) {
  30         checkMatrix(m, n, c, ldc)
  31         switch side {
  32         case blas.Left:
  33                 checkVector(m, v, 1)
  34                 if m > 10 && len(work) < n {
  35                         panic(badWork)
  36                 }
  37         case blas.Right:
  38                 checkVector(n, v, 1)
  39                 if n > 10 && len(work) < m {
  40                         panic(badWork)
  41                 }
  42         default:
  43                 panic(badSide)
  44         }
  45
  46         if tau == 0 {
  47                 return
  48         }
  49
  50         if side == blas.Left {
  51                 // Form H * C, where H has order m.
  52                 switch m {
  53                 default: // Code for general m.
  54                         impl.Dlarf(side, m, n, v, 1, tau, c, ldc, work)
  55                         return
  56
  57                 case 0: // No-op for zero size matrix.
  58                         return
  59
  60                 case 1: // Special code for 1×1 Householder matrix.
  61                         t0 := 1 - tau*v[0]*v[0]
  62                         for j := 0; j < n; j++ {
  63                                 c[j] *= t0
  64                         }
  65                         return
  66
  67                 case 2: // Special code for 2×2 Householder matrix.
  68                         v0 := v[0]
  69                         t0 := tau * v0
  70                         v1 := v[1]
  71                         t1 := tau * v1
  72                         for j := 0; j < n; j++ {
  73                                 sum := v0*c[j] + v1*c[ldc+j]
  74                                 c[j] -= sum * t0
  75                                 c[ldc+j] -= sum * t1
  76                         }
  77                         return
  78
  79                 case 3: // Special code for 3×3 Householder matrix.
  80                         v0 := v[0]
  81                         t0 := tau * v0
  82                         v1 := v[1]
  83                         t1 := tau * v1
  84                         v2 := v[2]
  85                         t2 := tau * v2
  86                         for j := 0; j < n; j++ {
  87                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j]
  88                                 c[j] -= sum * t0
  89                                 c[ldc+j] -= sum * t1
  90                                 c[2*ldc+j] -= sum * t2
  91                         }
  92                         return
  93
  94                 case 4: // Special code for 4×4 Householder matrix.
  95                         v0 := v[0]
  96                         t0 := tau * v0
  97                         v1 := v[1]
  98                         t1 := tau * v1
  99                         v2 := v[2]
 100                         t2 := tau * v2
 101                         v3 := v[3]
 102                         t3 := tau * v3
 103                         for j := 0; j < n; j++ {
 104                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j]
 105                                 c[j] -= sum * t0
 106                                 c[ldc+j] -= sum * t1
 107                                 c[2*ldc+j] -= sum * t2
 108                                 c[3*ldc+j] -= sum * t3
 109                         }
 110                         return
 111
 112                 case 5: // Special code for 5×5 Householder matrix.
 113                         v0 := v[0]
 114                         t0 := tau * v0
 115                         v1 := v[1]
 116                         t1 := tau * v1
 117                         v2 := v[2]
 118                         t2 := tau * v2
 119                         v3 := v[3]
 120                         t3 := tau * v3
 121                         v4 := v[4]
 122                         t4 := tau * v4
 123                         for j := 0; j < n; j++ {
 124                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j]
 125                                 c[j] -= sum * t0
 126                                 c[ldc+j] -= sum * t1
 127                                 c[2*ldc+j] -= sum * t2
 128                                 c[3*ldc+j] -= sum * t3
 129                                 c[4*ldc+j] -= sum * t4
 130                         }
 131                         return
 132
 133                 case 6: // Special code for 6×6 Householder matrix.
 134                         v0 := v[0]
 135                         t0 := tau * v0
 136                         v1 := v[1]
 137                         t1 := tau * v1
 138                         v2 := v[2]
 139                         t2 := tau * v2
 140                         v3 := v[3]
 141                         t3 := tau * v3
 142                         v4 := v[4]
 143                         t4 := tau * v4
 144                         v5 := v[5]
 145                         t5 := tau * v5
 146                         for j := 0; j < n; j++ {
 147                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
 148                                         v5*c[5*ldc+j]
 149                                 c[j] -= sum * t0
 150                                 c[ldc+j] -= sum * t1
 151                                 c[2*ldc+j] -= sum * t2
 152                                 c[3*ldc+j] -= sum * t3
 153                                 c[4*ldc+j] -= sum * t4
 154                                 c[5*ldc+j] -= sum * t5
 155                         }
 156                         return
 157
 158                 case 7: // Special code for 7×7 Householder matrix.
 159                         v0 := v[0]
 160                         t0 := tau * v0
 161                         v1 := v[1]
 162                         t1 := tau * v1
 163                         v2 := v[2]
 164                         t2 := tau * v2
 165                         v3 := v[3]
 166                         t3 := tau * v3
 167                         v4 := v[4]
 168                         t4 := tau * v4
 169                         v5 := v[5]
 170                         t5 := tau * v5
 171                         v6 := v[6]
 172                         t6 := tau * v6
 173                         for j := 0; j < n; j++ {
 174                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
 175                                         v5*c[5*ldc+j] + v6*c[6*ldc+j]
 176                                 c[j] -= sum * t0
 177                                 c[ldc+j] -= sum * t1
 178                                 c[2*ldc+j] -= sum * t2
 179                                 c[3*ldc+j] -= sum * t3
 180                                 c[4*ldc+j] -= sum * t4
 181                                 c[5*ldc+j] -= sum * t5
 182                                 c[6*ldc+j] -= sum * t6
 183                         }
 184                         return
 185
 186                 case 8: // Special code for 8×8 Householder matrix.
 187                         v0 := v[0]
 188                         t0 := tau * v0
 189                         v1 := v[1]
 190                         t1 := tau * v1
 191                         v2 := v[2]
 192                         t2 := tau * v2
 193                         v3 := v[3]
 194                         t3 := tau * v3
 195                         v4 := v[4]
 196                         t4 := tau * v4
 197                         v5 := v[5]
 198                         t5 := tau * v5
 199                         v6 := v[6]
 200                         t6 := tau * v6
 201                         v7 := v[7]
 202                         t7 := tau * v7
 203                         for j := 0; j < n; j++ {
 204                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
 205                                         v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j]
 206                                 c[j] -= sum * t0
 207                                 c[ldc+j] -= sum * t1
 208                                 c[2*ldc+j] -= sum * t2
 209                                 c[3*ldc+j] -= sum * t3
 210                                 c[4*ldc+j] -= sum * t4
 211                                 c[5*ldc+j] -= sum * t5
 212                                 c[6*ldc+j] -= sum * t6
 213                                 c[7*ldc+j] -= sum * t7
 214                         }
 215                         return
 216
 217                 case 9: // Special code for 9×9 Householder matrix.
 218                         v0 := v[0]
 219                         t0 := tau * v0
 220                         v1 := v[1]
 221                         t1 := tau * v1
 222                         v2 := v[2]
 223                         t2 := tau * v2
 224                         v3 := v[3]
 225                         t3 := tau * v3
 226                         v4 := v[4]
 227                         t4 := tau * v4
 228                         v5 := v[5]
 229                         t5 := tau * v5
 230                         v6 := v[6]
 231                         t6 := tau * v6
 232                         v7 := v[7]
 233                         t7 := tau * v7
 234                         v8 := v[8]
 235                         t8 := tau * v8
 236                         for j := 0; j < n; j++ {
 237                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
 238                                         v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j] + v8*c[8*ldc+j]
 239                                 c[j] -= sum * t0
 240                                 c[ldc+j] -= sum * t1
 241                                 c[2*ldc+j] -= sum * t2
 242                                 c[3*ldc+j] -= sum * t3
 243                                 c[4*ldc+j] -= sum * t4
 244                                 c[5*ldc+j] -= sum * t5
 245                                 c[6*ldc+j] -= sum * t6
 246                                 c[7*ldc+j] -= sum * t7
 247                                 c[8*ldc+j] -= sum * t8
 248                         }
 249                         return
 250
 251                 case 10: // Special code for 10×10 Householder matrix.
 252                         v0 := v[0]
 253                         t0 := tau * v0
 254                         v1 := v[1]
 255                         t1 := tau * v1
 256                         v2 := v[2]
 257                         t2 := tau * v2
 258                         v3 := v[3]
 259                         t3 := tau * v3
 260                         v4 := v[4]
 261                         t4 := tau * v4
 262                         v5 := v[5]
 263                         t5 := tau * v5
 264                         v6 := v[6]
 265                         t6 := tau * v6
 266                         v7 := v[7]
 267                         t7 := tau * v7
 268                         v8 := v[8]
 269                         t8 := tau * v8
 270                         v9 := v[9]
 271                         t9 := tau * v9
 272                         for j := 0; j < n; j++ {
 273                                 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
 274                                         v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j] + v8*c[8*ldc+j] + v9*c[9*ldc+j]
 275                                 c[j] -= sum * t0
 276                                 c[ldc+j] -= sum * t1
 277                                 c[2*ldc+j] -= sum * t2
 278                                 c[3*ldc+j] -= sum * t3
 279                                 c[4*ldc+j] -= sum * t4
 280                                 c[5*ldc+j] -= sum * t5
 281                                 c[6*ldc+j] -= sum * t6
 282                                 c[7*ldc+j] -= sum * t7
 283                                 c[8*ldc+j] -= sum * t8
 284                                 c[9*ldc+j] -= sum * t9
 285                         }
 286                         return
 287                 }
 288         }
 289
 290         // Form C * H, where H has order n.
 291         switch n {
 292         default: // Code for general n.
 293                 impl.Dlarf(side, m, n, v, 1, tau, c, ldc, work)
 294                 return
 295
 296         case 0: // No-op for zero size matrix.
 297                 return
 298
 299         case 1: // Special code for 1×1 Householder matrix.
 300                 t0 := 1 - tau*v[0]*v[0]
 301                 for j := 0; j < m; j++ {
 302                         c[j*ldc] *= t0
 303                 }
 304                 return
 305
 306         case 2: // Special code for 2×2 Householder matrix.
 307                 v0 := v[0]
 308                 t0 := tau * v0
 309                 v1 := v[1]
 310                 t1 := tau * v1
 311                 for j := 0; j < m; j++ {
 312                         cs := c[j*ldc:]
 313                         sum := v0*cs[0] + v1*cs[1]
 314                         cs[0] -= sum * t0
 315                         cs[1] -= sum * t1
 316                 }
 317                 return
 318
 319         case 3: // Special code for 3×3 Householder matrix.
 320                 v0 := v[0]
 321                 t0 := tau * v0
 322                 v1 := v[1]
 323                 t1 := tau * v1
 324                 v2 := v[2]
 325                 t2 := tau * v2
 326                 for j := 0; j < m; j++ {
 327                         cs := c[j*ldc:]
 328                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2]
 329                         cs[0] -= sum * t0
 330                         cs[1] -= sum * t1
 331                         cs[2] -= sum * t2
 332                 }
 333                 return
 334
 335         case 4: // Special code for 4×4 Householder matrix.
 336                 v0 := v[0]
 337                 t0 := tau * v0
 338                 v1 := v[1]
 339                 t1 := tau * v1
 340                 v2 := v[2]
 341                 t2 := tau * v2
 342                 v3 := v[3]
 343                 t3 := tau * v3
 344                 for j := 0; j < m; j++ {
 345                         cs := c[j*ldc:]
 346                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3]
 347                         cs[0] -= sum * t0
 348                         cs[1] -= sum * t1
 349                         cs[2] -= sum * t2
 350                         cs[3] -= sum * t3
 351                 }
 352                 return
 353
 354         case 5: // Special code for 5×5 Householder matrix.
 355                 v0 := v[0]
 356                 t0 := tau * v0
 357                 v1 := v[1]
 358                 t1 := tau * v1
 359                 v2 := v[2]
 360                 t2 := tau * v2
 361                 v3 := v[3]
 362                 t3 := tau * v3
 363                 v4 := v[4]
 364                 t4 := tau * v4
 365                 for j := 0; j < m; j++ {
 366                         cs := c[j*ldc:]
 367                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4]
 368                         cs[0] -= sum * t0
 369                         cs[1] -= sum * t1
 370                         cs[2] -= sum * t2
 371                         cs[3] -= sum * t3
 372                         cs[4] -= sum * t4
 373                 }
 374                 return
 375
 376         case 6: // Special code for 6×6 Householder matrix.
 377                 v0 := v[0]
 378                 t0 := tau * v0
 379                 v1 := v[1]
 380                 t1 := tau * v1
 381                 v2 := v[2]
 382                 t2 := tau * v2
 383                 v3 := v[3]
 384                 t3 := tau * v3
 385                 v4 := v[4]
 386                 t4 := tau * v4
 387                 v5 := v[5]
 388                 t5 := tau * v5
 389                 for j := 0; j < m; j++ {
 390                         cs := c[j*ldc:]
 391                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] + v5*cs[5]
 392                         cs[0] -= sum * t0
 393                         cs[1] -= sum * t1
 394                         cs[2] -= sum * t2
 395                         cs[3] -= sum * t3
 396                         cs[4] -= sum * t4
 397                         cs[5] -= sum * t5
 398                 }
 399                 return
 400
 401         case 7: // Special code for 7×7 Householder matrix.
 402                 v0 := v[0]
 403                 t0 := tau * v0
 404                 v1 := v[1]
 405                 t1 := tau * v1
 406                 v2 := v[2]
 407                 t2 := tau * v2
 408                 v3 := v[3]
 409                 t3 := tau * v3
 410                 v4 := v[4]
 411                 t4 := tau * v4
 412                 v5 := v[5]
 413                 t5 := tau * v5
 414                 v6 := v[6]
 415                 t6 := tau * v6
 416                 for j := 0; j < m; j++ {
 417                         cs := c[j*ldc:]
 418                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
 419                                 v5*cs[5] + v6*cs[6]
 420                         cs[0] -= sum * t0
 421                         cs[1] -= sum * t1
 422                         cs[2] -= sum * t2
 423                         cs[3] -= sum * t3
 424                         cs[4] -= sum * t4
 425                         cs[5] -= sum * t5
 426                         cs[6] -= sum * t6
 427                 }
 428                 return
 429
 430         case 8: // Special code for 8×8 Householder matrix.
 431                 v0 := v[0]
 432                 t0 := tau * v0
 433                 v1 := v[1]
 434                 t1 := tau * v1
 435                 v2 := v[2]
 436                 t2 := tau * v2
 437                 v3 := v[3]
 438                 t3 := tau * v3
 439                 v4 := v[4]
 440                 t4 := tau * v4
 441                 v5 := v[5]
 442                 t5 := tau * v5
 443                 v6 := v[6]
 444                 t6 := tau * v6
 445                 v7 := v[7]
 446                 t7 := tau * v7
 447                 for j := 0; j < m; j++ {
 448                         cs := c[j*ldc:]
 449                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
 450                                 v5*cs[5] + v6*cs[6] + v7*cs[7]
 451                         cs[0] -= sum * t0
 452                         cs[1] -= sum * t1
 453                         cs[2] -= sum * t2
 454                         cs[3] -= sum * t3
 455                         cs[4] -= sum * t4
 456                         cs[5] -= sum * t5
 457                         cs[6] -= sum * t6
 458                         cs[7] -= sum * t7
 459                 }
 460                 return
 461
 462         case 9: // Special code for 9×9 Householder matrix.
 463                 v0 := v[0]
 464                 t0 := tau * v0
 465                 v1 := v[1]
 466                 t1 := tau * v1
 467                 v2 := v[2]
 468                 t2 := tau * v2
 469                 v3 := v[3]
 470                 t3 := tau * v3
 471                 v4 := v[4]
 472                 t4 := tau * v4
 473                 v5 := v[5]
 474                 t5 := tau * v5
 475                 v6 := v[6]
 476                 t6 := tau * v6
 477                 v7 := v[7]
 478                 t7 := tau * v7
 479                 v8 := v[8]
 480                 t8 := tau * v8
 481                 for j := 0; j < m; j++ {
 482                         cs := c[j*ldc:]
 483                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
 484                                 v5*cs[5] + v6*cs[6] + v7*cs[7] + v8*cs[8]
 485                         cs[0] -= sum * t0
 486                         cs[1] -= sum * t1
 487                         cs[2] -= sum * t2
 488                         cs[3] -= sum * t3
 489                         cs[4] -= sum * t4
 490                         cs[5] -= sum * t5
 491                         cs[6] -= sum * t6
 492                         cs[7] -= sum * t7
 493                         cs[8] -= sum * t8
 494                 }
 495                 return
 496
 497         case 10: // Special code for 10×10 Householder matrix.
 498                 v0 := v[0]
 499                 t0 := tau * v0
 500                 v1 := v[1]
 501                 t1 := tau * v1
 502                 v2 := v[2]
 503                 t2 := tau * v2
 504                 v3 := v[3]
 505                 t3 := tau * v3
 506                 v4 := v[4]
 507                 t4 := tau * v4
 508                 v5 := v[5]
 509                 t5 := tau * v5
 510                 v6 := v[6]
 511                 t6 := tau * v6
 512                 v7 := v[7]
 513                 t7 := tau * v7
 514                 v8 := v[8]
 515                 t8 := tau * v8
 516                 v9 := v[9]
 517                 t9 := tau * v9
 518                 for j := 0; j < m; j++ {
 519                         cs := c[j*ldc:]
 520                         sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
 521                                 v5*cs[5] + v6*cs[6] + v7*cs[7] + v8*cs[8] + v9*cs[9]
 522                         cs[0] -= sum * t0
 523                         cs[1] -= sum * t1
 524                         cs[2] -= sum * t2
 525                         cs[3] -= sum * t3
 526                         cs[4] -= sum * t4
 527                         cs[5] -= sum * t5
 528                         cs[6] -= sum * t6
 529                         cs[7] -= sum * t7
 530                         cs[8] -= sum * t8
 531                         cs[9] -= sum * t9
 532                 }
 533                 return
 534         }
 535 }