gonum.org/v1/gonum@v0.14.0/blas/gonum/level2cmplx64.go (about)

     1  // Code generated by "go generate gonum.org/v1/gonum/blas/gonum”; DO NOT EDIT.
     2  
     3  // Copyright ©2017 The Gonum Authors. All rights reserved.
     4  // Use of this source code is governed by a BSD-style
     5  // license that can be found in the LICENSE file.
     6  
     7  package gonum
     8  
     9  import (
    10  	cmplx "gonum.org/v1/gonum/internal/cmplx64"
    11  
    12  	"gonum.org/v1/gonum/blas"
    13  	"gonum.org/v1/gonum/internal/asm/c64"
    14  )
    15  
    16  var _ blas.Complex64Level2 = Implementation{}
    17  
    18  // Cgbmv performs one of the matrix-vector operations
    19  //
    20  //	y = alpha * A * x + beta * y   if trans = blas.NoTrans
    21  //	y = alpha * Aᵀ * x + beta * y  if trans = blas.Trans
    22  //	y = alpha * Aᴴ * x + beta * y  if trans = blas.ConjTrans
    23  //
    24  // where alpha and beta are scalars, x and y are vectors, and A is an m×n band matrix
    25  // with kL sub-diagonals and kU super-diagonals.
    26  //
    27  // Complex64 implementations are autogenerated and not directly tested.
    28  func (Implementation) Cgbmv(trans blas.Transpose, m, n, kL, kU int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) {
    29  	switch trans {
    30  	default:
    31  		panic(badTranspose)
    32  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
    33  	}
    34  	if m < 0 {
    35  		panic(mLT0)
    36  	}
    37  	if n < 0 {
    38  		panic(nLT0)
    39  	}
    40  	if kL < 0 {
    41  		panic(kLLT0)
    42  	}
    43  	if kU < 0 {
    44  		panic(kULT0)
    45  	}
    46  	if lda < kL+kU+1 {
    47  		panic(badLdA)
    48  	}
    49  	if incX == 0 {
    50  		panic(zeroIncX)
    51  	}
    52  	if incY == 0 {
    53  		panic(zeroIncY)
    54  	}
    55  
    56  	// Quick return if possible.
    57  	if m == 0 || n == 0 {
    58  		return
    59  	}
    60  
    61  	// For zero matrix size the following slice length checks are trivially satisfied.
    62  	if len(a) < lda*(min(m, n+kL)-1)+kL+kU+1 {
    63  		panic(shortA)
    64  	}
    65  	var lenX, lenY int
    66  	if trans == blas.NoTrans {
    67  		lenX, lenY = n, m
    68  	} else {
    69  		lenX, lenY = m, n
    70  	}
    71  	if (incX > 0 && len(x) <= (lenX-1)*incX) || (incX < 0 && len(x) <= (1-lenX)*incX) {
    72  		panic(shortX)
    73  	}
    74  	if (incY > 0 && len(y) <= (lenY-1)*incY) || (incY < 0 && len(y) <= (1-lenY)*incY) {
    75  		panic(shortY)
    76  	}
    77  
    78  	// Quick return if possible.
    79  	if alpha == 0 && beta == 1 {
    80  		return
    81  	}
    82  
    83  	var kx int
    84  	if incX < 0 {
    85  		kx = (1 - lenX) * incX
    86  	}
    87  	var ky int
    88  	if incY < 0 {
    89  		ky = (1 - lenY) * incY
    90  	}
    91  
    92  	// Form y = beta*y.
    93  	if beta != 1 {
    94  		if incY == 1 {
    95  			if beta == 0 {
    96  				for i := range y[:lenY] {
    97  					y[i] = 0
    98  				}
    99  			} else {
   100  				c64.ScalUnitary(beta, y[:lenY])
   101  			}
   102  		} else {
   103  			iy := ky
   104  			if beta == 0 {
   105  				for i := 0; i < lenY; i++ {
   106  					y[iy] = 0
   107  					iy += incY
   108  				}
   109  			} else {
   110  				if incY > 0 {
   111  					c64.ScalInc(beta, y, uintptr(lenY), uintptr(incY))
   112  				} else {
   113  					c64.ScalInc(beta, y, uintptr(lenY), uintptr(-incY))
   114  				}
   115  			}
   116  		}
   117  	}
   118  
   119  	nRow := min(m, n+kL)
   120  	nCol := kL + 1 + kU
   121  	switch trans {
   122  	case blas.NoTrans:
   123  		iy := ky
   124  		if incX == 1 {
   125  			for i := 0; i < nRow; i++ {
   126  				l := max(0, kL-i)
   127  				u := min(nCol, n+kL-i)
   128  				aRow := a[i*lda+l : i*lda+u]
   129  				off := max(0, i-kL)
   130  				xtmp := x[off : off+u-l]
   131  				var sum complex64
   132  				for j, v := range aRow {
   133  					sum += xtmp[j] * v
   134  				}
   135  				y[iy] += alpha * sum
   136  				iy += incY
   137  			}
   138  		} else {
   139  			for i := 0; i < nRow; i++ {
   140  				l := max(0, kL-i)
   141  				u := min(nCol, n+kL-i)
   142  				aRow := a[i*lda+l : i*lda+u]
   143  				off := max(0, i-kL) * incX
   144  				jx := kx
   145  				var sum complex64
   146  				for _, v := range aRow {
   147  					sum += x[off+jx] * v
   148  					jx += incX
   149  				}
   150  				y[iy] += alpha * sum
   151  				iy += incY
   152  			}
   153  		}
   154  	case blas.Trans:
   155  		if incX == 1 {
   156  			for i := 0; i < nRow; i++ {
   157  				l := max(0, kL-i)
   158  				u := min(nCol, n+kL-i)
   159  				aRow := a[i*lda+l : i*lda+u]
   160  				off := max(0, i-kL) * incY
   161  				alphaxi := alpha * x[i]
   162  				jy := ky
   163  				for _, v := range aRow {
   164  					y[off+jy] += alphaxi * v
   165  					jy += incY
   166  				}
   167  			}
   168  		} else {
   169  			ix := kx
   170  			for i := 0; i < nRow; i++ {
   171  				l := max(0, kL-i)
   172  				u := min(nCol, n+kL-i)
   173  				aRow := a[i*lda+l : i*lda+u]
   174  				off := max(0, i-kL) * incY
   175  				alphaxi := alpha * x[ix]
   176  				jy := ky
   177  				for _, v := range aRow {
   178  					y[off+jy] += alphaxi * v
   179  					jy += incY
   180  				}
   181  				ix += incX
   182  			}
   183  		}
   184  	case blas.ConjTrans:
   185  		if incX == 1 {
   186  			for i := 0; i < nRow; i++ {
   187  				l := max(0, kL-i)
   188  				u := min(nCol, n+kL-i)
   189  				aRow := a[i*lda+l : i*lda+u]
   190  				off := max(0, i-kL) * incY
   191  				alphaxi := alpha * x[i]
   192  				jy := ky
   193  				for _, v := range aRow {
   194  					y[off+jy] += alphaxi * cmplx.Conj(v)
   195  					jy += incY
   196  				}
   197  			}
   198  		} else {
   199  			ix := kx
   200  			for i := 0; i < nRow; i++ {
   201  				l := max(0, kL-i)
   202  				u := min(nCol, n+kL-i)
   203  				aRow := a[i*lda+l : i*lda+u]
   204  				off := max(0, i-kL) * incY
   205  				alphaxi := alpha * x[ix]
   206  				jy := ky
   207  				for _, v := range aRow {
   208  					y[off+jy] += alphaxi * cmplx.Conj(v)
   209  					jy += incY
   210  				}
   211  				ix += incX
   212  			}
   213  		}
   214  	}
   215  }
   216  
   217  // Cgemv performs one of the matrix-vector operations
   218  //
   219  //	y = alpha * A * x + beta * y   if trans = blas.NoTrans
   220  //	y = alpha * Aᵀ * x + beta * y  if trans = blas.Trans
   221  //	y = alpha * Aᴴ * x + beta * y  if trans = blas.ConjTrans
   222  //
   223  // where alpha and beta are scalars, x and y are vectors, and A is an m×n dense matrix.
   224  //
   225  // Complex64 implementations are autogenerated and not directly tested.
   226  func (Implementation) Cgemv(trans blas.Transpose, m, n int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) {
   227  	switch trans {
   228  	default:
   229  		panic(badTranspose)
   230  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
   231  	}
   232  	if m < 0 {
   233  		panic(mLT0)
   234  	}
   235  	if n < 0 {
   236  		panic(nLT0)
   237  	}
   238  	if lda < max(1, n) {
   239  		panic(badLdA)
   240  	}
   241  	if incX == 0 {
   242  		panic(zeroIncX)
   243  	}
   244  	if incY == 0 {
   245  		panic(zeroIncY)
   246  	}
   247  
   248  	// Quick return if possible.
   249  	if m == 0 || n == 0 {
   250  		return
   251  	}
   252  
   253  	// For zero matrix size the following slice length checks are trivially satisfied.
   254  	var lenX, lenY int
   255  	if trans == blas.NoTrans {
   256  		lenX = n
   257  		lenY = m
   258  	} else {
   259  		lenX = m
   260  		lenY = n
   261  	}
   262  	if len(a) < lda*(m-1)+n {
   263  		panic(shortA)
   264  	}
   265  	if (incX > 0 && len(x) <= (lenX-1)*incX) || (incX < 0 && len(x) <= (1-lenX)*incX) {
   266  		panic(shortX)
   267  	}
   268  	if (incY > 0 && len(y) <= (lenY-1)*incY) || (incY < 0 && len(y) <= (1-lenY)*incY) {
   269  		panic(shortY)
   270  	}
   271  
   272  	// Quick return if possible.
   273  	if alpha == 0 && beta == 1 {
   274  		return
   275  	}
   276  
   277  	var kx int
   278  	if incX < 0 {
   279  		kx = (1 - lenX) * incX
   280  	}
   281  	var ky int
   282  	if incY < 0 {
   283  		ky = (1 - lenY) * incY
   284  	}
   285  
   286  	// Form y = beta*y.
   287  	if beta != 1 {
   288  		if incY == 1 {
   289  			if beta == 0 {
   290  				for i := range y[:lenY] {
   291  					y[i] = 0
   292  				}
   293  			} else {
   294  				c64.ScalUnitary(beta, y[:lenY])
   295  			}
   296  		} else {
   297  			iy := ky
   298  			if beta == 0 {
   299  				for i := 0; i < lenY; i++ {
   300  					y[iy] = 0
   301  					iy += incY
   302  				}
   303  			} else {
   304  				if incY > 0 {
   305  					c64.ScalInc(beta, y, uintptr(lenY), uintptr(incY))
   306  				} else {
   307  					c64.ScalInc(beta, y, uintptr(lenY), uintptr(-incY))
   308  				}
   309  			}
   310  		}
   311  	}
   312  
   313  	if alpha == 0 {
   314  		return
   315  	}
   316  
   317  	switch trans {
   318  	default:
   319  		// Form y = alpha*A*x + y.
   320  		iy := ky
   321  		if incX == 1 {
   322  			for i := 0; i < m; i++ {
   323  				y[iy] += alpha * c64.DotuUnitary(a[i*lda:i*lda+n], x[:n])
   324  				iy += incY
   325  			}
   326  			return
   327  		}
   328  		for i := 0; i < m; i++ {
   329  			y[iy] += alpha * c64.DotuInc(a[i*lda:i*lda+n], x, uintptr(n), 1, uintptr(incX), 0, uintptr(kx))
   330  			iy += incY
   331  		}
   332  		return
   333  
   334  	case blas.Trans:
   335  		// Form y = alpha*Aᵀ*x + y.
   336  		ix := kx
   337  		if incY == 1 {
   338  			for i := 0; i < m; i++ {
   339  				c64.AxpyUnitary(alpha*x[ix], a[i*lda:i*lda+n], y[:n])
   340  				ix += incX
   341  			}
   342  			return
   343  		}
   344  		for i := 0; i < m; i++ {
   345  			c64.AxpyInc(alpha*x[ix], a[i*lda:i*lda+n], y, uintptr(n), 1, uintptr(incY), 0, uintptr(ky))
   346  			ix += incX
   347  		}
   348  		return
   349  
   350  	case blas.ConjTrans:
   351  		// Form y = alpha*Aᴴ*x + y.
   352  		ix := kx
   353  		if incY == 1 {
   354  			for i := 0; i < m; i++ {
   355  				tmp := alpha * x[ix]
   356  				for j := 0; j < n; j++ {
   357  					y[j] += tmp * cmplx.Conj(a[i*lda+j])
   358  				}
   359  				ix += incX
   360  			}
   361  			return
   362  		}
   363  		for i := 0; i < m; i++ {
   364  			tmp := alpha * x[ix]
   365  			jy := ky
   366  			for j := 0; j < n; j++ {
   367  				y[jy] += tmp * cmplx.Conj(a[i*lda+j])
   368  				jy += incY
   369  			}
   370  			ix += incX
   371  		}
   372  		return
   373  	}
   374  }
   375  
   376  // Cgerc performs the rank-one operation
   377  //
   378  //	A += alpha * x * yᴴ
   379  //
   380  // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector,
   381  // and y is an n element vector.
   382  //
   383  // Complex64 implementations are autogenerated and not directly tested.
   384  func (Implementation) Cgerc(m, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, a []complex64, lda int) {
   385  	if m < 0 {
   386  		panic(mLT0)
   387  	}
   388  	if n < 0 {
   389  		panic(nLT0)
   390  	}
   391  	if lda < max(1, n) {
   392  		panic(badLdA)
   393  	}
   394  	if incX == 0 {
   395  		panic(zeroIncX)
   396  	}
   397  	if incY == 0 {
   398  		panic(zeroIncY)
   399  	}
   400  
   401  	// Quick return if possible.
   402  	if m == 0 || n == 0 {
   403  		return
   404  	}
   405  
   406  	// For zero matrix size the following slice length checks are trivially satisfied.
   407  	if (incX > 0 && len(x) <= (m-1)*incX) || (incX < 0 && len(x) <= (1-m)*incX) {
   408  		panic(shortX)
   409  	}
   410  	if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) {
   411  		panic(shortY)
   412  	}
   413  	if len(a) < lda*(m-1)+n {
   414  		panic(shortA)
   415  	}
   416  
   417  	// Quick return if possible.
   418  	if alpha == 0 {
   419  		return
   420  	}
   421  
   422  	var kx, jy int
   423  	if incX < 0 {
   424  		kx = (1 - m) * incX
   425  	}
   426  	if incY < 0 {
   427  		jy = (1 - n) * incY
   428  	}
   429  	for j := 0; j < n; j++ {
   430  		if y[jy] != 0 {
   431  			tmp := alpha * cmplx.Conj(y[jy])
   432  			c64.AxpyInc(tmp, x, a[j:], uintptr(m), uintptr(incX), uintptr(lda), uintptr(kx), 0)
   433  		}
   434  		jy += incY
   435  	}
   436  }
   437  
   438  // Cgeru performs the rank-one operation
   439  //
   440  //	A += alpha * x * yᵀ
   441  //
   442  // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector,
   443  // and y is an n element vector.
   444  //
   445  // Complex64 implementations are autogenerated and not directly tested.
   446  func (Implementation) Cgeru(m, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, a []complex64, lda int) {
   447  	if m < 0 {
   448  		panic(mLT0)
   449  	}
   450  	if n < 0 {
   451  		panic(nLT0)
   452  	}
   453  	if lda < max(1, n) {
   454  		panic(badLdA)
   455  	}
   456  	if incX == 0 {
   457  		panic(zeroIncX)
   458  	}
   459  	if incY == 0 {
   460  		panic(zeroIncY)
   461  	}
   462  
   463  	// Quick return if possible.
   464  	if m == 0 || n == 0 {
   465  		return
   466  	}
   467  
   468  	// For zero matrix size the following slice length checks are trivially satisfied.
   469  	if (incX > 0 && len(x) <= (m-1)*incX) || (incX < 0 && len(x) <= (1-m)*incX) {
   470  		panic(shortX)
   471  	}
   472  	if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) {
   473  		panic(shortY)
   474  	}
   475  	if len(a) < lda*(m-1)+n {
   476  		panic(shortA)
   477  	}
   478  
   479  	// Quick return if possible.
   480  	if alpha == 0 {
   481  		return
   482  	}
   483  
   484  	var kx int
   485  	if incX < 0 {
   486  		kx = (1 - m) * incX
   487  	}
   488  	if incY == 1 {
   489  		for i := 0; i < m; i++ {
   490  			if x[kx] != 0 {
   491  				tmp := alpha * x[kx]
   492  				c64.AxpyUnitary(tmp, y[:n], a[i*lda:i*lda+n])
   493  			}
   494  			kx += incX
   495  		}
   496  		return
   497  	}
   498  	var jy int
   499  	if incY < 0 {
   500  		jy = (1 - n) * incY
   501  	}
   502  	for i := 0; i < m; i++ {
   503  		if x[kx] != 0 {
   504  			tmp := alpha * x[kx]
   505  			c64.AxpyInc(tmp, y, a[i*lda:i*lda+n], uintptr(n), uintptr(incY), 1, uintptr(jy), 0)
   506  		}
   507  		kx += incX
   508  	}
   509  }
   510  
   511  // Chbmv performs the matrix-vector operation
   512  //
   513  //	y = alpha * A * x + beta * y
   514  //
   515  // where alpha and beta are scalars, x and y are vectors, and A is an n×n
   516  // Hermitian band matrix with k super-diagonals. The imaginary parts of
   517  // the diagonal elements of A are ignored and assumed to be zero.
   518  //
   519  // Complex64 implementations are autogenerated and not directly tested.
   520  func (Implementation) Chbmv(uplo blas.Uplo, n, k int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) {
   521  	switch uplo {
   522  	default:
   523  		panic(badUplo)
   524  	case blas.Upper, blas.Lower:
   525  	}
   526  	if n < 0 {
   527  		panic(nLT0)
   528  	}
   529  	if k < 0 {
   530  		panic(kLT0)
   531  	}
   532  	if lda < k+1 {
   533  		panic(badLdA)
   534  	}
   535  	if incX == 0 {
   536  		panic(zeroIncX)
   537  	}
   538  	if incY == 0 {
   539  		panic(zeroIncY)
   540  	}
   541  
   542  	// Quick return if possible.
   543  	if n == 0 {
   544  		return
   545  	}
   546  
   547  	// For zero matrix size the following slice length checks are trivially satisfied.
   548  	if len(a) < lda*(n-1)+k+1 {
   549  		panic(shortA)
   550  	}
   551  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
   552  		panic(shortX)
   553  	}
   554  	if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) {
   555  		panic(shortY)
   556  	}
   557  
   558  	// Quick return if possible.
   559  	if alpha == 0 && beta == 1 {
   560  		return
   561  	}
   562  
   563  	// Set up the start indices in X and Y.
   564  	var kx int
   565  	if incX < 0 {
   566  		kx = (1 - n) * incX
   567  	}
   568  	var ky int
   569  	if incY < 0 {
   570  		ky = (1 - n) * incY
   571  	}
   572  
   573  	// Form y = beta*y.
   574  	if beta != 1 {
   575  		if incY == 1 {
   576  			if beta == 0 {
   577  				for i := range y[:n] {
   578  					y[i] = 0
   579  				}
   580  			} else {
   581  				for i, v := range y[:n] {
   582  					y[i] = beta * v
   583  				}
   584  			}
   585  		} else {
   586  			iy := ky
   587  			if beta == 0 {
   588  				for i := 0; i < n; i++ {
   589  					y[iy] = 0
   590  					iy += incY
   591  				}
   592  			} else {
   593  				for i := 0; i < n; i++ {
   594  					y[iy] = beta * y[iy]
   595  					iy += incY
   596  				}
   597  			}
   598  		}
   599  	}
   600  
   601  	if alpha == 0 {
   602  		return
   603  	}
   604  
   605  	// The elements of A are accessed sequentially with one pass through a.
   606  	switch uplo {
   607  	case blas.Upper:
   608  		iy := ky
   609  		if incX == 1 {
   610  			for i := 0; i < n; i++ {
   611  				aRow := a[i*lda:]
   612  				alphaxi := alpha * x[i]
   613  				sum := alphaxi * complex(real(aRow[0]), 0)
   614  				u := min(k+1, n-i)
   615  				jy := incY
   616  				for j := 1; j < u; j++ {
   617  					v := aRow[j]
   618  					sum += alpha * x[i+j] * v
   619  					y[iy+jy] += alphaxi * cmplx.Conj(v)
   620  					jy += incY
   621  				}
   622  				y[iy] += sum
   623  				iy += incY
   624  			}
   625  		} else {
   626  			ix := kx
   627  			for i := 0; i < n; i++ {
   628  				aRow := a[i*lda:]
   629  				alphaxi := alpha * x[ix]
   630  				sum := alphaxi * complex(real(aRow[0]), 0)
   631  				u := min(k+1, n-i)
   632  				jx := incX
   633  				jy := incY
   634  				for j := 1; j < u; j++ {
   635  					v := aRow[j]
   636  					sum += alpha * x[ix+jx] * v
   637  					y[iy+jy] += alphaxi * cmplx.Conj(v)
   638  					jx += incX
   639  					jy += incY
   640  				}
   641  				y[iy] += sum
   642  				ix += incX
   643  				iy += incY
   644  			}
   645  		}
   646  	case blas.Lower:
   647  		iy := ky
   648  		if incX == 1 {
   649  			for i := 0; i < n; i++ {
   650  				l := max(0, k-i)
   651  				alphaxi := alpha * x[i]
   652  				jy := l * incY
   653  				aRow := a[i*lda:]
   654  				for j := l; j < k; j++ {
   655  					v := aRow[j]
   656  					y[iy] += alpha * v * x[i-k+j]
   657  					y[iy-k*incY+jy] += alphaxi * cmplx.Conj(v)
   658  					jy += incY
   659  				}
   660  				y[iy] += alphaxi * complex(real(aRow[k]), 0)
   661  				iy += incY
   662  			}
   663  		} else {
   664  			ix := kx
   665  			for i := 0; i < n; i++ {
   666  				l := max(0, k-i)
   667  				alphaxi := alpha * x[ix]
   668  				jx := l * incX
   669  				jy := l * incY
   670  				aRow := a[i*lda:]
   671  				for j := l; j < k; j++ {
   672  					v := aRow[j]
   673  					y[iy] += alpha * v * x[ix-k*incX+jx]
   674  					y[iy-k*incY+jy] += alphaxi * cmplx.Conj(v)
   675  					jx += incX
   676  					jy += incY
   677  				}
   678  				y[iy] += alphaxi * complex(real(aRow[k]), 0)
   679  				ix += incX
   680  				iy += incY
   681  			}
   682  		}
   683  	}
   684  }
   685  
   686  // Chemv performs the matrix-vector operation
   687  //
   688  //	y = alpha * A * x + beta * y
   689  //
   690  // where alpha and beta are scalars, x and y are vectors, and A is an n×n
   691  // Hermitian matrix. The imaginary parts of the diagonal elements of A are
   692  // ignored and assumed to be zero.
   693  //
   694  // Complex64 implementations are autogenerated and not directly tested.
   695  func (Implementation) Chemv(uplo blas.Uplo, n int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) {
   696  	switch uplo {
   697  	default:
   698  		panic(badUplo)
   699  	case blas.Upper, blas.Lower:
   700  	}
   701  	if n < 0 {
   702  		panic(nLT0)
   703  	}
   704  	if lda < max(1, n) {
   705  		panic(badLdA)
   706  	}
   707  	if incX == 0 {
   708  		panic(zeroIncX)
   709  	}
   710  	if incY == 0 {
   711  		panic(zeroIncY)
   712  	}
   713  
   714  	// Quick return if possible.
   715  	if n == 0 {
   716  		return
   717  	}
   718  
   719  	// For zero matrix size the following slice length checks are trivially satisfied.
   720  	if len(a) < lda*(n-1)+n {
   721  		panic(shortA)
   722  	}
   723  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
   724  		panic(shortX)
   725  	}
   726  	if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) {
   727  		panic(shortY)
   728  	}
   729  
   730  	// Quick return if possible.
   731  	if alpha == 0 && beta == 1 {
   732  		return
   733  	}
   734  
   735  	// Set up the start indices in X and Y.
   736  	var kx int
   737  	if incX < 0 {
   738  		kx = (1 - n) * incX
   739  	}
   740  	var ky int
   741  	if incY < 0 {
   742  		ky = (1 - n) * incY
   743  	}
   744  
   745  	// Form y = beta*y.
   746  	if beta != 1 {
   747  		if incY == 1 {
   748  			if beta == 0 {
   749  				for i := range y[:n] {
   750  					y[i] = 0
   751  				}
   752  			} else {
   753  				for i, v := range y[:n] {
   754  					y[i] = beta * v
   755  				}
   756  			}
   757  		} else {
   758  			iy := ky
   759  			if beta == 0 {
   760  				for i := 0; i < n; i++ {
   761  					y[iy] = 0
   762  					iy += incY
   763  				}
   764  			} else {
   765  				for i := 0; i < n; i++ {
   766  					y[iy] = beta * y[iy]
   767  					iy += incY
   768  				}
   769  			}
   770  		}
   771  	}
   772  
   773  	if alpha == 0 {
   774  		return
   775  	}
   776  
   777  	// The elements of A are accessed sequentially with one pass through
   778  	// the triangular part of A.
   779  
   780  	if uplo == blas.Upper {
   781  		// Form y when A is stored in upper triangle.
   782  		if incX == 1 && incY == 1 {
   783  			for i := 0; i < n; i++ {
   784  				tmp1 := alpha * x[i]
   785  				var tmp2 complex64
   786  				for j := i + 1; j < n; j++ {
   787  					y[j] += tmp1 * cmplx.Conj(a[i*lda+j])
   788  					tmp2 += a[i*lda+j] * x[j]
   789  				}
   790  				aii := complex(real(a[i*lda+i]), 0)
   791  				y[i] += tmp1*aii + alpha*tmp2
   792  			}
   793  		} else {
   794  			ix := kx
   795  			iy := ky
   796  			for i := 0; i < n; i++ {
   797  				tmp1 := alpha * x[ix]
   798  				var tmp2 complex64
   799  				jx := ix
   800  				jy := iy
   801  				for j := i + 1; j < n; j++ {
   802  					jx += incX
   803  					jy += incY
   804  					y[jy] += tmp1 * cmplx.Conj(a[i*lda+j])
   805  					tmp2 += a[i*lda+j] * x[jx]
   806  				}
   807  				aii := complex(real(a[i*lda+i]), 0)
   808  				y[iy] += tmp1*aii + alpha*tmp2
   809  				ix += incX
   810  				iy += incY
   811  			}
   812  		}
   813  		return
   814  	}
   815  
   816  	// Form y when A is stored in lower triangle.
   817  	if incX == 1 && incY == 1 {
   818  		for i := 0; i < n; i++ {
   819  			tmp1 := alpha * x[i]
   820  			var tmp2 complex64
   821  			for j := 0; j < i; j++ {
   822  				y[j] += tmp1 * cmplx.Conj(a[i*lda+j])
   823  				tmp2 += a[i*lda+j] * x[j]
   824  			}
   825  			aii := complex(real(a[i*lda+i]), 0)
   826  			y[i] += tmp1*aii + alpha*tmp2
   827  		}
   828  	} else {
   829  		ix := kx
   830  		iy := ky
   831  		for i := 0; i < n; i++ {
   832  			tmp1 := alpha * x[ix]
   833  			var tmp2 complex64
   834  			jx := kx
   835  			jy := ky
   836  			for j := 0; j < i; j++ {
   837  				y[jy] += tmp1 * cmplx.Conj(a[i*lda+j])
   838  				tmp2 += a[i*lda+j] * x[jx]
   839  				jx += incX
   840  				jy += incY
   841  			}
   842  			aii := complex(real(a[i*lda+i]), 0)
   843  			y[iy] += tmp1*aii + alpha*tmp2
   844  			ix += incX
   845  			iy += incY
   846  		}
   847  	}
   848  }
   849  
   850  // Cher performs the Hermitian rank-one operation
   851  //
   852  //	A += alpha * x * xᴴ
   853  //
   854  // where A is an n×n Hermitian matrix, alpha is a real scalar, and x is an n
   855  // element vector. On entry, the imaginary parts of the diagonal elements of A
   856  // are ignored and assumed to be zero, on return they will be set to zero.
   857  //
   858  // Complex64 implementations are autogenerated and not directly tested.
   859  func (Implementation) Cher(uplo blas.Uplo, n int, alpha float32, x []complex64, incX int, a []complex64, lda int) {
   860  	switch uplo {
   861  	default:
   862  		panic(badUplo)
   863  	case blas.Upper, blas.Lower:
   864  	}
   865  	if n < 0 {
   866  		panic(nLT0)
   867  	}
   868  	if lda < max(1, n) {
   869  		panic(badLdA)
   870  	}
   871  	if incX == 0 {
   872  		panic(zeroIncX)
   873  	}
   874  
   875  	// Quick return if possible.
   876  	if n == 0 {
   877  		return
   878  	}
   879  
   880  	// For zero matrix size the following slice length checks are trivially satisfied.
   881  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
   882  		panic(shortX)
   883  	}
   884  	if len(a) < lda*(n-1)+n {
   885  		panic(shortA)
   886  	}
   887  
   888  	// Quick return if possible.
   889  	if alpha == 0 {
   890  		return
   891  	}
   892  
   893  	var kx int
   894  	if incX < 0 {
   895  		kx = (1 - n) * incX
   896  	}
   897  	if uplo == blas.Upper {
   898  		if incX == 1 {
   899  			for i := 0; i < n; i++ {
   900  				if x[i] != 0 {
   901  					tmp := complex(alpha*real(x[i]), alpha*imag(x[i]))
   902  					aii := real(a[i*lda+i])
   903  					xtmp := real(tmp * cmplx.Conj(x[i]))
   904  					a[i*lda+i] = complex(aii+xtmp, 0)
   905  					for j := i + 1; j < n; j++ {
   906  						a[i*lda+j] += tmp * cmplx.Conj(x[j])
   907  					}
   908  				} else {
   909  					aii := real(a[i*lda+i])
   910  					a[i*lda+i] = complex(aii, 0)
   911  				}
   912  			}
   913  			return
   914  		}
   915  
   916  		ix := kx
   917  		for i := 0; i < n; i++ {
   918  			if x[ix] != 0 {
   919  				tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix]))
   920  				aii := real(a[i*lda+i])
   921  				xtmp := real(tmp * cmplx.Conj(x[ix]))
   922  				a[i*lda+i] = complex(aii+xtmp, 0)
   923  				jx := ix + incX
   924  				for j := i + 1; j < n; j++ {
   925  					a[i*lda+j] += tmp * cmplx.Conj(x[jx])
   926  					jx += incX
   927  				}
   928  			} else {
   929  				aii := real(a[i*lda+i])
   930  				a[i*lda+i] = complex(aii, 0)
   931  			}
   932  			ix += incX
   933  		}
   934  		return
   935  	}
   936  
   937  	if incX == 1 {
   938  		for i := 0; i < n; i++ {
   939  			if x[i] != 0 {
   940  				tmp := complex(alpha*real(x[i]), alpha*imag(x[i]))
   941  				for j := 0; j < i; j++ {
   942  					a[i*lda+j] += tmp * cmplx.Conj(x[j])
   943  				}
   944  				aii := real(a[i*lda+i])
   945  				xtmp := real(tmp * cmplx.Conj(x[i]))
   946  				a[i*lda+i] = complex(aii+xtmp, 0)
   947  			} else {
   948  				aii := real(a[i*lda+i])
   949  				a[i*lda+i] = complex(aii, 0)
   950  			}
   951  		}
   952  		return
   953  	}
   954  
   955  	ix := kx
   956  	for i := 0; i < n; i++ {
   957  		if x[ix] != 0 {
   958  			tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix]))
   959  			jx := kx
   960  			for j := 0; j < i; j++ {
   961  				a[i*lda+j] += tmp * cmplx.Conj(x[jx])
   962  				jx += incX
   963  			}
   964  			aii := real(a[i*lda+i])
   965  			xtmp := real(tmp * cmplx.Conj(x[ix]))
   966  			a[i*lda+i] = complex(aii+xtmp, 0)
   967  
   968  		} else {
   969  			aii := real(a[i*lda+i])
   970  			a[i*lda+i] = complex(aii, 0)
   971  		}
   972  		ix += incX
   973  	}
   974  }
   975  
   976  // Cher2 performs the Hermitian rank-two operation
   977  //
   978  //	A += alpha * x * yᴴ + conj(alpha) * y * xᴴ
   979  //
   980  // where alpha is a scalar, x and y are n element vectors and A is an n×n
   981  // Hermitian matrix. On entry, the imaginary parts of the diagonal elements are
   982  // ignored and assumed to be zero. On return they will be set to zero.
   983  //
   984  // Complex64 implementations are autogenerated and not directly tested.
   985  func (Implementation) Cher2(uplo blas.Uplo, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, a []complex64, lda int) {
   986  	switch uplo {
   987  	default:
   988  		panic(badUplo)
   989  	case blas.Upper, blas.Lower:
   990  	}
   991  	if n < 0 {
   992  		panic(nLT0)
   993  	}
   994  	if lda < max(1, n) {
   995  		panic(badLdA)
   996  	}
   997  	if incX == 0 {
   998  		panic(zeroIncX)
   999  	}
  1000  	if incY == 0 {
  1001  		panic(zeroIncY)
  1002  	}
  1003  
  1004  	// Quick return if possible.
  1005  	if n == 0 {
  1006  		return
  1007  	}
  1008  
  1009  	// For zero matrix size the following slice length checks are trivially satisfied.
  1010  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  1011  		panic(shortX)
  1012  	}
  1013  	if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) {
  1014  		panic(shortY)
  1015  	}
  1016  	if len(a) < lda*(n-1)+n {
  1017  		panic(shortA)
  1018  	}
  1019  
  1020  	// Quick return if possible.
  1021  	if alpha == 0 {
  1022  		return
  1023  	}
  1024  
  1025  	var kx, ky int
  1026  	var ix, iy int
  1027  	if incX != 1 || incY != 1 {
  1028  		if incX < 0 {
  1029  			kx = (1 - n) * incX
  1030  		}
  1031  		if incY < 0 {
  1032  			ky = (1 - n) * incY
  1033  		}
  1034  		ix = kx
  1035  		iy = ky
  1036  	}
  1037  	if uplo == blas.Upper {
  1038  		if incX == 1 && incY == 1 {
  1039  			for i := 0; i < n; i++ {
  1040  				if x[i] != 0 || y[i] != 0 {
  1041  					tmp1 := alpha * x[i]
  1042  					tmp2 := cmplx.Conj(alpha) * y[i]
  1043  					aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
  1044  					a[i*lda+i] = complex(aii, 0)
  1045  					for j := i + 1; j < n; j++ {
  1046  						a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
  1047  					}
  1048  				} else {
  1049  					aii := real(a[i*lda+i])
  1050  					a[i*lda+i] = complex(aii, 0)
  1051  				}
  1052  			}
  1053  			return
  1054  		}
  1055  		for i := 0; i < n; i++ {
  1056  			if x[ix] != 0 || y[iy] != 0 {
  1057  				tmp1 := alpha * x[ix]
  1058  				tmp2 := cmplx.Conj(alpha) * y[iy]
  1059  				aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
  1060  				a[i*lda+i] = complex(aii, 0)
  1061  				jx := ix + incX
  1062  				jy := iy + incY
  1063  				for j := i + 1; j < n; j++ {
  1064  					a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
  1065  					jx += incX
  1066  					jy += incY
  1067  				}
  1068  			} else {
  1069  				aii := real(a[i*lda+i])
  1070  				a[i*lda+i] = complex(aii, 0)
  1071  			}
  1072  			ix += incX
  1073  			iy += incY
  1074  		}
  1075  		return
  1076  	}
  1077  
  1078  	if incX == 1 && incY == 1 {
  1079  		for i := 0; i < n; i++ {
  1080  			if x[i] != 0 || y[i] != 0 {
  1081  				tmp1 := alpha * x[i]
  1082  				tmp2 := cmplx.Conj(alpha) * y[i]
  1083  				for j := 0; j < i; j++ {
  1084  					a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
  1085  				}
  1086  				aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
  1087  				a[i*lda+i] = complex(aii, 0)
  1088  			} else {
  1089  				aii := real(a[i*lda+i])
  1090  				a[i*lda+i] = complex(aii, 0)
  1091  			}
  1092  		}
  1093  		return
  1094  	}
  1095  	for i := 0; i < n; i++ {
  1096  		if x[ix] != 0 || y[iy] != 0 {
  1097  			tmp1 := alpha * x[ix]
  1098  			tmp2 := cmplx.Conj(alpha) * y[iy]
  1099  			jx := kx
  1100  			jy := ky
  1101  			for j := 0; j < i; j++ {
  1102  				a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
  1103  				jx += incX
  1104  				jy += incY
  1105  			}
  1106  			aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
  1107  			a[i*lda+i] = complex(aii, 0)
  1108  		} else {
  1109  			aii := real(a[i*lda+i])
  1110  			a[i*lda+i] = complex(aii, 0)
  1111  		}
  1112  		ix += incX
  1113  		iy += incY
  1114  	}
  1115  }
  1116  
  1117  // Chpmv performs the matrix-vector operation
  1118  //
  1119  //	y = alpha * A * x + beta * y
  1120  //
  1121  // where alpha and beta are scalars, x and y are vectors, and A is an n×n
  1122  // Hermitian matrix in packed form. The imaginary parts of the diagonal
  1123  // elements of A are ignored and assumed to be zero.
  1124  //
  1125  // Complex64 implementations are autogenerated and not directly tested.
  1126  func (Implementation) Chpmv(uplo blas.Uplo, n int, alpha complex64, ap []complex64, x []complex64, incX int, beta complex64, y []complex64, incY int) {
  1127  	switch uplo {
  1128  	default:
  1129  		panic(badUplo)
  1130  	case blas.Upper, blas.Lower:
  1131  	}
  1132  	if n < 0 {
  1133  		panic(nLT0)
  1134  	}
  1135  	if incX == 0 {
  1136  		panic(zeroIncX)
  1137  	}
  1138  	if incY == 0 {
  1139  		panic(zeroIncY)
  1140  	}
  1141  
  1142  	// Quick return if possible.
  1143  	if n == 0 {
  1144  		return
  1145  	}
  1146  
  1147  	// For zero matrix size the following slice length checks are trivially satisfied.
  1148  	if len(ap) < n*(n+1)/2 {
  1149  		panic(shortAP)
  1150  	}
  1151  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  1152  		panic(shortX)
  1153  	}
  1154  	if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) {
  1155  		panic(shortY)
  1156  	}
  1157  
  1158  	// Quick return if possible.
  1159  	if alpha == 0 && beta == 1 {
  1160  		return
  1161  	}
  1162  
  1163  	// Set up the start indices in X and Y.
  1164  	var kx int
  1165  	if incX < 0 {
  1166  		kx = (1 - n) * incX
  1167  	}
  1168  	var ky int
  1169  	if incY < 0 {
  1170  		ky = (1 - n) * incY
  1171  	}
  1172  
  1173  	// Form y = beta*y.
  1174  	if beta != 1 {
  1175  		if incY == 1 {
  1176  			if beta == 0 {
  1177  				for i := range y[:n] {
  1178  					y[i] = 0
  1179  				}
  1180  			} else {
  1181  				for i, v := range y[:n] {
  1182  					y[i] = beta * v
  1183  				}
  1184  			}
  1185  		} else {
  1186  			iy := ky
  1187  			if beta == 0 {
  1188  				for i := 0; i < n; i++ {
  1189  					y[iy] = 0
  1190  					iy += incY
  1191  				}
  1192  			} else {
  1193  				for i := 0; i < n; i++ {
  1194  					y[iy] *= beta
  1195  					iy += incY
  1196  				}
  1197  			}
  1198  		}
  1199  	}
  1200  
  1201  	if alpha == 0 {
  1202  		return
  1203  	}
  1204  
  1205  	// The elements of A are accessed sequentially with one pass through ap.
  1206  
  1207  	var kk int
  1208  	if uplo == blas.Upper {
  1209  		// Form y when ap contains the upper triangle.
  1210  		// Here, kk points to the current diagonal element in ap.
  1211  		if incX == 1 && incY == 1 {
  1212  			for i := 0; i < n; i++ {
  1213  				tmp1 := alpha * x[i]
  1214  				y[i] += tmp1 * complex(real(ap[kk]), 0)
  1215  				var tmp2 complex64
  1216  				k := kk + 1
  1217  				for j := i + 1; j < n; j++ {
  1218  					y[j] += tmp1 * cmplx.Conj(ap[k])
  1219  					tmp2 += ap[k] * x[j]
  1220  					k++
  1221  				}
  1222  				y[i] += alpha * tmp2
  1223  				kk += n - i
  1224  			}
  1225  		} else {
  1226  			ix := kx
  1227  			iy := ky
  1228  			for i := 0; i < n; i++ {
  1229  				tmp1 := alpha * x[ix]
  1230  				y[iy] += tmp1 * complex(real(ap[kk]), 0)
  1231  				var tmp2 complex64
  1232  				jx := ix
  1233  				jy := iy
  1234  				for k := kk + 1; k < kk+n-i; k++ {
  1235  					jx += incX
  1236  					jy += incY
  1237  					y[jy] += tmp1 * cmplx.Conj(ap[k])
  1238  					tmp2 += ap[k] * x[jx]
  1239  				}
  1240  				y[iy] += alpha * tmp2
  1241  				ix += incX
  1242  				iy += incY
  1243  				kk += n - i
  1244  			}
  1245  		}
  1246  		return
  1247  	}
  1248  
  1249  	// Form y when ap contains the lower triangle.
  1250  	// Here, kk points to the beginning of current row in ap.
  1251  	if incX == 1 && incY == 1 {
  1252  		for i := 0; i < n; i++ {
  1253  			tmp1 := alpha * x[i]
  1254  			var tmp2 complex64
  1255  			k := kk
  1256  			for j := 0; j < i; j++ {
  1257  				y[j] += tmp1 * cmplx.Conj(ap[k])
  1258  				tmp2 += ap[k] * x[j]
  1259  				k++
  1260  			}
  1261  			aii := complex(real(ap[kk+i]), 0)
  1262  			y[i] += tmp1*aii + alpha*tmp2
  1263  			kk += i + 1
  1264  		}
  1265  	} else {
  1266  		ix := kx
  1267  		iy := ky
  1268  		for i := 0; i < n; i++ {
  1269  			tmp1 := alpha * x[ix]
  1270  			var tmp2 complex64
  1271  			jx := kx
  1272  			jy := ky
  1273  			for k := kk; k < kk+i; k++ {
  1274  				y[jy] += tmp1 * cmplx.Conj(ap[k])
  1275  				tmp2 += ap[k] * x[jx]
  1276  				jx += incX
  1277  				jy += incY
  1278  			}
  1279  			aii := complex(real(ap[kk+i]), 0)
  1280  			y[iy] += tmp1*aii + alpha*tmp2
  1281  			ix += incX
  1282  			iy += incY
  1283  			kk += i + 1
  1284  		}
  1285  	}
  1286  }
  1287  
  1288  // Chpr performs the Hermitian rank-1 operation
  1289  //
  1290  //	A += alpha * x * xᴴ
  1291  //
  1292  // where alpha is a real scalar, x is a vector, and A is an n×n hermitian matrix
  1293  // in packed form. On entry, the imaginary parts of the diagonal elements are
  1294  // assumed to be zero, and on return they are set to zero.
  1295  //
  1296  // Complex64 implementations are autogenerated and not directly tested.
  1297  func (Implementation) Chpr(uplo blas.Uplo, n int, alpha float32, x []complex64, incX int, ap []complex64) {
  1298  	switch uplo {
  1299  	default:
  1300  		panic(badUplo)
  1301  	case blas.Upper, blas.Lower:
  1302  	}
  1303  	if n < 0 {
  1304  		panic(nLT0)
  1305  	}
  1306  	if incX == 0 {
  1307  		panic(zeroIncX)
  1308  	}
  1309  
  1310  	// Quick return if possible.
  1311  	if n == 0 {
  1312  		return
  1313  	}
  1314  
  1315  	// For zero matrix size the following slice length checks are trivially satisfied.
  1316  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  1317  		panic(shortX)
  1318  	}
  1319  	if len(ap) < n*(n+1)/2 {
  1320  		panic(shortAP)
  1321  	}
  1322  
  1323  	// Quick return if possible.
  1324  	if alpha == 0 {
  1325  		return
  1326  	}
  1327  
  1328  	// Set up start index in X.
  1329  	var kx int
  1330  	if incX < 0 {
  1331  		kx = (1 - n) * incX
  1332  	}
  1333  
  1334  	// The elements of A are accessed sequentially with one pass through ap.
  1335  
  1336  	var kk int
  1337  	if uplo == blas.Upper {
  1338  		// Form A when upper triangle is stored in AP.
  1339  		// Here, kk points to the current diagonal element in ap.
  1340  		if incX == 1 {
  1341  			for i := 0; i < n; i++ {
  1342  				xi := x[i]
  1343  				if xi != 0 {
  1344  					aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi)
  1345  					ap[kk] = complex(aii, 0)
  1346  
  1347  					tmp := complex(alpha, 0) * xi
  1348  					a := ap[kk+1 : kk+n-i]
  1349  					x := x[i+1 : n]
  1350  					for j, v := range x {
  1351  						a[j] += tmp * cmplx.Conj(v)
  1352  					}
  1353  				} else {
  1354  					ap[kk] = complex(real(ap[kk]), 0)
  1355  				}
  1356  				kk += n - i
  1357  			}
  1358  		} else {
  1359  			ix := kx
  1360  			for i := 0; i < n; i++ {
  1361  				xi := x[ix]
  1362  				if xi != 0 {
  1363  					aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi)
  1364  					ap[kk] = complex(aii, 0)
  1365  
  1366  					tmp := complex(alpha, 0) * xi
  1367  					jx := ix + incX
  1368  					a := ap[kk+1 : kk+n-i]
  1369  					for k := range a {
  1370  						a[k] += tmp * cmplx.Conj(x[jx])
  1371  						jx += incX
  1372  					}
  1373  				} else {
  1374  					ap[kk] = complex(real(ap[kk]), 0)
  1375  				}
  1376  				ix += incX
  1377  				kk += n - i
  1378  			}
  1379  		}
  1380  		return
  1381  	}
  1382  
  1383  	// Form A when lower triangle is stored in AP.
  1384  	// Here, kk points to the beginning of current row in ap.
  1385  	if incX == 1 {
  1386  		for i := 0; i < n; i++ {
  1387  			xi := x[i]
  1388  			if xi != 0 {
  1389  				tmp := complex(alpha, 0) * xi
  1390  				a := ap[kk : kk+i]
  1391  				for j, v := range x[:i] {
  1392  					a[j] += tmp * cmplx.Conj(v)
  1393  				}
  1394  
  1395  				aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi)
  1396  				ap[kk+i] = complex(aii, 0)
  1397  			} else {
  1398  				ap[kk+i] = complex(real(ap[kk+i]), 0)
  1399  			}
  1400  			kk += i + 1
  1401  		}
  1402  	} else {
  1403  		ix := kx
  1404  		for i := 0; i < n; i++ {
  1405  			xi := x[ix]
  1406  			if xi != 0 {
  1407  				tmp := complex(alpha, 0) * xi
  1408  				a := ap[kk : kk+i]
  1409  				jx := kx
  1410  				for k := range a {
  1411  					a[k] += tmp * cmplx.Conj(x[jx])
  1412  					jx += incX
  1413  				}
  1414  
  1415  				aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi)
  1416  				ap[kk+i] = complex(aii, 0)
  1417  			} else {
  1418  				ap[kk+i] = complex(real(ap[kk+i]), 0)
  1419  			}
  1420  			ix += incX
  1421  			kk += i + 1
  1422  		}
  1423  	}
  1424  }
  1425  
  1426  // Chpr2 performs the Hermitian rank-2 operation
  1427  //
  1428  //	A += alpha * x * yᴴ + conj(alpha) * y * xᴴ
  1429  //
  1430  // where alpha is a complex scalar, x and y are n element vectors, and A is an
  1431  // n×n Hermitian matrix, supplied in packed form. On entry, the imaginary parts
  1432  // of the diagonal elements are assumed to be zero, and on return they are set to zero.
  1433  //
  1434  // Complex64 implementations are autogenerated and not directly tested.
  1435  func (Implementation) Chpr2(uplo blas.Uplo, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, ap []complex64) {
  1436  	switch uplo {
  1437  	default:
  1438  		panic(badUplo)
  1439  	case blas.Upper, blas.Lower:
  1440  	}
  1441  	if n < 0 {
  1442  		panic(nLT0)
  1443  	}
  1444  	if incX == 0 {
  1445  		panic(zeroIncX)
  1446  	}
  1447  	if incY == 0 {
  1448  		panic(zeroIncY)
  1449  	}
  1450  
  1451  	// Quick return if possible.
  1452  	if n == 0 {
  1453  		return
  1454  	}
  1455  
  1456  	// For zero matrix size the following slice length checks are trivially satisfied.
  1457  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  1458  		panic(shortX)
  1459  	}
  1460  	if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) {
  1461  		panic(shortY)
  1462  	}
  1463  	if len(ap) < n*(n+1)/2 {
  1464  		panic(shortAP)
  1465  	}
  1466  
  1467  	// Quick return if possible.
  1468  	if alpha == 0 {
  1469  		return
  1470  	}
  1471  
  1472  	// Set up start indices in X and Y.
  1473  	var kx int
  1474  	if incX < 0 {
  1475  		kx = (1 - n) * incX
  1476  	}
  1477  	var ky int
  1478  	if incY < 0 {
  1479  		ky = (1 - n) * incY
  1480  	}
  1481  
  1482  	// The elements of A are accessed sequentially with one pass through ap.
  1483  
  1484  	var kk int
  1485  	if uplo == blas.Upper {
  1486  		// Form A when upper triangle is stored in AP.
  1487  		// Here, kk points to the current diagonal element in ap.
  1488  		if incX == 1 && incY == 1 {
  1489  			for i := 0; i < n; i++ {
  1490  				if x[i] != 0 || y[i] != 0 {
  1491  					tmp1 := alpha * x[i]
  1492  					tmp2 := cmplx.Conj(alpha) * y[i]
  1493  					aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
  1494  					ap[kk] = complex(aii, 0)
  1495  					k := kk + 1
  1496  					for j := i + 1; j < n; j++ {
  1497  						ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
  1498  						k++
  1499  					}
  1500  				} else {
  1501  					ap[kk] = complex(real(ap[kk]), 0)
  1502  				}
  1503  				kk += n - i
  1504  			}
  1505  		} else {
  1506  			ix := kx
  1507  			iy := ky
  1508  			for i := 0; i < n; i++ {
  1509  				if x[ix] != 0 || y[iy] != 0 {
  1510  					tmp1 := alpha * x[ix]
  1511  					tmp2 := cmplx.Conj(alpha) * y[iy]
  1512  					aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
  1513  					ap[kk] = complex(aii, 0)
  1514  					jx := ix + incX
  1515  					jy := iy + incY
  1516  					for k := kk + 1; k < kk+n-i; k++ {
  1517  						ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
  1518  						jx += incX
  1519  						jy += incY
  1520  					}
  1521  				} else {
  1522  					ap[kk] = complex(real(ap[kk]), 0)
  1523  				}
  1524  				ix += incX
  1525  				iy += incY
  1526  				kk += n - i
  1527  			}
  1528  		}
  1529  		return
  1530  	}
  1531  
  1532  	// Form A when lower triangle is stored in AP.
  1533  	// Here, kk points to the beginning of current row in ap.
  1534  	if incX == 1 && incY == 1 {
  1535  		for i := 0; i < n; i++ {
  1536  			if x[i] != 0 || y[i] != 0 {
  1537  				tmp1 := alpha * x[i]
  1538  				tmp2 := cmplx.Conj(alpha) * y[i]
  1539  				k := kk
  1540  				for j := 0; j < i; j++ {
  1541  					ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
  1542  					k++
  1543  				}
  1544  				aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
  1545  				ap[kk+i] = complex(aii, 0)
  1546  			} else {
  1547  				ap[kk+i] = complex(real(ap[kk+i]), 0)
  1548  			}
  1549  			kk += i + 1
  1550  		}
  1551  	} else {
  1552  		ix := kx
  1553  		iy := ky
  1554  		for i := 0; i < n; i++ {
  1555  			if x[ix] != 0 || y[iy] != 0 {
  1556  				tmp1 := alpha * x[ix]
  1557  				tmp2 := cmplx.Conj(alpha) * y[iy]
  1558  				jx := kx
  1559  				jy := ky
  1560  				for k := kk; k < kk+i; k++ {
  1561  					ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
  1562  					jx += incX
  1563  					jy += incY
  1564  				}
  1565  				aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
  1566  				ap[kk+i] = complex(aii, 0)
  1567  			} else {
  1568  				ap[kk+i] = complex(real(ap[kk+i]), 0)
  1569  			}
  1570  			ix += incX
  1571  			iy += incY
  1572  			kk += i + 1
  1573  		}
  1574  	}
  1575  }
  1576  
  1577  // Ctbmv performs one of the matrix-vector operations
  1578  //
  1579  //	x = A * x   if trans = blas.NoTrans
  1580  //	x = Aᵀ * x  if trans = blas.Trans
  1581  //	x = Aᴴ * x  if trans = blas.ConjTrans
  1582  //
  1583  // where x is an n element vector and A is an n×n triangular band matrix, with
  1584  // (k+1) diagonals.
  1585  //
  1586  // Complex64 implementations are autogenerated and not directly tested.
  1587  func (Implementation) Ctbmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, k int, a []complex64, lda int, x []complex64, incX int) {
  1588  	switch trans {
  1589  	default:
  1590  		panic(badTranspose)
  1591  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
  1592  	}
  1593  	switch uplo {
  1594  	default:
  1595  		panic(badUplo)
  1596  	case blas.Upper, blas.Lower:
  1597  	}
  1598  	switch diag {
  1599  	default:
  1600  		panic(badDiag)
  1601  	case blas.NonUnit, blas.Unit:
  1602  	}
  1603  	if n < 0 {
  1604  		panic(nLT0)
  1605  	}
  1606  	if k < 0 {
  1607  		panic(kLT0)
  1608  	}
  1609  	if lda < k+1 {
  1610  		panic(badLdA)
  1611  	}
  1612  	if incX == 0 {
  1613  		panic(zeroIncX)
  1614  	}
  1615  
  1616  	// Quick return if possible.
  1617  	if n == 0 {
  1618  		return
  1619  	}
  1620  
  1621  	// For zero matrix size the following slice length checks are trivially satisfied.
  1622  	if len(a) < lda*(n-1)+k+1 {
  1623  		panic(shortA)
  1624  	}
  1625  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  1626  		panic(shortX)
  1627  	}
  1628  
  1629  	// Set up start index in X.
  1630  	var kx int
  1631  	if incX < 0 {
  1632  		kx = (1 - n) * incX
  1633  	}
  1634  
  1635  	switch trans {
  1636  	case blas.NoTrans:
  1637  		if uplo == blas.Upper {
  1638  			if incX == 1 {
  1639  				for i := 0; i < n; i++ {
  1640  					xi := x[i]
  1641  					if diag == blas.NonUnit {
  1642  						xi *= a[i*lda]
  1643  					}
  1644  					kk := min(k, n-i-1)
  1645  					for j, aij := range a[i*lda+1 : i*lda+kk+1] {
  1646  						xi += x[i+j+1] * aij
  1647  					}
  1648  					x[i] = xi
  1649  				}
  1650  			} else {
  1651  				ix := kx
  1652  				for i := 0; i < n; i++ {
  1653  					xi := x[ix]
  1654  					if diag == blas.NonUnit {
  1655  						xi *= a[i*lda]
  1656  					}
  1657  					kk := min(k, n-i-1)
  1658  					jx := ix + incX
  1659  					for _, aij := range a[i*lda+1 : i*lda+kk+1] {
  1660  						xi += x[jx] * aij
  1661  						jx += incX
  1662  					}
  1663  					x[ix] = xi
  1664  					ix += incX
  1665  				}
  1666  			}
  1667  		} else {
  1668  			if incX == 1 {
  1669  				for i := n - 1; i >= 0; i-- {
  1670  					xi := x[i]
  1671  					if diag == blas.NonUnit {
  1672  						xi *= a[i*lda+k]
  1673  					}
  1674  					kk := min(k, i)
  1675  					for j, aij := range a[i*lda+k-kk : i*lda+k] {
  1676  						xi += x[i-kk+j] * aij
  1677  					}
  1678  					x[i] = xi
  1679  				}
  1680  			} else {
  1681  				ix := kx + (n-1)*incX
  1682  				for i := n - 1; i >= 0; i-- {
  1683  					xi := x[ix]
  1684  					if diag == blas.NonUnit {
  1685  						xi *= a[i*lda+k]
  1686  					}
  1687  					kk := min(k, i)
  1688  					jx := ix - kk*incX
  1689  					for _, aij := range a[i*lda+k-kk : i*lda+k] {
  1690  						xi += x[jx] * aij
  1691  						jx += incX
  1692  					}
  1693  					x[ix] = xi
  1694  					ix -= incX
  1695  				}
  1696  			}
  1697  		}
  1698  	case blas.Trans:
  1699  		if uplo == blas.Upper {
  1700  			if incX == 1 {
  1701  				for i := n - 1; i >= 0; i-- {
  1702  					kk := min(k, n-i-1)
  1703  					xi := x[i]
  1704  					for j, aij := range a[i*lda+1 : i*lda+kk+1] {
  1705  						x[i+j+1] += xi * aij
  1706  					}
  1707  					if diag == blas.NonUnit {
  1708  						x[i] *= a[i*lda]
  1709  					}
  1710  				}
  1711  			} else {
  1712  				ix := kx + (n-1)*incX
  1713  				for i := n - 1; i >= 0; i-- {
  1714  					kk := min(k, n-i-1)
  1715  					jx := ix + incX
  1716  					xi := x[ix]
  1717  					for _, aij := range a[i*lda+1 : i*lda+kk+1] {
  1718  						x[jx] += xi * aij
  1719  						jx += incX
  1720  					}
  1721  					if diag == blas.NonUnit {
  1722  						x[ix] *= a[i*lda]
  1723  					}
  1724  					ix -= incX
  1725  				}
  1726  			}
  1727  		} else {
  1728  			if incX == 1 {
  1729  				for i := 0; i < n; i++ {
  1730  					kk := min(k, i)
  1731  					xi := x[i]
  1732  					for j, aij := range a[i*lda+k-kk : i*lda+k] {
  1733  						x[i-kk+j] += xi * aij
  1734  					}
  1735  					if diag == blas.NonUnit {
  1736  						x[i] *= a[i*lda+k]
  1737  					}
  1738  				}
  1739  			} else {
  1740  				ix := kx
  1741  				for i := 0; i < n; i++ {
  1742  					kk := min(k, i)
  1743  					jx := ix - kk*incX
  1744  					xi := x[ix]
  1745  					for _, aij := range a[i*lda+k-kk : i*lda+k] {
  1746  						x[jx] += xi * aij
  1747  						jx += incX
  1748  					}
  1749  					if diag == blas.NonUnit {
  1750  						x[ix] *= a[i*lda+k]
  1751  					}
  1752  					ix += incX
  1753  				}
  1754  			}
  1755  		}
  1756  	case blas.ConjTrans:
  1757  		if uplo == blas.Upper {
  1758  			if incX == 1 {
  1759  				for i := n - 1; i >= 0; i-- {
  1760  					kk := min(k, n-i-1)
  1761  					xi := x[i]
  1762  					for j, aij := range a[i*lda+1 : i*lda+kk+1] {
  1763  						x[i+j+1] += xi * cmplx.Conj(aij)
  1764  					}
  1765  					if diag == blas.NonUnit {
  1766  						x[i] *= cmplx.Conj(a[i*lda])
  1767  					}
  1768  				}
  1769  			} else {
  1770  				ix := kx + (n-1)*incX
  1771  				for i := n - 1; i >= 0; i-- {
  1772  					kk := min(k, n-i-1)
  1773  					jx := ix + incX
  1774  					xi := x[ix]
  1775  					for _, aij := range a[i*lda+1 : i*lda+kk+1] {
  1776  						x[jx] += xi * cmplx.Conj(aij)
  1777  						jx += incX
  1778  					}
  1779  					if diag == blas.NonUnit {
  1780  						x[ix] *= cmplx.Conj(a[i*lda])
  1781  					}
  1782  					ix -= incX
  1783  				}
  1784  			}
  1785  		} else {
  1786  			if incX == 1 {
  1787  				for i := 0; i < n; i++ {
  1788  					kk := min(k, i)
  1789  					xi := x[i]
  1790  					for j, aij := range a[i*lda+k-kk : i*lda+k] {
  1791  						x[i-kk+j] += xi * cmplx.Conj(aij)
  1792  					}
  1793  					if diag == blas.NonUnit {
  1794  						x[i] *= cmplx.Conj(a[i*lda+k])
  1795  					}
  1796  				}
  1797  			} else {
  1798  				ix := kx
  1799  				for i := 0; i < n; i++ {
  1800  					kk := min(k, i)
  1801  					jx := ix - kk*incX
  1802  					xi := x[ix]
  1803  					for _, aij := range a[i*lda+k-kk : i*lda+k] {
  1804  						x[jx] += xi * cmplx.Conj(aij)
  1805  						jx += incX
  1806  					}
  1807  					if diag == blas.NonUnit {
  1808  						x[ix] *= cmplx.Conj(a[i*lda+k])
  1809  					}
  1810  					ix += incX
  1811  				}
  1812  			}
  1813  		}
  1814  	}
  1815  }
  1816  
  1817  // Ctbsv solves one of the systems of equations
  1818  //
  1819  //	A * x = b   if trans == blas.NoTrans
  1820  //	Aᵀ * x = b  if trans == blas.Trans
  1821  //	Aᴴ * x = b  if trans == blas.ConjTrans
  1822  //
  1823  // where b and x are n element vectors and A is an n×n triangular band matrix
  1824  // with (k+1) diagonals.
  1825  //
  1826  // On entry, x contains the values of b, and the solution is
  1827  // stored in-place into x.
  1828  //
  1829  // No test for singularity or near-singularity is included in this
  1830  // routine. Such tests must be performed before calling this routine.
  1831  //
  1832  // Complex64 implementations are autogenerated and not directly tested.
  1833  func (Implementation) Ctbsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, k int, a []complex64, lda int, x []complex64, incX int) {
  1834  	switch trans {
  1835  	default:
  1836  		panic(badTranspose)
  1837  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
  1838  	}
  1839  	switch uplo {
  1840  	default:
  1841  		panic(badUplo)
  1842  	case blas.Upper, blas.Lower:
  1843  	}
  1844  	switch diag {
  1845  	default:
  1846  		panic(badDiag)
  1847  	case blas.NonUnit, blas.Unit:
  1848  	}
  1849  	if n < 0 {
  1850  		panic(nLT0)
  1851  	}
  1852  	if k < 0 {
  1853  		panic(kLT0)
  1854  	}
  1855  	if lda < k+1 {
  1856  		panic(badLdA)
  1857  	}
  1858  	if incX == 0 {
  1859  		panic(zeroIncX)
  1860  	}
  1861  
  1862  	// Quick return if possible.
  1863  	if n == 0 {
  1864  		return
  1865  	}
  1866  
  1867  	// For zero matrix size the following slice length checks are trivially satisfied.
  1868  	if len(a) < lda*(n-1)+k+1 {
  1869  		panic(shortA)
  1870  	}
  1871  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  1872  		panic(shortX)
  1873  	}
  1874  
  1875  	// Set up start index in X.
  1876  	var kx int
  1877  	if incX < 0 {
  1878  		kx = (1 - n) * incX
  1879  	}
  1880  
  1881  	switch trans {
  1882  	case blas.NoTrans:
  1883  		if uplo == blas.Upper {
  1884  			if incX == 1 {
  1885  				for i := n - 1; i >= 0; i-- {
  1886  					kk := min(k, n-i-1)
  1887  					var sum complex64
  1888  					for j, aij := range a[i*lda+1 : i*lda+kk+1] {
  1889  						sum += x[i+1+j] * aij
  1890  					}
  1891  					x[i] -= sum
  1892  					if diag == blas.NonUnit {
  1893  						x[i] /= a[i*lda]
  1894  					}
  1895  				}
  1896  			} else {
  1897  				ix := kx + (n-1)*incX
  1898  				for i := n - 1; i >= 0; i-- {
  1899  					kk := min(k, n-i-1)
  1900  					var sum complex64
  1901  					jx := ix + incX
  1902  					for _, aij := range a[i*lda+1 : i*lda+kk+1] {
  1903  						sum += x[jx] * aij
  1904  						jx += incX
  1905  					}
  1906  					x[ix] -= sum
  1907  					if diag == blas.NonUnit {
  1908  						x[ix] /= a[i*lda]
  1909  					}
  1910  					ix -= incX
  1911  				}
  1912  			}
  1913  		} else {
  1914  			if incX == 1 {
  1915  				for i := 0; i < n; i++ {
  1916  					kk := min(k, i)
  1917  					var sum complex64
  1918  					for j, aij := range a[i*lda+k-kk : i*lda+k] {
  1919  						sum += x[i-kk+j] * aij
  1920  					}
  1921  					x[i] -= sum
  1922  					if diag == blas.NonUnit {
  1923  						x[i] /= a[i*lda+k]
  1924  					}
  1925  				}
  1926  			} else {
  1927  				ix := kx
  1928  				for i := 0; i < n; i++ {
  1929  					kk := min(k, i)
  1930  					var sum complex64
  1931  					jx := ix - kk*incX
  1932  					for _, aij := range a[i*lda+k-kk : i*lda+k] {
  1933  						sum += x[jx] * aij
  1934  						jx += incX
  1935  					}
  1936  					x[ix] -= sum
  1937  					if diag == blas.NonUnit {
  1938  						x[ix] /= a[i*lda+k]
  1939  					}
  1940  					ix += incX
  1941  				}
  1942  			}
  1943  		}
  1944  	case blas.Trans:
  1945  		if uplo == blas.Upper {
  1946  			if incX == 1 {
  1947  				for i := 0; i < n; i++ {
  1948  					if diag == blas.NonUnit {
  1949  						x[i] /= a[i*lda]
  1950  					}
  1951  					kk := min(k, n-i-1)
  1952  					xi := x[i]
  1953  					for j, aij := range a[i*lda+1 : i*lda+kk+1] {
  1954  						x[i+1+j] -= xi * aij
  1955  					}
  1956  				}
  1957  			} else {
  1958  				ix := kx
  1959  				for i := 0; i < n; i++ {
  1960  					if diag == blas.NonUnit {
  1961  						x[ix] /= a[i*lda]
  1962  					}
  1963  					kk := min(k, n-i-1)
  1964  					xi := x[ix]
  1965  					jx := ix + incX
  1966  					for _, aij := range a[i*lda+1 : i*lda+kk+1] {
  1967  						x[jx] -= xi * aij
  1968  						jx += incX
  1969  					}
  1970  					ix += incX
  1971  				}
  1972  			}
  1973  		} else {
  1974  			if incX == 1 {
  1975  				for i := n - 1; i >= 0; i-- {
  1976  					if diag == blas.NonUnit {
  1977  						x[i] /= a[i*lda+k]
  1978  					}
  1979  					kk := min(k, i)
  1980  					xi := x[i]
  1981  					for j, aij := range a[i*lda+k-kk : i*lda+k] {
  1982  						x[i-kk+j] -= xi * aij
  1983  					}
  1984  				}
  1985  			} else {
  1986  				ix := kx + (n-1)*incX
  1987  				for i := n - 1; i >= 0; i-- {
  1988  					if diag == blas.NonUnit {
  1989  						x[ix] /= a[i*lda+k]
  1990  					}
  1991  					kk := min(k, i)
  1992  					xi := x[ix]
  1993  					jx := ix - kk*incX
  1994  					for _, aij := range a[i*lda+k-kk : i*lda+k] {
  1995  						x[jx] -= xi * aij
  1996  						jx += incX
  1997  					}
  1998  					ix -= incX
  1999  				}
  2000  			}
  2001  		}
  2002  	case blas.ConjTrans:
  2003  		if uplo == blas.Upper {
  2004  			if incX == 1 {
  2005  				for i := 0; i < n; i++ {
  2006  					if diag == blas.NonUnit {
  2007  						x[i] /= cmplx.Conj(a[i*lda])
  2008  					}
  2009  					kk := min(k, n-i-1)
  2010  					xi := x[i]
  2011  					for j, aij := range a[i*lda+1 : i*lda+kk+1] {
  2012  						x[i+1+j] -= xi * cmplx.Conj(aij)
  2013  					}
  2014  				}
  2015  			} else {
  2016  				ix := kx
  2017  				for i := 0; i < n; i++ {
  2018  					if diag == blas.NonUnit {
  2019  						x[ix] /= cmplx.Conj(a[i*lda])
  2020  					}
  2021  					kk := min(k, n-i-1)
  2022  					xi := x[ix]
  2023  					jx := ix + incX
  2024  					for _, aij := range a[i*lda+1 : i*lda+kk+1] {
  2025  						x[jx] -= xi * cmplx.Conj(aij)
  2026  						jx += incX
  2027  					}
  2028  					ix += incX
  2029  				}
  2030  			}
  2031  		} else {
  2032  			if incX == 1 {
  2033  				for i := n - 1; i >= 0; i-- {
  2034  					if diag == blas.NonUnit {
  2035  						x[i] /= cmplx.Conj(a[i*lda+k])
  2036  					}
  2037  					kk := min(k, i)
  2038  					xi := x[i]
  2039  					for j, aij := range a[i*lda+k-kk : i*lda+k] {
  2040  						x[i-kk+j] -= xi * cmplx.Conj(aij)
  2041  					}
  2042  				}
  2043  			} else {
  2044  				ix := kx + (n-1)*incX
  2045  				for i := n - 1; i >= 0; i-- {
  2046  					if diag == blas.NonUnit {
  2047  						x[ix] /= cmplx.Conj(a[i*lda+k])
  2048  					}
  2049  					kk := min(k, i)
  2050  					xi := x[ix]
  2051  					jx := ix - kk*incX
  2052  					for _, aij := range a[i*lda+k-kk : i*lda+k] {
  2053  						x[jx] -= xi * cmplx.Conj(aij)
  2054  						jx += incX
  2055  					}
  2056  					ix -= incX
  2057  				}
  2058  			}
  2059  		}
  2060  	}
  2061  }
  2062  
  2063  // Ctpmv performs one of the matrix-vector operations
  2064  //
  2065  //	x = A * x   if trans = blas.NoTrans
  2066  //	x = Aᵀ * x  if trans = blas.Trans
  2067  //	x = Aᴴ * x  if trans = blas.ConjTrans
  2068  //
  2069  // where x is an n element vector and A is an n×n triangular matrix, supplied in
  2070  // packed form.
  2071  //
  2072  // Complex64 implementations are autogenerated and not directly tested.
  2073  func (Implementation) Ctpmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex64, x []complex64, incX int) {
  2074  	switch uplo {
  2075  	default:
  2076  		panic(badUplo)
  2077  	case blas.Upper, blas.Lower:
  2078  	}
  2079  	switch trans {
  2080  	default:
  2081  		panic(badTranspose)
  2082  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
  2083  	}
  2084  	switch diag {
  2085  	default:
  2086  		panic(badDiag)
  2087  	case blas.NonUnit, blas.Unit:
  2088  	}
  2089  	if n < 0 {
  2090  		panic(nLT0)
  2091  	}
  2092  	if incX == 0 {
  2093  		panic(zeroIncX)
  2094  	}
  2095  
  2096  	// Quick return if possible.
  2097  	if n == 0 {
  2098  		return
  2099  	}
  2100  
  2101  	// For zero matrix size the following slice length checks are trivially satisfied.
  2102  	if len(ap) < n*(n+1)/2 {
  2103  		panic(shortAP)
  2104  	}
  2105  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  2106  		panic(shortX)
  2107  	}
  2108  
  2109  	// Set up start index in X.
  2110  	var kx int
  2111  	if incX < 0 {
  2112  		kx = (1 - n) * incX
  2113  	}
  2114  
  2115  	// The elements of A are accessed sequentially with one pass through A.
  2116  
  2117  	if trans == blas.NoTrans {
  2118  		// Form x = A*x.
  2119  		if uplo == blas.Upper {
  2120  			// kk points to the current diagonal element in ap.
  2121  			kk := 0
  2122  			if incX == 1 {
  2123  				x = x[:n]
  2124  				for i := range x {
  2125  					if diag == blas.NonUnit {
  2126  						x[i] *= ap[kk]
  2127  					}
  2128  					if n-i-1 > 0 {
  2129  						x[i] += c64.DotuUnitary(ap[kk+1:kk+n-i], x[i+1:])
  2130  					}
  2131  					kk += n - i
  2132  				}
  2133  			} else {
  2134  				ix := kx
  2135  				for i := 0; i < n; i++ {
  2136  					if diag == blas.NonUnit {
  2137  						x[ix] *= ap[kk]
  2138  					}
  2139  					if n-i-1 > 0 {
  2140  						x[ix] += c64.DotuInc(ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
  2141  					}
  2142  					ix += incX
  2143  					kk += n - i
  2144  				}
  2145  			}
  2146  		} else {
  2147  			// kk points to the beginning of current row in ap.
  2148  			kk := n*(n+1)/2 - n
  2149  			if incX == 1 {
  2150  				for i := n - 1; i >= 0; i-- {
  2151  					if diag == blas.NonUnit {
  2152  						x[i] *= ap[kk+i]
  2153  					}
  2154  					if i > 0 {
  2155  						x[i] += c64.DotuUnitary(ap[kk:kk+i], x[:i])
  2156  					}
  2157  					kk -= i
  2158  				}
  2159  			} else {
  2160  				ix := kx + (n-1)*incX
  2161  				for i := n - 1; i >= 0; i-- {
  2162  					if diag == blas.NonUnit {
  2163  						x[ix] *= ap[kk+i]
  2164  					}
  2165  					if i > 0 {
  2166  						x[ix] += c64.DotuInc(ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
  2167  					}
  2168  					ix -= incX
  2169  					kk -= i
  2170  				}
  2171  			}
  2172  		}
  2173  		return
  2174  	}
  2175  
  2176  	if trans == blas.Trans {
  2177  		// Form x = Aᵀ*x.
  2178  		if uplo == blas.Upper {
  2179  			// kk points to the current diagonal element in ap.
  2180  			kk := n*(n+1)/2 - 1
  2181  			if incX == 1 {
  2182  				for i := n - 1; i >= 0; i-- {
  2183  					xi := x[i]
  2184  					if diag == blas.NonUnit {
  2185  						x[i] *= ap[kk]
  2186  					}
  2187  					if n-i-1 > 0 {
  2188  						c64.AxpyUnitary(xi, ap[kk+1:kk+n-i], x[i+1:n])
  2189  					}
  2190  					kk -= n - i + 1
  2191  				}
  2192  			} else {
  2193  				ix := kx + (n-1)*incX
  2194  				for i := n - 1; i >= 0; i-- {
  2195  					xi := x[ix]
  2196  					if diag == blas.NonUnit {
  2197  						x[ix] *= ap[kk]
  2198  					}
  2199  					if n-i-1 > 0 {
  2200  						c64.AxpyInc(xi, ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
  2201  					}
  2202  					ix -= incX
  2203  					kk -= n - i + 1
  2204  				}
  2205  			}
  2206  		} else {
  2207  			// kk points to the beginning of current row in ap.
  2208  			kk := 0
  2209  			if incX == 1 {
  2210  				x = x[:n]
  2211  				for i := range x {
  2212  					if i > 0 {
  2213  						c64.AxpyUnitary(x[i], ap[kk:kk+i], x[:i])
  2214  					}
  2215  					if diag == blas.NonUnit {
  2216  						x[i] *= ap[kk+i]
  2217  					}
  2218  					kk += i + 1
  2219  				}
  2220  			} else {
  2221  				ix := kx
  2222  				for i := 0; i < n; i++ {
  2223  					if i > 0 {
  2224  						c64.AxpyInc(x[ix], ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
  2225  					}
  2226  					if diag == blas.NonUnit {
  2227  						x[ix] *= ap[kk+i]
  2228  					}
  2229  					ix += incX
  2230  					kk += i + 1
  2231  				}
  2232  			}
  2233  		}
  2234  		return
  2235  	}
  2236  
  2237  	// Form x = Aᴴ*x.
  2238  	if uplo == blas.Upper {
  2239  		// kk points to the current diagonal element in ap.
  2240  		kk := n*(n+1)/2 - 1
  2241  		if incX == 1 {
  2242  			for i := n - 1; i >= 0; i-- {
  2243  				xi := x[i]
  2244  				if diag == blas.NonUnit {
  2245  					x[i] *= cmplx.Conj(ap[kk])
  2246  				}
  2247  				k := kk + 1
  2248  				for j := i + 1; j < n; j++ {
  2249  					x[j] += xi * cmplx.Conj(ap[k])
  2250  					k++
  2251  				}
  2252  				kk -= n - i + 1
  2253  			}
  2254  		} else {
  2255  			ix := kx + (n-1)*incX
  2256  			for i := n - 1; i >= 0; i-- {
  2257  				xi := x[ix]
  2258  				if diag == blas.NonUnit {
  2259  					x[ix] *= cmplx.Conj(ap[kk])
  2260  				}
  2261  				jx := ix + incX
  2262  				k := kk + 1
  2263  				for j := i + 1; j < n; j++ {
  2264  					x[jx] += xi * cmplx.Conj(ap[k])
  2265  					jx += incX
  2266  					k++
  2267  				}
  2268  				ix -= incX
  2269  				kk -= n - i + 1
  2270  			}
  2271  		}
  2272  	} else {
  2273  		// kk points to the beginning of current row in ap.
  2274  		kk := 0
  2275  		if incX == 1 {
  2276  			x = x[:n]
  2277  			for i, xi := range x {
  2278  				for j := 0; j < i; j++ {
  2279  					x[j] += xi * cmplx.Conj(ap[kk+j])
  2280  				}
  2281  				if diag == blas.NonUnit {
  2282  					x[i] *= cmplx.Conj(ap[kk+i])
  2283  				}
  2284  				kk += i + 1
  2285  			}
  2286  		} else {
  2287  			ix := kx
  2288  			for i := 0; i < n; i++ {
  2289  				xi := x[ix]
  2290  				jx := kx
  2291  				for j := 0; j < i; j++ {
  2292  					x[jx] += xi * cmplx.Conj(ap[kk+j])
  2293  					jx += incX
  2294  				}
  2295  				if diag == blas.NonUnit {
  2296  					x[ix] *= cmplx.Conj(ap[kk+i])
  2297  				}
  2298  				ix += incX
  2299  				kk += i + 1
  2300  			}
  2301  		}
  2302  	}
  2303  }
  2304  
  2305  // Ctpsv solves one of the systems of equations
  2306  //
  2307  //	A * x = b   if trans == blas.NoTrans
  2308  //	Aᵀ * x = b  if trans == blas.Trans
  2309  //	Aᴴ * x = b  if trans == blas.ConjTrans
  2310  //
  2311  // where b and x are n element vectors and A is an n×n triangular matrix in
  2312  // packed form.
  2313  //
  2314  // On entry, x contains the values of b, and the solution is
  2315  // stored in-place into x.
  2316  //
  2317  // No test for singularity or near-singularity is included in this
  2318  // routine. Such tests must be performed before calling this routine.
  2319  //
  2320  // Complex64 implementations are autogenerated and not directly tested.
  2321  func (Implementation) Ctpsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex64, x []complex64, incX int) {
  2322  	switch uplo {
  2323  	default:
  2324  		panic(badUplo)
  2325  	case blas.Upper, blas.Lower:
  2326  	}
  2327  	switch trans {
  2328  	default:
  2329  		panic(badTranspose)
  2330  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
  2331  	}
  2332  	switch diag {
  2333  	default:
  2334  		panic(badDiag)
  2335  	case blas.NonUnit, blas.Unit:
  2336  	}
  2337  	if n < 0 {
  2338  		panic(nLT0)
  2339  	}
  2340  	if incX == 0 {
  2341  		panic(zeroIncX)
  2342  	}
  2343  
  2344  	// Quick return if possible.
  2345  	if n == 0 {
  2346  		return
  2347  	}
  2348  
  2349  	// For zero matrix size the following slice length checks are trivially satisfied.
  2350  	if len(ap) < n*(n+1)/2 {
  2351  		panic(shortAP)
  2352  	}
  2353  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  2354  		panic(shortX)
  2355  	}
  2356  
  2357  	// Set up start index in X.
  2358  	var kx int
  2359  	if incX < 0 {
  2360  		kx = (1 - n) * incX
  2361  	}
  2362  
  2363  	// The elements of A are accessed sequentially with one pass through ap.
  2364  
  2365  	if trans == blas.NoTrans {
  2366  		// Form x = inv(A)*x.
  2367  		if uplo == blas.Upper {
  2368  			kk := n*(n+1)/2 - 1
  2369  			if incX == 1 {
  2370  				for i := n - 1; i >= 0; i-- {
  2371  					aii := ap[kk]
  2372  					if n-i-1 > 0 {
  2373  						x[i] -= c64.DotuUnitary(x[i+1:n], ap[kk+1:kk+n-i])
  2374  					}
  2375  					if diag == blas.NonUnit {
  2376  						x[i] /= aii
  2377  					}
  2378  					kk -= n - i + 1
  2379  				}
  2380  			} else {
  2381  				ix := kx + (n-1)*incX
  2382  				for i := n - 1; i >= 0; i-- {
  2383  					aii := ap[kk]
  2384  					if n-i-1 > 0 {
  2385  						x[ix] -= c64.DotuInc(x, ap[kk+1:kk+n-i], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0)
  2386  					}
  2387  					if diag == blas.NonUnit {
  2388  						x[ix] /= aii
  2389  					}
  2390  					ix -= incX
  2391  					kk -= n - i + 1
  2392  				}
  2393  			}
  2394  		} else {
  2395  			kk := 0
  2396  			if incX == 1 {
  2397  				for i := 0; i < n; i++ {
  2398  					if i > 0 {
  2399  						x[i] -= c64.DotuUnitary(x[:i], ap[kk:kk+i])
  2400  					}
  2401  					if diag == blas.NonUnit {
  2402  						x[i] /= ap[kk+i]
  2403  					}
  2404  					kk += i + 1
  2405  				}
  2406  			} else {
  2407  				ix := kx
  2408  				for i := 0; i < n; i++ {
  2409  					if i > 0 {
  2410  						x[ix] -= c64.DotuInc(x, ap[kk:kk+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0)
  2411  					}
  2412  					if diag == blas.NonUnit {
  2413  						x[ix] /= ap[kk+i]
  2414  					}
  2415  					ix += incX
  2416  					kk += i + 1
  2417  				}
  2418  			}
  2419  		}
  2420  		return
  2421  	}
  2422  
  2423  	if trans == blas.Trans {
  2424  		// Form x = inv(Aᵀ)*x.
  2425  		if uplo == blas.Upper {
  2426  			kk := 0
  2427  			if incX == 1 {
  2428  				for j := 0; j < n; j++ {
  2429  					if diag == blas.NonUnit {
  2430  						x[j] /= ap[kk]
  2431  					}
  2432  					if n-j-1 > 0 {
  2433  						c64.AxpyUnitary(-x[j], ap[kk+1:kk+n-j], x[j+1:n])
  2434  					}
  2435  					kk += n - j
  2436  				}
  2437  			} else {
  2438  				jx := kx
  2439  				for j := 0; j < n; j++ {
  2440  					if diag == blas.NonUnit {
  2441  						x[jx] /= ap[kk]
  2442  					}
  2443  					if n-j-1 > 0 {
  2444  						c64.AxpyInc(-x[jx], ap[kk+1:kk+n-j], x, uintptr(n-j-1), 1, uintptr(incX), 0, uintptr(jx+incX))
  2445  					}
  2446  					jx += incX
  2447  					kk += n - j
  2448  				}
  2449  			}
  2450  		} else {
  2451  			kk := n*(n+1)/2 - n
  2452  			if incX == 1 {
  2453  				for j := n - 1; j >= 0; j-- {
  2454  					if diag == blas.NonUnit {
  2455  						x[j] /= ap[kk+j]
  2456  					}
  2457  					if j > 0 {
  2458  						c64.AxpyUnitary(-x[j], ap[kk:kk+j], x[:j])
  2459  					}
  2460  					kk -= j
  2461  				}
  2462  			} else {
  2463  				jx := kx + (n-1)*incX
  2464  				for j := n - 1; j >= 0; j-- {
  2465  					if diag == blas.NonUnit {
  2466  						x[jx] /= ap[kk+j]
  2467  					}
  2468  					if j > 0 {
  2469  						c64.AxpyInc(-x[jx], ap[kk:kk+j], x, uintptr(j), 1, uintptr(incX), 0, uintptr(kx))
  2470  					}
  2471  					jx -= incX
  2472  					kk -= j
  2473  				}
  2474  			}
  2475  		}
  2476  		return
  2477  	}
  2478  
  2479  	// Form x = inv(Aᴴ)*x.
  2480  	if uplo == blas.Upper {
  2481  		kk := 0
  2482  		if incX == 1 {
  2483  			for j := 0; j < n; j++ {
  2484  				if diag == blas.NonUnit {
  2485  					x[j] /= cmplx.Conj(ap[kk])
  2486  				}
  2487  				xj := x[j]
  2488  				k := kk + 1
  2489  				for i := j + 1; i < n; i++ {
  2490  					x[i] -= xj * cmplx.Conj(ap[k])
  2491  					k++
  2492  				}
  2493  				kk += n - j
  2494  			}
  2495  		} else {
  2496  			jx := kx
  2497  			for j := 0; j < n; j++ {
  2498  				if diag == blas.NonUnit {
  2499  					x[jx] /= cmplx.Conj(ap[kk])
  2500  				}
  2501  				xj := x[jx]
  2502  				ix := jx + incX
  2503  				k := kk + 1
  2504  				for i := j + 1; i < n; i++ {
  2505  					x[ix] -= xj * cmplx.Conj(ap[k])
  2506  					ix += incX
  2507  					k++
  2508  				}
  2509  				jx += incX
  2510  				kk += n - j
  2511  			}
  2512  		}
  2513  	} else {
  2514  		kk := n*(n+1)/2 - n
  2515  		if incX == 1 {
  2516  			for j := n - 1; j >= 0; j-- {
  2517  				if diag == blas.NonUnit {
  2518  					x[j] /= cmplx.Conj(ap[kk+j])
  2519  				}
  2520  				xj := x[j]
  2521  				for i := 0; i < j; i++ {
  2522  					x[i] -= xj * cmplx.Conj(ap[kk+i])
  2523  				}
  2524  				kk -= j
  2525  			}
  2526  		} else {
  2527  			jx := kx + (n-1)*incX
  2528  			for j := n - 1; j >= 0; j-- {
  2529  				if diag == blas.NonUnit {
  2530  					x[jx] /= cmplx.Conj(ap[kk+j])
  2531  				}
  2532  				xj := x[jx]
  2533  				ix := kx
  2534  				for i := 0; i < j; i++ {
  2535  					x[ix] -= xj * cmplx.Conj(ap[kk+i])
  2536  					ix += incX
  2537  				}
  2538  				jx -= incX
  2539  				kk -= j
  2540  			}
  2541  		}
  2542  	}
  2543  }
  2544  
  2545  // Ctrmv performs one of the matrix-vector operations
  2546  //
  2547  //	x = A * x   if trans = blas.NoTrans
  2548  //	x = Aᵀ * x  if trans = blas.Trans
  2549  //	x = Aᴴ * x  if trans = blas.ConjTrans
  2550  //
  2551  // where x is a vector, and A is an n×n triangular matrix.
  2552  //
  2553  // Complex64 implementations are autogenerated and not directly tested.
  2554  func (Implementation) Ctrmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex64, lda int, x []complex64, incX int) {
  2555  	switch trans {
  2556  	default:
  2557  		panic(badTranspose)
  2558  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
  2559  	}
  2560  	switch uplo {
  2561  	default:
  2562  		panic(badUplo)
  2563  	case blas.Upper, blas.Lower:
  2564  	}
  2565  	switch diag {
  2566  	default:
  2567  		panic(badDiag)
  2568  	case blas.NonUnit, blas.Unit:
  2569  	}
  2570  	if n < 0 {
  2571  		panic(nLT0)
  2572  	}
  2573  	if lda < max(1, n) {
  2574  		panic(badLdA)
  2575  	}
  2576  	if incX == 0 {
  2577  		panic(zeroIncX)
  2578  	}
  2579  
  2580  	// Quick return if possible.
  2581  	if n == 0 {
  2582  		return
  2583  	}
  2584  
  2585  	// For zero matrix size the following slice length checks are trivially satisfied.
  2586  	if len(a) < lda*(n-1)+n {
  2587  		panic(shortA)
  2588  	}
  2589  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  2590  		panic(shortX)
  2591  	}
  2592  
  2593  	// Set up start index in X.
  2594  	var kx int
  2595  	if incX < 0 {
  2596  		kx = (1 - n) * incX
  2597  	}
  2598  
  2599  	// The elements of A are accessed sequentially with one pass through A.
  2600  
  2601  	if trans == blas.NoTrans {
  2602  		// Form x = A*x.
  2603  		if uplo == blas.Upper {
  2604  			if incX == 1 {
  2605  				for i := 0; i < n; i++ {
  2606  					if diag == blas.NonUnit {
  2607  						x[i] *= a[i*lda+i]
  2608  					}
  2609  					if n-i-1 > 0 {
  2610  						x[i] += c64.DotuUnitary(a[i*lda+i+1:i*lda+n], x[i+1:n])
  2611  					}
  2612  				}
  2613  			} else {
  2614  				ix := kx
  2615  				for i := 0; i < n; i++ {
  2616  					if diag == blas.NonUnit {
  2617  						x[ix] *= a[i*lda+i]
  2618  					}
  2619  					if n-i-1 > 0 {
  2620  						x[ix] += c64.DotuInc(a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
  2621  					}
  2622  					ix += incX
  2623  				}
  2624  			}
  2625  		} else {
  2626  			if incX == 1 {
  2627  				for i := n - 1; i >= 0; i-- {
  2628  					if diag == blas.NonUnit {
  2629  						x[i] *= a[i*lda+i]
  2630  					}
  2631  					if i > 0 {
  2632  						x[i] += c64.DotuUnitary(a[i*lda:i*lda+i], x[:i])
  2633  					}
  2634  				}
  2635  			} else {
  2636  				ix := kx + (n-1)*incX
  2637  				for i := n - 1; i >= 0; i-- {
  2638  					if diag == blas.NonUnit {
  2639  						x[ix] *= a[i*lda+i]
  2640  					}
  2641  					if i > 0 {
  2642  						x[ix] += c64.DotuInc(a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
  2643  					}
  2644  					ix -= incX
  2645  				}
  2646  			}
  2647  		}
  2648  		return
  2649  	}
  2650  
  2651  	if trans == blas.Trans {
  2652  		// Form x = Aᵀ*x.
  2653  		if uplo == blas.Upper {
  2654  			if incX == 1 {
  2655  				for i := n - 1; i >= 0; i-- {
  2656  					xi := x[i]
  2657  					if diag == blas.NonUnit {
  2658  						x[i] *= a[i*lda+i]
  2659  					}
  2660  					if n-i-1 > 0 {
  2661  						c64.AxpyUnitary(xi, a[i*lda+i+1:i*lda+n], x[i+1:n])
  2662  					}
  2663  				}
  2664  			} else {
  2665  				ix := kx + (n-1)*incX
  2666  				for i := n - 1; i >= 0; i-- {
  2667  					xi := x[ix]
  2668  					if diag == blas.NonUnit {
  2669  						x[ix] *= a[i*lda+i]
  2670  					}
  2671  					if n-i-1 > 0 {
  2672  						c64.AxpyInc(xi, a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
  2673  					}
  2674  					ix -= incX
  2675  				}
  2676  			}
  2677  		} else {
  2678  			if incX == 1 {
  2679  				for i := 0; i < n; i++ {
  2680  					if i > 0 {
  2681  						c64.AxpyUnitary(x[i], a[i*lda:i*lda+i], x[:i])
  2682  					}
  2683  					if diag == blas.NonUnit {
  2684  						x[i] *= a[i*lda+i]
  2685  					}
  2686  				}
  2687  			} else {
  2688  				ix := kx
  2689  				for i := 0; i < n; i++ {
  2690  					if i > 0 {
  2691  						c64.AxpyInc(x[ix], a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
  2692  					}
  2693  					if diag == blas.NonUnit {
  2694  						x[ix] *= a[i*lda+i]
  2695  					}
  2696  					ix += incX
  2697  				}
  2698  			}
  2699  		}
  2700  		return
  2701  	}
  2702  
  2703  	// Form x = Aᴴ*x.
  2704  	if uplo == blas.Upper {
  2705  		if incX == 1 {
  2706  			for i := n - 1; i >= 0; i-- {
  2707  				xi := x[i]
  2708  				if diag == blas.NonUnit {
  2709  					x[i] *= cmplx.Conj(a[i*lda+i])
  2710  				}
  2711  				for j := i + 1; j < n; j++ {
  2712  					x[j] += xi * cmplx.Conj(a[i*lda+j])
  2713  				}
  2714  			}
  2715  		} else {
  2716  			ix := kx + (n-1)*incX
  2717  			for i := n - 1; i >= 0; i-- {
  2718  				xi := x[ix]
  2719  				if diag == blas.NonUnit {
  2720  					x[ix] *= cmplx.Conj(a[i*lda+i])
  2721  				}
  2722  				jx := ix + incX
  2723  				for j := i + 1; j < n; j++ {
  2724  					x[jx] += xi * cmplx.Conj(a[i*lda+j])
  2725  					jx += incX
  2726  				}
  2727  				ix -= incX
  2728  			}
  2729  		}
  2730  	} else {
  2731  		if incX == 1 {
  2732  			for i := 0; i < n; i++ {
  2733  				for j := 0; j < i; j++ {
  2734  					x[j] += x[i] * cmplx.Conj(a[i*lda+j])
  2735  				}
  2736  				if diag == blas.NonUnit {
  2737  					x[i] *= cmplx.Conj(a[i*lda+i])
  2738  				}
  2739  			}
  2740  		} else {
  2741  			ix := kx
  2742  			for i := 0; i < n; i++ {
  2743  				jx := kx
  2744  				for j := 0; j < i; j++ {
  2745  					x[jx] += x[ix] * cmplx.Conj(a[i*lda+j])
  2746  					jx += incX
  2747  				}
  2748  				if diag == blas.NonUnit {
  2749  					x[ix] *= cmplx.Conj(a[i*lda+i])
  2750  				}
  2751  				ix += incX
  2752  			}
  2753  		}
  2754  	}
  2755  }
  2756  
  2757  // Ctrsv solves one of the systems of equations
  2758  //
  2759  //	A * x = b   if trans == blas.NoTrans
  2760  //	Aᵀ * x = b  if trans == blas.Trans
  2761  //	Aᴴ * x = b  if trans == blas.ConjTrans
  2762  //
  2763  // where b and x are n element vectors and A is an n×n triangular matrix.
  2764  //
  2765  // On entry, x contains the values of b, and the solution is
  2766  // stored in-place into x.
  2767  //
  2768  // No test for singularity or near-singularity is included in this
  2769  // routine. Such tests must be performed before calling this routine.
  2770  //
  2771  // Complex64 implementations are autogenerated and not directly tested.
  2772  func (Implementation) Ctrsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex64, lda int, x []complex64, incX int) {
  2773  	switch trans {
  2774  	default:
  2775  		panic(badTranspose)
  2776  	case blas.NoTrans, blas.Trans, blas.ConjTrans:
  2777  	}
  2778  	switch uplo {
  2779  	default:
  2780  		panic(badUplo)
  2781  	case blas.Upper, blas.Lower:
  2782  	}
  2783  	switch diag {
  2784  	default:
  2785  		panic(badDiag)
  2786  	case blas.NonUnit, blas.Unit:
  2787  	}
  2788  	if n < 0 {
  2789  		panic(nLT0)
  2790  	}
  2791  	if lda < max(1, n) {
  2792  		panic(badLdA)
  2793  	}
  2794  	if incX == 0 {
  2795  		panic(zeroIncX)
  2796  	}
  2797  
  2798  	// Quick return if possible.
  2799  	if n == 0 {
  2800  		return
  2801  	}
  2802  
  2803  	// For zero matrix size the following slice length checks are trivially satisfied.
  2804  	if len(a) < lda*(n-1)+n {
  2805  		panic(shortA)
  2806  	}
  2807  	if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) {
  2808  		panic(shortX)
  2809  	}
  2810  
  2811  	// Set up start index in X.
  2812  	var kx int
  2813  	if incX < 0 {
  2814  		kx = (1 - n) * incX
  2815  	}
  2816  
  2817  	// The elements of A are accessed sequentially with one pass through A.
  2818  
  2819  	if trans == blas.NoTrans {
  2820  		// Form x = inv(A)*x.
  2821  		if uplo == blas.Upper {
  2822  			if incX == 1 {
  2823  				for i := n - 1; i >= 0; i-- {
  2824  					aii := a[i*lda+i]
  2825  					if n-i-1 > 0 {
  2826  						x[i] -= c64.DotuUnitary(x[i+1:n], a[i*lda+i+1:i*lda+n])
  2827  					}
  2828  					if diag == blas.NonUnit {
  2829  						x[i] /= aii
  2830  					}
  2831  				}
  2832  			} else {
  2833  				ix := kx + (n-1)*incX
  2834  				for i := n - 1; i >= 0; i-- {
  2835  					aii := a[i*lda+i]
  2836  					if n-i-1 > 0 {
  2837  						x[ix] -= c64.DotuInc(x, a[i*lda+i+1:i*lda+n], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0)
  2838  					}
  2839  					if diag == blas.NonUnit {
  2840  						x[ix] /= aii
  2841  					}
  2842  					ix -= incX
  2843  				}
  2844  			}
  2845  		} else {
  2846  			if incX == 1 {
  2847  				for i := 0; i < n; i++ {
  2848  					if i > 0 {
  2849  						x[i] -= c64.DotuUnitary(x[:i], a[i*lda:i*lda+i])
  2850  					}
  2851  					if diag == blas.NonUnit {
  2852  						x[i] /= a[i*lda+i]
  2853  					}
  2854  				}
  2855  			} else {
  2856  				ix := kx
  2857  				for i := 0; i < n; i++ {
  2858  					if i > 0 {
  2859  						x[ix] -= c64.DotuInc(x, a[i*lda:i*lda+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0)
  2860  					}
  2861  					if diag == blas.NonUnit {
  2862  						x[ix] /= a[i*lda+i]
  2863  					}
  2864  					ix += incX
  2865  				}
  2866  			}
  2867  		}
  2868  		return
  2869  	}
  2870  
  2871  	if trans == blas.Trans {
  2872  		// Form x = inv(Aᵀ)*x.
  2873  		if uplo == blas.Upper {
  2874  			if incX == 1 {
  2875  				for j := 0; j < n; j++ {
  2876  					if diag == blas.NonUnit {
  2877  						x[j] /= a[j*lda+j]
  2878  					}
  2879  					if n-j-1 > 0 {
  2880  						c64.AxpyUnitary(-x[j], a[j*lda+j+1:j*lda+n], x[j+1:n])
  2881  					}
  2882  				}
  2883  			} else {
  2884  				jx := kx
  2885  				for j := 0; j < n; j++ {
  2886  					if diag == blas.NonUnit {
  2887  						x[jx] /= a[j*lda+j]
  2888  					}
  2889  					if n-j-1 > 0 {
  2890  						c64.AxpyInc(-x[jx], a[j*lda+j+1:j*lda+n], x, uintptr(n-j-1), 1, uintptr(incX), 0, uintptr(jx+incX))
  2891  					}
  2892  					jx += incX
  2893  				}
  2894  			}
  2895  		} else {
  2896  			if incX == 1 {
  2897  				for j := n - 1; j >= 0; j-- {
  2898  					if diag == blas.NonUnit {
  2899  						x[j] /= a[j*lda+j]
  2900  					}
  2901  					xj := x[j]
  2902  					if j > 0 {
  2903  						c64.AxpyUnitary(-xj, a[j*lda:j*lda+j], x[:j])
  2904  					}
  2905  				}
  2906  			} else {
  2907  				jx := kx + (n-1)*incX
  2908  				for j := n - 1; j >= 0; j-- {
  2909  					if diag == blas.NonUnit {
  2910  						x[jx] /= a[j*lda+j]
  2911  					}
  2912  					if j > 0 {
  2913  						c64.AxpyInc(-x[jx], a[j*lda:j*lda+j], x, uintptr(j), 1, uintptr(incX), 0, uintptr(kx))
  2914  					}
  2915  					jx -= incX
  2916  				}
  2917  			}
  2918  		}
  2919  		return
  2920  	}
  2921  
  2922  	// Form x = inv(Aᴴ)*x.
  2923  	if uplo == blas.Upper {
  2924  		if incX == 1 {
  2925  			for j := 0; j < n; j++ {
  2926  				if diag == blas.NonUnit {
  2927  					x[j] /= cmplx.Conj(a[j*lda+j])
  2928  				}
  2929  				xj := x[j]
  2930  				for i := j + 1; i < n; i++ {
  2931  					x[i] -= xj * cmplx.Conj(a[j*lda+i])
  2932  				}
  2933  			}
  2934  		} else {
  2935  			jx := kx
  2936  			for j := 0; j < n; j++ {
  2937  				if diag == blas.NonUnit {
  2938  					x[jx] /= cmplx.Conj(a[j*lda+j])
  2939  				}
  2940  				xj := x[jx]
  2941  				ix := jx + incX
  2942  				for i := j + 1; i < n; i++ {
  2943  					x[ix] -= xj * cmplx.Conj(a[j*lda+i])
  2944  					ix += incX
  2945  				}
  2946  				jx += incX
  2947  			}
  2948  		}
  2949  	} else {
  2950  		if incX == 1 {
  2951  			for j := n - 1; j >= 0; j-- {
  2952  				if diag == blas.NonUnit {
  2953  					x[j] /= cmplx.Conj(a[j*lda+j])
  2954  				}
  2955  				xj := x[j]
  2956  				for i := 0; i < j; i++ {
  2957  					x[i] -= xj * cmplx.Conj(a[j*lda+i])
  2958  				}
  2959  			}
  2960  		} else {
  2961  			jx := kx + (n-1)*incX
  2962  			for j := n - 1; j >= 0; j-- {
  2963  				if diag == blas.NonUnit {
  2964  					x[jx] /= cmplx.Conj(a[j*lda+j])
  2965  				}
  2966  				xj := x[jx]
  2967  				ix := kx
  2968  				for i := 0; i < j; i++ {
  2969  					x[ix] -= xj * cmplx.Conj(a[j*lda+i])
  2970  					ix += incX
  2971  				}
  2972  				jx -= incX
  2973  			}
  2974  		}
  2975  	}
  2976  }