aboutsummaryrefslogtreecommitdiff
path: root/src/libmath/blas2body
diff options
context:
space:
mode:
Diffstat (limited to 'src/libmath/blas2body')
-rw-r--r--src/libmath/blas2body256
1 files changed, 256 insertions, 0 deletions
diff --git a/src/libmath/blas2body b/src/libmath/blas2body
new file mode 100644
index 0000000..45baf67
--- /dev/null
+++ b/src/libmath/blas2body
@@ -0,0 +1,256 @@
+/* general matrix multiply */
+error
+func(gemv)(uint flag, INT nrow, INT ncol, FLOAT a, FLOAT *m, INT incm, FLOAT *x, INT incx, FLOAT b, FLOAT *y, INT incy)
+{
+ INT r, c, nr, nc;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX, INCY) row[I] = m + (r+I)*incm, reg[I] = 0;
+#define TERM(J, I, INCX, INCY) row[I][c+J] * x[INCX*(c+J)]
+#define KERN(I, INCX, INCY, LENGTH) reg[I] += EXPAND(LENGTH,0,TERM,+,I,INCX,INCY);
+#define FINI(I, INCX, INCY) y[INCY*(r+I)] = b*y[INCY*(r+I)] + a*reg[I];
+
+ if (!flag) {
+ BODY_RECT();
+ } else {
+ func(scale)(ncol, b, y, incy);
+#undef KERN
+#undef FINI
+#undef INIT
+#undef TERM
+#define INIT(I, INCX, INCY) row[I] = m + (r+I)*incm, reg[I] = a*x[INCX*(r+I)];
+#define TERM(J, I, INCX, INCY) row[I][c+J] * reg[I]
+#define KERN(I, INCX, INCY, LENGTH) y[INCY*(c+I)] += EXPAND(LENGTH,0,TERM,+,I,INCX,INCY);
+#define FINI(I, INCX, INCY)
+ BODY_RECT();
+ }
+
+ return 0;
+#undef INIT
+#undef TERM
+#undef KERN
+#undef FINI
+}
+
+/* symmetric matrix vector multiply (different layouts) */
+void
+func(symv)(uint upper, INT n, FLOAT a, FLOAT *m, INT incm, FLOAT *x, INT incx, FLOAT b, FLOAT *y, INT incy)
+{
+ INT r, c, nr, nc, i;
+ FLOAT *row[UNROW], reg1[UNROW], reg2[UNROW];
+
+#define INIT(I, INCX, INCY) row[I] = m + (r XX I)*incm, reg1[I] = 0, reg2[I] = 0;
+#define TERM1(J, I, INCX, INCY) row[I][c XX J]*x[INCX*(c XX J)]
+#define TERM2(J, I, INCX, INCY) row[I][c XX J]*x[INCX*((n-c-1) XX J)]
+#define KERN(I, INCX, INCY, REPEAT) reg1[I] += EXPAND(REPEAT,0,TERM1,+,I,INCX,INCY), \
+ reg2[I] += EXPAND(REPEAT,0,TERM2,+,I,INCX,INCY);
+#define FINI(I, INCX, INCY) y[INCY*(r+I)] += a*(reg1[I] + row[I][r]*x[INCX*r]), \
+ y[INCY*(n-r-1+I)] += a*reg2[I];
+
+ func(scale)(n, b, y, incy);
+#define XX +
+ if (!upper) {
+ BODY_LOTRI_XY();
+ } else {
+#undef XX
+#define XX -
+ BODY_UPTRI_XY();
+ }
+#undef XX
+
+#undef INIT
+}
+
+void
+func(spmv)(uint upper, INT n, FLOAT a, FLOAT *m, FLOAT *x, INT incx, FLOAT b, FLOAT *y, INT incy)
+{
+ INT r, c, nr, nc, i;
+ FLOAT *row[UNROW], reg1[UNROW], reg2[UNROW];
+
+#define INIT(I, INCX, INCY) row[I] = m + ((r+I)*(r+I+1))*(1/2), reg1[I] = 0, reg2[I] = 0;
+
+ func(scale)(n, b, y, incy);
+#define XX +
+ if (!upper) {
+ BODY_LOTRI_XY();
+ } else {
+#undef XX
+#undef INIT
+#define INIT(I, INCX, INCY) row[I] = m + ((2*n-r-I)*(r+I+1))*(1/2), reg1[I] = 0, reg2[I] = 0;
+#define XX -
+ BODY_UPTRI_XY();
+ }
+#undef XX
+
+#undef TERM
+#undef INIT
+#undef KERN
+#undef FINI
+}
+
+/* triangular multiply (different layouts) */
+void
+func(trmv)(uint upper, INT n, FLOAT *m, INT incm, FLOAT *x, INT incx)
+{
+ INT r, c, nr, nc, i;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX) row[I] = m + (r XX I)*incm, reg[I] = 0;
+#define TERM(J, I, INCX) row[I][c XX J]*x[INCX*(c XX J)]
+#define KERN(I, INCX, REPEAT) reg[I] += EXPAND(REPEAT,0,TERM,+,I,INCX);
+#define FINI(I, INCX) x[INCX*(r XX I)] = row[I][r XX I]*x[INCX*(r XX I)] + reg[I];
+
+#define XX +
+ if (!upper) {
+ BODY_LOTRI();
+ } else {
+#undef XX
+#define XX -
+ BODY_UPTRI();
+ }
+#undef XX
+
+#undef INIT
+}
+
+void
+func(tpmv)(uint upper, INT n, FLOAT *m, FLOAT *x, INT incx)
+{
+ INT r, c, nr, nc, i;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX) row[I] = m + ((r+I)*(r+I+1))*(1/2), reg[I] = 0;
+
+#define XX +
+ if (!upper) {
+ BODY_LOTRI();
+ } else {
+#undef XX
+#undef INIT
+#define INIT(I, INCX) row[I] = m + ((2*n-r-I)*(r+I+1))*(1/2), reg[I] = 0;
+#define XX -
+ BODY_UPTRI();
+ }
+#undef XX
+
+#undef INIT
+#undef TERM
+#undef KERN
+#undef FINI
+}
+
+/* triangular solve (different layouts) */
+void
+func(trsv)(uint upper, INT n, FLOAT *m, INT incm, FLOAT *x, INT incx)
+{
+ INT r, c, nr, nc, i;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX) row[I] = m + (r XX I)*incm, reg[I] = 0;
+#define TERM(J, I, INCX) row[I][c XX J]*x[c XX J]
+#define KERN(I, INCX, REPEAT) reg[I] += EXPAND(REPEAT,0,TERM,+,I,INCX);
+#define SOLN(J, I, INCX) reg[J] += row[I][r XX I]*x[INCX*(r XX I)]
+#define FINI(I, INCX) x[INCX*(r XX I)] = reg[I] / row[I][r XX I]; EXPAND_TRI(UNROW,INC(I),SOLN,;,I,INCX);
+
+#define XX +
+ if (!upper) {
+ BODY_LOTRI();
+ } else {
+#undef XX
+#define XX -
+ BODY_UPTRI();
+ }
+#undef XX
+
+#undef INIT
+}
+
+void
+func(tpsv)(uint upper, INT n, FLOAT *m, FLOAT *x, INT incx)
+{
+ INT r, c, nr, nc, i;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX) row[I] = m + ((r+I)*(r+I+1))*(1/2), reg[I] = 0;
+
+#define XX +
+ if (!upper) {
+ BODY_LOTRI();
+ } else {
+#undef XX
+#undef INIT
+#define INIT(I, INCX) row[I] = m + ((2*n-r-I)*(r+I+1))*(1/2), reg[I] = 0;
+#define XX -
+ BODY_UPTRI();
+ }
+#undef XX
+
+#undef INIT
+#undef TERM
+#undef KERN
+#undef SOLN
+#undef FINI
+}
+
+/* rank 1 update */
+void
+func(ger)(INT nrow, INT ncol, FLOAT a, FLOAT *x, INT incx, FLOAT *y, INT incy, FLOAT *m, INT incm)
+{
+ INT r, c, nr, nc;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX, INCY) row[I] = m + (r+I)*incm, reg[I] = a*x[INCX*(r+I)];
+#define TERM(J, I, INCX, INCY) row[I][c+J] += reg[I] * y[INCY*(c+J)]
+#define KERN(I, INCX, INCY, LENGTH) EXPAND(LENGTH,0,TERM,;,I,INCX, INCY);
+#define FINI(I, ...)
+
+ BODY_RECT();
+
+#undef INIT
+#undef TERM
+#undef KERN
+#undef FINI
+}
+
+/* symmetric rank 1 update (different memory layouts) */
+void
+func(syr)(uint upper, INT n, FLOAT a, FLOAT *x, INT incx, FLOAT *m, INT incm)
+{
+ INT r, c, nr, nc;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX) row[I] = m + (r XX I)*incm, reg[I] = a*x[INCX*(r XX I)];
+#define TERM(J, I, INCX) row[I][c XX J] += reg[I] * x[INCX*(c XX J)]
+#define KERN(I, INCX, LENGTH) EXPAND(LENGTH,0,TERM,;,I,INCX);
+#define FINI(I, ...)
+
+#define XX +
+ if (!upper) {
+ BODY_LOTRI();
+ } else {
+#undef XX
+#define XX -
+ BODY_UPTRI();
+ }
+#undef XX
+
+#undef INIT
+}
+
+void
+func(spr)(uint upper, INT n, FLOAT a, FLOAT *x, INT incx, FLOAT *m)
+{
+ INT r, c, nr, nc;
+ FLOAT *row[UNROW], reg[UNROW];
+#define INIT(I, INCX) row[I] = m + ((r+I)*(r+I+1))*(1/2), reg[I] = 0;
+
+#define XX +
+ if (!upper) {
+ BODY_LOTRI();
+ } else {
+#undef XX
+#undef INIT
+#define INIT(I, INCX) row[I] = m + ((2*n-r-I)*(r+I+1))*(1/2), reg[I] = 0;
+#define XX -
+ BODY_UPTRI();
+ }
+#undef XX
+
+#undef INIT
+#undef TERM
+#undef KERN
+#undef FINI
+}
generated by cgit v1.2.1 (git 2.18.0) at 2024-05-21 08:23:58 +0000