level3float32.go 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877
  1. // Code generated by "go generate gonum.org/v1/gonum/blas/gonum”; DO NOT EDIT.
  2. // Copyright ©2014 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. "gonum.org/v1/gonum/blas"
  8. "gonum.org/v1/gonum/internal/asm/f32"
  9. )
  10. var _ blas.Float32Level3 = Implementation{}
  11. // Strsm solves one of the matrix equations
  12. // A * X = alpha * B if tA == blas.NoTrans and side == blas.Left
  13. // Aᵀ * X = alpha * B if tA == blas.Trans or blas.ConjTrans, and side == blas.Left
  14. // X * A = alpha * B if tA == blas.NoTrans and side == blas.Right
  15. // X * Aᵀ = alpha * B if tA == blas.Trans or blas.ConjTrans, and side == blas.Right
  16. // where A is an n×n or m×m triangular matrix, X and B are m×n matrices, and alpha is a
  17. // scalar.
  18. //
  19. // At entry to the function, X contains the values of B, and the result is
  20. // stored in-place into X.
  21. //
  22. // No check is made that A is invertible.
  23. //
  24. // Float32 implementations are autogenerated and not directly tested.
  25. func (Implementation) Strsm(s blas.Side, ul blas.Uplo, tA blas.Transpose, d blas.Diag, m, n int, alpha float32, a []float32, lda int, b []float32, ldb int) {
  26. if s != blas.Left && s != blas.Right {
  27. panic(badSide)
  28. }
  29. if ul != blas.Lower && ul != blas.Upper {
  30. panic(badUplo)
  31. }
  32. if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
  33. panic(badTranspose)
  34. }
  35. if d != blas.NonUnit && d != blas.Unit {
  36. panic(badDiag)
  37. }
  38. if m < 0 {
  39. panic(mLT0)
  40. }
  41. if n < 0 {
  42. panic(nLT0)
  43. }
  44. k := n
  45. if s == blas.Left {
  46. k = m
  47. }
  48. if lda < max(1, k) {
  49. panic(badLdA)
  50. }
  51. if ldb < max(1, n) {
  52. panic(badLdB)
  53. }
  54. // Quick return if possible.
  55. if m == 0 || n == 0 {
  56. return
  57. }
  58. // For zero matrix size the following slice length checks are trivially satisfied.
  59. if len(a) < lda*(k-1)+k {
  60. panic(shortA)
  61. }
  62. if len(b) < ldb*(m-1)+n {
  63. panic(shortB)
  64. }
  65. if alpha == 0 {
  66. for i := 0; i < m; i++ {
  67. btmp := b[i*ldb : i*ldb+n]
  68. for j := range btmp {
  69. btmp[j] = 0
  70. }
  71. }
  72. return
  73. }
  74. nonUnit := d == blas.NonUnit
  75. if s == blas.Left {
  76. if tA == blas.NoTrans {
  77. if ul == blas.Upper {
  78. for i := m - 1; i >= 0; i-- {
  79. btmp := b[i*ldb : i*ldb+n]
  80. if alpha != 1 {
  81. f32.ScalUnitary(alpha, btmp)
  82. }
  83. for ka, va := range a[i*lda+i+1 : i*lda+m] {
  84. if va != 0 {
  85. k := ka + i + 1
  86. f32.AxpyUnitary(-va, b[k*ldb:k*ldb+n], btmp)
  87. }
  88. }
  89. if nonUnit {
  90. tmp := 1 / a[i*lda+i]
  91. f32.ScalUnitary(tmp, btmp)
  92. }
  93. }
  94. return
  95. }
  96. for i := 0; i < m; i++ {
  97. btmp := b[i*ldb : i*ldb+n]
  98. if alpha != 1 {
  99. f32.ScalUnitary(alpha, btmp)
  100. }
  101. for k, va := range a[i*lda : i*lda+i] {
  102. if va != 0 {
  103. f32.AxpyUnitary(-va, b[k*ldb:k*ldb+n], btmp)
  104. }
  105. }
  106. if nonUnit {
  107. tmp := 1 / a[i*lda+i]
  108. f32.ScalUnitary(tmp, btmp)
  109. }
  110. }
  111. return
  112. }
  113. // Cases where a is transposed
  114. if ul == blas.Upper {
  115. for k := 0; k < m; k++ {
  116. btmpk := b[k*ldb : k*ldb+n]
  117. if nonUnit {
  118. tmp := 1 / a[k*lda+k]
  119. f32.ScalUnitary(tmp, btmpk)
  120. }
  121. for ia, va := range a[k*lda+k+1 : k*lda+m] {
  122. if va != 0 {
  123. i := ia + k + 1
  124. f32.AxpyUnitary(-va, btmpk, b[i*ldb:i*ldb+n])
  125. }
  126. }
  127. if alpha != 1 {
  128. f32.ScalUnitary(alpha, btmpk)
  129. }
  130. }
  131. return
  132. }
  133. for k := m - 1; k >= 0; k-- {
  134. btmpk := b[k*ldb : k*ldb+n]
  135. if nonUnit {
  136. tmp := 1 / a[k*lda+k]
  137. f32.ScalUnitary(tmp, btmpk)
  138. }
  139. for i, va := range a[k*lda : k*lda+k] {
  140. if va != 0 {
  141. f32.AxpyUnitary(-va, btmpk, b[i*ldb:i*ldb+n])
  142. }
  143. }
  144. if alpha != 1 {
  145. f32.ScalUnitary(alpha, btmpk)
  146. }
  147. }
  148. return
  149. }
  150. // Cases where a is to the right of X.
  151. if tA == blas.NoTrans {
  152. if ul == blas.Upper {
  153. for i := 0; i < m; i++ {
  154. btmp := b[i*ldb : i*ldb+n]
  155. if alpha != 1 {
  156. f32.ScalUnitary(alpha, btmp)
  157. }
  158. for k, vb := range btmp {
  159. if vb == 0 {
  160. continue
  161. }
  162. if nonUnit {
  163. btmp[k] /= a[k*lda+k]
  164. }
  165. f32.AxpyUnitary(-btmp[k], a[k*lda+k+1:k*lda+n], btmp[k+1:n])
  166. }
  167. }
  168. return
  169. }
  170. for i := 0; i < m; i++ {
  171. btmp := b[i*ldb : i*ldb+n]
  172. if alpha != 1 {
  173. f32.ScalUnitary(alpha, btmp)
  174. }
  175. for k := n - 1; k >= 0; k-- {
  176. if btmp[k] == 0 {
  177. continue
  178. }
  179. if nonUnit {
  180. btmp[k] /= a[k*lda+k]
  181. }
  182. f32.AxpyUnitary(-btmp[k], a[k*lda:k*lda+k], btmp[:k])
  183. }
  184. }
  185. return
  186. }
  187. // Cases where a is transposed.
  188. if ul == blas.Upper {
  189. for i := 0; i < m; i++ {
  190. btmp := b[i*ldb : i*ldb+n]
  191. for j := n - 1; j >= 0; j-- {
  192. tmp := alpha*btmp[j] - f32.DotUnitary(a[j*lda+j+1:j*lda+n], btmp[j+1:])
  193. if nonUnit {
  194. tmp /= a[j*lda+j]
  195. }
  196. btmp[j] = tmp
  197. }
  198. }
  199. return
  200. }
  201. for i := 0; i < m; i++ {
  202. btmp := b[i*ldb : i*ldb+n]
  203. for j := 0; j < n; j++ {
  204. tmp := alpha*btmp[j] - f32.DotUnitary(a[j*lda:j*lda+j], btmp[:j])
  205. if nonUnit {
  206. tmp /= a[j*lda+j]
  207. }
  208. btmp[j] = tmp
  209. }
  210. }
  211. }
  212. // Ssymm performs one of the matrix-matrix operations
  213. // C = alpha * A * B + beta * C if side == blas.Left
  214. // C = alpha * B * A + beta * C if side == blas.Right
  215. // where A is an n×n or m×m symmetric matrix, B and C are m×n matrices, and alpha
  216. // is a scalar.
  217. //
  218. // Float32 implementations are autogenerated and not directly tested.
  219. func (Implementation) Ssymm(s blas.Side, ul blas.Uplo, m, n int, alpha float32, a []float32, lda int, b []float32, ldb int, beta float32, c []float32, ldc int) {
  220. if s != blas.Right && s != blas.Left {
  221. panic(badSide)
  222. }
  223. if ul != blas.Lower && ul != blas.Upper {
  224. panic(badUplo)
  225. }
  226. if m < 0 {
  227. panic(mLT0)
  228. }
  229. if n < 0 {
  230. panic(nLT0)
  231. }
  232. k := n
  233. if s == blas.Left {
  234. k = m
  235. }
  236. if lda < max(1, k) {
  237. panic(badLdA)
  238. }
  239. if ldb < max(1, n) {
  240. panic(badLdB)
  241. }
  242. if ldc < max(1, n) {
  243. panic(badLdC)
  244. }
  245. // Quick return if possible.
  246. if m == 0 || n == 0 {
  247. return
  248. }
  249. // For zero matrix size the following slice length checks are trivially satisfied.
  250. if len(a) < lda*(k-1)+k {
  251. panic(shortA)
  252. }
  253. if len(b) < ldb*(m-1)+n {
  254. panic(shortB)
  255. }
  256. if len(c) < ldc*(m-1)+n {
  257. panic(shortC)
  258. }
  259. // Quick return if possible.
  260. if alpha == 0 && beta == 1 {
  261. return
  262. }
  263. if alpha == 0 {
  264. if beta == 0 {
  265. for i := 0; i < m; i++ {
  266. ctmp := c[i*ldc : i*ldc+n]
  267. for j := range ctmp {
  268. ctmp[j] = 0
  269. }
  270. }
  271. return
  272. }
  273. for i := 0; i < m; i++ {
  274. ctmp := c[i*ldc : i*ldc+n]
  275. for j := 0; j < n; j++ {
  276. ctmp[j] *= beta
  277. }
  278. }
  279. return
  280. }
  281. isUpper := ul == blas.Upper
  282. if s == blas.Left {
  283. for i := 0; i < m; i++ {
  284. atmp := alpha * a[i*lda+i]
  285. btmp := b[i*ldb : i*ldb+n]
  286. ctmp := c[i*ldc : i*ldc+n]
  287. for j, v := range btmp {
  288. ctmp[j] *= beta
  289. ctmp[j] += atmp * v
  290. }
  291. for k := 0; k < i; k++ {
  292. var atmp float32
  293. if isUpper {
  294. atmp = a[k*lda+i]
  295. } else {
  296. atmp = a[i*lda+k]
  297. }
  298. atmp *= alpha
  299. f32.AxpyUnitary(atmp, b[k*ldb:k*ldb+n], ctmp)
  300. }
  301. for k := i + 1; k < m; k++ {
  302. var atmp float32
  303. if isUpper {
  304. atmp = a[i*lda+k]
  305. } else {
  306. atmp = a[k*lda+i]
  307. }
  308. atmp *= alpha
  309. f32.AxpyUnitary(atmp, b[k*ldb:k*ldb+n], ctmp)
  310. }
  311. }
  312. return
  313. }
  314. if isUpper {
  315. for i := 0; i < m; i++ {
  316. for j := n - 1; j >= 0; j-- {
  317. tmp := alpha * b[i*ldb+j]
  318. var tmp2 float32
  319. atmp := a[j*lda+j+1 : j*lda+n]
  320. btmp := b[i*ldb+j+1 : i*ldb+n]
  321. ctmp := c[i*ldc+j+1 : i*ldc+n]
  322. for k, v := range atmp {
  323. ctmp[k] += tmp * v
  324. tmp2 += btmp[k] * v
  325. }
  326. c[i*ldc+j] *= beta
  327. c[i*ldc+j] += tmp*a[j*lda+j] + alpha*tmp2
  328. }
  329. }
  330. return
  331. }
  332. for i := 0; i < m; i++ {
  333. for j := 0; j < n; j++ {
  334. tmp := alpha * b[i*ldb+j]
  335. var tmp2 float32
  336. atmp := a[j*lda : j*lda+j]
  337. btmp := b[i*ldb : i*ldb+j]
  338. ctmp := c[i*ldc : i*ldc+j]
  339. for k, v := range atmp {
  340. ctmp[k] += tmp * v
  341. tmp2 += btmp[k] * v
  342. }
  343. c[i*ldc+j] *= beta
  344. c[i*ldc+j] += tmp*a[j*lda+j] + alpha*tmp2
  345. }
  346. }
  347. }
  348. // Ssyrk performs one of the symmetric rank-k operations
  349. // C = alpha * A * Aᵀ + beta * C if tA == blas.NoTrans
  350. // C = alpha * Aᵀ * A + beta * C if tA == blas.Trans or tA == blas.ConjTrans
  351. // where A is an n×k or k×n matrix, C is an n×n symmetric matrix, and alpha and
  352. // beta are scalars.
  353. //
  354. // Float32 implementations are autogenerated and not directly tested.
  355. func (Implementation) Ssyrk(ul blas.Uplo, tA blas.Transpose, n, k int, alpha float32, a []float32, lda int, beta float32, c []float32, ldc int) {
  356. if ul != blas.Lower && ul != blas.Upper {
  357. panic(badUplo)
  358. }
  359. if tA != blas.Trans && tA != blas.NoTrans && tA != blas.ConjTrans {
  360. panic(badTranspose)
  361. }
  362. if n < 0 {
  363. panic(nLT0)
  364. }
  365. if k < 0 {
  366. panic(kLT0)
  367. }
  368. row, col := k, n
  369. if tA == blas.NoTrans {
  370. row, col = n, k
  371. }
  372. if lda < max(1, col) {
  373. panic(badLdA)
  374. }
  375. if ldc < max(1, n) {
  376. panic(badLdC)
  377. }
  378. // Quick return if possible.
  379. if n == 0 {
  380. return
  381. }
  382. // For zero matrix size the following slice length checks are trivially satisfied.
  383. if len(a) < lda*(row-1)+col {
  384. panic(shortA)
  385. }
  386. if len(c) < ldc*(n-1)+n {
  387. panic(shortC)
  388. }
  389. if alpha == 0 {
  390. if beta == 0 {
  391. if ul == blas.Upper {
  392. for i := 0; i < n; i++ {
  393. ctmp := c[i*ldc+i : i*ldc+n]
  394. for j := range ctmp {
  395. ctmp[j] = 0
  396. }
  397. }
  398. return
  399. }
  400. for i := 0; i < n; i++ {
  401. ctmp := c[i*ldc : i*ldc+i+1]
  402. for j := range ctmp {
  403. ctmp[j] = 0
  404. }
  405. }
  406. return
  407. }
  408. if ul == blas.Upper {
  409. for i := 0; i < n; i++ {
  410. ctmp := c[i*ldc+i : i*ldc+n]
  411. for j := range ctmp {
  412. ctmp[j] *= beta
  413. }
  414. }
  415. return
  416. }
  417. for i := 0; i < n; i++ {
  418. ctmp := c[i*ldc : i*ldc+i+1]
  419. for j := range ctmp {
  420. ctmp[j] *= beta
  421. }
  422. }
  423. return
  424. }
  425. if tA == blas.NoTrans {
  426. if ul == blas.Upper {
  427. for i := 0; i < n; i++ {
  428. ctmp := c[i*ldc+i : i*ldc+n]
  429. atmp := a[i*lda : i*lda+k]
  430. if beta == 0 {
  431. for jc := range ctmp {
  432. j := jc + i
  433. ctmp[jc] = alpha * f32.DotUnitary(atmp, a[j*lda:j*lda+k])
  434. }
  435. } else {
  436. for jc, vc := range ctmp {
  437. j := jc + i
  438. ctmp[jc] = vc*beta + alpha*f32.DotUnitary(atmp, a[j*lda:j*lda+k])
  439. }
  440. }
  441. }
  442. return
  443. }
  444. for i := 0; i < n; i++ {
  445. ctmp := c[i*ldc : i*ldc+i+1]
  446. atmp := a[i*lda : i*lda+k]
  447. if beta == 0 {
  448. for j := range ctmp {
  449. ctmp[j] = alpha * f32.DotUnitary(a[j*lda:j*lda+k], atmp)
  450. }
  451. } else {
  452. for j, vc := range ctmp {
  453. ctmp[j] = vc*beta + alpha*f32.DotUnitary(a[j*lda:j*lda+k], atmp)
  454. }
  455. }
  456. }
  457. return
  458. }
  459. // Cases where a is transposed.
  460. if ul == blas.Upper {
  461. for i := 0; i < n; i++ {
  462. ctmp := c[i*ldc+i : i*ldc+n]
  463. if beta == 0 {
  464. for j := range ctmp {
  465. ctmp[j] = 0
  466. }
  467. } else if beta != 1 {
  468. for j := range ctmp {
  469. ctmp[j] *= beta
  470. }
  471. }
  472. for l := 0; l < k; l++ {
  473. tmp := alpha * a[l*lda+i]
  474. if tmp != 0 {
  475. f32.AxpyUnitary(tmp, a[l*lda+i:l*lda+n], ctmp)
  476. }
  477. }
  478. }
  479. return
  480. }
  481. for i := 0; i < n; i++ {
  482. ctmp := c[i*ldc : i*ldc+i+1]
  483. if beta != 1 {
  484. for j := range ctmp {
  485. ctmp[j] *= beta
  486. }
  487. }
  488. for l := 0; l < k; l++ {
  489. tmp := alpha * a[l*lda+i]
  490. if tmp != 0 {
  491. f32.AxpyUnitary(tmp, a[l*lda:l*lda+i+1], ctmp)
  492. }
  493. }
  494. }
  495. }
  496. // Ssyr2k performs one of the symmetric rank 2k operations
  497. // C = alpha * A * Bᵀ + alpha * B * Aᵀ + beta * C if tA == blas.NoTrans
  498. // C = alpha * Aᵀ * B + alpha * Bᵀ * A + beta * C if tA == blas.Trans or tA == blas.ConjTrans
  499. // where A and B are n×k or k×n matrices, C is an n×n symmetric matrix, and
  500. // alpha and beta are scalars.
  501. //
  502. // Float32 implementations are autogenerated and not directly tested.
  503. func (Implementation) Ssyr2k(ul blas.Uplo, tA blas.Transpose, n, k int, alpha float32, a []float32, lda int, b []float32, ldb int, beta float32, c []float32, ldc int) {
  504. if ul != blas.Lower && ul != blas.Upper {
  505. panic(badUplo)
  506. }
  507. if tA != blas.Trans && tA != blas.NoTrans && tA != blas.ConjTrans {
  508. panic(badTranspose)
  509. }
  510. if n < 0 {
  511. panic(nLT0)
  512. }
  513. if k < 0 {
  514. panic(kLT0)
  515. }
  516. row, col := k, n
  517. if tA == blas.NoTrans {
  518. row, col = n, k
  519. }
  520. if lda < max(1, col) {
  521. panic(badLdA)
  522. }
  523. if ldb < max(1, col) {
  524. panic(badLdB)
  525. }
  526. if ldc < max(1, n) {
  527. panic(badLdC)
  528. }
  529. // Quick return if possible.
  530. if n == 0 {
  531. return
  532. }
  533. // For zero matrix size the following slice length checks are trivially satisfied.
  534. if len(a) < lda*(row-1)+col {
  535. panic(shortA)
  536. }
  537. if len(b) < ldb*(row-1)+col {
  538. panic(shortB)
  539. }
  540. if len(c) < ldc*(n-1)+n {
  541. panic(shortC)
  542. }
  543. if alpha == 0 {
  544. if beta == 0 {
  545. if ul == blas.Upper {
  546. for i := 0; i < n; i++ {
  547. ctmp := c[i*ldc+i : i*ldc+n]
  548. for j := range ctmp {
  549. ctmp[j] = 0
  550. }
  551. }
  552. return
  553. }
  554. for i := 0; i < n; i++ {
  555. ctmp := c[i*ldc : i*ldc+i+1]
  556. for j := range ctmp {
  557. ctmp[j] = 0
  558. }
  559. }
  560. return
  561. }
  562. if ul == blas.Upper {
  563. for i := 0; i < n; i++ {
  564. ctmp := c[i*ldc+i : i*ldc+n]
  565. for j := range ctmp {
  566. ctmp[j] *= beta
  567. }
  568. }
  569. return
  570. }
  571. for i := 0; i < n; i++ {
  572. ctmp := c[i*ldc : i*ldc+i+1]
  573. for j := range ctmp {
  574. ctmp[j] *= beta
  575. }
  576. }
  577. return
  578. }
  579. if tA == blas.NoTrans {
  580. if ul == blas.Upper {
  581. for i := 0; i < n; i++ {
  582. atmp := a[i*lda : i*lda+k]
  583. btmp := b[i*ldb : i*ldb+k]
  584. ctmp := c[i*ldc+i : i*ldc+n]
  585. for jc := range ctmp {
  586. j := i + jc
  587. var tmp1, tmp2 float32
  588. binner := b[j*ldb : j*ldb+k]
  589. for l, v := range a[j*lda : j*lda+k] {
  590. tmp1 += v * btmp[l]
  591. tmp2 += atmp[l] * binner[l]
  592. }
  593. ctmp[jc] *= beta
  594. ctmp[jc] += alpha * (tmp1 + tmp2)
  595. }
  596. }
  597. return
  598. }
  599. for i := 0; i < n; i++ {
  600. atmp := a[i*lda : i*lda+k]
  601. btmp := b[i*ldb : i*ldb+k]
  602. ctmp := c[i*ldc : i*ldc+i+1]
  603. for j := 0; j <= i; j++ {
  604. var tmp1, tmp2 float32
  605. binner := b[j*ldb : j*ldb+k]
  606. for l, v := range a[j*lda : j*lda+k] {
  607. tmp1 += v * btmp[l]
  608. tmp2 += atmp[l] * binner[l]
  609. }
  610. ctmp[j] *= beta
  611. ctmp[j] += alpha * (tmp1 + tmp2)
  612. }
  613. }
  614. return
  615. }
  616. if ul == blas.Upper {
  617. for i := 0; i < n; i++ {
  618. ctmp := c[i*ldc+i : i*ldc+n]
  619. if beta != 1 {
  620. for j := range ctmp {
  621. ctmp[j] *= beta
  622. }
  623. }
  624. for l := 0; l < k; l++ {
  625. tmp1 := alpha * b[l*ldb+i]
  626. tmp2 := alpha * a[l*lda+i]
  627. btmp := b[l*ldb+i : l*ldb+n]
  628. if tmp1 != 0 || tmp2 != 0 {
  629. for j, v := range a[l*lda+i : l*lda+n] {
  630. ctmp[j] += v*tmp1 + btmp[j]*tmp2
  631. }
  632. }
  633. }
  634. }
  635. return
  636. }
  637. for i := 0; i < n; i++ {
  638. ctmp := c[i*ldc : i*ldc+i+1]
  639. if beta != 1 {
  640. for j := range ctmp {
  641. ctmp[j] *= beta
  642. }
  643. }
  644. for l := 0; l < k; l++ {
  645. tmp1 := alpha * b[l*ldb+i]
  646. tmp2 := alpha * a[l*lda+i]
  647. btmp := b[l*ldb : l*ldb+i+1]
  648. if tmp1 != 0 || tmp2 != 0 {
  649. for j, v := range a[l*lda : l*lda+i+1] {
  650. ctmp[j] += v*tmp1 + btmp[j]*tmp2
  651. }
  652. }
  653. }
  654. }
  655. }
  656. // Strmm performs one of the matrix-matrix operations
  657. // B = alpha * A * B if tA == blas.NoTrans and side == blas.Left
  658. // B = alpha * Aᵀ * B if tA == blas.Trans or blas.ConjTrans, and side == blas.Left
  659. // B = alpha * B * A if tA == blas.NoTrans and side == blas.Right
  660. // B = alpha * B * Aᵀ if tA == blas.Trans or blas.ConjTrans, and side == blas.Right
  661. // where A is an n×n or m×m triangular matrix, B is an m×n matrix, and alpha is a scalar.
  662. //
  663. // Float32 implementations are autogenerated and not directly tested.
  664. func (Implementation) Strmm(s blas.Side, ul blas.Uplo, tA blas.Transpose, d blas.Diag, m, n int, alpha float32, a []float32, lda int, b []float32, ldb int) {
  665. if s != blas.Left && s != blas.Right {
  666. panic(badSide)
  667. }
  668. if ul != blas.Lower && ul != blas.Upper {
  669. panic(badUplo)
  670. }
  671. if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
  672. panic(badTranspose)
  673. }
  674. if d != blas.NonUnit && d != blas.Unit {
  675. panic(badDiag)
  676. }
  677. if m < 0 {
  678. panic(mLT0)
  679. }
  680. if n < 0 {
  681. panic(nLT0)
  682. }
  683. k := n
  684. if s == blas.Left {
  685. k = m
  686. }
  687. if lda < max(1, k) {
  688. panic(badLdA)
  689. }
  690. if ldb < max(1, n) {
  691. panic(badLdB)
  692. }
  693. // Quick return if possible.
  694. if m == 0 || n == 0 {
  695. return
  696. }
  697. // For zero matrix size the following slice length checks are trivially satisfied.
  698. if len(a) < lda*(k-1)+k {
  699. panic(shortA)
  700. }
  701. if len(b) < ldb*(m-1)+n {
  702. panic(shortB)
  703. }
  704. if alpha == 0 {
  705. for i := 0; i < m; i++ {
  706. btmp := b[i*ldb : i*ldb+n]
  707. for j := range btmp {
  708. btmp[j] = 0
  709. }
  710. }
  711. return
  712. }
  713. nonUnit := d == blas.NonUnit
  714. if s == blas.Left {
  715. if tA == blas.NoTrans {
  716. if ul == blas.Upper {
  717. for i := 0; i < m; i++ {
  718. tmp := alpha
  719. if nonUnit {
  720. tmp *= a[i*lda+i]
  721. }
  722. btmp := b[i*ldb : i*ldb+n]
  723. f32.ScalUnitary(tmp, btmp)
  724. for ka, va := range a[i*lda+i+1 : i*lda+m] {
  725. k := ka + i + 1
  726. if va != 0 {
  727. f32.AxpyUnitary(alpha*va, b[k*ldb:k*ldb+n], btmp)
  728. }
  729. }
  730. }
  731. return
  732. }
  733. for i := m - 1; i >= 0; i-- {
  734. tmp := alpha
  735. if nonUnit {
  736. tmp *= a[i*lda+i]
  737. }
  738. btmp := b[i*ldb : i*ldb+n]
  739. f32.ScalUnitary(tmp, btmp)
  740. for k, va := range a[i*lda : i*lda+i] {
  741. if va != 0 {
  742. f32.AxpyUnitary(alpha*va, b[k*ldb:k*ldb+n], btmp)
  743. }
  744. }
  745. }
  746. return
  747. }
  748. // Cases where a is transposed.
  749. if ul == blas.Upper {
  750. for k := m - 1; k >= 0; k-- {
  751. btmpk := b[k*ldb : k*ldb+n]
  752. for ia, va := range a[k*lda+k+1 : k*lda+m] {
  753. i := ia + k + 1
  754. btmp := b[i*ldb : i*ldb+n]
  755. if va != 0 {
  756. f32.AxpyUnitary(alpha*va, btmpk, btmp)
  757. }
  758. }
  759. tmp := alpha
  760. if nonUnit {
  761. tmp *= a[k*lda+k]
  762. }
  763. if tmp != 1 {
  764. f32.ScalUnitary(tmp, btmpk)
  765. }
  766. }
  767. return
  768. }
  769. for k := 0; k < m; k++ {
  770. btmpk := b[k*ldb : k*ldb+n]
  771. for i, va := range a[k*lda : k*lda+k] {
  772. btmp := b[i*ldb : i*ldb+n]
  773. if va != 0 {
  774. f32.AxpyUnitary(alpha*va, btmpk, btmp)
  775. }
  776. }
  777. tmp := alpha
  778. if nonUnit {
  779. tmp *= a[k*lda+k]
  780. }
  781. if tmp != 1 {
  782. f32.ScalUnitary(tmp, btmpk)
  783. }
  784. }
  785. return
  786. }
  787. // Cases where a is on the right
  788. if tA == blas.NoTrans {
  789. if ul == blas.Upper {
  790. for i := 0; i < m; i++ {
  791. btmp := b[i*ldb : i*ldb+n]
  792. for k := n - 1; k >= 0; k-- {
  793. tmp := alpha * btmp[k]
  794. if tmp == 0 {
  795. continue
  796. }
  797. btmp[k] = tmp
  798. if nonUnit {
  799. btmp[k] *= a[k*lda+k]
  800. }
  801. f32.AxpyUnitary(tmp, a[k*lda+k+1:k*lda+n], btmp[k+1:n])
  802. }
  803. }
  804. return
  805. }
  806. for i := 0; i < m; i++ {
  807. btmp := b[i*ldb : i*ldb+n]
  808. for k := 0; k < n; k++ {
  809. tmp := alpha * btmp[k]
  810. if tmp == 0 {
  811. continue
  812. }
  813. btmp[k] = tmp
  814. if nonUnit {
  815. btmp[k] *= a[k*lda+k]
  816. }
  817. f32.AxpyUnitary(tmp, a[k*lda:k*lda+k], btmp[:k])
  818. }
  819. }
  820. return
  821. }
  822. // Cases where a is transposed.
  823. if ul == blas.Upper {
  824. for i := 0; i < m; i++ {
  825. btmp := b[i*ldb : i*ldb+n]
  826. for j, vb := range btmp {
  827. tmp := vb
  828. if nonUnit {
  829. tmp *= a[j*lda+j]
  830. }
  831. tmp += f32.DotUnitary(a[j*lda+j+1:j*lda+n], btmp[j+1:n])
  832. btmp[j] = alpha * tmp
  833. }
  834. }
  835. return
  836. }
  837. for i := 0; i < m; i++ {
  838. btmp := b[i*ldb : i*ldb+n]
  839. for j := n - 1; j >= 0; j-- {
  840. tmp := btmp[j]
  841. if nonUnit {
  842. tmp *= a[j*lda+j]
  843. }
  844. tmp += f32.DotUnitary(a[j*lda:j*lda+j], btmp[:j])
  845. btmp[j] = alpha * tmp
  846. }
  847. }
  848. }