level2cmplx64.go 62 KB

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